otbPixelWiseBlockMatchingImageFilter.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 __otbPixelWiseBlockMatchingImageFilter_h
#define __otbPixelWiseBlockMatchingImageFilter_h


#include "itkImageToImageFilter.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"
#include "otbImage.h"

namespace otb
{

namespace Functor
{
template <class TInputImage, class TOutputMetricImage>
ITK_EXPORT class SSDBlockMatching
{
public:
  typedef itk::ConstNeighborhoodIterator<TInputImage> ConstNeighborhoodIteratorType;
  typedef typename TOutputMetricImage::ValueType      MetricValueType;

  // Implement the SSD operator
  inline MetricValueType operator()(ConstNeighborhoodIteratorType & a, ConstNeighborhoodIteratorType & b) const
  {
    MetricValueType ssd = 0;
    
    // For some reason, iterators do not work on neighborhoods
    for(unsigned int i = 0; i<a.Size(); ++i)
      {
      ssd += (a.GetPixel(i)-b.GetPixel(i))*(a.GetPixel(i)-b.GetPixel(i));
      }

    return ssd;
  }
};
template <class TInputImage, class TOutputMetricImage>
ITK_EXPORT class NCCBlockMatching
{
public:
  typedef itk::ConstNeighborhoodIterator<TInputImage> ConstNeighborhoodIteratorType;
  typedef typename TOutputMetricImage::ValueType      MetricValueType;

  // Implement the NCC operator
  inline MetricValueType operator()(ConstNeighborhoodIteratorType & a, ConstNeighborhoodIteratorType & b) const
  {
    //MetricValueType meanA(0),meanB(0), sigmaA(0), sigmaB(0), cov(0), ncc(0);
    double  meanA=0.0;
    double  meanB=0.0;
    double  sigmaA=0.0;
    double  sigmaB=0.0;
    double  cov=0.0;
    double  ncc=0.0;
    double  valueA;
    double  valueB;
    double size=a.Size();
    // For some reason, iterators do not work on neighborhoods
    for(unsigned int i = 0; i<size; ++i)
      {
      meanA+=static_cast<double>(a.GetPixel(i));
      meanB+=static_cast<double>(b.GetPixel(i));
      }
    
    // Compute mean
    meanA/=size;
    meanB/=size;

    for(unsigned int i = 0; i<size; ++i)
      {
      valueA=static_cast<double>(a.GetPixel(i));
      valueB=static_cast<double>(b.GetPixel(i));
      cov+=(valueA-meanA)*(valueB-meanB);
      sigmaA+=(valueA-meanA)*(valueA-meanA);
      sigmaB+=(valueB-meanB)*(valueB-meanB);
      }
    
    cov/=size-1;
    sigmaA/=size-1;
    sigmaB/=size-1;
    sigmaA = vcl_sqrt(sigmaA);
    sigmaB = vcl_sqrt(sigmaB);

    if(sigmaA > 1e-20 && sigmaB > 1e-20)
      {
      ncc = vcl_abs(cov)/(sigmaA*sigmaB);
      }
    else
      {
      ncc = 0;
      }

    return static_cast<MetricValueType>(ncc);
  }
};

template <class TInputImage, class TOutputMetricImage>
ITK_EXPORT class LPBlockMatching
{
public:
  typedef itk::ConstNeighborhoodIterator<TInputImage> ConstNeighborhoodIteratorType;
  typedef typename TOutputMetricImage::ValueType      MetricValueType;
  
  LPBlockMatching(): m_P(1)
    {
    }
  
  void SetP(double p)
    {
    if (p > 0.0)
      {
      m_P = p;
      }
    else
      {
      m_P = 1.0;
      }
    }
  
  // Implement the Lp metric
  inline MetricValueType operator()(ConstNeighborhoodIteratorType & a, ConstNeighborhoodIteratorType & b) const
  {
    MetricValueType score(0);
    
    // For some reason, iterators do not work on neighborhoods
    for(unsigned int i = 0; i<a.Size(); ++i)
      {
      score += vcl_pow( vcl_abs(static_cast<double>(a.GetPixel(i)-b.GetPixel(i))) , m_P);
      }
    
    return score;
  }
  
private:
  
  double m_P;
};

} // End Namespace Functor

template <class TInputImage, class TOutputMetricImage, class TOutputDisparityImage = TOutputMetricImage, class TMaskImage = otb::Image<unsigned char>,
          class TBlockMatchingFunctor = Functor::SSDBlockMatching<TInputImage,TOutputMetricImage> >
class ITK_EXPORT PixelWiseBlockMatchingImageFilter :
    public itk::ImageToImageFilter<TInputImage,TOutputDisparityImage>
{
public:
  typedef PixelWiseBlockMatchingImageFilter       Self;
  typedef itk::ImageToImageFilter<TInputImage,
                                  TOutputDisparityImage>    Superclass;
  typedef itk::SmartPointer<Self>                           Pointer;
  typedef itk::SmartPointer<const Self>                     ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(PixelWiseBlockMatchingImageFilter, ImageToImageFilter);

  typedef TInputImage                                       InputImageType;
  typedef TOutputMetricImage                                OutputMetricImageType;
  typedef TOutputDisparityImage                             OutputDisparityImageType;
  typedef TMaskImage                                        InputMaskImageType;
  typedef TBlockMatchingFunctor                             BlockMatchingFunctorType;

  typedef typename InputImageType::SizeType                 SizeType;
  typedef typename InputImageType::IndexType                IndexType;
  typedef typename InputImageType::RegionType               RegionType;
  
  typedef typename TOutputMetricImage::ValueType            MetricValueType;
  
  typedef itk::ConstNeighborhoodIterator<TInputImage>       ConstNeighborhoodIteratorType;

  void SetLeftInput( const TInputImage * image);

  void SetRightInput( const TInputImage * image);

  void SetLeftMaskInput(const TMaskImage * image);
  
  void SetRightMaskInput(const TMaskImage * image);

  const TInputImage * GetLeftInput() const;
  const TInputImage * GetRightInput() const;
  const TMaskImage  * GetLeftMaskInput() const;
  const TMaskImage  * GetRightMaskInput() const;

  const TOutputMetricImage * GetMetricOutput() const;
  TOutputMetricImage * GetMetricOutput();

  const TOutputDisparityImage * GetHorizontalDisparityOutput() const;
  TOutputDisparityImage * GetHorizontalDisparityOutput();

  const TOutputDisparityImage * GetVerticalDisparityOutput() const;
  TOutputDisparityImage * GetVerticalDisparityOutput();


  void SetRadius(unsigned int radius)
  {
    m_Radius.Fill(radius);
  }

  itkSetMacro(Radius, SizeType);
  itkGetConstReferenceMacro(Radius, SizeType);
  
  /*** Set/Get the minimum disparity to explore */
  itkSetMacro(MinimumHorizontalDisparity,int);
  itkGetConstReferenceMacro(MinimumHorizontalDisparity,int);

  /*** Set/Get the maximum disparity to explore */
  itkSetMacro(MaximumHorizontalDisparity,int);
  itkGetConstReferenceMacro(MaximumHorizontalDisparity,int);

  /*** Set/Get the minimum disparity to explore */
  itkSetMacro(MinimumVerticalDisparity,int);
  itkGetConstReferenceMacro(MinimumVerticalDisparity,int);

  /*** Set/Get the maximum disparity to explore */
  itkSetMacro(MaximumVerticalDisparity,int);
  itkGetConstReferenceMacro(MaximumVerticalDisparity,int);

  itkSetMacro(Minimize, bool);
  itkGetConstReferenceMacro(Minimize,bool);
  itkBooleanMacro(Minimize);
  
  itkSetMacro(ExplorationRadius, SizeType);
  itkGetConstReferenceMacro(ExplorationRadius, SizeType);
  
  itkSetMacro(InitHorizontalDisparity,int);
  itkGetConstReferenceMacro(InitHorizontalDisparity,int);
  
  itkSetMacro(InitVerticalDisparity,int);
  itkGetConstReferenceMacro(InitVerticalDisparity,int);
  
  BlockMatchingFunctorType &  GetFunctor()
  {
    return m_Functor;
  }

  const BlockMatchingFunctorType &  GetFunctor() const
  {
    return m_Functor;
  }

  void SetHorizontalDisparityInput( const TOutputDisparityImage * hfield);

  void SetVerticalDisparityInput( const TOutputDisparityImage * vfield);
  
  const TOutputDisparityImage * GetHorizontalDisparityInput() const;
  const TOutputDisparityImage * GetVerticalDisparityInput() const;

protected:
  PixelWiseBlockMatchingImageFilter();

  virtual ~PixelWiseBlockMatchingImageFilter();

  virtual void GenerateInputRequestedRegion();

  virtual void BeforeThreadedGenerateData();
  
  virtual void ThreadedGenerateData(const RegionType & outputRegionForThread, int threadId);
  
private:
  PixelWiseBlockMatchingImageFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemeFnted
 
  SizeType                      m_Radius;

  int                           m_MinimumHorizontalDisparity;

  int                           m_MaximumHorizontalDisparity;

  int                           m_MinimumVerticalDisparity;

  int                           m_MaximumVerticalDisparity;

  bool                          m_Minimize;
  
  SizeType                      m_ExplorationRadius;

  BlockMatchingFunctorType      m_Functor;
  
  int                           m_InitHorizontalDisparity;
  
  int                           m_InitVerticalDisparity;
};
} // end namespace otb

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

#endif