From 2a2b639787dbb17bb7d0f6ad8c5c36572b8ca0d1 Mon Sep 17 00:00:00 2001
From: Gregoire Mercier <gregoire.mercier@telecom-bretagne.eu>
Date: Wed, 2 Mar 2011 19:06:35 +0100
Subject: [PATCH] ADD: FastICA
---
Code/Hyperspectral/otbFastICAImageFilter.h | 196 +++++++++++++++
Code/Hyperspectral/otbFastICAImageFilter.txx | 250 +++++++++++++++++++
2 files changed, 446 insertions(+)
create mode 100644 Code/Hyperspectral/otbFastICAImageFilter.h
create mode 100644 Code/Hyperspectral/otbFastICAImageFilter.txx
diff --git a/Code/Hyperspectral/otbFastICAImageFilter.h b/Code/Hyperspectral/otbFastICAImageFilter.h
new file mode 100644
index 0000000000..15553eb7be
--- /dev/null
+++ b/Code/Hyperspectral/otbFastICAImageFilter.h
@@ -0,0 +1,196 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbFastICAImageFilter_h
+#define __otbFastICAImageFilter_h
+
+#include "otbMacro.h"
+#include "itkImageToImageFilter.h"
+#include "otbPCAImageFilter.h"
+#include "otbMatrixMultiplyImageFilter.h"
+
+namespace otb
+{
+/** \class FastICAImageFilter
+ * \brief Performs a Independent Component Analysis
+ *
+ * The internal structure of this filter is a filter-to-filter like structure.
+ * The estimation of the covariance matrix has persistent capabilities...
+ *
+ * \sa PCAImageFilter
+ */
+template <class TInputImage, class TOutputImage, Transform::TransformDirection TDirectionOfTransformation >
+class ITK_EXPORT FastICAImageFilter
+ : public itk::ImageToImageFilter<TInputImage, TOutputImage>
+{
+public:
+ /** Standard typedefs */
+ typedef FastICAImageFilter Self;
+ typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Type macro */
+ itkNewMacro(Self);
+
+ /** Creation through object factory macro */
+ itkTypeMacro(FastICAImageFilter, ImageToImageFilter);
+
+ /** Dimension */
+ itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+ itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);
+
+ typedef Transform::TransformDirection TransformDirectionEnumType;
+ itkStaticConstMacro(DirectionOfTransformation,TransformDirectionEnumType,TDirectionOfTransformation);
+
+ /** typedefs */
+ typedef TInputImage InputImageType;
+ typedef TOutputImage OutputImageType;
+
+ typedef PCAImageFilter< InputImageType, OutputImageType, TDirectionOfTransformation >
+ PCAFilterType;
+ typedef typename PCAFilterType::Pointer PCAFilterPointerType;
+
+ typedef typename PCAFilterType::RealType RealType;
+ typedef typename PCAFilterType::VectorType VectorType;
+ typedef typename PCAFilterType::MatrixType MatrixType;
+ typedef typename MatrixType::InternalMatrixType InternalMatrixType;
+ typedef typename InternalMatrixType::element_type MatrixElementType;
+
+ typedef MatrixMultiplyImageFilter< TInputImage, TOutputImage, RealType > TransformFilterType;
+ typedef typename TransformFilterType::Pointer TransformFilterPointerType;
+
+ typedef double (*ContrastFunctionType) ( double );
+
+
+ /**
+ * Set/Get the number of required largest principal components.
+ */
+ itkGetMacro(NumberOfPrincipalComponentsRequired,unsigned int);
+ itkSetMacro(NumberOfPrincipalComponentsRequired,unsigned int);
+
+ itkGetConstMacro(PCAFilter,PCAFilterType *);
+ itkGetMacro(PCAFilter,PCAFilterType *);
+ itkSetMacro(PCAFilter,PCAFilterType *);
+
+ itkGetConstMacro(TransformFilter,TransformFilterType *);
+ itkGetMacro(TransformFilter,TransformFilterType *);
+ itkSetMacro(TransformFilter,TransformFilterType *);
+
+ VectorType GetMeanValues () const
+ {
+ return this->GetPCAFilter()->GetMeanValues();
+ }
+ void SetMeanValues ( const VectorType & vec )
+ {
+ m_PCAFilter->SetMeanValues(vec);
+ }
+
+ VectorType GetStdDevValues ( ) const
+ {
+ return this->GetPCAFilter()->GetStdDevValues();
+ }
+ void SetStdDevValues ( const VectorType & vec )
+ {
+ m_PCAFilter->SetStdDevValues(vec);
+ }
+
+ MatrixType GetPCATransformationMatrix () const
+ {
+ return this->GetPCAFilter()->GetTransformationMatrix();
+ }
+ void SetPCACovarianceMatrix ( const MatrixType & mat, bool isForward = true )
+ {
+ m_PCAFilter->SetTransformationMatrix(mat,isForward);
+ }
+
+ itkGetConstMacro(TransformationMatrix,MatrixType);
+ itkGetMacro(TransformationMatrix,MatrixType);
+ void SetTransformationMatrix ( const MatrixType & mat, bool isForward = true )
+ {
+ m_IsTransformationForward = isForward;
+ m_GivenTransformationMatrix = true;
+ m_TransformationMatrix = mat;
+ this->Modified();
+ }
+
+ itkGetMacro(MaximumOfIterations,unsigned int);
+ itkSetMacro(MaximumOfIterations,unsigned int);
+
+ itkGetMacro(ConvergenceThreshold,double);
+ itkSetMacro(ConvergenceThreshold,double);
+
+ itkGetMacro(ContrastFunction,ContrastFunctionType);
+
+protected:
+ FastICAImageFilter ();
+ virtual ~FastICAImageFilter() { }
+
+ /** GenerateOutputInformation
+ * Propagate vector length info and modify if needed
+ * NumberOfPrincipalComponentsRequired
+ *
+ * In REVERSE mode, the covariance matrix or the transformation matrix
+ * (which may not be square) has to be given,
+ * otherwize, GenerateOutputInformation throws an itk::ExceptionObject
+ */
+ virtual void GenerateOutputInformation();
+
+ /** GenerateData
+ * Through a filter of filter structure
+ */
+ virtual void GenerateData ();
+
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+ /** Internal methods */
+ virtual void ForwardGenerateData();
+ virtual void ReverseGenerateData();
+
+ /** this is the specifical part of FastICA */
+ virtual void GenerateTransformationMatrix();
+
+ unsigned int m_NumberOfPrincipalComponentsRequired;
+
+ /** Transformation matrix refers to the ICA step (not PCA) */
+ bool m_GivenTransformationMatrix;
+ bool m_IsTransformationForward;
+ MatrixType m_TransformationMatrix;
+
+ /** FastICA parameters */
+ unsigned int m_MaximumOfIterations; // def is 50
+ double m_ConvergenceThreshold; // def is 1e-4
+ ContrastFunctionType m_ConstrastFunction; // see g() function in the biblio. Def is tanh
+
+ PCAFilterPointerType m_PCAFilter;
+ TransformFilterPointerType m_TransformFilter;
+
+private:
+ FastICAImageFilter ( const Self & );
+ void operator= ( const Self & );
+}; // end of class
+
+} // end of namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbFastICAImageFilter.txx"
+#endif
+
+
+#endif
+
+
diff --git a/Code/Hyperspectral/otbFastICAImageFilter.txx b/Code/Hyperspectral/otbFastICAImageFilter.txx
new file mode 100644
index 0000000000..05ff275418
--- /dev/null
+++ b/Code/Hyperspectral/otbFastICAImageFilter.txx
@@ -0,0 +1,250 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbFastICAImageFilter_txx
+#define __otbFastICAImageFilter_txx
+#include "otbFastICAImageFilter.h"
+
+#include "otbMacro.h"
+
+#include "itkExceptionObject.h"
+#include "itkNumericTraits.h"
+
+#include <vnl/vnl_matrix.h>
+#include <vnl/algo/vnl_matrix_inverse.h>
+#include <vnl/algo/vnl_generalized_eigensystem.h>
+
+namespace otb
+{
+
+template < class TInputImage, class TOutputImage,
+ Transform::TransformDirection TDirectionOfTransformation >
+FastICAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
+::FastICAImageFilter ()
+{
+ this->SetNumberOfRequiredInputs(1);
+
+ m_NumberOfPrincipalComponentsRequired = 0;
+
+ m_GivenTransformationMatrix = false;
+ m_IsTransformationForward = true;
+
+ m_MaximumOfIterations = 50;
+ m_ConvergenceThreshold = 1E-4;
+ m_ConstrastFunction = &vcl_tanh;
+
+ m_PCAFilter = PCAFilterType::New();
+ m_PCAFilter->UseNormalization();
+
+ m_TransformFilter = TransformFilterType::New();
+}
+
+template < class TInputImage, class TOutputImage,
+ Transform::TransformDirection TDirectionOfTransformation >
+void
+FastICAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
+::GenerateOutputInformation()
+// throw itk::ExceptionObject
+{
+ Superclass::GenerateOutputInformation();
+
+ switch ( DirectionOfTransformation )
+ {
+ case Transform::FORWARD:
+ {
+ if ( m_NumberOfPrincipalComponentsRequired == 0
+ || m_NumberOfPrincipalComponentsRequired
+ > this->GetInput()->GetNumberOfComponentsPerPixel() )
+ {
+ m_NumberOfPrincipalComponentsRequired =
+ this->GetInput()->GetNumberOfComponentsPerPixel();
+ }
+
+ this->GetOutput()->SetNumberOfComponentsPerPixel(
+ m_NumberOfPrincipalComponentsRequired );
+ break;
+ }
+ case Transform::INVERSE:
+ {
+ unsigned int theOutputDimension = 0;
+ if ( m_GivenTransformationMatrix )
+ {
+ theOutputDimension = m_TransformationMatrix.Rows() >= m_TransformationMatrix.Cols() ?
+ m_TransformationMatrix.Rows() : m_TransformationMatrix.Cols();
+ }
+ else if ( m_GivenCovarianceMatrix )
+ {
+ theOutputDimension = m_CovarianceMatrix.Rows() >= m_CovarianceMatrix.Cols() ?
+ m_CovarianceMatrix.Rows() : m_CovarianceMatrix.Cols();
+ }
+ else
+ {
+ throw itk::ExceptionObject(__FILE__, __LINE__,
+ "Mixture matrix is required to know the output size",
+ ITK_LOCATION);
+ }
+
+ this->GetOutput()->SetNumberOfComponentsPerPixel( theOutputDimension );
+
+ break;
+ }
+ default:
+ throw itk::ExceptionObject(__FILE__, __LINE__,
+ "Class should be templeted with FORWARD or INVERSE only...",
+ ITK_LOCATION );
+ }
+}
+
+template < class TInputImage, class TOutputImage,
+ Transform::TransformDirection TDirectionOfTransformation >
+void
+FastICAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
+::GenerateData ()
+{
+ switch ( DirectionOfTransformation )
+ {
+ case Transform::FORWARD:
+ return ForwardGenerateData();
+ case Transform::INVERSE:
+ return ReverseGenerateData();
+ default:
+ throw itk::ExceptionObject(__FILE__, __LINE__,
+ "Class should be templated with FORWARD or INVERSE only...",
+ ITK_LOCATION );
+ }
+}
+
+template < class TInputImage, class TOutputImage,
+ Transform::TransformDirection TDirectionOfTransformation >
+void
+FastICAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
+::ForwardGenerateData ()
+{
+ typename InputImageType::Pointer inputImgPtr
+ = const_cast<InputImageType*>( this->GetInput() );
+
+ m_PCAFilter->SetInput( inputImgPtr );
+ m_PCAFilter->Update();
+
+ if ( !m_GivenTransformationMatrix )
+ {
+ GenerateTransformationMatrix();
+ }
+ else if ( !m_IsTransformationForward )
+ {
+ // prevent from multiple inversion in the pipelines
+ m_IsTransformationForward = true;
+ vnl_svd< MatrixElementType > invertor ( m_TransformationMatrix.GetVnlMatrix() );
+ m_TransformationMatrix = invertor.pinverse();
+ }
+
+ if ( m_TransformationMatrix.GetVnlMatrix().empty() )
+ {
+ throw itk::ExceptionObject( __FILE__, __LINE__,
+ "Empty transformation matrix",
+ ITK_LOCATION);
+ }
+
+ m_TransformFilter->SetInput( m_PCAFilter->GetOutput() );
+ m_TransformFilter->SetMatrix( m_TransformationMatrix.GetVnlMatrix() );
+ m_TransformFilter->GraftOutput( this->GetOutput() );
+ m_TransformFilter->Update();
+
+ this->GraftOutput( m_TransformFilter->GetOutput() );
+}
+
+template < class TInputImage, class TOutputImage,
+ Transform::TransformDirection TDirectionOfTransformation >
+void
+FastICAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
+::ReverseGenerateData ()
+{
+ if ( !m_GivenTransformationMatrix )
+ {
+ throw itk::ExceptionObject( __FILE__, __LINE__,
+ "No Transformation matrix given",
+ ITK_LOCATION );
+ }
+
+ if ( m_TransformationMatrix.GetVnlMatrix().empty() )
+ {
+ throw itk::ExceptionObject( __FILE__, __LINE__,
+ "Empty transformation matrix",
+ ITK_LOCATION);
+ }
+
+ if ( m_IsTransformationForward )
+ {
+ // prevent from multiple inversion in the pipelines
+ m_IsTransformationForward = false;
+ vnl_svd< MatrixElementType > invertor ( m_TransformationMatrix.GetVnlMatrix() );
+ m_TransformationMatrix = invertor.pinverse();
+ }
+
+ m_TransformFilter->SetInput( this->GetInput() );
+ m_TransformFilter->SetMatrix( m_TransformationMatrix.GetVnlMatrix() );
+
+ /*
+ * PCA filter may throw exception if
+ * the mean, stdDev and transformation matrix
+ * have not been given at this point
+ */
+ m_PCAFilter->SetInput( m_TransformFilter->GetOutput() );
+ m_PCAFilter->GraftOutput( this->GetOutput() );
+ m_PCAFilter->Update();
+
+ this->GraftOutput( m_PCAFilter->GetOutput() );
+}
+
+template < class TInputImage, class TOutputImage,
+ Transform::TransformDirection TDirectionOfTransformation >
+void
+FastICAImageFilter< TInputImage, TOutputImage, TDirectionOfTransformation >
+::GenerateTransformationMatrix ()
+{
+ double convergence = itk::NumericTraits<double>::Max();
+ unsigned int iteration = 0;
+
+ unsigned int size = this->GetInput()->GetNumberOfComponentsPerPixel();
+
+ // transformation matrix
+ InternalMatrixType W ( size, size, vnl_matrix_identity );
+
+ while ( iteration++ < GetMaximumOfIterations()
+ && convergence > GetConvergenceThreshold() )
+ {
+ otbMsgDebugMacro( "Iteration " << iteration << " / " << GetMaximumOfIterations()
+ << ", MSE = " << convergence );
+
+ InternalMatrixType W_old ( W );
+
+ // TODO le premier coup ne sert a rien
+ TransformFilterPointerType transformer = TransformFilterType::New();
+ transformer->SetIntput( GetPCAFilter()->GetOutput() );
+ transformer->SetMatrix( W );
+ transformer->Update();
+
+ // Faire un image to image filter...
+
+ } // end of while loop
+}
+
+} // end of namespace otb
+
+#endif
+
+
--
GitLab