diff --git a/Code/Hyperspectral/otbPCAImageFilter.h b/Code/Hyperspectral/otbPCAImageFilter.h index 0758c8059dda517793ae524d6c43dd27d69cfef4..314302c85a715ca01d7722b7c090979f1b761c61 100644 --- a/Code/Hyperspectral/otbPCAImageFilter.h +++ b/Code/Hyperspectral/otbPCAImageFilter.h @@ -78,6 +78,8 @@ public: typedef typename CovarianceEstimatorFilterType::RealPixelType VectorType; typedef typename CovarianceEstimatorFilterType::MatrixObjectType MatrixObjectType; typedef typename MatrixObjectType::ComponentType MatrixType; + typedef typename MatrixType::InternalMatrixType InternalMatrixType; + typedef typename InternalMatrixType::element_type MatrixElementType; typedef MatrixMultiplyImageFilter< TInputImage, TOutputImage, RealType > TransformFilterType; typedef typename TransformFilterType::Pointer TransformFilterPointerType; @@ -180,7 +182,7 @@ protected: virtual void ForwardGenerateData(); virtual void ReverseGenerateData(); - void GetTransformationMatrixFromCovarianceMatrix(); + void GenerateTransformationMatrix(); /** Internal attributes */ unsigned int m_NumberOfPrincipalComponentsRequired; diff --git a/Code/Hyperspectral/otbPCAImageFilter.txx b/Code/Hyperspectral/otbPCAImageFilter.txx index 6ce15760d007d9189f64c3934afedea7ad8ebd04..9756736305b1de85b039993c2c54ab7aa5121f42 100644 --- a/Code/Hyperspectral/otbPCAImageFilter.txx +++ b/Code/Hyperspectral/otbPCAImageFilter.txx @@ -188,7 +188,7 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation > m_Transformer->SetInput( inputImgPtr ); } - GetTransformationMatrixFromCovarianceMatrix(); + GenerateTransformationMatrix(); } else if ( !m_IsTransformationMatrixForward ) { @@ -225,7 +225,7 @@ PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation > ITK_LOCATION ); } - GetTransformationMatrixFromCovarianceMatrix(); + GenerateTransformationMatrix(); m_TransformationMatrix = m_TransformationMatrix.GetTranspose(); } else if ( m_IsTransformationMatrixForward ) @@ -291,8 +291,13 @@ template < class TInputImage, class TOutputImage, Transform::TransformDirection TDirectionOfTransformation > void PCAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation > -::GetTransformationMatrixFromCovarianceMatrix () +::GenerateTransformationMatrix () { +#if 0 + /* + * Old stuff but did work ! + */ + MatrixType Id ( m_CovarianceMatrix ); Id.SetIdentity(); @@ -331,6 +336,42 @@ 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] ); +#else + vnl_svd< MatrixElementType > solver ( m_CovarianceMatrix.GetVnlMatrix() ); + InternalMatrixType transf = solver.U(); + + InternalMatrixType valP = solver.W(); + /* We used normalized PCA */ + for ( unsigned int i = 0; i < valP.rows(); i++ ) + { + if ( valP(i,i) == 0. ) + valP(i,i) = 1. / vcl_sqrt( vcl_abs( valP(i,i) ) ); + else + otbMsgDebugMacro( << "ValP(" << i << ") null : " << valP(i,i) ); + //throw itk::ExceptionObject( __FILE__, __LINE__, + // "Null Eigen value !!", ITK_LOCATION ); + } + transf = valP * transf.transpose(); + + 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,i) ); + + std::cerr << "Matrice de transformation \n"; + std::cerr << m_TransformationMatrix << "\n"; + + std::cerr << "Vecteurs propres\n"; + std::cerr << m_EigenValues << "\n"; + + std::cerr << "Matrice W\n" << solver.W(); +#endif }