diff --git a/Code/BasicFilters/otbLocalGradientVectorImageFilter.h b/Code/BasicFilters/otbLocalGradientVectorImageFilter.h
index f58e523bdf79c3f104aec816003845e4480a3cc4..a1c371e2d8c2a13b0b9702ec7553edb3efe026bd 100644
--- a/Code/BasicFilters/otbLocalGradientVectorImageFilter.h
+++ b/Code/BasicFilters/otbLocalGradientVectorImageFilter.h
@@ -47,9 +47,9 @@ public:
for ( unsigned int i = 0; i < length; i++ )
{
output[i] = static_cast<typename TOutput::ValueType>(
- input.GetPixel(4)[i]
- - input.GetPixel(1)[i] / 2.
- - input.GetPixel(3)[i] / 2. );
+ input.GetCenterPixel()[i]
+ - input.GetPixel(5)[i] / 2.
+ - input.GetPixel(7)[i] / 2. );
}
return output;
}
diff --git a/Code/Hyperspectral/otbMNFImageFilter.txx b/Code/Hyperspectral/otbMNFImageFilter.txx
index 8b250542a53ce0ed7070c661ca30cd0e6ad2cf9a..3998bd90a094ca215954894c3fae98e1c1da2156 100644
--- a/Code/Hyperspectral/otbMNFImageFilter.txx
+++ b/Code/Hyperspectral/otbMNFImageFilter.txx
@@ -154,6 +154,16 @@ MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransf
m_Normalizer->SetUseStdDev( false );
m_Normalizer->SetInput( inputImgPtr );
+ m_Normalizer->Update();
+
+ if ( !m_GivenMeanValues )
+ m_MeanValues = m_Normalizer->GetFunctor().GetMean();
+
+ if ( m_UseNormalization )
+ {
+ if ( !m_GivenStdDevValues )
+ m_StdDevValues = m_Normalizer->GetFunctor().GetStdDev();
+ }
if ( !m_GivenTransformationMatrix )
{
@@ -205,15 +215,6 @@ MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransf
this->GraftOutput( m_Transformer->GetOutput() );
- /** Once the Normalizer has been updated */
- if ( !m_GivenMeanValues )
- m_MeanValues = m_Normalizer->GetFunctor().GetMean();
-
- if ( m_UseNormalization )
- {
- if ( !m_GivenStdDevValues )
- m_StdDevValues = m_Normalizer->GetFunctor().GetStdDev();
- }
}
template <class TInputImage, class TOutputImage,
@@ -258,7 +259,7 @@ MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransf
if ( m_TransformationMatrix.Rows() == m_TransformationMatrix.Cols() )
{
m_TransformationMatrix = vnl_matrix_inverse< MatrixElementType >
- ( m_TransformationMatrix.GetTranspose() );
+ ( m_TransformationMatrix.GetVnlMatrix() );
}
else
{
@@ -329,21 +330,9 @@ MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransf
InternalMatrixType An = m_NoiseCovarianceMatrix.GetVnlMatrix();
InternalMatrixType W = An * Ax_inv;
- std::cerr << "Matrice C \n";
- W.print( std::cerr );
-
- InternalMatrixType transf;
- vnl_vector<double> vectValP;
-
- vnl_symmetric_eigensystem_compute( W, transf, vectValP );
-
- std::cerr << "Matrice de transf\n";
- transf.print( std::cerr );
-
- std::cerr << "VectPropres\n";
- for ( unsigned int i = 0; i < vectValP.size(); i++ )
- std::cerr << vectValP[i] << " ";
- std::cerr << "\n";
+ vnl_svd< MatrixElementType > solver ( m_CovarianceMatrix.GetVnlMatrix() );
+ InternalMatrixType transf = solver.U();
+ InternalMatrixType valP = solver.W();
InternalMatrixType normMat //= transf.transpose() * An * transf;
= transf.transpose() * Ax_inv * transf;
@@ -356,7 +345,6 @@ MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransf
}
transf.inplace_transpose();
- transf.fliplr();
if ( m_NumberOfPrincipalComponentsRequired
!= this->GetInput()->GetNumberOfComponentsPerPixel() )
@@ -364,14 +352,9 @@ MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransf
else
m_TransformationMatrix = transf;
- std::cerr << "Matrice de projection \n" << m_TransformationMatrix << "\n";
-
m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; i++ )
- m_EigenValues[i]
- = static_cast< RealType >( vectValP[m_NumberOfPrincipalComponentsRequired-1-i] );
-
- std::cerr << "Valeurs propores \n" << m_EigenValues << "\n";
+ m_EigenValues[i] = static_cast< RealType >( valP(i,i) );
}
template <class TInputImage, class TOutputImage,
diff --git a/Code/Hyperspectral/otbNAPCAImageFilter.h b/Code/Hyperspectral/otbNAPCAImageFilter.h
index ea11b034341cce07b1174329bafa4354f83cc050..27c3c08f0f3aa7a86e54c22d9fbb130c4304c68a 100644
--- a/Code/Hyperspectral/otbNAPCAImageFilter.h
+++ b/Code/Hyperspectral/otbNAPCAImageFilter.h
@@ -37,7 +37,7 @@ namespace otb {
* for a better scalability...
*
* This class is very similar to the MNFImageFilter one since only the projection
- * matrix definition (through GetTransformationMatrixFromCovarianceMatrix function)
+ * matrix definition (through GenerateTransformationMatrix function)
* differs.
*
* TODO? Utiliser une 2e entree pour donner directement une image de bruit ??
@@ -92,7 +92,7 @@ protected:
virtual ~NAPCAImageFilter () { }
/** Specific functionality of NAPCA */
- virtual void GetTransformationMatrixFromCovarianceMatrix();
+ virtual void GenerateTransformationMatrix();
}; // end of class
} // end of namespace otb
diff --git a/Code/Hyperspectral/otbNAPCAImageFilter.txx b/Code/Hyperspectral/otbNAPCAImageFilter.txx
index eb91bc27fdcc24e704a7f0cbbe52b4c8bf967923..42ac1bbc64a5eab86efa14f2864e71fd7b9b9982 100644
--- a/Code/Hyperspectral/otbNAPCAImageFilter.txx
+++ b/Code/Hyperspectral/otbNAPCAImageFilter.txx
@@ -33,22 +33,24 @@ template <class TInputImage, class TOutputImage,
Transform::TransformDirection TDirectionOfTransformation >
void
NAPCAImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransformation >
-::GetTransformationMatrixFromCovarianceMatrix ()
+::GenerateTransformationMatrix ()
{
InternalMatrixType An = this->GetNoiseCovarianceMatrix().GetVnlMatrix();
- InternalMatrixType Ax = this->GetCovarianceMatrix().GetVnlMatrix();
+ InternalMatrixType Fn;
+ vnl_vector<double> vectValPn;
+ vnl_symmetric_eigensystem_compute( An, Fn, vectValPn );
+ Fn.inplace_transpose();
- vnl_svd< MatrixElementType > An_solver ( An );
- InternalMatrixType U_cholesky = An_solver.U();
- InternalMatrixType U = vnl_matrix_inverse< MatrixElementType > ( U_cholesky );
- InternalMatrixType C = U.transpose() * Ax * U;
+ InternalMatrixType Ax = this->GetCovarianceMatrix().GetVnlMatrix();
+ InternalMatrixType Aadj = Fn.transpose() * Ax * Fn;
- InternalMatrixType Id ( C.rows(), C.cols(), vnl_matrix_identity );
- vnl_generalized_eigensystem solver ( C, Id );
+ InternalMatrixType Fadj;
+ vnl_vector<double> vectValPadj;
+ vnl_symmetric_eigensystem_compute( Aadj, Fadj, vectValPadj );
- InternalMatrixType transf = solver.V;
- transf *= U.transpose();
+ InternalMatrixType transf = Fn * Fadj;
transf.fliplr();
+ transf.inplace_transpose();
if ( this->GetNumberOfPrincipalComponentsRequired()
!= this->GetInput()->GetNumberOfComponentsPerPixel() )
@@ -56,12 +58,10 @@ NAPCAImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTran
else
this->m_TransformationMatrix = transf;
- vnl_vector< double > valP = solver.D.diagonal();
- valP.flip();
-
this->m_EigenValues.SetSize( this->GetNumberOfPrincipalComponentsRequired() );
for ( unsigned int i = 0; i < this->GetNumberOfPrincipalComponentsRequired(); i++ )
- this->m_EigenValues[i] = static_cast< RealType >( valP[i] );
+ this->m_EigenValues[this->GetNumberOfPrincipalComponentsRequired()-1-i]
+ = static_cast< RealType >( vectValPadj[i] );
}
diff --git a/Code/Hyperspectral/otbPCAImageFilter.txx b/Code/Hyperspectral/otbPCAImageFilter.txx
index d0cf826eeddcdbdf5727a67ff2b86fbe87966f72..4f19ebf1ac7114c4f21e12e3413ccc604e8e1bbf 100644
--- a/Code/Hyperspectral/otbPCAImageFilter.txx
+++ b/Code/Hyperspectral/otbPCAImageFilter.txx
@@ -341,10 +341,6 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; i++ )
m_EigenValues[i] = static_cast< RealType >( valP[i] );
#else
- //vnl_svd< MatrixElementType > solver ( m_CovarianceMatrix.GetVnlMatrix() );
- //InternalMatrixType transf = solver.U();
- //InternalMatrixType valP = solver.W();
-
InternalMatrixType transf;
vnl_vector<double> vectValP;
vnl_symmetric_eigensystem_compute( m_CovarianceMatrix.GetVnlMatrix(), transf, vectValP );
@@ -372,6 +368,7 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
}
}
transf = valP * transf.transpose();
+ transf.flipud();
if ( m_NumberOfPrincipalComponentsRequired
!= this->GetInput()->GetNumberOfComponentsPerPixel() )
@@ -379,14 +376,10 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
else
m_TransformationMatrix = transf;
- // std::cerr << "Matrice de projection \n" << m_TransformationMatrix << "\n";
-
m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; i++ )
m_EigenValues[i] = static_cast< RealType >( valP(i,i) );
- // std::cerr << "Valeurs propores \n" << m_EigenValues << "\n";
-
#endif
}