From 0a275268db617b8a4ffae4f7e802e412500581c6 Mon Sep 17 00:00:00 2001
From: Christophe Palmann <christophe.palmann@c-s.fr>
Date: Fri, 9 Oct 2015 16:45:35 +0200
Subject: [PATCH] ENH: use target vectors for ReciprocalCoherencyMatrix
ReciprocalCovarianceMatrix calculation
---
...nclairToReciprocalCoherencyMatrixFunctor.h | 34 +++++++++++--------
...clairToReciprocalCovarianceMatrixFunctor.h | 22 ++++++++----
2 files changed, 36 insertions(+), 20 deletions(-)
diff --git a/Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCoherencyMatrixFunctor.h b/Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCoherencyMatrixFunctor.h
index e1d6743c6f..b3754f8d88 100644
--- a/Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCoherencyMatrixFunctor.h
+++ b/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);
}
diff --git a/Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCovarianceMatrixFunctor.h b/Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCovarianceMatrixFunctor.h
index 8ad714abef..632033676a 100644
--- a/Modules/Filtering/Polarimetry/include/otbSinclairToReciprocalCovarianceMatrixFunctor.h
+++ b/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);
}
--
GitLab