diff --git a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx b/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx index 24f5951f6bf8afb75168f47686f24d995e367f45..3e5cbf8e7a4ad2ecbac6182a3b42cf13e56c9141 100644 --- a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx +++ b/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx @@ -130,6 +130,10 @@ private: AddParameter(ParameterType_OutputFilename, "method.pca.outeigenvalues", "Output file containing eigenvalues (txt format)"); SetParameterDescription("method.pca.outeigenvalues", "Output file containing eigenvalues (txt format)."); MandatoryOff("method.pca.outeigenvalues"); + AddParameter(ParameterType_Bool, "method.pca.whiten", "Perform pca whitening"); + SetParameterDescription("method.pca.whiten", "Perform whitening and ensure uncorrelated outputs with unit component wise variances"); + SetParameterInt("method.pca.whiten", 1); + MandatoryOff("method.pca.whiten"); AddChoice("method.napca", "NA-PCA"); SetParameterDescription("method.napca", "Noise Adjusted Principal Component Analysis."); @@ -257,6 +261,7 @@ private: filter->SetInput(GetParameterFloatVectorImage("in")); filter->SetNumberOfPrincipalComponentsRequired(nbComp); + filter->SetWhitening(GetParameterInt("method.pca.whiten")); // Center AND reduce the input data. if (normalize) diff --git a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h index 3368163846bb09339b0fe29eeed9981043f66ac5..b90faae53679d0d806cc3886d10d72bdb465e762 100644 --- a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h +++ b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h @@ -97,6 +97,9 @@ public: */ itkSetMacro(NumberOfPrincipalComponentsRequired, unsigned int); itkGetMacro(NumberOfPrincipalComponentsRequired, unsigned int); + + itkSetMacro(Whitening, bool); + itkGetMacro(Whitening, bool); itkGetMacro(CovarianceEstimator, CovarianceEstimatorFilterType *); itkGetMacro(Transformer, TransformFilterType *); @@ -197,6 +200,7 @@ protected: bool m_GivenCovarianceMatrix; bool m_GivenTransformationMatrix; bool m_IsTransformationMatrixForward; + bool m_Whitening; VectorType m_MeanValues; VectorType m_StdDevValues; diff --git a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx index 3c2487119ce13d2c04d66a2866de545ee67a6d6b..51633d1eafda79562947a590da19baaf9a33b7e3 100644 --- a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx +++ b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx @@ -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; @@ -181,7 +181,7 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation > { 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))); + m_CovarianceMatrix(r,c) = cov(r,c) / std::sqrt(cov(r,r)*cov(c,c)); } } } @@ -376,19 +376,23 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation > m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired ); for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; ++i ) - m_EigenValues[m_NumberOfPrincipalComponentsRequired - 1 - i] = static_cast< RealType >( valP(i, 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 {