otbKullbackLeiblerDistanceImageFilter.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.

  Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
  See GETCopyright.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 __otbKullbackLeiblerDistanceImageFilter_h
#define __otbKullbackLeiblerDistanceImageFilter_h

#include "itkVariableLengthVector.h"
#include "otbBinaryFunctorNeighborhoodImageFilter.h"

namespace otb
{

template <class TInput>
class CumulantsForEdgeworth
{
public:
  CumulantsForEdgeworth (const TInput& input);
  CumulantsForEdgeworth (const itk::Image<typename TInput::ImageType::PixelType, 1> * input);
  virtual ~CumulantsForEdgeworth () {}

  template <class TInput2>
  double Divergence(const CumulantsForEdgeworth<TInput2>& cumulants);

  inline double  GetMean()  const
  {
    return this->fMean;
  }
  inline double  GetVariance()  const
  {
    return this->fVariance;
  }
  inline double  GetSkewness()  const
  {
    return this->fSkewness;
  }
  inline double  GetKurtosis()  const
  {
    return this->fKurtosis;
  }
  inline bool IsDataAvailable() const
  {
    return this->fDataAvailable;
  }

protected:

  void  MakeSumAndMoments(const TInput& input);
  void  MakeSumAndMoments(const itk::Image<typename TInput::ImageType::PixelType, 1> * input);
  void MakeCumulants();

  double fSum0, fSum1, fSum2, fSum3, fSum4;
  double fMu1, fMu2, fMu3, fMu4;
  double fMean, fVariance, fSkewness, fKurtosis;

  bool fDataAvailable;
};

namespace Functor
{
template <class TInput1, class TInput2, class TOutput>
class KullbackLeiblerDistance
{
public:
  KullbackLeiblerDistance () {}
  virtual ~KullbackLeiblerDistance () {}
  TOutput operator ()(const TInput1& it1, const TInput2& it2)
  {
    CumulantsForEdgeworth<TInput1> cum1 (it1);
    if (!cum1.IsDataAvailable()) return static_cast<TOutput>(0.);

    CumulantsForEdgeworth<TInput2> cum2 (it2);
    if (!cum2.IsDataAvailable()) return static_cast<TOutput>(0.);

    return static_cast<TOutput> (cum1.Divergence(cum2)
                                 + cum2.Divergence(cum1));
  }
};

} // Functor

template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT KullbackLeiblerDistanceImageFilter :
  public otb::BinaryFunctorNeighborhoodImageFilter<
      TInputImage1, TInputImage2, TOutputImage,
      Functor::KullbackLeiblerDistance<
          typename itk::ConstNeighborhoodIterator<TInputImage1>,
          typename itk::ConstNeighborhoodIterator<TInputImage2>,
          typename TOutputImage::PixelType> >
{
public:
  typedef KullbackLeiblerDistanceImageFilter Self;
  typedef otb::BinaryFunctorNeighborhoodImageFilter<
      TInputImage1, TInputImage2, TOutputImage,
      Functor::KullbackLeiblerDistance<
          typename itk::ConstNeighborhoodIterator<TInputImage1>,
          typename itk::ConstNeighborhoodIterator<TInputImage2>,
          typename TOutputImage::PixelType>
      >  Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

  itkNewMacro(Self);

protected:
  KullbackLeiblerDistanceImageFilter() {}
  virtual ~KullbackLeiblerDistanceImageFilter() {}

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

};

} // namespace otb

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

#endif