ADD: add DotProductImageFilter and ProjectiveProjectionImageFilter

Code/BasicFilters/otbDotProductImageFilter.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 __otbDotProductImageFilter_h
+#define __otbDotProductImageFilter_h
+
+#include "otbUnaryFunctorImageFilter.h"
+
+namespace otb
+{
+
+
+namespace Functor {
+
+/** \class DotProductFunctor
+ *
+ * \brief Apply a matrix multiplication to a VariableLengthVector
+ *
+ */
+template<class TInput, class TOutput>
+class DotProductFunctor
+{
+public:
+  typedef TInput     InputType;
+  typedef TOutput    OutputType;
+
+  DotProductFunctor() {}
+  virtual ~DotProductFunctor() {}
+
+  bool operator !=(const DotProductFunctor& other) const
+  {
+    return false;
+  }
+
+  bool operator ==(const DotProductFunctor& other) const
+  {
+    return !(*this != other);
+  }
+
+  const InputType& GetVector()
+  {
+    return m_Vector;
+  }
+
+  void SetVector(const InputType& m)
+  {
+    m_Vector = m;
+  }
+
+  OutputType operator ()(const InputType& in)
+  {
+    OutputType result = 0;
+    for(unsigned int i = 0; i < in.Size(); i++)
+      {
+      result += in[i] * m_Vector[i];
+      }
+    return result;
+  }
+
+private:
+  InputType m_Vector;
+};
+}
+
+/** \class DotProductImageFilter
+ *
+ * \brief Applies pixel-wise dot product to a VectorImage
+ *
+ * Given a vector, this filter outputs the dot product of each pixel of a multiband image
+ *
+ * \ingroup Streamed
+ * \ingroup Threaded
+ */
+template <class TInputImage, class TOutputImage>
+class ITK_EXPORT DotProductImageFilter :
+  public itk::UnaryFunctorImageFilter<TInputImage, TOutputImage,
+      Functor::DotProductFunctor<typename TInputImage::PixelType,
+          typename TOutputImage::PixelType> >
+{
+public:
+  /** Standard class typedefs. */
+  typedef DotProductImageFilter Self;
+  typedef itk::UnaryFunctorImageFilter
+     <TInputImage,
+      TOutputImage,
+      Functor::DotProductFunctor<
+          typename TInputImage::PixelType,
+          typename TOutputImage::PixelType>
+     >                                 Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  typedef Functor::DotProductFunctor<
+      typename TInputImage::PixelType,
+      typename TOutputImage::PixelType> FunctorType;
+  typedef typename FunctorType::MatrixType MatrixType;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(DotProductImageFilter, itk::UnaryFunctorImageFilter);
+
+  /** Pixel types. */
+  typedef typename TInputImage::PixelType  InputPixelType;
+  typedef typename TOutputImage::PixelType OutputPixelType;
+
+  const InputPixelType& GetVector()
+  {
+    return this->GetFunctor().GetVector();
+  }
+
+  void SetVector(const InputPixelType& p)
+  {
+    this->GetFunctor().SetVector(p);
+    this->Modified();
+  }
+
+protected:
+  DotProductImageFilter();
+
+  virtual ~DotProductImageFilter() {}
+
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+  DotProductImageFilter(const Self &); //purposely not implemented
+  void operator =(const Self&); //purposely not implemented
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbDotProductImageFilter.txx"
+#endif
+
+#endif

Code/BasicFilters/otbDotProductImageFilter.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 __otbDotProductImageFilter_txx
+#define __otbDotProductImageFilter_txx
+
+#include "otbDotProductImageFilter.h"
+
+namespace otb
+{
+
+/**
+ *
+ */
+template <class TInputImage, class TOutputImage>
+DotProductImageFilter<TInputImage, TOutputImage>
+::DotProductImageFilter()
+{
+}
+
+template <class TInputImage, class TOutputImage>
+void
+DotProductImageFilter<TInputImage, TOutputImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+}
+
+} // end namespace
+
+#endif

Code/BasicFilters/otbProjectiveProjectionImageFilter.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 __otbProjectiveProjectionImageFilter_h
+#define __otbProjectiveProjectionImageFilter_h
+
+#include "otbUnaryFunctorImageFilter.h"
+#include "vnl/algo/vnl_svd.h"
+#include <boost/shared_ptr.hpp>
+
+namespace otb
+{
+
+
+namespace Functor {
+
+/** \class ProjectiveProjectionFunctor
+ *
+ * \brief TODO
+ *
+ */
+template<class TInput, class TOutput, class TPrecision>
+class ProjectiveProjectionFunctor
+{
+public:
+  typedef TInput     InputType;
+  typedef TOutput    OutputType;
+  typedef TPrecision PrecisionType;
+
+  ProjectiveProjectionFunctor() : m_OutputSize(0) {}
+  virtual ~ProjectiveProjectionFunctor() {}
+
+  unsigned int GetOutputSize()
+  {
+    return m_OutputSize;
+  }
+
+  const InputType& GetProjectionDirection()
+  {
+    return m_ProjectionDirection;
+  }
+
+  void SetProjectionDirection(const InputType& p)
+  {
+    m_ProjectionDirection = p;
+    m_OutputSize = m_ProjectionDirection.Size();
+  }
+
+  bool operator !=(const ProjectiveProjectionFunctor& other) const
+  {
+    return false;
+  }
+
+  bool operator ==(const ProjectiveProjectionFunctor& other) const
+  {
+    return !(*this != other);
+  }
+
+  OutputType operator ()(const InputType& in)
+  {
+    PrecisionType dotProduct = 0;
+    for (unsigned int i = 0; i < in.Size(); i++)
+      {
+      dotProduct += in[i] * m_ProjectionDirection[i];
+      }
+
+    OutputType projected(in.Size());
+    for (unsigned int j = 0; j < in.Size(); j++)
+      {
+      projected[j] = in[j] / dotProduct;
+      }
+
+    return projected;
+  }
+
+private:
+  unsigned int   m_OutputSize;
+  InputType      m_ProjectionDirection;
+};
+}
+
+/** \class ProjectiveProjectionImageFilter
+ *
+ * \brief Applies a projective projection to each pixel of an image
+ *
+ * \ingroup Streamed
+ * \ingroup Threaded
+ */
+template <class TInputImage, class TOutputImage, class TPrecision>
+class ITK_EXPORT ProjectiveProjectionImageFilter :
+  public otb::UnaryFunctorImageFilter<TInputImage, TOutputImage,
+      Functor::ProjectiveProjectionFunctor<typename TInputImage::PixelType,
+          typename TOutputImage::PixelType, TPrecision> >
+{
+public:
+  /** Standard class typedefs. */
+  typedef ProjectiveProjectionImageFilter Self;
+  typedef otb::UnaryFunctorImageFilter
+     <TInputImage,
+      TOutputImage,
+      Functor::ProjectiveProjectionFunctor<
+          typename TInputImage::PixelType,
+          typename TOutputImage::PixelType,
+          TPrecision>
+     >                                 Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  typedef Functor::ProjectiveProjectionFunctor<
+      typename TInputImage::PixelType,
+      typename TOutputImage::PixelType,
+      TPrecision> FunctorType;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(ProjectiveProjectionImageFilter, otb::UnaryFunctorImageFilter);
+
+  /** Pixel types. */
+  typedef typename TInputImage::PixelType  InputPixelType;
+  typedef typename TOutputImage::PixelType OutputPixelType;
+
+  const InputPixelType& GetProjectionDirection()
+  {
+    return this->GetFunctor().GetProjectionDirection();
+  }
+
+  void SetProjectionDirection(const InputPixelType& p)
+  {
+    this->GetFunctor().SetProjectionDirection(p);
+    this->Modified();
+  }
+
+protected:
+  ProjectiveProjectionImageFilter();
+
+  virtual ~ProjectiveProjectionImageFilter() {}
+
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  virtual void GenerateOutputInformation()
+  {
+    Superclass::GenerateOutputInformation();
+    typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
+    typename Superclass::InputImagePointer  inputPtr = this->GetInput();
+    outputPtr->SetNumberOfComponentsPerPixel( // propagate vector length info
+      this->GetInput()->GetNumberOfComponentsPerPixel());
+  }
+
+private:
+  ProjectiveProjectionImageFilter(const Self &); //purposely not implemented
+
+  void operator =(const Self&); //purposely not implemented
+
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbProjectiveProjectionImageFilter.txx"
+#endif
+
+#endif
+

Code/BasicFilters/otbProjectiveProjectionImageFilter.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 __otbProjectiveProjectionImageFilter_txx
+#define __otbProjectiveProjectionImageFilter_txx
+
+#include "otbProjectiveProjectionImageFilter.h"
+
+namespace otb
+{
+
+/**
+ *
+ */
+template <class TInputImage, class TOutputImage, class TPrecision>
+ProjectiveProjectionImageFilter<TInputImage, TOutputImage, TPrecision>
+::ProjectiveProjectionImageFilter()
+{
+}
+
+template <class TInputImage, class TOutputImage, class TPrecision>
+void
+ProjectiveProjectionImageFilter<TInputImage, TOutputImage, TPrecision>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+}
+
+} // end namespace
+
+#endif

Testing/CMakeLists.txt

@@ -11,6 +11,7 @@ SET( otbHyperTests1_SRC
   otbMatrixMultiplyImageFilter.cxx
   otbUnConstrainedLeastSquareImageFilter.cxx
   otbVectorImageToMatrixImageFilter.cxx
+  otbProjectiveProjection.cxx
 )
 
 ADD_EXECUTABLE(otbHyperTests1 ${otbHyperTests1_SRC} otbHyperTests1.cxx)
@@ -115,4 +116,10 @@ ADD_TEST(hyEigenvalueLikelihoodMaximization2
          ${DATA}/SYNTHETIC/hsi_cube.tif
          ${TEMP}/synthetic_elm.txt)
 
-         
\ No newline at end of file
+ADD_TEST(testProjectiveProjection
+         ${TESTEXE_DIR}/otbHyperTests1
+         otbProjectiveProjectionTestHighSNR
+         ${DATA}/SYNTHETIC/hsi_cube.tif
+         8
+         ${TEMP}/hsi_cube_projectiveprojected_8.tif )
+

Testing/otbHyperTests1.cxx

@@ -36,4 +36,5 @@ void RegisterTests()
   REGISTER_TEST(otbVectorImageToMatrixTest);
   REGISTER_TEST(vahineVCA);
   REGISTER_TEST(vahineElmVCA);
+  REGISTER_TEST(otbProjectiveProjectionTestHighSNR);
 }

Testing/otbProjectiveProjection.cxx

+/*=========================================================================
+
+  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.
+
+=========================================================================*/
+#include "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+
+#include "otbDotProductImageFilter.h"
+#include "otbProjectiveProjectionImageFilter.h"
+#include "otbMatrixMultiplyImageFilter.h"
+#include "otbVectorImageToMatrixImageFilter.h"
+#include "otbStreamingStatisticsImageFilter.h"
+#include "otbStreamingStatisticsVectorImageFilter2.h"
+
+const unsigned int Dimension = 2;
+typedef double PixelType;
+typedef double PrecisionType;
+
+typedef otb::Image<PixelType, Dimension> ImageType;
+typedef otb::VectorImage<PixelType, Dimension> VectorImageType;
+typedef otb::ImageFileReader<VectorImageType> ReaderType;
+typedef otb::ProjectiveProjectionImageFilter<VectorImageType,VectorImageType,PrecisionType> ProjectiveProjectionImageFilterType;
+typedef otb::DotProductImageFilter<VectorImageType,ImageType> DotProductImageFilterType;
+typedef otb::MatrixMultiplyImageFilter<VectorImageType,VectorImageType,PrecisionType> MatrixMultiplyImageFilterType;
+typedef otb::VectorImageToMatrixImageFilter<VectorImageType> VectorImageToMatrixImageFilterType;
+typedef otb::ImageFileWriter<VectorImageType> WriterType;
+typedef otb::StreamingStatisticsVectorImageFilter2<VectorImageType> StreamingStatisticsVectorImageFilterType;
+typedef otb::StreamingStatisticsImageFilter<ImageType> StreamingStatisticsImageFilterType;
+
+typedef StreamingStatisticsVectorImageFilterType::MatrixType MatrixType;
+
+int otbProjectiveProjectionTestHighSNR(int argc, char * argv[])
+{
+  const char * inputImage = argv[1];
+  const unsigned int nbEndmembers = atoi(argv[2]);
+  const char * outputImage = argv[3];
+
+
+  ReaderType::Pointer readerImage = ReaderType::New();
+  readerImage->SetFileName(inputImage);
+
+  std::cout << "Computing image stats" << std::endl;
+  StreamingStatisticsVectorImageFilterType::Pointer statsInput = \
+      StreamingStatisticsVectorImageFilterType::New();
+
+  statsInput->SetInput(readerImage->GetOutput());
+  statsInput->Update();
+
+  std::cout << "Computing SVD of correlation matrix" << std::endl;
+  // Take the correlation matrix
+  vnl_matrix<PrecisionType> R = statsInput->GetCorrelation().GetVnlMatrix();
+
+  // Apply SVD
+  vnl_svd<PrecisionType>    svd(R);
+  vnl_matrix<PrecisionType> U = svd.U();
+  std::cout << "U : "<< U.rows() << " " << U.cols() << std::endl;
+  vnl_matrix<PrecisionType> Ud = U.get_n_columns(0, nbEndmembers).transpose();
+  std::cout << "Ud : "<< Ud.rows() << " " << Ud.cols() << std::endl;
+
+  std::cout << "Apply dimensionnality reduction" << std::endl;
+  // Xd = Ud.'*M;
+  MatrixMultiplyImageFilterType::Pointer mulUd = MatrixMultiplyImageFilterType::New();
+  mulUd->SetInput(readerImage->GetOutput());
+  mulUd->SetMatrix(Ud);
+  mulUd->UpdateOutputInformation();
+
+  VectorImageType::Pointer Xd = mulUd->GetOutput();
+  std::cout << "Xd = " << Xd << std::endl;
+
+/*
+  std::cout << "Write Xd output" << std::endl;
+  WriterType::Pointer writer = WriterType::New();
+  writer->SetFileName(outputImage);
+  wrtier->SetInput(Xd);
+  writer->Update();
+*/
+
+  // Compute mean(Xd)
+  std::cout << "Compute mean(Xd)" << std::endl;
+  StreamingStatisticsVectorImageFilterType::Pointer statsXd = \
+      StreamingStatisticsVectorImageFilterType::New();
+  statsXd->SetInput(Xd);
+  statsXd->Update();
+  VectorImageType::PixelType Xdmean = statsXd->GetMean();
+  std::cout << "mean(Xd) = " << Xdmean << std::endl;
+  std::cout << "mean(Xd) size = " << Xdmean.Size() << std::endl;
+
+  // Compute Xd ./ repmat( sum( Xd .* repmat(u,[1 N]) ) ,[d 1]);
+  // -> divides each pixel component by the dot product <Xd(i,j), mean(Xd)>
+  std::cout << "Compute projective projection" << std::endl;
+  ProjectiveProjectionImageFilterType::Pointer proj = ProjectiveProjectionImageFilterType::New();
+  proj->SetInput(Xd);
+  proj->SetProjectionDirection(Xdmean);
+
+  std::cout << "Write output" << std::endl;
+  WriterType::Pointer writer = WriterType::New();
+  writer->SetFileName(outputImage);
+  writer->SetInput(proj->GetOutput());
+  writer->Update();
+
+  return EXIT_SUCCESS;
+}