Orfeo Toolbox  4.0
otbReciprocalHAlphaImageFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 
19 #ifndef __ReciprocalHAlphaImageFilter_h
20 #define __ReciprocalHAlphaImageFilter_h
21 
23 #include "otbMath.h"
24 #include "vnl/algo/vnl_complex_eigensystem.h"
25 #include <algorithm>
26 
27 namespace otb
28  {
29 
30 namespace Functor {
31 
56 template< class TInput, class TOutput>
58 {
59 public:
60  typedef typename std::complex<double> ComplexType;
61  typedef vnl_matrix<ComplexType> VNLMatrixType;
62  typedef vnl_vector<ComplexType> VNLVectorType;
63  typedef vnl_vector<double> VNLDoubleVectorType;
64  typedef std::vector<double> VectorType;
65  typedef typename TOutput::ValueType OutputValueType;
66 
67 
68  inline TOutput operator()( const TInput & Coherency ) const
69  {
70  TOutput result;
71  result.SetSize(m_NumberOfComponentsPerPixel);
72 
73  const double T0 = static_cast<double>(Coherency[0].real());
74  const double T1 = static_cast<double>(Coherency[3].real());
75  const double T2 = static_cast<double>(Coherency[5].real());
76 
77  VNLMatrixType vnlMat(3, 3, 0.);
78  vnlMat[0][0] = ComplexType(T0, 0.);
79  vnlMat[0][1] = std::conj(ComplexType(Coherency[1]));
80  vnlMat[0][2] = std::conj(ComplexType(Coherency[2]));
81  vnlMat[1][0] = ComplexType(Coherency[1]);
82  vnlMat[1][1] = ComplexType(T1, 0.);
83  vnlMat[1][2] = std::conj(ComplexType(Coherency[4]));
84  vnlMat[2][0] = ComplexType(Coherency[2]);
85  vnlMat[2][1] = ComplexType(Coherency[4]);
86  vnlMat[2][2] = ComplexType(T2, 0.);
87 
88  // Only compute the left symetry to respect the previous Hermitian Analisys code
89  vnl_complex_eigensystem syst(vnlMat, false, true);
90  const VNLMatrixType eigenVectors( syst.L );
91  const VNLVectorType eigenValues(syst.W);
92 
93  // Entropy estimation
94  double totalEigenValues(0.0);
95  double p[3];
96  double entropy;
97  double alpha;
98  double anisotropy;
99 
100  // Sort eigen values in decreasing order
101  VectorType sortedRealEigenValues(3, eigenValues[0].real());
102  sortedRealEigenValues[1] = eigenValues[1].real();
103  sortedRealEigenValues[2] = eigenValues[2].real();
104  std::sort( sortedRealEigenValues.begin(), sortedRealEigenValues.end() );
105  std::reverse( sortedRealEigenValues.begin(), sortedRealEigenValues.end() );
106 
107  // Extract the first component of each the eigen vector sorted by eigen value decrease order
108  VNLVectorType sortedGreaterEigenVector(3, eigenVectors[0][0]);
109  for(unsigned int i=0; i<3; ++i)
110  {
111  if( vcl_abs( eigenValues[1].real()-sortedRealEigenValues[i] ) < m_Epsilon )
112  {
113  sortedGreaterEigenVector[i] = eigenVectors[1][0];
114  }
115  else if( vcl_abs( eigenValues[2].real()-sortedRealEigenValues[i] ) < m_Epsilon )
116  {
117  sortedGreaterEigenVector[i] = eigenVectors[2][0];
118  }
119  }
120 
121  totalEigenValues = sortedRealEigenValues[0] + sortedRealEigenValues[1] + sortedRealEigenValues[2];
122  if (totalEigenValues <m_Epsilon)
123  {
124  totalEigenValues = m_Epsilon;
125  }
126 
127  for (unsigned int k = 0; k < 3; ++k)
128  {
129  p[k] = std::max(sortedRealEigenValues[k], 0.) / totalEigenValues;
130  }
131 
132  if ( (p[0] < m_Epsilon) || (p[1] < m_Epsilon) || (p[2] < m_Epsilon) )
133  {
134  entropy =0.0;
135  }
136  else
137  {
138  entropy = p[0]*log(p[0]) + p[1]*log(p[1]) + p[2]*log(p[2]);
139  entropy /= -log(3.);
140  }
141 
142  // alpha estimation
143  double val0, val1, val2;
144  double a0, a1, a2;
145 
146  for(unsigned int k = 0; k < 3; ++k)
147  {
148  p[k] = sortedRealEigenValues[k] / totalEigenValues;
149 
150  if (p[k] < 0.)
151  {
152  p[k] = 0.;
153  }
154 
155  if (p[k] > 1.)
156  {
157  p[k] = 1.;
158  }
159  }
160 
161  val0 = std::abs(sortedGreaterEigenVector[0]);
162  a0=acos(vcl_abs(val0)) * CONST_180_PI;
163 
164  val1 = std::abs(sortedGreaterEigenVector[1]);
165  a1=acos(vcl_abs(val1)) * CONST_180_PI;
166 
167  val2= std::abs(sortedGreaterEigenVector[2]);
168  a2=acos(vcl_abs(val2)) * CONST_180_PI;
169 
170  alpha=p[0]*a0 + p[1]*a1 + p[2]*a2;
171 
172  if (alpha>90)
173  {
174  alpha=0.;
175  }
176 
177  // Anisotropy estimation
178  anisotropy=(sortedRealEigenValues[1] - sortedRealEigenValues[2])/(sortedRealEigenValues[1] + sortedRealEigenValues[2] + m_Epsilon);
179 
180  result[0] = static_cast<OutputValueType>(entropy);
181  result[1] = static_cast<OutputValueType>(alpha);
182  result[2] = static_cast<OutputValueType>(anisotropy);
183 
184  return result;
185  }
186 
187  unsigned int GetOutputSize()
188  {
190  }
191 
194 
197 
198 private:
199  itkStaticConstMacro(m_NumberOfComponentsPerPixel, unsigned int, 3);
200  const double m_Epsilon;
201 };
202 }
203 
204 
214 template <class TInputImage, class TOutputImage>
215 class ITK_EXPORT ReciprocalHAlphaImageFilter :
216  public otb::UnaryFunctorImageFilter<TInputImage, TOutputImage, Functor::ReciprocalHAlphaFunctor<
217  typename TInputImage::PixelType, typename TOutputImage::PixelType> >
218 {
219 public:
222  typedef typename Functor::ReciprocalHAlphaFunctor<
223  typename TInputImage::PixelType, typename TOutputImage::PixelType> FunctionType;
227 
229  itkNewMacro(Self);
230 
233 
234 protected:
237 
238 private:
239  ReciprocalHAlphaImageFilter(const Self&); //purposely not implemented
240  void operator=(const Self&); //purposely not implemented
241 
242 };
243 
244 } // end namespace otb
245 
246 #endif

Generated at Sat Mar 8 2014 16:15:18 for Orfeo Toolbox with doxygen 1.8.3.1