Commit 6e96b931 authored by Cédric Traizet's avatar Cédric Traizet

Merge branch 'fix_dimensionality_reduction' into 'develop'

Fix PCA in Dimensionality reduction

See merge request !585
parents e9c43047 26a78e67
Pipeline #2568 passed with stages
in 6 minutes and 46 seconds
......@@ -35,7 +35,7 @@ otb_test_application(NAME apTvFEDimensionalityReductionPCA
-out ${TEMP}/apTvChDimensionalityReductionPCA.tif
-method pca
VALID --compare-image 0.025
${BASELINE}/bfTvPCAImageFilter3.tif
${BASELINE}/apTvChDimensionalityReductionPCA.tif
${TEMP}/apTvChDimensionalityReductionPCA.tif)
#-------------------------------------------------------------------------------
......
......@@ -145,6 +145,19 @@ public:
m_IsTransformationMatrixForward = isForward;
}
void SetStatisticsUserIgnoredValue ( RealType value )
{
/** User ignored value for the normalizer */
m_Normalizer->GetCovarianceEstimator()->SetUserIgnoredValue(value);
m_Normalizer->GetCovarianceEstimator()->SetIgnoreUserDefinedValue(true);
/** User ignored value for the covariance estimator */
m_CovarianceEstimator->SetUserIgnoredValue(value);
m_CovarianceEstimator->SetIgnoreUserDefinedValue(true);
/** User ignored value for the noise covariance estimator */
m_NoiseCovarianceEstimator->SetUserIgnoredValue(value);
m_NoiseCovarianceEstimator->SetIgnoreUserDefinedValue(true);
}
itkGetConstMacro(EigenValues, VectorType);
protected:
......
......@@ -376,54 +376,6 @@ MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransf
InternalMatrixType transf = Rn_inv * U;
#if 0
InternalMatrixType Ax = m_CovarianceMatrix.GetVnlMatrix();
InternalMatrixType Ax_inv = vnl_matrix_inverse< MatrixElementType > ( Ax );
InternalMatrixType An = m_NoiseCovarianceMatrix.GetVnlMatrix();
InternalMatrixType W = An * Ax_inv;
#endif
#if 0
InternalMatrixType transf;
vnl_vector< MatrixElementType > vectValP;
vnl_symmetric_eigensystem_compute( W, 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];
#endif
#if 0
vnl_svd< MatrixElementType > solver ( W );
InternalMatrixType transf = solver.U();
InternalMatrixType valP = solver.W();
#endif
#if 0
MatrixType Id ( m_CovarianceMatrix );
Id.SetIdentity();
vnl_generalized_eigensystem solver ( W, Id.GetVnlMatrix() );
InternalMatrixType transf = solver.V;
transf.fliplr();
InternalMatrixType valP = solver.D;
valP.fliplr();
valP.flipud();
#endif
/*
InternalMatrixType normMat = transf.transpose() * Ax * transf;
for ( unsigned int i = 0; i < transf.rows(); ++i )
{
double norm = 1. / std::sqrt( normMat.get(i, i) );
for ( unsigned int j = 0; j < transf.cols(); ++j )
transf.put( j, i, transf.get(j, i) * norm );
}
*/
transf.inplace_transpose();
if ( m_NumberOfPrincipalComponentsRequired
......
......@@ -98,6 +98,9 @@ public:
itkSetMacro(NumberOfPrincipalComponentsRequired, unsigned int);
itkGetMacro(NumberOfPrincipalComponentsRequired, unsigned int);
itkSetMacro(Whitening, bool);
itkGetMacro(Whitening, bool);
itkGetMacro(CovarianceEstimator, CovarianceEstimatorFilterType *);
itkGetMacro(Transformer, TransformFilterType *);
......@@ -158,7 +161,17 @@ public:
m_StdDevValues = vec;
this->Modified();
}
void SetStatisticsUserIgnoredValue ( RealType value )
{
/** User ignored value for the normalizer */
m_Normalizer->GetCovarianceEstimator()->SetUserIgnoredValue(value);
m_Normalizer->GetCovarianceEstimator()->SetIgnoreUserDefinedValue(true);
/** User ignored value for the covariance estimator */
m_CovarianceEstimator->SetUserIgnoredValue(value);
m_CovarianceEstimator->SetIgnoreUserDefinedValue(true);
}
protected:
PCAImageFilter();
~PCAImageFilter() override { }
......@@ -197,6 +210,7 @@ protected:
bool m_GivenCovarianceMatrix;
bool m_GivenTransformationMatrix;
bool m_IsTransformationMatrixForward;
bool m_Whitening;
VectorType m_MeanValues;
VectorType m_StdDevValues;
......
......@@ -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;
......@@ -154,14 +154,43 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
if ( !m_GivenMeanValues )
{
m_MeanValues = m_Normalizer->GetCovarianceEstimator()->GetMean();
// Set User mean value so the filter won't recompute the stats
m_Normalizer->SetMean( m_MeanValues );
if ( !m_GivenStdDevValues )
{
m_StdDevValues = m_Normalizer->GetFunctor().GetStdDev();
}
if ( m_UseVarianceForNormalization )
m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCorrelation();
{
// 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)*cov(c,c));
}
}
}
else
{
m_Normalizer->SetUseStdDev(false);
m_CovarianceMatrix = m_Normalizer->GetCovarianceEstimator()->GetCovariance();
}
}
else
{
......@@ -337,50 +366,6 @@ void
PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
::GenerateTransformationMatrix ()
{
#if 0
/*
* Old stuff but did work !
*/
MatrixType Id ( m_CovarianceMatrix );
Id.SetIdentity();
typename MatrixType::InternalMatrixType A = m_CovarianceMatrix.GetVnlMatrix();
typename MatrixType::InternalMatrixType I = Id.GetVnlMatrix();
vnl_generalized_eigensystem solver ( A, I );
typename MatrixType::InternalMatrixType transf = solver.V;
transf.fliplr();
transf.inplace_transpose();
vnl_vector< double > valP = solver.D.diagonal();
valP.flip();
/*
* We used normalized PCA
*/
for ( unsigned int i = 0; i < valP.size(); ++i )
{
if ( valP[i] != 0. )
valP[i] = 1. / std::sqrt( std::abs( valP[i] ) );
else
throw itk::ExceptionObject( __FILE__, __LINE__,
"Null Eigen value !!", ITK_LOCATION );
}
valP.post_multiply( transf );
if ( m_NumberOfPrincipalComponentsRequired
!= this->GetInput()->GetNumberOfComponentsPerPixel() )
m_TransformationMatrix = transf.get_n_rows( 0, m_NumberOfPrincipalComponentsRequired );
else
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;
vnl_symmetric_eigensystem_compute( m_CovarianceMatrix.GetVnlMatrix(), transf, vectValP );
......@@ -389,17 +374,25 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
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 ( 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
{
......@@ -415,12 +408,6 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
m_TransformationMatrix = transf.get_n_rows( 0, m_NumberOfPrincipalComponentsRequired );
else
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, i) );
#endif
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment