otbMDMDNMFImageFilter.h

Go to the documentation of this file.

/*=========================================================================

  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 __otbMDMDNMFImageFilter_h
#define __otbMDMDNMFImageFilter_h

#include "itkImageToImageFilter.h"
#include "itkVariableSizeMatrix.h"
#include "otbVectorImageToMatrixImageFilter.h"

#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector.h"

namespace otb
{
template <class TInputImage, class TOutputImage>
class ITK_EXPORT MDMDNMFImageFilter:
    public itk::ImageToImageFilter<TInputImage, TOutputImage>
{
public:

  typedef MDMDNMFImageFilter                                    Self;
  typedef itk::ImageToImageFilter< TInputImage, TOutputImage >  Superclass;
  typedef itk::SmartPointer<Self>                               Pointer;
  typedef itk::SmartPointer<const Self>                         ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(MDMDNMFImageFilter, ImageToImageFilter);

  typedef          TInputImage                   InputImageType;
  typedef typename InputImageType::Pointer       InputPointerType;
  typedef typename InputImageType::ConstPointer  InputConstPointerType;
  typedef typename InputImageType::IndexType     InputIndexType;
  typedef typename InputImageType::SizeType      InputSizeType;

  typedef          TOutputImage                  OutputImageType;
  typedef typename OutputImageType::Pointer      OutputPointerType;
  typedef typename OutputImageType::IndexType    OutputIndexType;
  typedef typename OutputImageType::OffsetType   OutputOffsetType;
  typedef typename OutputImageType::SizeType     OutputSizeType;
  typedef typename OutputImageType::RegionType   OutputImageRegionType;

  /* Matrix type */
  typedef double                                                       PrecisionType;
  typedef vnl_matrix<PrecisionType>                                    MatrixType;
  typedef typename otb::VectorImageToMatrixImageFilter<InputImageType> VectorImageToMatrixImageFilterType;
  
  /* vector type */
  typedef vnl_vector<double> VectorType;

  MatrixType GetEndmembersMatrix() const
  {
    return m_Endmembers;
  }
  void SetEndmembersMatrix(const MatrixType& m)
  {
    m_Endmembers = m;
  }

  itkSetMacro(MaxIter, unsigned int);
  itkGetMacro(MaxIter, unsigned int);
  
  itkSetMacro(CritStopValue, double);
  itkGetMacro(CritStopValue, double);
  
  itkSetMacro(Delt, double);
  itkGetMacro(Delt, double);
  
  itkSetMacro(LambdD, double);
  itkGetMacro(LambdD, double);
  
  itkSetMacro(LambdS, double);
  itkGetMacro(LambdS, double);

protected:
  MDMDNMFImageFilter();
  virtual ~MDMDNMFImageFilter() {}
  void PrintSelf(std::ostream& os, itk::Indent indent) const;

  virtual void GenerateInputRequestedRegion();

  virtual void GenerateOutputInformation();

  virtual void GenerateData();


private:
  MDMDNMFImageFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

  static void AddOneRowOfOnes(const MatrixType & m, MatrixType & M);
  
  static double Criterion(const MatrixType & X,
                   const MatrixType & A,
                   const MatrixType & S,
                   const double &delt,
                   const double &lambdS,
                   const double &lambdD);

  static void EvalGradS(const MatrixType &X,
                 const MatrixType &A,
                 const MatrixType &S,
                 const double &delt,
                 const double &lambdS,
                 MatrixType & gradS);

  static void EvalGradA(const MatrixType &X,
                 const MatrixType &A,
                 const MatrixType &S,
                 const double &delt,
                 const double &lambdD,
                 MatrixType &gradA);

  static double Call(const MatrixType & variMat,
              const MatrixType & fixedMat,
              const MatrixType & X,
              const double & delt,
              const double & lambdS,
              const double & lambdD, const bool isDirectEvalDirection);
  
  static void ProjGradOneStep(
    const MatrixType & X,
    const MatrixType & fixedMat,
    const MatrixType & gradVariMat,
    const double & sig,
    const double & betinit,
    const double & delt,
    const double & lambdS,
    const double & lambdD,
    MatrixType & variMat,
    double & alph, const bool isDirectEvalDirection );

  static void SetNegativeCoefficientsToZero(MatrixType & M);
  
  static MatrixType TermByTermMatrixProduct(const MatrixType & M1, const MatrixType & M2);

  static double SumMatrixElements(const MatrixType & M);
  
  static bool ArmijoTest(const double & sig,
                  const MatrixType variMat,
                  const MatrixType & newVariMat,
                  const double & evalf,
                  const double & newEvalf,
                  const MatrixType & gradVariMat,
                  const double & alph);
  
  // attributes
  MatrixType   m_Endmembers;
  unsigned int m_MaxIter;
  double       m_CritStopValue;

  // Optimized function parameters
  double m_Delt;
  double m_LambdD;
  double m_LambdS;
};

} // end namespace otb

#ifndef OTB_MANUAL_INSTANTIATION
#include "otbMDMDNMFImageFilter.txx"
#endif

#endif