From f4a1b37f665e981b2a491adf48a5520b47f6b360 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=A9dric=20Traizet?= <cedric.traizet@c-s.fr>
Date: Tue, 3 Sep 2019 17:16:51 +0200
Subject: [PATCH] BUG: center the data by default and use the correct
correlation matrix when the data is reduced
---
.../app/otbDimensionalityReduction.cxx | 24 ++++++++----
.../include/otbPCAImageFilter.hxx | 38 +++++++++++++++----
2 files changed, 48 insertions(+), 14 deletions(-)
diff --git a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx b/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
index c287154cab..59f31167f6 100644
--- a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
+++ b/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)
{
diff --git a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx
index a0ad9ddb75..24a31cc6b7 100644
--- a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx
+++ b/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 )
{
--
GitLab