diff --git a/Data/Baseline/OTB/Images/apTvChDimensionalityReductionPCA.tif b/Data/Baseline/OTB/Images/apTvChDimensionalityReductionPCA.tif
new file mode 100644
index 0000000000000000000000000000000000000000..c694726c2717746b6b9dbc678830a5f2723bcd0a
--- /dev/null
+++ b/Data/Baseline/OTB/Images/apTvChDimensionalityReductionPCA.tif
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c12364434ffad8319d964d79a467f6eee094c2472a72234bb6b86a9673bec437
+size 180416
diff --git a/Data/Baseline/OTB/Images/bfTvPCAImageFilter3Norm.tif b/Data/Baseline/OTB/Images/bfTvPCAImageFilter3Norm.tif
index 466bd1d886ae000064b0a0208177d079c97f9cf1..2d918a952e15a63f4bf4e6b09e38e533473f8ef9 100644
--- a/Data/Baseline/OTB/Images/bfTvPCAImageFilter3Norm.tif
+++ b/Data/Baseline/OTB/Images/bfTvPCAImageFilter3Norm.tif
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:1cd38ed0a735e6e7e332237b02982d955a44759565c5892381e508a9859f1628
-size 450331
+oid sha256:7f42429c4f3375ba587c64078a8ebfd08f8b864d104fd14efb8bd9285f074b06
+size 360816
diff --git a/Data/Baseline/OTB/Images/bfTvPCAImageFilter4InvNorm.tif b/Data/Baseline/OTB/Images/bfTvPCAImageFilter4InvNorm.tif
index f01499d2690a5498128e939f3b2b499c5c97f10c..884d6b1670fb03a5d108316d832b9af17a3822d2 100644
--- a/Data/Baseline/OTB/Images/bfTvPCAImageFilter4InvNorm.tif
+++ b/Data/Baseline/OTB/Images/bfTvPCAImageFilter4InvNorm.tif
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:3f8a7ac9f4e604573cb693afe0208adfd09af409be9446c42626453574de6993
-size 429742
+oid sha256:5bfcdfb81cff472ceb447e39fc5a6718224805de34106ecb477535fec8506b84
+size 360816
diff --git a/Data/Baseline/OTB/Images/bfTvPCAImageFilter4Norm.tif b/Data/Baseline/OTB/Images/bfTvPCAImageFilter4Norm.tif
index 253afa7ec9d0b3fda9468201a9d08eff61c6430d..6b36f055b2d5745b5e6056b34a052d85f6503dbe 100644
--- a/Data/Baseline/OTB/Images/bfTvPCAImageFilter4Norm.tif
+++ b/Data/Baseline/OTB/Images/bfTvPCAImageFilter4Norm.tif
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:b071be1eae39fa5a77f9c0851cfa5b02e459e29b351a13da995db7174c5d4412
-size 181939
+oid sha256:538cd26eb6c8a6da6fff2799f090e6d5ef28cad9aca0d0aa813c35c33507b24d
+size 144332
diff --git a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx b/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
index cecb7e9045780d5011ed840afb4a4b679469c352..7ceb4f578efe22251d448e753e602430eed83d5f 100644
--- a/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
+++ b/Modules/Applications/AppDimensionalityReduction/app/otbDimensionalityReduction.cxx
@@ -90,31 +90,33 @@ private:
SetName("DimensionalityReduction");
SetDescription("Perform Dimension reduction of the input image.");
SetDocLongDescription("Performs dimensionality reduction on input image. PCA,NA-PCA,MAF,ICA methods are available. It is also possible to compute the inverse transform to reconstruct the image and to optionally export the transformation matrix to a text file.");
- SetDocLimitations("This application does not provide the inverse transform and the transformation matrix export for the MAF.");
+ SetDocLimitations("This application does not provide the inverse transform and the transformation matrix export for the MAF. the background value option is not supported for MAF and ICA.");
SetDocAuthors("OTB-Team");
SetDocSeeAlso(
"\"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.");
AddParameter(ParameterType_OutputImage, "out", "Output Image");
SetParameterDescription("out", "output image. Components are ordered by decreasing eigenvalues.");
MandatoryOff("out");
- AddParameter(ParameterType_Group, "rescale", "Rescale Output");
+
+ AddParameter(ParameterType_Choice, "rescale", "Rescale Output");
+ SetParameterDescription("rescale", "Enable rescaling of the reduced output image.");
MandatoryOff("rescale");
- // AddChoice("rescale.no","No rescale");
- // AddChoice("rescale.minmax","rescale to min max value");
+ AddChoice("rescale.no", "No rescale");
+ AddChoice("rescale.minmax", "rescale to min max value");
- AddParameter(ParameterType_Float, "rescale.outmin", "Output min value");
- AddParameter(ParameterType_Float, "rescale.outmax", "Output max value");
- SetDefaultParameterFloat("rescale.outmin", 0.0);
- SetParameterDescription("rescale.outmin", "Minimum value of the output image.");
- SetDefaultParameterFloat("rescale.outmax", 255.0);
- SetParameterDescription("rescale.outmax", "Maximum value of the output image.");
+ AddParameter(ParameterType_Float, "rescale.minmax.outmin", "Output min value");
+ AddParameter(ParameterType_Float, "rescale.minmax.outmax", "Output max value");
+ SetDefaultParameterFloat("rescale.minmax.outmin", 0.0);
+ SetParameterDescription("rescale.minmax.outmin", "Minimum value of the output image.");
+ SetDefaultParameterFloat("rescale.minmax.outmax", 255.0);
+ SetParameterDescription("rescale.minmax.outmax", "Maximum value of the output image.");
AddParameter(ParameterType_OutputImage, "outinv", " Inverse Output Image");
SetParameterDescription("outinv", "reconstruct output image.");
@@ -125,6 +127,14 @@ private:
AddChoice("method.pca", "PCA");
SetParameterDescription("method.pca", "Principal Component Analysis.");
+ 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.");
AddParameter(ParameterType_Int, "method.napca.radiusx", "Set the x radius of the sliding window");
@@ -136,6 +146,7 @@ private:
AddChoice("method.maf", "MAF");
SetParameterDescription("method.maf", "Maximum Autocorrelation Factor.");
+
AddChoice("method.ica", "ICA");
SetParameterDescription("method.ica", "Independent Component Analysis using a stabilized fixed point FastICA algorithm.");
AddParameter(ParameterType_Int, "method.ica.iter", "number of iterations");
@@ -164,14 +175,17 @@ private:
MandatoryOff("nbcomp");
SetMinimumParameterIntValue("nbcomp", 0);
- AddParameter(ParameterType_Bool, "normalize", "Normalize");
- SetParameterDescription("normalize", "Center and reduce data before Dimensionality reduction.");
+ AddParameter(ParameterType_Bool, "normalize", "Center and reduce data");
+ SetParameterDescription("normalize", "Center and reduce data before Dimensionality reduction (if this parameter is set to false, the data will be centered but not reduced.");
AddParameter(ParameterType_OutputFilename, "outmatrix", "Transformation matrix output (text format)");
SetParameterDescription("outmatrix", "Filename to store the transformation matrix (csv format)");
MandatoryOff("outmatrix");
DisableParameter("outmatrix");
+ AddParameter(ParameterType_Float, "bv", "Background Value");
+ SetParameterDescription( "bv", "Background value to ignore in computation of the transformation matrix. Note that all pixels will still be processed when applying the transformation." );
+ MandatoryOff("bv");
AddRAMParameter();
@@ -237,74 +251,111 @@ private:
// Get Parameters
int nbComp = GetParameterInt("nbcomp");
bool normalize = GetParameterInt("normalize");
- bool rescale = IsParameterEnabled("rescale");
-
bool invTransform = HasValue("outinv") && IsParameterEnabled("outinv");
switch (GetParameterInt("method"))
{
// 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;
- otbAppLogDEBUG( << "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->SetWhitening(GetParameterInt("method.pca.whiten"));
- filter->SetInput(GetParameterFloatVectorImage("in"));
- filter->SetNumberOfPrincipalComponentsRequired(nbComp);
- filter->SetUseNormalization(normalize);
- m_ForwardFilter->GetOutput()->UpdateOutputInformation();
-
- if (invTransform)
- {
- invFilter->SetInput(m_ForwardFilter->GetOutput());
+ // 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);
+ }
+
+ if( HasValue( "bv" ) )
+ {
+ filter->SetStatisticsUserIgnoredValue(GetParameterFloat("bv"));
+ }
+
+ m_ForwardFilter->GetOutput()->UpdateOutputInformation();
+
+ // Write eigenvalues
+ std::ofstream outFile;
+ outFile.open(this->GetParameterString("method.pca.outeigenvalues"));
+ if (outFile.is_open())
+ {
+ outFile << std::fixed;
+ outFile.precision(10);
+
+ outFile << filter->GetEigenValues() << std::endl;
+ ;
+ outFile.close();
+ }
+
+ if (invTransform)
+ {
+ invFilter->SetInput(m_ForwardFilter->GetOutput());
+ // Data has been centered and reduced by the forward filter
+ if (normalize)
{
otbAppLogINFO( << "Normalization MeanValue:"<<filter->GetMeanValues() );
invFilter->SetMeanValues(filter->GetMeanValues());
- // By default normalization by std dev is deactivated in
- //forward filter, and GetStdDevValues() returns an empty
- //vector, which confuses the invFilter.
- //invFilter->SetStdDevValues(filter->GetStdDevValues());
+ invFilter->SetStdDevValues(filter->GetStdDevValues());
+ }
+ // Data has been centered by the forward filter
+ else
+ {
+ invFilter->SetMeanValues(filter->GetMeanValues());
}
- invFilter->SetTransformationMatrix(filter->GetTransformationMatrix());
- m_TransformationMatrix = invFilter->GetTransformationMatrix();
+ invFilter->SetTransformationMatrix(filter->GetTransformationMatrix());
+ m_TransformationMatrix = invFilter->GetTransformationMatrix();
}
m_TransformationMatrix = filter->GetTransformationMatrix();
+ otbAppLogINFO("PCA transform has been computed.");
break;
}
case 1:
{
- otbAppLogDEBUG( << "NA-PCA Algorithm ");
+ otbAppLogINFO("Using the NA-PCA Algorithm ");
- // NA-PCA
+ // NA-PCA
- unsigned int radiusX = static_cast<unsigned int> (GetParameterInt("method.napca.radiusx"));
- unsigned int radiusY = static_cast<unsigned int> (GetParameterInt("method.napca.radiusy"));
+ unsigned int radiusX = static_cast<unsigned int>(GetParameterInt("method.napca.radiusx"));
+ unsigned int radiusY = static_cast<unsigned int>(GetParameterInt("method.napca.radiusy"));
- // Noise filtering
- NoiseFilterType::RadiusType radius = { { radiusX, radiusY } };
+ // Noise filtering
+ NoiseFilterType::RadiusType radius = {{radiusX, radiusY}};
- NAPCAForwardFilterType::Pointer filter = NAPCAForwardFilterType::New();
- m_ForwardFilter = filter;
- NAPCAInverseFilterType::Pointer invFilter = NAPCAInverseFilterType::New();
- m_InverseFilter = invFilter;
+ NAPCAForwardFilterType::Pointer filter = NAPCAForwardFilterType::New();
+ m_ForwardFilter = filter;
+ NAPCAInverseFilterType::Pointer invFilter = NAPCAInverseFilterType::New();
+ m_InverseFilter = invFilter;
- filter->SetInput(GetParameterFloatVectorImage("in"));
- filter->SetNumberOfPrincipalComponentsRequired(nbComp);
- filter->SetUseNormalization(normalize);
- filter->GetNoiseImageFilter()->SetRadius(radius);
+ filter->SetInput(GetParameterFloatVectorImage("in"));
+ filter->SetNumberOfPrincipalComponentsRequired(nbComp);
+ filter->SetUseNormalization(normalize);
+ filter->GetNoiseImageFilter()->SetRadius(radius);
- m_ForwardFilter->GetOutput()->UpdateOutputInformation();
-
- if (invTransform)
+ if( HasValue( "bv" ) )
+ {
+ filter->SetStatisticsUserIgnoredValue(GetParameterFloat("bv"));
+ }
+
+ m_ForwardFilter->GetOutput()->UpdateOutputInformation();
+
+ if (invTransform)
{
otbAppLogDEBUG( << "Compute Inverse Transform");
invFilter->SetInput(m_ForwardFilter->GetOutput());
@@ -322,37 +373,37 @@ private:
}
m_TransformationMatrix = filter->GetTransformationMatrix();
-
+ otbAppLogINFO("NA-PCA transform has been computed.");
break;
}
case 2:
{
- otbAppLogDEBUG( << "MAF Algorithm ");
- MAFForwardFilterType::Pointer filter = MAFForwardFilterType::New();
- m_ForwardFilter = filter;
- filter->SetInput(GetParameterFloatVectorImage("in"));
- otbAppLogINFO( << "V :"<<std::endl<<filter->GetV()<<"Auto-Correlation :"<<std::endl <<filter->GetAutoCorrelation() );
+ otbAppLogINFO("Using the MAF Algorithm ");
+ MAFForwardFilterType::Pointer filter = MAFForwardFilterType::New();
+ m_ForwardFilter = filter;
+ filter->SetInput(GetParameterFloatVectorImage("in"));
+ otbAppLogINFO(<< "V :" << std::endl << filter->GetV() << "Auto-Correlation :" << std::endl << filter->GetAutoCorrelation());
- break;
+ break;
}
case 3:
{
- otbAppLogDEBUG( << "Fast ICA Algorithm ");
+ otbAppLogINFO("Using the fast ICA Algorithm ");
- unsigned int nbIterations = static_cast<unsigned int> (GetParameterInt("method.ica.iter"));
- double mu = static_cast<double> (GetParameterFloat("method.ica.mu"));
+ unsigned int nbIterations = static_cast<unsigned int>(GetParameterInt("method.ica.iter"));
+ double mu = static_cast<double>(GetParameterFloat("method.ica.mu"));
- ICAForwardFilterType::Pointer filter = ICAForwardFilterType::New();
- m_ForwardFilter = filter;
- ICAInverseFilterType::Pointer invFilter = ICAInverseFilterType::New();
- m_InverseFilter = invFilter;
+ ICAForwardFilterType::Pointer filter = ICAForwardFilterType::New();
+ m_ForwardFilter = filter;
+ ICAInverseFilterType::Pointer invFilter = ICAInverseFilterType::New();
+ m_InverseFilter = invFilter;
- filter->SetInput(GetParameterFloatVectorImage("in"));
- filter->SetNumberOfPrincipalComponentsRequired(nbComp);
- filter->SetNumberOfIterations(nbIterations);
- filter->SetMu(mu);
-
- switch (GetParameterInt("method.ica.g"))
+ filter->SetInput(GetParameterFloatVectorImage("in"));
+ filter->SetNumberOfPrincipalComponentsRequired(nbComp);
+ filter->SetNumberOfIterations(nbIterations);
+ filter->SetMu(mu);
+
+ switch (GetParameterInt("method.ica.g"))
{
// tanh
case 0:
@@ -406,10 +457,10 @@ private:
invFilter->SetPCATransformationMatrix(filter->GetPCATransformationMatrix());
invFilter->SetTransformationMatrix(filter->GetTransformationMatrix());
}
+ otbAppLogINFO("ICA transform has been computed.");
+ m_TransformationMatrix = filter->GetTransformationMatrix();
- m_TransformationMatrix = filter->GetTransformationMatrix();
-
- break;
+ break;
}
default:
@@ -443,52 +494,54 @@ private:
//Write transformation matrix
std::ofstream outFile;
outFile.open(this->GetParameterString("outmatrix"));
- outFile << std::fixed;
- outFile.precision(10);
+ if (outFile.is_open())
+ {
+ outFile << std::fixed;
+ outFile.precision(10);
- outFile << m_TransformationMatrix;
- outFile.close();
+ outFile << m_TransformationMatrix;
+ outFile.close();
+ }
}
}
- if (!rescale)
+ if (GetParameterString("rescale") == "no")
{
SetParameterOutputImage("out", m_ForwardFilter->GetOutput());
}
else
{
- otbAppLogDEBUG( << "Rescaling " )
- otbAppLogDEBUG( << "Starting Min/Max computation" )
+ otbAppLogINFO("Starting Min/Max computation for rescaling");
- m_MinMaxFilter = MinMaxFilterType::New();
- m_MinMaxFilter->SetInput(m_ForwardFilter->GetOutput());
- //m_MinMaxFilter->GetStreamer()->SetNumberOfLinesStrippedStreaming(50);
- m_MinMaxFilter->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
+ m_MinMaxFilter = MinMaxFilterType::New();
+ m_MinMaxFilter->SetInput(m_ForwardFilter->GetOutput());
+ m_MinMaxFilter->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
- AddProcess(m_MinMaxFilter->GetStreamer(), "Min/Max computing");
- m_MinMaxFilter->Update();
+ AddProcess(m_MinMaxFilter->GetStreamer(), "Min/Max computing");
+ m_MinMaxFilter->Update();
- otbAppLogDEBUG( << "Min/Max computation done : min=" << m_MinMaxFilter->GetMinimum()
- << " max=" << m_MinMaxFilter->GetMaximum() )
+ otbAppLogINFO(<< "Min/Max computation done : min=" << m_MinMaxFilter->GetMinimum() << " max=" << m_MinMaxFilter->GetMaximum())
- FloatVectorImageType::PixelType inMin, inMax;
+ FloatVectorImageType::PixelType inMin,
+ inMax;
- m_RescaleFilter = RescaleImageFilterType::New();
- m_RescaleFilter->SetInput(m_ForwardFilter->GetOutput());
- m_RescaleFilter->SetInputMinimum(m_MinMaxFilter->GetMinimum());
- m_RescaleFilter->SetInputMaximum(m_MinMaxFilter->GetMaximum());
+ m_RescaleFilter = RescaleImageFilterType::New();
+ m_RescaleFilter->SetInput(m_ForwardFilter->GetOutput());
+ m_RescaleFilter->SetAutomaticInputMinMaxComputation(false);
+ m_RescaleFilter->SetInputMinimum(m_MinMaxFilter->GetMinimum());
+ m_RescaleFilter->SetInputMaximum(m_MinMaxFilter->GetMaximum());
- FloatVectorImageType::PixelType outMin, outMax;
- outMin.SetSize(m_ForwardFilter->GetOutput()->GetNumberOfComponentsPerPixel());
- outMax.SetSize(m_ForwardFilter->GetOutput()->GetNumberOfComponentsPerPixel());
- outMin.Fill(GetParameterFloat("rescale.outmin"));
- outMax.Fill(GetParameterFloat("rescale.outmax"));
+ FloatVectorImageType::PixelType outMin, outMax;
+ outMin.SetSize(m_ForwardFilter->GetOutput()->GetNumberOfComponentsPerPixel());
+ outMax.SetSize(m_ForwardFilter->GetOutput()->GetNumberOfComponentsPerPixel());
+ outMin.Fill(GetParameterFloat("rescale.minmax.outmin"));
+ outMax.Fill(GetParameterFloat("rescale.minmax.outmax"));
- m_RescaleFilter->SetOutputMinimum(outMin);
- m_RescaleFilter->SetOutputMaximum(outMax);
- m_RescaleFilter->UpdateOutputInformation();
+ m_RescaleFilter->SetOutputMinimum(outMin);
+ m_RescaleFilter->SetOutputMaximum(outMax);
+ m_RescaleFilter->UpdateOutputInformation();
- SetParameterOutputImage("out", m_RescaleFilter->GetOutput());
+ SetParameterOutputImage("out", m_RescaleFilter->GetOutput());
}
diff --git a/Modules/Applications/AppDimensionalityReduction/test/CMakeLists.txt b/Modules/Applications/AppDimensionalityReduction/test/CMakeLists.txt
index e4037adcb0a388f34ae1ba4cc49a70200bab4d08..82d8b013b9e79e65ba95b9d872b40cfb5d0d061a 100644
--- a/Modules/Applications/AppDimensionalityReduction/test/CMakeLists.txt
+++ b/Modules/Applications/AppDimensionalityReduction/test/CMakeLists.txt
@@ -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)
#-------------------------------------------------------------------------------
diff --git a/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.h b/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.h
index 9d35211893a0519a220620a151dd070721ec4744..94b9b74d9c5eee286ece13ac16222114564da46b 100644
--- a/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.h
+++ b/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.h
@@ -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:
diff --git a/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.hxx b/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.hxx
index d120bbde8ed76bf3efef89be4b33afc55828d5c6..68579a81cf7b9b4e5f6dd5679e36309879b0d3a5 100644
--- a/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.hxx
+++ b/Modules/Filtering/DimensionalityReduction/include/otbMNFImageFilter.hxx
@@ -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
diff --git a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h
index 3368163846bb09339b0fe29eeed9981043f66ac5..6dc5702172e5f2a2240e06e1082e37b9027cdda0 100644
--- a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h
+++ b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.h
@@ -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;
diff --git a/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx b/Modules/Filtering/DimensionalityReduction/include/otbPCAImageFilter.hxx
index 6e49b249ec8ffd5a589e0f900cee9be387aef40c..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;
@@ -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
}