BUG: center the data by default and use the correct correlation matrix when the data is reduced

Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx

@@ -96,7 +96,7 @@ private:
                   "\"Kernel maximum autocorrelation factor and minimum noise fraction transformations,\" IEEE Transactions on Image Processing, vol. 20, no. 3, pp. 612-624, (2011)");
 
     AddDocTag(Tags::Filter);
-	AddDocTag(Tags::DimensionReduction);
+    AddDocTag(Tags::DimensionReduction);
 
     AddParameter(ParameterType_InputImage, "in", "Input Image");
     SetParameterDescription("in", "The input image to apply dimensionality reduction.");
@@ -250,20 +250,31 @@ private:
       // PCA Algorithm
       case 0:
         {
-
         otbAppLogINFO("Using the PCA Algorithm ");
         PCAForwardFilterType::Pointer filter = PCAForwardFilterType::New();
         m_ForwardFilter = filter;
         PCAInverseFilterType::Pointer invFilter = PCAInverseFilterType::New();
         m_InverseFilter = invFilter;
 
-
         filter->SetInput(GetParameterFloatVectorImage("in"));
         filter->SetNumberOfPrincipalComponentsRequired(nbComp);
-        filter->SetUseNormalization(normalize);        
+        
+        // Center AND reduce the input data.
+        if (normalize)
+        {
+          filter->SetUseNormalization(true);
+          filter->SetUseVarianceForNormalization(true);
+        }
+        // Only center the input data.
+        else
+        {
+          filter->SetUseNormalization(true);
+          filter->SetUseVarianceForNormalization(false);
+        }
+        
         m_ForwardFilter->GetOutput()->UpdateOutputInformation();
         
-        //Write transformation matrix
+        //Write eigenvalues
         std::ofstream outFile;
         outFile.open(this->GetParameterString("method.pca.outeigenvalues"));
         if (outFile.is_open())
@@ -275,9 +286,8 @@ private:
           outFile.close();
           }
         
-      
         if (invTransform)
-          {  
+          {
           invFilter->SetInput(m_ForwardFilter->GetOutput());
           if (normalize)
             {

Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx

@@ -160,13 +160,37 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
           if ( !m_GivenStdDevValues )
           {
             m_StdDevValues = m_Normalizer->GetFunctor().GetStdDev();
+            
+            
+          }
+          if ( m_UseVarianceForNormalization )
+          {
             // Set User std value so the filter won't recompute the stats
             m_Normalizer->SetStdDev( m_StdDevValues );
+            
+            // Compute the correlation matrix, note that GetCovarianceEstimator()->GetCorrelation()
+            // would give us the matrix with component E[XY], which is not what we want., we want
+            // the matrix defined by its component (E[XY]-E[X]E[Y])/(sigmaX*sigmaY)
+            
+            m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
+            
+            auto cov = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
+            auto numberOfComponent = this->GetInput()->GetNumberOfComponentsPerPixel();
+            
+            for (unsigned int r = 0; r < numberOfComponent; ++r)
+            {
+              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)));
+              }
+            }
           }
-          if ( m_UseVarianceForNormalization )
-            m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCorrelation();
           else
+          {
+            m_Normalizer->SetUseStdDev(false);
             m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
+          }
+        
         }
         else
         {
@@ -362,6 +386,9 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
   vnl_vector< double > valP = solver.D.diagonal();
   valP.flip();
 
+  m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
+  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
+    m_EigenValues[i] = static_cast< RealType >( valP[i] );
   /*
    * We used normalized PCA
    */
@@ -382,9 +409,6 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
     m_TransformationMatrix = transf;
 
 
-  m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
-  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
-    m_EigenValues[i] = static_cast< RealType >( valP[i] );
 #else
   InternalMatrixType transf;
   vnl_vector<double> vectValP;
@@ -396,8 +420,8 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
 
   m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
   for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i )
-    m_EigenValues[i] = static_cast< RealType >( valP(i, i) );
-
+    m_EigenValues[m_NumberOfPrincipalComponentsRequired - 1 - i] = static_cast< RealType >( valP(i, i) );
+  
   /* We used normalized PCA */
   for ( unsigned int i = 0; i < valP.rows(); ++i )
   {