From 88c2b29d765b75e3f5dce962c1aff1621ff26963 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=A9dric=20Traizet?= <cedric.traizet@c-s.fr>
Date: Wed, 4 Sep 2019 17:17:33 +0200
Subject: [PATCH] ENH: add parameter to disable whitening in pca

---
 .../app/otbDimensionalityReduction.cxx             |  5 +++++
 .../include/otbPCAImageFilter.h                    |  4 ++++
 .../include/otbPCAImageFilter.hxx                  | 14 +++++++++-----
 3 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx b/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
index 24f5951f6b..3e5cbf8e7a 100644
--- a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
+++ b/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
@@ -130,6 +130,10 @@ private:
     AddParameter(ParameterType_OutputFilename, "method.pca.outeigenvalues", "Output file containing eigenvalues (txt format)");
     SetParameterDescription("method.pca.outeigenvalues", "Output file containing eigenvalues (txt format).");
     MandatoryOff("method.pca.outeigenvalues");
+    AddParameter(ParameterType_Bool, "method.pca.whiten", "Perform pca whitening");
+    SetParameterDescription("method.pca.whiten", "Perform whitening and ensure uncorrelated outputs with unit component wise variances");
+    SetParameterInt("method.pca.whiten", 1);
+    MandatoryOff("method.pca.whiten");
     
     AddChoice("method.napca", "NA-PCA");
     SetParameterDescription("method.napca", "Noise Adjusted Principal Component Analysis.");
@@ -257,6 +261,7 @@ private:
 
         filter->SetInput(GetParameterFloatVectorImage("in"));
         filter->SetNumberOfPrincipalComponentsRequired(nbComp);
+        filter->SetWhitening(GetParameterInt("method.pca.whiten"));
         
         // Center AND reduce the input data.
         if (normalize)
diff --git a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h
index 3368163846..b90faae536 100644
--- a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h
+++ b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h
@@ -97,6 +97,9 @@ public:
    */
   itkSetMacro(NumberOfPrincipalComponentsRequired, unsigned int);
   itkGetMacro(NumberOfPrincipalComponentsRequired, unsigned int);
+  
+  itkSetMacro(Whitening, bool);
+  itkGetMacro(Whitening, bool);
 
   itkGetMacro(CovarianceEstimator, CovarianceEstimatorFilterType *);
   itkGetMacro(Transformer, TransformFilterType *);
@@ -197,6 +200,7 @@ protected:
   bool         m_GivenCovarianceMatrix;
   bool         m_GivenTransformationMatrix;
   bool         m_IsTransformationMatrixForward;
+  bool         m_Whitening;
 
   VectorType m_MeanValues;
   VectorType m_StdDevValues;
diff --git a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx
index 3c2487119c..51633d1eaf 100644
--- a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx
+++ b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx
@@ -39,7 +39,7 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
   this->SetNumberOfRequiredInputs(1);
 
   m_NumberOfPrincipalComponentsRequired = 0;
-
+  m_Whitening = true;
   m_UseNormalization = false;
   m_UseVarianceForNormalization = false;
   m_GivenMeanValues = false;
@@ -181,7 +181,7 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
             {
               for (unsigned int c = 0; c < numberOfComponent; ++c)
               {
-                m_CovarianceMatrix(r,c) = cov(r,c) / (std::sqrt(cov(r,r)) *std::sqrt( cov(c,c)));
+                m_CovarianceMatrix(r,c) = cov(r,c) / std::sqrt(cov(r,r)*cov(c,c));
               }
             }
           }
@@ -376,19 +376,23 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
 
   m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
   for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
-    m_EigenValues[m_NumberOfPrincipalComponentsRequired - 1 - i] = static_cast< RealType >( valP(i, i) );
+    m_EigenValues[m_NumberOfPrincipalComponentsRequired - 1 - i] = static_cast< RealType >( vectValP[i] );
   
   /* We used normalized PCA */
   for ( unsigned int i = 0; i < valP.rows(); ++i )
   {
     if (  valP(i, i) > 0. )
     {
-      valP(i, i) = 1. / std::sqrt( valP(i, i) );
+      if (m_Whitening)
+        valP(i, i) = 1. / std::sqrt( valP(i, i) );
     }
     else if ( valP(i, i) < 0. )
     {
       otbMsgDebugMacro( << "ValP(" << i << ") neg : " << valP(i, i) << " taking abs value" );
-      valP(i, i) = 1. / std::sqrt( std::abs( valP(i, i) ) );
+      if (m_Whitening)
+        valP(i, i) = 1. / std::sqrt( std::abs( valP(i, i) ) );
+      else
+        valP(i, i) = std::abs( valP(i, i));
     }
     else
     {
-- 
GitLab