otbPixelWiseBlockMatchingImageFilter.h

Go to the documentation of this file.

/*
 * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
 
#ifndef otbPixelWiseBlockMatchingImageFilter_h
#define otbPixelWiseBlockMatchingImageFilter_h
 
 
#include "itkImageToImageFilter.h"
#include "itkConstNeighborhoodIterator.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 SSDDivMeanBlockMatching
{
public:
  typedef itk::ConstNeighborhoodIterator<TInputImage> ConstNeighborhoodIteratorType;
  typedef typename TOutputMetricImage::ValueType      MetricValueType;
 
  // Implement the SSD DivMean operator
  inline MetricValueType operator()(ConstNeighborhoodIteratorType& a, ConstNeighborhoodIteratorType& b) const
  {
    MetricValueType ssd   = 0;
    MetricValueType meana = 0;
    MetricValueType meanb = 0;
 
    for (unsigned int i = 0; i < a.Size(); ++i)
    {
      meana += a.GetPixel(i);
      meanb += b.GetPixel(i);
    }
    meana /= a.Size();
    meanb /= b.Size();
 
    // For some reason, iterators do not work on neighborhoods
    for (unsigned int i = 0; i < a.Size(); ++i)
    {
      ssd += (a.GetPixel(i) / meana - b.GetPixel(i) / meanb) * (a.GetPixel(i) / meana - b.GetPixel(i) / meanb);
    }
 
    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 = std::sqrt(sigmaA);
    sigmaB = std::sqrt(sigmaB);
 
    if (sigmaA > 1e-20 && sigmaB > 1e-20)
    {
      ncc = std::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 += std::pow(std::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 InputImageType::SpacingType SpacingType;
  typedef typename InputImageType::PointType   PointType;
 
  typedef typename TOutputMetricImage::ValueType MetricValueType;
 
  typedef typename OutputDisparityImageType::PixelType DisparityPixelType;
 
  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);
 
  itkSetMacro(MinimumHorizontalDisparity, int);
  itkGetConstReferenceMacro(MinimumHorizontalDisparity, int);
 
  itkSetMacro(MaximumHorizontalDisparity, int);
  itkGetConstReferenceMacro(MaximumHorizontalDisparity, int);
 
  itkSetMacro(MinimumVerticalDisparity, int);
  itkGetConstReferenceMacro(MinimumVerticalDisparity, int);
 
  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;
 
  itkSetMacro(Step, unsigned int);
  itkGetMacro(Step, unsigned int);
 
  itkSetMacro(GridIndex, IndexType);
  itkGetConstReferenceMacro(GridIndex, IndexType);
 
  static RegionType ConvertFullToSubsampledRegion(RegionType full, unsigned int step, IndexType index);
 
  static RegionType ConvertSubsampledToFullRegion(RegionType sub, unsigned int step, IndexType index);
 
protected:
  PixelWiseBlockMatchingImageFilter();
 
  ~PixelWiseBlockMatchingImageFilter() override;
 
  void GenerateOutputInformation() override;
 
  void GenerateInputRequestedRegion() override;
 
  void BeforeThreadedGenerateData() override;
 
  void ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId) override;
 
private:
  PixelWiseBlockMatchingImageFilter(const Self&) = delete;
  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;
 
  unsigned int m_Step;
 
  IndexType m_GridIndex;
};
} // end namespace otb
 
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbPixelWiseBlockMatchingImageFilter.hxx"
#endif
 
#endif