ADD:NAPCA functionality

Code/Hyperspectral/otbNAPCAImageFilter.h

+/*=========================================================================
+
+  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 __otbNAPCAImageFilter_h
+#define __otbNAPCAImageFilter_h
+
+#include "otbMacro.h"
+#include "itkImageToImageFilter.h"
+#include "otbStreamingStatisticsVectorImageFilter2.h"
+#include "otbMatrixMultiplyImageFilter.h"
+#include "otbNormalizeVectorImageFilter.h"
+#include "otbMNFImageFilter.h"
+
+namespace otb {
+
+/** \class NAPCAImageFilter
+ * \brief Performs a Noise Adjusted PCA analysis of a vector image.
+ *
+ * The internal structure of this filter is a filter-to-filter like structure. 
+ * The estimation of the covariance matrix has persistent capabilities...
+ *
+ * The high pass filter which has to be used for the noise estimation is templated
+ * for a better scalability...
+ *
+ * This class is very similar to the MNFImageFilter one since only the projection 
+ * matrix definition (through GetTransformationMatrixFromCovarianceMatrix function)
+ * differs.
+ *
+ * TODO? Utiliser une 2e entree pour donner directement une image de bruit ??
+ *
+ * \sa otbStreamingStatisticsVectorImageFilter2
+ * \sa MNFImageFilter
+ */
+template <class TInputImage, class TOutputImage, 
+            class TNoiseImageFilter, 
+            Transform::TransformDirection TDirectionOfTransformation >
+class ITK_EXPORT NAPCAImageFilter 
+  : public MNFImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransformation >
+{
+public:
+  typedef NAPCAImageFilter Self;
+  typedef MNFImageFilter< TInputImage, TOutputImage, 
+            TNoiseImageFilter, TDirectionOfTransformation > Superclass;
+  typedef itk::SmartPointer<Self>                            Pointer;
+  typedef itk::SmartPointer<const Self>                      ConstPointer;
+
+  /** Type macro */
+  itkNewMacro(Self);
+
+  /** Creation through object factory macro */
+  itkTypeMacro(NAPCAImageFilter,MNFImageFilter);
+
+  /** Template parameters typedefs */
+  typedef typename Superclass::InputImageType InputImageType;
+  typedef typename Superclass::OutputImageType OutputImageType;
+
+  typedef typename Superclass::CovarianceEstimatorFilterType CovarianceEstimatorFilterType;
+  typedef typename Superclass::CovarianceEstimatorFilterPointerType CovarianceEstimatorFilterPointerType;
+
+  typedef typename Superclass::RealType RealType;
+  typedef typename Superclass::VectorType VectorType;
+  typedef typename Superclass::MatrixObjectType MatrixObjectType;
+  typedef typename Superclass::MatrixType MatrixType;
+  typedef typename Superclass::InternalMatrixType InternalMatrixType;
+  typedef typename Superclass::MatrixElementType MatrixElementType;
+
+  typedef typename Superclass::TransformFilterType TransformFilterType;
+  typedef typename Superclass::TransformFilterPointerType TransformFilterPointerType;
+
+  typedef typename Superclass::NoiseImageFilterType NoiseImageFilterType;
+  typedef typename Superclass::NoiseImageFilterPointerType NoiseImageFilterPointerType;
+
+  typedef typename Superclass::NormalizeFilterType NormalizeFilterType;
+  typedef typename Superclass::NormalizeFilterPointerType NormalizeFilterPointerType;
+
+protected:
+  NAPCAImageFilter() { }
+  virtual ~NAPCAImageFilter () { }
+
+  /** Specific functionality of NAPCA */
+  virtual void GetTransformationMatrixFromCovarianceMatrix();
+}; // end of class
+
+} // end of namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbNAPCAImageFilter.txx"
+#endif
+
+#endif
+

Code/Hyperspectral/otbNAPCAImageFilter.txx

+/*=========================================================================
+
+  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 __otbNAPCAImageFilter_txx
+#define __otbNAPCAImageFilter_txx
+#include "otbNAPCAImageFilter.h"
+
+#include "itkExceptionObject.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, 
+            class TNoiseImageFilter, 
+            Transform::TransformDirection TDirectionOfTransformation >
+void
+NAPCAImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransformation >
+::GetTransformationMatrixFromCovarianceMatrix ()
+{
+  InternalMatrixType An = m_NoiseCovarianceMatrix.GetVnlMatrix();
+  InternalMatrixType Ax = m_CovarianceMatrix.GetVnlMatrix();
+
+  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 Id ( C.rows(), C.cols(), vnl_matrix_identity );
+  vnl_generalized_eigensystem solver ( C, I );
+
+  InternalMatrixType transf = solver.V;
+  transf *= U.transpose();
+  transf.fliplr();
+
+  if ( m_NumberOfPrincipalComponentsRequired 
+      != this->GetInput()->GetNumberOfComponentsPerPixel() )
+    m_TransformationMatrix = transf.get_n_rows( 0, m_NumberOfPrincipalComponentsRequired );
+  else
+    m_TransformationMatrix = transf;
+
+  vnl_vector< double > valP = solver.D.diagonal();
+  valP.flip();
+
+  m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired );
+  for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; i++ )
+    m_EigenValues[i] = static_cast< RealType >( valP[i] );
+}
+
+
+} // end of namespace otb
+
+#endif
+
+