OTB  6.7.0
Orfeo Toolbox
otbReciprocalHAlphaImageFilter.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 
22 #ifndef otbReciprocalHAlphaImageFilter_h
23 #define otbReciprocalHAlphaImageFilter_h
24 
25 #include "otbMath.h"
26 #include "vnl/algo/vnl_complex_eigensystem.h"
27 #include <algorithm>
28 #include <vector>
29 #include <complex>
30 
31 #include "otbFunctorImageFilter.h"
32 
33 namespace otb
34  {
35 
36 namespace Functor {
37 
63 template< class TInput, class TOutput>
65 {
66 public:
67  typedef typename std::complex<double> ComplexType;
68  typedef vnl_matrix<ComplexType> VNLMatrixType;
69  typedef vnl_vector<ComplexType> VNLVectorType;
70  typedef vnl_vector<double> VNLDoubleVectorType;
71  typedef std::vector<double> VectorType;
72  typedef typename TOutput::ValueType OutputValueType;
73 
74 
75  inline void operator()(TOutput& result, const TInput& Coherency) const
76  {
77  const double T0 = static_cast<double>(Coherency[0].real());
78  const double T1 = static_cast<double>(Coherency[3].real());
79  const double T2 = static_cast<double>(Coherency[5].real());
80 
81  VNLMatrixType vnlMat(3, 3, 0.);
82  vnlMat[0][0] = ComplexType(T0, 0.);
83  vnlMat[0][1] = ComplexType(Coherency[1]);
84  vnlMat[0][2] = ComplexType(Coherency[2]);
85  vnlMat[1][0] = std::conj(ComplexType(Coherency[1]));
86  vnlMat[1][1] = ComplexType(T1, 0.);
87  vnlMat[1][2] = ComplexType(Coherency[4]);
88  vnlMat[2][0] = std::conj(ComplexType(Coherency[2]));
89  vnlMat[2][1] = std::conj(ComplexType(Coherency[4]));
90  vnlMat[2][2] = ComplexType(T2, 0.);
91 
92  // Only compute the left symmetry to respect the previous Hermitian Analisys code
93  vnl_complex_eigensystem syst(vnlMat, false, true);
94  const VNLMatrixType eigenVectors( syst.L );
95  const VNLVectorType eigenValues(syst.W);
96 
97  // Entropy estimation
98  double totalEigenValues(0.0);
99  double p[3];
100  double plog[3];
101  double entropy;
102  double alpha;
103  double anisotropy;
104 
105  // Sort eigen values in decreasing order
106  VectorType sortedRealEigenValues(3, eigenValues[0].real());
107  sortedRealEigenValues[1] = eigenValues[1].real();
108  sortedRealEigenValues[2] = eigenValues[2].real();
109  std::sort( sortedRealEigenValues.begin(), sortedRealEigenValues.end() );
110  std::reverse( sortedRealEigenValues.begin(), sortedRealEigenValues.end() );
111 
112  // Extract the first component of each the eigen vector sorted by eigen value decrease order
113  VNLVectorType sortedGreaterEigenVector(3, eigenVectors[0][0]);
114  for(unsigned int i=0; i<3; ++i)
115  {
116  if( std::abs( eigenValues[1].real()-sortedRealEigenValues[i] ) < m_Epsilon )
117  {
118  sortedGreaterEigenVector[i] = eigenVectors[1][0];
119  }
120  else if( std::abs( eigenValues[2].real()-sortedRealEigenValues[i] ) < m_Epsilon )
121  {
122  sortedGreaterEigenVector[i] = eigenVectors[2][0];
123  }
124  }
125 
126  totalEigenValues = 0.0;
127  for (unsigned int k = 0; k < 3; ++k)
128  {
129  sortedRealEigenValues[k] = std::max(sortedRealEigenValues[k], 0.);
130  totalEigenValues += sortedRealEigenValues[k];
131  }
132 
133 
134  for (unsigned int k = 0; k < 3; ++k)
135  {
136  p[k] = sortedRealEigenValues[k] / totalEigenValues;
137 
138  if (p[k]<m_Epsilon) //n=log(n)-->0 when n-->0
139  plog[k]=0.0;
140  else
141  plog[k]=-p[k]*log(p[k])/log(3.0);
142 
143  }
144 
145  entropy = 0.0;
146  for (unsigned int k = 0; k < 3; ++k)
147  entropy += plog[k];
148 
149  // alpha estimation
150  double val0, val1, val2;
151  double a0, a1, a2;
152 
153  val0 = std::abs(sortedGreaterEigenVector[0]);
154  a0=acos(std::abs(val0)) * CONST_180_PI;
155 
156  val1 = std::abs(sortedGreaterEigenVector[1]);
157  a1=acos(std::abs(val1)) * CONST_180_PI;
158 
159  val2= std::abs(sortedGreaterEigenVector[2]);
160  a2=acos(std::abs(val2)) * CONST_180_PI;
161 
162  alpha=p[0]*a0 + p[1]*a1 + p[2]*a2;
163 
164  // Anisotropy estimation
165  anisotropy=(sortedRealEigenValues[1] - sortedRealEigenValues[2])/(sortedRealEigenValues[1] + sortedRealEigenValues[2] + m_Epsilon);
166 
167  result[0] = static_cast<OutputValueType>(entropy);
168  result[1] = static_cast<OutputValueType>(alpha);
169  result[2] = static_cast<OutputValueType>(anisotropy);
170  }
171 
172  constexpr size_t OutputSize(...) const
173  {
174  // Size of the result (entropy, alpha, anisotropy)
175  return 3;
176  }
177 
178  private:
179  static constexpr double m_Epsilon = 1e-6;
180 };
181 } // namespace Functor
182 
195 template <typename TInputImage, typename TOutputImage>
197 } // end namespace otb
198 
199 #endif
void operator()(TOutput &result, const TInput &Coherency) const
A generic functor filter templated by its functor.
constexpr double CONST_180_PI
Definition: otbMath.h:56