Merge branch 'fix_pca' into 'develop'

FIX PCA transformation matrix computation

See merge request orfeotoolbox/otb!641

Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx

@@ -325,36 +325,30 @@ void PCAImageFilter<TInputImage, TOutputImage, TDirectionOfTransformation>::Gene
   vnl_vector<double> vectValP;
   vnl_symmetric_eigensystem_compute(m_CovarianceMatrix.GetVnlMatrix(), transf, vectValP);
 
-  InternalMatrixType valP(vectValP.size(), vectValP.size(), vnl_matrix_null);
-  for (unsigned int i = 0; i < vectValP.size(); ++i)
-    valP(i, i) = vectValP[i];
 
   m_EigenValues.SetSize(m_NumberOfPrincipalComponentsRequired);
   for (unsigned int i                                            = 0; i < m_NumberOfPrincipalComponentsRequired; ++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 (m_Whitening)
   {
-    if (valP(i, i) > 0.)
-    {
-      if (m_Whitening)
-        valP(i, i) = 1. / std::sqrt(valP(i, i));
-    }
-    else if (valP(i, i) < 0.)
+    InternalMatrixType valP(vectValP.size(), vectValP.size(), vnl_matrix_null);
+    for (unsigned int i = 0; i < vectValP.size(); ++i)
+    valP(i, i) = vectValP[i];
+
+    for (unsigned int i = 0; i < valP.rows(); ++i)
     {
-      otbMsgDebugMacro(<< "ValP(" << i << ") neg : " << valP(i, i) << " taking abs value");
-      if (m_Whitening)
-        valP(i, i) = 1. / std::sqrt(std::abs(valP(i, i)));
+      if (valP(i,i) != 0.0)
+        valP(i,i) = 1.0 / std::sqrt(std::abs(valP(i,i)));
       else
-        valP(i, i) = std::abs(valP(i, i));
-    }
-    else
-    {
-      throw itk::ExceptionObject(__FILE__, __LINE__, "Null Eigen value !!", ITK_LOCATION);
+        throw itk::ExceptionObject(__FILE__, __LINE__, "Null Eigen value !!", ITK_LOCATION);
     }
+    transf = valP * transf.transpose();
   }
-  transf = valP * transf.transpose();
+  else {
+      transf = transf.transpose();
+  }
+
   transf.flipud();
 
   if (m_NumberOfPrincipalComponentsRequired != this->GetInput()->GetNumberOfComponentsPerPixel())

Modules/Filtering/DimensionalityReduction/test/CMakeLists.txt

@@ -301,7 +301,8 @@ otb_add_test(NAME bfTvPCAImageFilter2 COMMAND otbDimensionalityReductionTestDriv
   ${TEMP}/bfTvPCAImageFilter2.tif
   ${TEMP}/bfTvPCAImageFilter2Inv.tif
   true
-  0)
+  0
+  true)
 
 otb_add_test(NAME bfTvPCAImageFilter4 COMMAND otbDimensionalityReductionTestDriver
   --compare-n-images ${EPSILON_7} 2
@@ -314,7 +315,8 @@ otb_add_test(NAME bfTvPCAImageFilter4 COMMAND otbDimensionalityReductionTestDriv
   ${TEMP}/bfTvPCAImageFilter4.tif
   ${TEMP}/bfTvPCAImageFilter4Inv.tif
   false
-  4)
+  4
+  true)
 
 otb_add_test(NAME bfTvPCAImageFilter4Norm COMMAND otbDimensionalityReductionTestDriver
   --compare-n-images ${EPSILON_7} 2
@@ -327,7 +329,8 @@ otb_add_test(NAME bfTvPCAImageFilter4Norm COMMAND otbDimensionalityReductionTest
   ${TEMP}/bfTvPCAImageFilter4Norm.tif
   ${TEMP}/bfTvPCAImageFilter4InvNorm.tif
   true
-  4)
+  4
+  true)
 
 otb_add_test(NAME bfTvPCAImageFilter3 COMMAND otbDimensionalityReductionTestDriver
   --compare-n-images ${EPSILON_7} 2
@@ -340,7 +343,22 @@ otb_add_test(NAME bfTvPCAImageFilter3 COMMAND otbDimensionalityReductionTestDriv
   ${TEMP}/bfTvPCAImageFilter3.tif
   ${TEMP}/bfTvPCAImageFilter3Inv.tif
   false
-  0)
+  0
+  true)
+
+otb_add_test(NAME bfTvPCAImageFilter5 COMMAND otbDimensionalityReductionTestDriver
+  --compare-n-images ${EPSILON_7} 2
+  ${BASELINE}/bfTvPCAImageFilter5.tif
+  ${TEMP}/bfTvPCAImageFilter5.tif
+  ${BASELINE}/bfTvPCAImageFilter5Inv.tif
+  ${TEMP}/bfTvPCAImageFilter5Inv.tif
+  otbPCAImageFilterTest
+  ${INPUTDATA}/cupriteSubHsi.tif
+  ${TEMP}/bfTvPCAImageFilter5.tif
+  ${TEMP}/bfTvPCAImageFilter5Inv.tif
+  false
+  0
+  false)
 
 otb_add_test(NAME bfTvPCAImageFilter1 COMMAND otbDimensionalityReductionTestDriver
   --compare-image ${EPSILON_7}
@@ -351,7 +369,8 @@ otb_add_test(NAME bfTvPCAImageFilter1 COMMAND otbDimensionalityReductionTestDriv
   ${TEMP}/bfTvPCAImageFilter1.tif
   ${TEMP}/bfTvPCAImageFilter1Inv.tif
   false
-  0)
+  0
+  true)
 
 
 otb_add_test(NAME bfTvPCAImageFilter3Norm COMMAND otbDimensionalityReductionTestDriver
@@ -365,7 +384,8 @@ otb_add_test(NAME bfTvPCAImageFilter3Norm COMMAND otbDimensionalityReductionTest
   ${TEMP}/bfTvPCAImageFilter3Norm.tif
   ${TEMP}/bfTvPCAImageFilter3InvNorm.tif
   true
-  0)
+  0
+  true)
 
 otb_add_test(NAME bfTvNAPCAImageFilter1 COMMAND otbDimensionalityReductionTestDriver
   --compare-image ${EPSILON_7}

Modules/Filtering/DimensionalityReduction/test/otbPCAImageFilter.cxx

@@ -33,6 +33,10 @@ int otbPCAImageFilterTest(int, char* argv[])
   bool normalization = false;
   if (std::string(argv[4]).compare("true") == 0)
     normalization = true;
+  
+  bool whitening = false;
+  if (std::string(argv[6]).compare("true") == 0)
+    whitening = true;
 
   // Main type definition
   const unsigned int Dimension = 2;
@@ -50,6 +54,7 @@ int otbPCAImageFilterTest(int, char* argv[])
   filter->SetInput(reader->GetOutput());
   filter->SetNumberOfPrincipalComponentsRequired(nbComponents);
   filter->SetUseNormalization(normalization);
+  filter->SetWhitening(whitening);
 
   typedef otb::CommandProgressUpdate<FilterType> CommandType;
   CommandType::Pointer                           observer = CommandType::New();