otbMeanShiftSmoothingImageFilter.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 otbMeanShiftSmoothingImageFilter_h
#define otbMeanShiftSmoothingImageFilter_h
 
#include "otbImage.h"
#include "otbVectorImage.h"
#include "itkImageToImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include <algorithm>
 
 
namespace otb
{
namespace Meanshift
{
 
template <typename T>
inline T simple_pow(T const& v, unsigned int p)
{
  T res = 1;
  for (unsigned int i = 0; i != p; ++i)
  {
    res *= v;
  }
  return res;
}
 
template <class TInputImage, class TOutputJointImage>
class SpatialRangeJointDomainTransform
{
public:
  typedef double RealType;
 
  SpatialRangeJointDomainTransform() : m_ImageDimension(0), m_NumberOfComponentsPerPixel(0), m_OutputSize(0)
  {
  }
 
  typename TOutputJointImage::PixelType operator()(const typename TInputImage::PixelType& inputPixel, const typename TInputImage::IndexType& index) const
  {
    typename TOutputJointImage::PixelType jointPixel(m_ImageDimension + m_NumberOfComponentsPerPixel);
 
    for (unsigned int comp = 0; comp < m_ImageDimension; comp++)
    {
      jointPixel[comp] = index[comp] + m_GlobalShift[comp];
    }
    for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++)
    {
      jointPixel[m_ImageDimension + comp] = inputPixel[comp];
    }
    return jointPixel;
  }
 
  void Initialize(unsigned int _ImageDimension, unsigned int numberOfComponentsPerPixel_, typename TInputImage::IndexType globalShift_)
  {
    m_ImageDimension             = _ImageDimension;
    m_NumberOfComponentsPerPixel = numberOfComponentsPerPixel_;
    m_OutputSize                 = m_ImageDimension + m_NumberOfComponentsPerPixel;
    m_GlobalShift                = globalShift_;
  }
 
  unsigned int GetOutputSize() const
  {
    return m_OutputSize;
  }
 
private:
  unsigned int                    m_ImageDimension;
  unsigned int                    m_NumberOfComponentsPerPixel;
  unsigned int                    m_OutputSize;
  typename TInputImage::IndexType m_GlobalShift;
};
 
class KernelUniform
{
public:
  typedef double RealType;
 
  // KernelUniform() {}
  // ~KernelUniform() {}
 
  RealType operator()(RealType x) const
  {
    return (x <= 1) ? 1.0 : 0.0;
  }
 
  RealType GetRadius(RealType bandwidth) const
  {
    return bandwidth;
  }
};
 
class KernelGaussian
{
public:
  typedef double RealType;
 
  // KernelGaussian() {}
  // ~KernelGaussian() {}
 
  RealType operator()(RealType x) const
  {
    return std::exp(-0.5 * x);
  }
 
  RealType GetRadius(RealType bandwidth) const
  {
    return 3.0 * bandwidth;
  }
};
 
template <typename TImage>
class FastImageRegionConstIterator : public itk::ImageRegionConstIterator<TImage>
{
public:
  typedef FastImageRegionConstIterator<TImage>  Self;
  typedef itk::ImageRegionConstIterator<TImage> Superclass;
 
  typedef typename Superclass::ImageType  ImageType;
  typedef typename Superclass::RegionType RegionType;
 
  typedef typename TImage::PixelType         PixelType;
  typedef typename TImage::InternalPixelType InternalPixelType;
 
  itkTypeMacro(FastImageRegionConstIterator, ImageRegionConstIterator);
  ;
 
  FastImageRegionConstIterator() : Superclass()
  {
  }
  FastImageRegionConstIterator(const ImageType* ptr, const RegionType& region) : Superclass(ptr, region)
  {
    m_NumberOfComponentsPerPixel = ptr->GetNumberOfComponentsPerPixel();
  }
 
  const InternalPixelType* GetPixelPointer() const
  {
    return this->m_Buffer + (this->m_Offset * m_NumberOfComponentsPerPixel);
  }
 
private:
  unsigned int m_NumberOfComponentsPerPixel;
};
 
#if 0 // disable bucket mode
 
template<class TImage>
class BucketImage
{
public:
  typedef TImage ImageType;
  typedef typename ImageType::ConstPointer ImageConstPointerType;
  typedef typename ImageType::PixelType PixelType;
  typedef typename ImageType::InternalPixelType InternalPixelType;
  typedef typename ImageType::RegionType RegionType;
  typedef typename ImageType::IndexType IndexType;
 
  typedef double RealType;
 
  static const unsigned int ImageDimension = ImageType::ImageDimension;
 
  typedef std::vector<typename ImageType::SizeType::SizeValueType> BucketImageSizeType;
  typedef std::vector<typename ImageType::IndexType::IndexValueType> BucketImageIndexType;
  //typedef std::vector<long> BucketImageIndexType;
 
 
  typedef const typename ImageType::InternalPixelType * ImageDataPointerType;
  typedef std::vector<ImageDataPointerType> BucketType;
  typedef std::vector<BucketType> BucketListType;
 
  BucketImage()
  {
  }
 
  BucketImage(ImageConstPointerType image, const RegionType & region, RealType spatialRadius, RealType rangeRadius,
              unsigned int spectralCoordinate) :
    m_Image(image), m_Region(region), m_SpatialRadius(spatialRadius), m_RangeRadius(rangeRadius),
        m_SpectralCoordinate(spectralCoordinate)
  {
 
    // Find max and min of the used spectral band
    itk::ImageRegionConstIterator<ImageType> inputIt(m_Image, m_Region);
    inputIt.GoToBegin();
    InternalPixelType minValue = inputIt.Get()[spectralCoordinate];
    InternalPixelType maxValue = minValue;
    ++inputIt;
    while (!inputIt.IsAtEnd())
      {
      const PixelType &p = inputIt.Get();
      minValue = std::min(minValue, p[m_SpectralCoordinate]);
      maxValue = std::max(maxValue, p[m_SpectralCoordinate]);
      ++inputIt;
      }
 
    m_MinValue = minValue;
    m_MaxValue = maxValue;
 
    // Compute bucket image dimensions. Note: empty buckets are at each border
    // to simplify image border issues
    m_DimensionVector.resize(ImageDimension + 1); // NB: pays for a 0-innit
    for (unsigned int dim = 0; dim < ImageDimension; ++dim)
      {
      m_DimensionVector[dim] = m_Region.GetSize()[dim] / m_SpatialRadius + 3;
      }
    m_DimensionVector[ImageDimension] = (unsigned int) ((maxValue - minValue) / m_RangeRadius) + 3;
 
    unsigned int numBuckets = m_DimensionVector[0];
    for (unsigned int dim = 1; dim <= ImageDimension; ++dim)
      numBuckets *= m_DimensionVector[dim];
 
    m_BucketList.resize(numBuckets);
    // Build buckets
    itk::ImageRegionConstIteratorWithIndex<ImageType> it(m_Image, m_Region);
    it.GoToBegin();
    // this iterator is only used to get the pixel data pointer
    FastImageRegionConstIterator<ImageType> fastIt(m_Image, m_Region);
    fastIt.GoToBegin();
 
    while (!it.IsAtEnd())
      {
      const IndexType & index = it.GetIndex();
      const PixelType & pixel = it.Get();
 
      // Find which bucket this pixel belongs to
      const BucketImageIndexType bucketIndex = GetBucketIndex(pixel, index);
 
      unsigned int bucketListIndex = BucketIndexToBucketListIndex(bucketIndex);
      assert(bucketListIndex < numBuckets);
      m_BucketList[bucketListIndex].push_back(fastIt.GetPixelPointer());
      ++it;
      ++fastIt;
      }
 
    // Prepare neighborhood offset vector
    // BucketImageIndexType zeroOffsetIndex(ImageDimension+1);
    std::vector<BucketImageIndexType> neighborsIndexList;
    neighborsIndexList.reserve(simple_pow(3, ImageDimension + 1));
    neighborsIndexList.resize(1, BucketImageIndexType(ImageDimension + 1)); // zeroOffsetIndex
    // neighborsIndexList.push_back(zeroOffsetIndex);
    for (unsigned dim = 0; dim <= ImageDimension; ++dim)
      {
      // take all neighbors already in the list and add their direct neighbor
      // along the current dim
      const unsigned int curSize = neighborsIndexList.size();
      for (unsigned int i = 0; i < curSize; ++i)
        {
        BucketImageIndexType index = neighborsIndexList[i];
        index[dim]--;
        neighborsIndexList.push_back(index);
        index[dim] += 2;
        neighborsIndexList.push_back(index);
        }
      }
    // Convert all neighbors n-dimensional indices to bucket list 1D indices
    const unsigned int neighborhoodOffsetVectorSize = neighborsIndexList.size();
    m_NeighborhoodOffsetVector.reserve(neighborhoodOffsetVectorSize);
    for (unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i)
      {
      const int listIndex = BucketIndexToBucketListIndex(neighborsIndexList[i]);
      m_NeighborhoodOffsetVector.push_back(listIndex);
      }
  }
 
  ~BucketImage()
  {
  }
 
  BucketImageIndexType GetBucketIndex(const PixelType & pixel, const IndexType & index)
  {
    BucketImageIndexType bucketIndex(ImageDimension + 1);
    for (unsigned int dim = 0; dim < ImageDimension; ++dim)
      {
      bucketIndex[dim] = (index[dim] - m_Region.GetIndex()[dim]) / m_SpatialRadius + 1;
      }
    bucketIndex[ImageDimension] = (pixel[m_SpectralCoordinate] - m_MinValue) / m_RangeRadius + 1;
    return bucketIndex;
  }
 
  int BucketIndexToBucketListIndex(const BucketImageIndexType & bucketIndex) const
  {
    int bucketListIndex = bucketIndex[0];
    for (unsigned int dim = 1; dim <= ImageDimension; ++dim)
      {
      bucketListIndex = bucketListIndex * m_DimensionVector[dim] + bucketIndex[dim];
      }
    return bucketListIndex;
  }
 
  std::vector<unsigned int> GetNeighborhoodBucketListIndices(int bucketIndex) const
  {
    const unsigned int neighborhoodOffsetVectorSize = m_NeighborhoodOffsetVector.size();
    std::vector<unsigned int> indices(neighborhoodOffsetVectorSize);
 
    for (unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i)
      {
      indices[i] = bucketIndex + m_NeighborhoodOffsetVector[i];
      }
    return indices;
  }
 
  /* Returns the list of pixels (actually pointer to pixel data) contained in a bucket */
  const BucketType & GetBucket(unsigned int index) const
  {
    return m_BucketList[index];
  }
 
  unsigned int GetNumberOfNeighborBuckets() const
  {
    return m_NeighborhoodOffsetVector.size();
  }
 
private:
  ImageConstPointerType m_Image;
 
  RegionType m_Region;
 
  RealType m_SpatialRadius;
 
  RealType m_RangeRadius;
 
  unsigned int m_SpectralCoordinate;
 
  InternalPixelType m_MinValue;
  InternalPixelType m_MaxValue;
 
  BucketListType m_BucketList;
 
  BucketImageSizeType m_DimensionVector;
 
  std::vector<int> m_NeighborhoodOffsetVector;
};
#endif
 
} // end namespace Meanshift
 
template <class TInputImage, class TOutputImage, class TKernel = Meanshift::KernelUniform,
          class TOutputIterationImage = otb::Image<unsigned int, TInputImage::ImageDimension>>
class ITK_EXPORT MeanShiftSmoothingImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>
{
public:
  typedef MeanShiftSmoothingImageFilter Self;
  typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
  typedef double                        RealType;
 
  itkTypeMacro(MeanShiftSmoothingImageFilter, ImageToImageFilter);
  ;
  itkNewMacro(Self);
  ;
 
  typedef TInputImage                             InputImageType;
  typedef typename InputImageType::Pointer        InputImagePointerType;
  typedef typename InputImageType::PixelType      InputPixelType;
  typedef typename InputImageType::IndexType      InputIndexType;
  typedef typename InputImageType::SizeType       InputSizeType;
  typedef typename InputImageType::IndexValueType InputIndexValueType;
  typedef typename InputImageType::PointType      PointType;
  typedef typename InputImageType::RegionType     RegionType;
  typedef typename InputImageType::SizeType       SizeType;
 
  typedef TOutputImage                         OutputImageType;
  typedef typename OutputImageType::Pointer    OutputImagePointerType;
  typedef typename OutputImageType::PixelType  OutputPixelType;
  typedef typename OutputImageType::RegionType OutputRegionType;
 
  typedef TOutputIterationImage OutputIterationImageType;
 
  typedef unsigned long LabelType;
  typedef otb::Image<LabelType, InputImageType::ImageDimension> OutputLabelImageType;
 
  typedef otb::VectorImage<RealType, InputImageType::ImageDimension> OutputSpatialImageType;
  typedef typename OutputSpatialImageType::Pointer   OutputSpatialImagePointerType;
  typedef typename OutputSpatialImageType::PixelType OutputSpatialPixelType;
 
  typedef TKernel KernelType;
 
  itkStaticConstMacro(ImageDimension, unsigned int, InputImageType::ImageDimension);
 
  typedef itk::VariableLengthVector<RealType> RealVector;
  typedef otb::VectorImage<RealType, InputImageType::ImageDimension> RealVectorImageType;
  typedef otb::Image<unsigned short, InputImageType::ImageDimension> ModeTableImageType;
 
  itkSetMacro(SpatialBandwidth, RealType);
  itkGetConstReferenceMacro(SpatialBandwidth, RealType);
 
  itkSetMacro(RangeBandwidth, RealType);
  itkGetConstReferenceMacro(RangeBandwidth, RealType);
 
  itkSetMacro(RangeBandwidthRamp, RealType);
  itkGetConstReferenceMacro(RangeBandwidthRamp, RealType);
 
  itkGetConstReferenceMacro(MaxIterationNumber, unsigned int);
  itkSetMacro(MaxIterationNumber, unsigned int);
 
  itkGetConstReferenceMacro(Threshold, double);
  itkSetMacro(Threshold, double);
 
  itkSetMacro(ModeSearch, bool);
  itkGetConstReferenceMacro(ModeSearch, bool);
 
#if 0
 
  itkSetMacro(BucketOptimization, bool);
  itkGetConstReferenceMacro(BucketOptimization, bool);
#endif
 
 
  itkSetMacro(GlobalShift, InputIndexType);
 
  const OutputSpatialImageType* GetSpatialOutput() const;
 
  const OutputImageType* GetRangeOutput() const;
 
  const OutputIterationImageType* GetIterationOutput() const;
 
  const OutputLabelImageType* GetLabelOutput() const;
 
  OutputSpatialImageType* GetSpatialOutput();
 
  OutputImageType* GetRangeOutput();
 
  OutputIterationImageType* GetIterationOutput();
 
  OutputLabelImageType* GetLabelOutput();
 
protected:
  void GenerateOutputInformation(void) override;
 
  void GenerateInputRequestedRegion() override;
 
  void BeforeThreadedGenerateData() override;
 
  void ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId) override;
 
  void AfterThreadedGenerateData() override;
 
  void AllocateOutputs() override;
 
  MeanShiftSmoothingImageFilter();
 
  ~MeanShiftSmoothingImageFilter() override;
 
  void PrintSelf(std::ostream& os, itk::Indent indent) const override;
 
  virtual void CalculateMeanShiftVector(const typename RealVectorImageType::Pointer inputImagePtr, const RealVector& jointPixel,
                                        const OutputRegionType& outputRegion, const RealVector& bandwidth, RealVector& meanShiftVector);
#if 0
  virtual void CalculateMeanShiftVectorBucket(const RealVector& jointPixel, RealVector& meanShiftVector);
#endif
 
private:
  MeanShiftSmoothingImageFilter(const Self&) = delete;
  void operator=(const Self&) = delete;
 
  RealType m_RangeBandwidth;
 
  RealType m_RangeBandwidthRamp;
 
  RealType m_SpatialBandwidth;
 
  InputSizeType m_SpatialRadius;
 
  double m_Threshold;
 
  unsigned int m_MaxIterationNumber;
 
  KernelType m_Kernel;
 
  unsigned int m_NumberOfComponentsPerPixel;
 
  typename RealVectorImageType::Pointer m_JointImage;
 
  typename ModeTableImageType::Pointer m_ModeTable;
 
  bool m_ModeSearch;
 
#if 0
 
  bool m_BucketOptimization;
#endif
 
  itk::VariableLengthVector<LabelType> m_NumLabels;
 
  unsigned int m_ThreadIdNumberOfBits;
 
#if 0
  typedef Meanshift::BucketImage<RealVectorImageType> BucketImageType;
  BucketImageType m_BucketImage;
#endif
 
  InputIndexType m_GlobalShift;
};
 
} // end namespace otb
 
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbMeanShiftSmoothingImageFilter.hxx"
#endif
 
#endif