itkPatchBasedDenoisingImageFilter.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  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.txt
 *
 *  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 __itkPatchBasedDenoisingImageFilter_h
#define __itkPatchBasedDenoisingImageFilter_h

#include "itkPatchBasedDenoisingBaseImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkVector.h"
#include "itkVectorImage.h"
#include "itkRGBPixel.h"
#include "itkRGBAPixel.h"
#include "itkDiffusionTensor3D.h"
#include "itkRegionConstrainedSubsampler.h"
#include "itkEnableIf.h"
#include "itkIsSame.h"

namespace itk
{

template <class TInputImage, class TOutputImage>
class ITK_EXPORT PatchBasedDenoisingImageFilter :
public PatchBasedDenoisingBaseImageFilter<TInputImage, TOutputImage>
{
 public:
  typedef PatchBasedDenoisingImageFilter                                Self;
  typedef PatchBasedDenoisingBaseImageFilter<TInputImage, TOutputImage> Superclass;
  typedef SmartPointer<Self>                                            Pointer;
  typedef SmartPointer<const Self>                                      ConstPointer;
  typedef typename Superclass::OutputImagePointer                       OutputImagePointer;

  itkNewMacro(Self);

  itkTypeMacro(PatchBasedDenoisingImageFilter, ImageToImageFilter);

  typedef typename Superclass::InputImageType   InputImageType;
  typedef typename Superclass::OutputImageType  OutputImageType;

  itkStaticConstMacro(ImageDimension, unsigned int,
                      Superclass::ImageDimension);

  typedef typename InputImageType::RegionType           InputImageRegionType;

  typedef ImageRegionIterator<OutputImageType>          OutputImageRegionIteratorType;
  typedef ImageRegionConstIterator<InputImageType>      InputImageRegionConstIteratorType;

  typedef typename Superclass::PixelType      PixelType;
  typedef typename Superclass::PixelValueType PixelValueType;

  typedef typename NumericTraits< PixelType >::RealType      RealType;
  typedef typename NumericTraits< PixelValueType >::RealType RealValueType;
  typedef Array<PixelValueType>                              PixelArrayType;
  typedef Array<RealValueType>                               RealArrayType;
  typedef Array<unsigned short>                              ShortArrayType;

  typedef typename Superclass::ListAdaptorType         ListAdaptorType;
  typedef typename Superclass::PatchRadiusType         PatchRadiusType;
  typedef typename Superclass::InputImagePatchIterator InputImagePatchIterator;
  typedef ListAdaptorType                              PatchSampleType;
  typedef typename Superclass::PatchWeightsType        PatchWeightsType;

  typedef itk::Statistics::RegionConstrainedSubsampler<
    PatchSampleType, InputImageRegionType >            BaseSamplerType;
  typedef typename BaseSamplerType::Pointer            BaseSamplerPointer;
  typedef typename BaseSamplerType::InstanceIdentifier InstanceIdentifier;

  struct ThreadDataStruct
  {
    ShortArrayType     validDerivatives;
    RealArrayType      entropyFirstDerivative;
    RealArrayType      entropySecondDerivative;
    ShortArrayType     validNorms;
    RealArrayType      minNorm;
    RealArrayType      maxNorm;
    BaseSamplerPointer sampler;
  };


  itkSetMacro(UseSmoothDiscPatchWeights, bool);
  itkBooleanMacro(UseSmoothDiscPatchWeights);
  itkGetConstMacro(UseSmoothDiscPatchWeights, bool);


  itkSetMacro(GaussianKernelSigma, RealArrayType);
  itkGetConstMacro(GaussianKernelSigma, RealArrayType);

  itkSetClampMacro(FractionPixelsForSigmaUpdate, double, 0.01, 1.0);
  itkGetConstReferenceMacro(FractionPixelsForSigmaUpdate, double);

  itkSetMacro(ComputeConditionalDerivatives, bool);
  itkBooleanMacro(ComputeConditionalDerivatives);
  itkGetConstMacro(ComputeConditionalDerivatives, bool);

  itkStaticConstMacro(MaxSigmaUpdateIterations, unsigned int,
                      20);

  itkSetClampMacro(SigmaMultiplicationFactor, double, 0.01, 100);
  itkGetConstReferenceMacro(SigmaMultiplicationFactor, double);

  void SetNoiseSigma(const RealType& sigma);
  itkGetConstMacro(NoiseSigma, RealType);


  itkSetObjectMacro(Sampler, BaseSamplerType);
  itkGetObjectMacro(Sampler, BaseSamplerType);

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

  virtual void AllocateUpdateBuffer();

  virtual void CopyInputToOutput();

  virtual void GenerateInputRequestedRegion();

  template< typename T>
    typename EnableIfC<
      IsSame<T, typename NumericTraits<T>::ValueType>::Value,
      T >::Type
  GetComponent(const T pix,
               unsigned int itkNotUsed( idx ) ) const
  {
    // The enable if idiom is used to overload this method for both
    // scalars and multi-component types. By exploiting that
    // NumericTraits' ValueType typedef (defines the per-element type
    // for multi-component types ) is different then the parameterize
    // type, the bracket operator is used only for multi-component
    // types.
    return pix;
  }
  template< typename T>
    typename DisableIfC<
      IsSame<T, typename NumericTraits<T>::ValueType>::Value,
       typename NumericTraits<T>::ValueType>::Type
  GetComponent(const T& pix,
               unsigned int idx ) const
  {
    return pix[idx];
  }

  template< typename T >
    void
  SetComponent( T &pix,
                unsigned int itkNotUsed( idx ),
                typename EnableIfC< IsSame<T, typename NumericTraits<T>::ValueType>::Value, RealValueType>::Type val) const
  {
    pix = val;
  }
  template< typename T >
    void
  SetComponent( T &pix,
                unsigned int idx,
                typename DisableIfC< IsSame<T, typename NumericTraits<T>::ValueType>::Value, RealValueType>::Type val) const
  {
    pix[idx] =  val;
  }


  void ComputeMinMax(const Image< DiffusionTensor3D<PixelValueType> , ImageDimension>* img)
  {
    if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
    {
      DispatchedRiemannianMinMax(img);
    }
    else
    {
      DispatchedArrayMinMax(img);
    }
  }

  template< typename TImageType>
    typename EnableIfC<
      IsSame<typename TImageType::PixelType, typename NumericTraits<typename TImageType::PixelType>::ValueType>::Value
      >::Type
    ComputeMinMax(const TImageType* img)
  {
    DispatchedMinMax(img);
  }
 template< typename TImageType>
    typename DisableIfC<
      IsSame<typename TImageType::PixelType, typename NumericTraits<typename TImageType::PixelType>::ValueType>::Value
      >::Type
   ComputeMinMax(const TImageType* img)
  {
    DispatchedArrayMinMax(img);
  }

  void ComputeDifferenceAndWeightedSquaredNorm(const DiffusionTensor3D<PixelValueType>& a,
                                               const DiffusionTensor3D<PixelValueType>& b,
                                               const RealArrayType& weight,
                                               RealType& diff, RealArrayType& norm)
  {
    if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
    {
      ComputeLogMapAndWeightedSquaredGeodesicDifference(a, b, weight, diff, norm);
    }
    else
    {
      ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(a, b, weight, diff, norm);
    }
  }
  template <class PixelT>
  void ComputeDifferenceAndWeightedSquaredNorm(const PixelT& a,
                                               const PixelT& b,
                                               const RealArrayType& weight,
                                               RealType& diff, RealArrayType& norm)
  {
    ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(a, b, weight, diff, norm);
  }

  RealType AddUpdate(const DiffusionTensor3D<RealValueType>& a,
                     const RealType& b)
  {
    if (this->m_ComponentSpace == Superclass::RIEMANNIAN)
    {
      return this->AddExponentialMapUpdate(a, b);
    }
    else
    {
      return this->AddEuclideanUpdate(a, b);
    }
  }
  template <class RealT>
  RealType AddUpdate(const RealT& a,
                     const RealType& b)
  {
    return this->AddEuclideanUpdate(a, b);
  }


  virtual void EnforceConstraints();

  virtual void Initialize();

  virtual void InitializeKernelSigma();

  virtual void InitializePatchWeights();
  virtual void InitializePatchWeightsSmoothDisc();

  virtual void InitializeIteration();

  virtual void ComputeKernelBandwidthUpdate(); // derived from base class; define here
  virtual ThreadDataStruct ThreadedComputeSigmaUpdate(const InputImageRegionType& regionToProcess,
                                                      const int itkNotUsed(threadId),
                                                      ThreadDataStruct threadData);
  virtual RealArrayType ResolveSigmaUpdate();


  virtual void ComputeImageUpdate();
  virtual ThreadDataStruct ThreadedComputeImageUpdate(const InputImageRegionType& regionToProcess,
                                                      const int threadId,
                                                      ThreadDataStruct threadData);
  virtual RealType ComputeGradientJointEntropy(InstanceIdentifier id,
                                               typename ListAdaptorType::Pointer& inList,
                                               BaseSamplerPointer& sampler);


  virtual void ApplyUpdate();
  virtual void ThreadedApplyUpdate(const InputImageRegionType& regionToProcess,
                                   const int itkNotUsed(threadId));


  virtual void PostProcessOutput();

  virtual void SetThreadData(int threadId, const ThreadDataStruct& data);
  virtual ThreadDataStruct GetThreadData(int threadId);

  std::vector<ThreadDataStruct> m_ThreadData;

  typename OutputImageType::Pointer m_UpdateBuffer;

  unsigned int     m_NumPixelComponents;
  unsigned int     m_NumIndependentComponents;
  unsigned int     m_TotalNumberPixels;
  //
  bool             m_UseSmoothDiscPatchWeights;
  //
  RealArrayType    m_GaussianKernelSigma;
  RealArrayType    m_IntensityRescaleInvFactor;
  PixelType        m_ZeroPixel;
  PixelArrayType   m_ImageMin;
  PixelArrayType   m_ImageMax;
  double           m_FractionPixelsForSigmaUpdate;
  bool             m_ComputeConditionalDerivatives;
  double           m_MinSigma;
  double           m_MinProbability;
  unsigned int     m_SigmaUpdateDecimationFactor;
  double           m_SigmaUpdateConvergenceTolerance;
  ShortArrayType   m_SigmaConverged;
  double           m_SigmaMultiplicationFactor;
  //
  RealType         m_NoiseSigma;
  RealType         m_NoiseSigmaSquared;
  bool             m_NoiseSigmaIsSet;
  //
  BaseSamplerPointer                m_Sampler;
  typename ListAdaptorType::Pointer m_SearchSpaceList;

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

  static ITK_THREAD_RETURN_TYPE ComputeSigmaUpdateThreaderCallback( void *arg );

  static ITK_THREAD_RETURN_TYPE ComputeImageUpdateThreaderCallback( void *arg );

  static ITK_THREAD_RETURN_TYPE ApplyUpdateThreaderCallback( void *arg );

  template <typename TInputImageType>
  void DispatchedMinMax(const TInputImageType* img);

  template <typename TInputImageType>
  void DispatchedArrayMinMax(const TInputImageType* img);

  template <typename TInputImageType>
  void DispatchedVectorMinMax(const TInputImageType* img);

    template <typename TInputImageType>
  void DispatchedRiemannianMinMax(const TInputImageType* img);
  static ITK_THREAD_RETURN_TYPE RiemannianMinMaxThreaderCallback( void *arg );

  ThreadDataStruct ThreadedRiemannianMinMax(const InputImageRegionType& regionToProcess,
                                            const int itkNotUsed(threadId),
                                            const InputImageType* img,
                                            ThreadDataStruct threadData);

  virtual void ResolveRiemannianMinMax();


  void ComputeSignedEuclideanDifferenceAndWeightedSquaredNorm(const PixelType& a, const PixelType& b,
                                                              const RealArrayType& weight,
                                                              RealType& diff, RealArrayType& norm);

  void ComputeLogMapAndWeightedSquaredGeodesicDifference(const DiffusionTensor3D<PixelValueType>& spdMatrixA,
                                                         const DiffusionTensor3D<PixelValueType>& spdMatrixB,
                                                         const RealArrayType& weight,
                                                         RealType& symMatrixLogMap, RealArrayType& geodesicDist);

  RealType AddEuclideanUpdate(const RealType& a, const RealType& b);
  RealType AddExponentialMapUpdate(const DiffusionTensor3D<RealValueType>& spdMatrix,
                                   const DiffusionTensor3D<RealValueType>& symMatrix);

  struct ThreadFilterStruct
  {
    PatchBasedDenoisingImageFilter *Filter;
    InputImageType                 *Img;
  };

 }; // end class PatchBasedDenoisingImageFilter

} // end namespace itk

// Define instantiation macro for this template
#define ITK_TEMPLATE_PatchBasedDenoisingImageFilter(_, EXPORT, x, y) namespace itk { \
    _(2(class EXPORT PatchBasedDenoisingImageFilter< ITK_TEMPLATE_2 x >)) | \
      namespace Templates { typedef PatchBasedDenoisingImageFilter< ITK_TEMPLATE_2 x > \
        PatchBasedDenoisingImageFilter##y; }                            \
  }

#if ITK_TEMPLATE_EXPLICIT
# include "Templates/itkPatchBasedDenoisingImageFilter+-.h"
#endif

#if ITK_TEMPLATE_TXX
# include "itkPatchBasedDenoisingImageFilter.hxx"
#endif

#endif