ENL: Use mean and variance in the PCA normalization

Code/Hyperspectral/otbPCAImageFilter.h

@@ -132,6 +132,18 @@ public:
     this->Modified();
   }
 
+  itkGetConstMacro(UseVarianceForNormalization,bool);
+  itkSetMacro(UseVarianceForNormalization,bool);
+
+  itkGetConstMacro(StdDevValues,VectorType);
+  void SetStdDevValues ( const VectorType & vec )
+  {
+    m_GivenStdDevValues = true;
+    m_UseVarianceForNormalization = true;
+    m_StdDevValues = vec;
+    this->Modified();
+  }
+
 protected:
   PCAImageFilter();
   virtual ~PCAImageFilter() { }
@@ -162,11 +174,14 @@ protected:
   /** Internal attributes */
   unsigned int m_NumberOfPrincipalComponentsRequired;
   bool m_GivenMeanValues;
+  bool m_GivenStdDevValues;
   bool m_GivenCovarianceMatrix;
   bool m_GivenTransformationMatrix;
   bool m_IsTransformationMatrixForward;
+  bool m_UseVarianceForNormalization;
 
   VectorType m_MeanValues;
+  VectorType m_StdDevValues;
   MatrixType m_CovarianceMatrix;
   VectorType m_EigenValues;
   MatrixType m_TransformationMatrix;

Code/Hyperspectral/otbPCAImageFilter.txx

@@ -36,9 +36,11 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
 
   m_NumberOfPrincipalComponentsRequired = 0;
   m_GivenMeanValues = false;
+  m_GivenStdDevValues = false;
   m_GivenCovarianceMatrix = false;
   m_GivenTransformationMatrix = false;
   m_IsTransformationMatrixForward = true;
+  m_UseVarianceForNormalization = false;
 
   m_CovarianceEstimator = CovarianceEstimatorFilterType::New();
   m_Transformer = TransformFilterType::New();
@@ -107,6 +109,7 @@ void
 PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
 ::GenerateData ()
 {
+  std::cerr << __PRETTY_FUNCTION__ << "\n";
   switch ( DirectionOfTransformation )
   {
     case Transform::FORWARD:
@@ -134,17 +137,27 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
     if ( !m_GivenCovarianceMatrix )
     {
       m_Normalizer->SetInput( inputImgPtr );
-      m_Normalizer->SetUseStdDev( false );
+      m_Normalizer->SetUseStdDev( m_UseVarianceForNormalization );
 
       if ( m_GivenMeanValues )
         m_Normalizer->SetMean( m_MeanValues );
+       
+      if ( m_GivenStdDevValues )
+        m_Normalizer->SetStdDev( m_StdDevValues );
       
       m_Normalizer->Update();
 
       if ( !m_GivenMeanValues )
       {
         m_MeanValues = m_Normalizer->GetCovarianceEstimator()->GetMean();
-        m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
+
+        if ( !m_GivenStdDevValues )
+          m_StdDevValues = m_Normalizer->GetFunctor().GetStdDev();
+
+        if ( m_UseVarianceForNormalization )
+          m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCorrelation();
+        else
+          m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
       }
       else
       {
@@ -203,7 +216,8 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
   }
   else if ( m_IsTransformationMatrixForward )
   {
-    std::cerr << "Transposition\n";
+    // prevents from multiple transpositions...
+    m_IsTransformationMatrixForward = false;
     m_TransformationMatrix = m_TransformationMatrix.GetTranspose();
   }
 
@@ -217,15 +231,35 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
   m_Transformer->SetInput( this->GetInput() );
   m_Transformer->SetMatrix( m_TransformationMatrix.GetVnlMatrix() );
 
-  if ( m_GivenMeanValues )
+  if ( m_GivenMeanValues || m_GivenStdDevValues )
   {
-    VectorType revMean ( m_MeanValues.Size() );
-    for ( unsigned int i = 0; i < m_MeanValues.Size(); i++ )
-      revMean[i] = -m_MeanValues[i];
-
     m_Normalizer->SetInput( m_Transformer->GetOutput() );
-    m_Normalizer->SetMean( revMean );
-    m_Normalizer->SetUseStdDev( false );
+
+    if ( m_GivenStdDevValues )
+    {
+      VectorType revStdDev ( m_StdDevValues.Size() );
+      for ( unsigned int i = 0; i < m_StdDevValues.Size(); i++ )
+        revStdDev[i] = 1. / m_StdDevValues[i];
+      m_Normalizer->SetStdDev( revStdDev );
+    }
+
+    if ( m_GivenMeanValues )
+    {
+      VectorType revMean ( m_MeanValues.Size() );
+      if ( m_GivenStdDevValues )
+      {
+        for ( unsigned int i = 0; i < m_MeanValues.Size(); i++ )
+          revMean[i] = - m_MeanValues[i] / m_StdDevValues[i];
+        m_Normalizer->SetUseStdDev( true );
+      }
+      else
+      {
+        for ( unsigned int i = 0; i < m_MeanValues.Size(); i++ )
+          revMean[i] = -m_MeanValues[i];
+        m_Normalizer->SetUseStdDev( false );
+      }
+      m_Normalizer->SetMean( revMean );
+    }
 
     m_Normalizer->GraftOutput( this->GetOutput() );
     m_Normalizer->Update();

Testing/otbPCAImageFilter.cxx

@@ -47,6 +47,7 @@ int otbPCAImageFilterTest ( int argc, char* argv[] )
   parser->AddInputImage();
   parser->AddOption( "--NumComponents", "Number of components to keep for output", "-n", 1, false );
   parser->AddOption( "--Inverse", "Performs also the inverse transformation (give the output name)", "-inv", 1, false );
+  parser->AddOption( "--NormalizeVariance", "center AND reduce data before PCA", "-norm", 0, false );
   parser->AddOutputImage();
 
   typedef otb::CommandLineArgumentParseResult ParserResultType;  
@@ -73,6 +74,7 @@ int otbPCAImageFilterTest ( int argc, char* argv[] )
   const char * outputImageName = parseResult->GetOutputImage().c_str();
   const unsigned int nbComponents = parseResult->IsOptionPresent("--NumComponents") ?
     parseResult->GetParameterUInt("--NumComponents") : 0;
+  const bool normalization = parseResult->IsOptionPresent("--NormalizeVariance");
 
   // Main type definition
   const unsigned int Dimension = 2;
@@ -89,6 +91,7 @@ int otbPCAImageFilterTest ( int argc, char* argv[] )
   FilterType::Pointer filter = FilterType::New();
   filter->SetInput( reader->GetOutput() );
   filter->SetNumberOfPrincipalComponentsRequired( nbComponents );
+  filter->SetUseVarianceForNormalization( normalization );
 
   typedef otb::CommandProgressUpdate< FilterType > CommandType;
   CommandType::Pointer observer = CommandType::New();
@@ -108,6 +111,8 @@ int otbPCAImageFilterTest ( int argc, char* argv[] )
     InvFilterType::Pointer invFilter = InvFilterType::New();
     invFilter->SetInput( filter->GetOutput() );
     invFilter->SetMeanValues( filter->GetMeanValues() );
+    if ( normalization )
+      invFilter->SetStdDevValues( filter->GetStdDevValues() );
     invFilter->SetTransformationMatrix( filter->GetTransformationMatrix() );
 
     typedef otb::CommandProgressUpdate< InvFilterType > CommandType2;