ENH: use target vectors for ReciprocalCoherencyMatrix ReciprocalCovarianceMatrix calculation

Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCoherencyMatrixFunctor.h

@@ -19,7 +19,8 @@
 #define __otbSinclairToReciprocalCoherencyMatrixFunctor_h
 
 #include "vcl_complex.h"
-#include "itkMacro.h"
+#include "otbMath.h"
+#include "vnl/vnl_matrix.h"
 
 namespace otb
 {
@@ -62,6 +63,7 @@ class SinclairToReciprocalCoherencyMatrixFunctor
 public:
   /** Some typedefs. */
   typedef typename std::complex <double>           ComplexType;
+  typedef vnl_matrix<ComplexType>       		   VNLMatrixType;
   typedef typename TOutput::ValueType              OutputValueType;
   
   itkStaticConstMacro(NumberOfComponentsPerPixel, unsigned int, 6);
@@ -75,19 +77,23 @@ public:
     const ComplexType S_hh = static_cast<ComplexType>(Shh);
     const ComplexType S_hv = static_cast<ComplexType>(Shv);
     const ComplexType S_vv = static_cast<ComplexType>(Svv);
-
-    const ComplexType HHPlusVV  = S_hh + S_vv;
-    const ComplexType HHMinusVV = S_hh - S_vv;
-    const ComplexType twoHV     = ComplexType( 2.0 ) * S_hv;
-
-    result[0] = static_cast<OutputValueType>( std::norm(HHPlusVV) );
-    result[1] = static_cast<OutputValueType>( HHPlusVV * vcl_conj(HHMinusVV) );
-    result[2] = static_cast<OutputValueType>( HHPlusVV * vcl_conj(twoHV) );
-    result[3] = static_cast<OutputValueType>( std::norm(HHMinusVV) );
-    result[4] = static_cast<OutputValueType>( HHMinusVV *vcl_conj(twoHV) );
-    result[5] = static_cast<OutputValueType>( std::norm(twoHV) );
-
-    result /= 2.0;
+   
+    
+    VNLMatrixType f3p(3, 1, 0.);
+    f3p[0][0]= (S_hh + S_vv) / ComplexType( std::sqrt(2.0) , 0.0);
+    f3p[1][0]= (S_hh - S_vv) / ComplexType( std::sqrt(2.0) , 0.0);
+    f3p[2][0]= ComplexType( std::sqrt(2.0) , 0.0) * S_hv;
+
+
+    VNLMatrixType res = f3p*f3p.conjugate_transpose();
+    
+    result[0] = static_cast<OutputValueType>( res[0][0] );
+    result[1] = static_cast<OutputValueType>( res[0][1] );
+    result[2] = static_cast<OutputValueType>( res[0][2] );
+    result[3] = static_cast<OutputValueType>( res[1][1] );
+    result[4] = static_cast<OutputValueType>( res[1][2] );
+    result[5] = static_cast<OutputValueType>( res[2][2] );
+    
 
     return (result);
   }

Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCovarianceMatrixFunctor.h

@@ -19,6 +19,8 @@
 #define __otbSinclairToReciprocalCovarianceMatrixFunctor_h
 
 #include "vcl_complex.h"
+#include "otbMath.h"
+#include "vnl/vnl_matrix.h"
 
 namespace otb
 {
@@ -61,6 +63,7 @@ class SinclairToReciprocalCovarianceMatrixFunctor
 public:
   /** Some typedefs. */
   typedef typename std::complex <double>           ComplexType;
+  typedef vnl_matrix<ComplexType>       		   VNLMatrixType;
   typedef typename TOutput::ValueType              OutputValueType;
   inline TOutput operator ()(const TInput1& Shh, const TInput2& Shv, const TInput3& Svv)
   {
@@ -72,12 +75,19 @@ public:
     const ComplexType S_hv = static_cast<ComplexType>(Shv);
     const ComplexType S_vv = static_cast<ComplexType>(Svv);
     
-    result[0] = static_cast<OutputValueType>( std::norm( S_hh ) );
-    result[1] = static_cast<OutputValueType>( vcl_sqrt(2.0)*S_hh*vcl_conj(S_hv) );
-    result[2] = static_cast<OutputValueType>( S_hh*vcl_conj(S_vv) );
-    result[3] = static_cast<OutputValueType>( 2.0*std::norm( S_hv ) );
-    result[4] = static_cast<OutputValueType>( vcl_sqrt(2.0)*S_hv*vcl_conj(S_vv) );
-    result[5] = static_cast<OutputValueType>( std::norm( S_vv ) );
+    VNLMatrixType f3l(3, 1, 0.);
+    f3l[0][0]=S_hh;
+    f3l[1][0]=ComplexType(std::sqrt(2.0),0.0)*S_hv;
+    f3l[2][0]=S_vv;
+    
+    VNLMatrixType res = f3l*f3l.conjugate_transpose();
+    
+    result[0] = static_cast<OutputValueType>( res[0][0] );
+    result[1] = static_cast<OutputValueType>( res[0][1] );
+    result[2] = static_cast<OutputValueType>( res[0][2] );
+    result[3] = static_cast<OutputValueType>( res[1][1] );
+    result[4] = static_cast<OutputValueType>( res[1][2] );
+    result[5] = static_cast<OutputValueType>( res[2][2] );
 
     return (result);
   }