itkBlockMatchingMethod.h

Go to the documentation of this file.

/*=============================================================================

  NifTK: A software platform for medical image computing.

  Copyright (c) University College London (UCL). All rights reserved.

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.

  See LICENSE.txt in the top level directory for details.

=============================================================================*/

#ifndef itkBlockMatchingMethod_h
#define itkBlockMatchingMethod_h


#include "itkMaskedImageRegistrationMethod.h"
#include <itkResampleImageFilter.h>
#include <itkArray.h>
#include <itkPointSet.h>
#include <itkRegionOfInterestImageFilter.h>
#include <itkNearestNeighborInterpolateImageFunction.h>
#include <itkSimilarityMeasure.h>
#include <itkDefaultStaticMeshTraits.h>
#include <itkPointSetToPointSetSingleValuedMetric.h>
#include <itkCovarianceSampleFilter.h>
#include <itkImageToListSampleAdaptor.h>
#include <itkGradientMagnitudeImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>

namespace itk
{
template <typename TImageType, class TScalarType>
class ITK_EXPORT BlockMatchingMethod 
: public MaskedImageRegistrationMethod<TImageType> 
{
public:

  typedef BlockMatchingMethod                                           Self;
  typedef MaskedImageRegistrationMethod<TImageType>                     Superclass;
  typedef SmartPointer<Self>                                            Pointer;
  typedef SmartPointer<const Self>                                      ConstPointer;
  
  itkNewMacro(Self);
  
  itkTypeMacro(BlockMatchingMethod, MaskedImageRegistrationMethod);

  itkStaticConstMacro(Dimension, unsigned int, TImageType::ImageDimension);

  typedef TImageType                                                    ImageType;
  typedef typename ImageType::SizeType                                  ImageSizeType;
  typedef typename ImageType::PixelType                                 ImagePixelType;
  typedef typename ImageType::SpacingType                               ImageSpacingType;
  typedef typename ImageType::RegionType                                ImageRegionType;
  typedef typename ImageType::IndexType                                 ImageIndexType;
  typedef ResampleImageFilter<TImageType, TImageType >                  ResampleFilterType;
  typedef typename ResampleFilterType::Pointer                          ResampleFilterPointer;
  typedef RegionOfInterestImageFilter<TImageType, TImageType>           RegionOfInterestFilterType;
  typedef typename RegionOfInterestFilterType::Pointer                  RegionOfInterestFilterPointer;
  
  // We either calculate the variance, based on the image, which will be whatever scalar type is passed in as template parameter
  typedef itk::Statistics::ImageToListSampleAdaptor< ImageType >                    ImageTypeListAdaptorType;
  typedef typename ImageTypeListAdaptorType::Pointer                                ImageTypeListAdaptorPointer;
  typedef itk::Statistics::CovarianceSampleFilter< ImageTypeListAdaptorType >         ImageTypeCovarianceSampleFilterType;
  typedef typename ImageTypeCovarianceSampleFilterType::Pointer                       ImageTypeCovarianceSampleFilterPointer;

  // Or we take the gradient magnitude image, so we should do it as a float.
  typedef Image<float, TImageType::ImageDimension>                                  GradientImageType;
  typedef GradientMagnitudeImageFilter<TImageType, GradientImageType>               GradientMagnitudeFilterType;
  typedef typename GradientMagnitudeFilterType::Pointer                             GradientMagnitudeFilterPointer;
  typedef itk::Statistics::ImageToListSampleAdaptor< GradientImageType >            GradientImageTypeListAdaptorType;
  typedef typename GradientImageTypeListAdaptorType::Pointer                        GradientImageTypeListAdaptorPointer;
  typedef itk::Statistics::CovarianceSampleFilter< GradientImageTypeListAdaptorType > GradientImageTypeCovarianceSampleFilterType;
  typedef typename GradientImageTypeCovarianceSampleFilterType::Pointer               GradientImageTypeCovarianceSampleFilterPointer;
  
  typedef SimilarityMeasure<TImageType, TImageType >                    SimilarityMeasureType;
  typedef SimilarityMeasureType*                                        SimilarityMeasurePointer;
  typedef NearestNeighborInterpolateImageFunction< TImageType, 
                                                   TScalarType>         DummyInterpolatorType;
  typedef typename DummyInterpolatorType::Pointer                       DummyInterpolatorPointer;
  typedef typename Superclass::ParametersType                           ParametersType;
  typedef PointSet< TScalarType, 
                    TImageType::ImageDimension,
                    DefaultStaticMeshTraits
                      < 
                        TScalarType, 
                        TImageType::ImageDimension, 
                        TImageType::ImageDimension,
                        TScalarType,
                        TScalarType
                      > 
                    >                                                   PointSetType;
  typedef typename PointSetType::PointType                              PointType;
  typedef typename PointSetType::PointsContainer::ConstIterator         PointIterator;
  typedef typename PointSetType::Pointer                                PointSetPointer;
  typedef typename PointSetType::PointsContainer                        PointsContainerType;
  typedef typename PointsContainerType::Pointer                         PointsContainerPointer;
  typedef PointSetToPointSetSingleValuedMetric<PointSetType, 
                                               PointSetType>            PointSetMetricType;
  typedef typename PointSetMetricType::Pointer                          PointSetMetricPointer;
  typedef typename Superclass::TransformType                            TransformType;
  typedef TransformType*                                                TransformPointer;
  typedef MinimumMaximumImageCalculator<TImageType>                     MinimumMaximumImageCalculatorType;
  typedef typename MinimumMaximumImageCalculatorType::Pointer           MinimumMaximumImageCalculatorPointer;
  
  itkSetObjectMacro(PointSetMetric, PointSetMetricType);
  itkGetObjectMacro(PointSetMetric, PointSetMetricType);
                                                 
  itkSetMacro(MaximumNumberOfIterationsRoundMainLoop, unsigned int);
  itkGetMacro(MaximumNumberOfIterationsRoundMainLoop, unsigned int);
  
  void SetBlockParameters(
    double blockSize,
    double blockHalfWidth,
    double blockSpacing,
    double blockSubSampling);
    
  itkGetMacro(BlockSize, double);

  itkGetMacro(BlockHalfWidth, double);

  itkGetMacro(BlockSpacing, double);

  itkGetMacro(BlockSubSampling, double);

  itkSetMacro(TransformedMovingImageFileName, std::string);
  itkGetMacro(TransformedMovingImageFileName, std::string);

  itkSetMacro(TransformedMovingImageFileExt, std::string);
  itkGetMacro(TransformedMovingImageFileExt, std::string);
  
  itkSetMacro(WriteTransformedMovingImage, bool);
  itkGetMacro(WriteTransformedMovingImage, bool);

  itkSetMacro(Epsilon, double);
  itkGetMacro(Epsilon, double);

  itkSetMacro(ParameterReductionFactor, double);
  itkGetMacro(ParameterReductionFactor, double);
  
  itkSetMacro(MinimumBlockSize, double);
  itkGetMacro(MinimumBlockSize, double);
  
  itkSetMacro(PercentageOfPointsToKeep, int);
  itkGetMacro(PercentageOfPointsToKeep, int);

  itkSetMacro(PercentageOfPointsInLeastTrimmedSquares, int);
  itkGetMacro(PercentageOfPointsInLeastTrimmedSquares, int);
  
  itkSetMacro(UseGradientMagnitudeVariance, bool);
  itkGetMacro(UseGradientMagnitudeVariance, bool);

  itkSetMacro(ScaleByMillimetres, bool);
  itkGetMacro(ScaleByMillimetres, bool);
  
  itkSetMacro(PointSetFileNameWithoutExtension, std::string);
  itkGetMacro(PointSetFileNameWithoutExtension, std::string);
  
  itkSetMacro(WritePointSet, bool);
  itkGetMacro(WritePointSet, bool);
  
  itkSetMacro(NoZero, bool);
  itkGetMacro(NoZero, bool);
  
  itkSetMacro(TransformedMovingImagePadValue, ImagePixelType);
  itkGetMacro(TransformedMovingImagePadValue, ImagePixelType);
  
protected:

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

  virtual void Initialize() throw (ExceptionObject);

  virtual void DoRegistration() throw (ExceptionObject);

  virtual void WritePointSet(const PointsContainerPointer& fixedPointContainer, 
      const PointsContainerPointer& movingPointContainer);
  
  virtual void GetPointCorrespondencies2D(
    ImageSizeType& size,
    ImageSizeType& bigN,
    ImageSizeType& bigOmega,
    ImageSizeType& bigDeltaOne,
    ImageSizeType& bigDeltaTwo,
    PointsContainerPointer& fixedPointContainer,
    PointsContainerPointer& movingPointContainer
    );
  
  virtual void GetPointCorrespondencies3D(
    ImageSizeType& size,
    ImageSizeType& bigN,
    ImageSizeType& bigOmega,
    ImageSizeType& bigDeltaOne,
    ImageSizeType& bigDeltaTwo,
    PointsContainerPointer& fixedPointContainer,
    PointsContainerPointer& movingPointContainer
    );

  virtual bool CheckEpsilon(ParametersType& previousParameters, ParametersType& currentParameters);
  
  virtual double CheckSinglePoint(ParametersType& previousParameters, ParametersType& currentParameters, ImageIndexType& index);
  
  virtual void TrimPoints(const TransformType* transform, 
      const PointsContainerType* fixedPoints,
      const PointsContainerType* movingPoints,
      PointsContainerType* trimmedFixedPoints,
      PointsContainerType* trimmedMovingPoints);
  
  class VarianceHeapDataType 
    {
      public:
        VarianceHeapDataType(const TScalarType& aKey, const ImageIndexType& anIndex)
          {
            this->key = aKey;
            this->index = anIndex;
          }
          
      void operator=(const VarianceHeapDataType& another)
        {
          this->key = another.key;
          this->index = another.index;
        }
        
      friend bool operator<(
        const VarianceHeapDataType& x, const VarianceHeapDataType& y) 
          {
            if(x.GetKey() < y.GetKey())
              return true;
            else
              return false;
          }

      void SetKey(const TScalarType& aKey) { this->key = aKey; }
      TScalarType GetKey() const { return this->key; }
      
      void SetIndex(const ImageIndexType& anIndex) { this->index = anIndex; }
      ImageIndexType GetIndex() const { return this->index; }
      
      private:
        TScalarType key;
        ImageIndexType index;  
    };
   
  typedef std::priority_queue<VarianceHeapDataType> VarianceHeap;
  
  class ResidualHeapDataType
    {

      public:
        ResidualHeapDataType(const TScalarType& aResidual, const PointType& aFixedPoint, const PointType& aMovingPoint)
          {
            this->residual = aResidual;
            this->fixed = aFixedPoint;
            this->moving = aMovingPoint;
          }
        
      void operator=(const ResidualHeapDataType& another)
        {
          this->residual = another.residual;
          this->fixed = another.fixed;
          this->moving = another.moving;
        }
      
      friend bool operator<(
        const ResidualHeapDataType& x, const ResidualHeapDataType& y) 
          {
            if(y.GetResidual() < x.GetResidual())
              return true;
            else
              return false;
          }

      void SetResidual(const TScalarType& aResidual) { this->residual = aResidual; }
      TScalarType GetResidual() const { return this->residual; }
    
      void SetFixed(const PointType& aFixedPoint) { this->fixed = aFixedPoint; }
      PointType GetFixed() const { return this->fixed; }

      void SetMoving(const PointType& aMovingPoint) { this->moving = aMovingPoint; }
      PointType GetMoving() const { return this->moving; }

      private:
        TScalarType residual;
        PointType fixed;
        PointType moving;
    };
  
  typedef std::priority_queue<ResidualHeapDataType> ResidualHeap;
  
private:
  
  BlockMatchingMethod(const Self&); // purposely not implemented
  void operator=(const Self&);      // purposely not implemented
  
  ResampleFilterPointer                        m_MovingImageResampler;
  
  RegionOfInterestFilterPointer                m_FixedImageRegionFilter;
  
  ImageTypeListAdaptorPointer                  m_FixedImageListAdaptor;
  ImageTypeCovarianceSampleFilterPointer       m_FixedImageCovarianceSampleFilter;

  GradientMagnitudeFilterPointer               m_GradientMagnitudeImageFilter;
  GradientImageTypeListAdaptorPointer          m_GradientMagnitudeListAdaptor;
  GradientImageTypeCovarianceSampleFilterPointer m_GradientMagnitudeCovarianceSampleFilter;
  
  DummyInterpolatorPointer                     m_DummyInterpolator;
 
  PointSetPointer                              m_FixedPointSet;
  PointsContainerPointer                       m_FixedPointSetContainer;
  
  PointSetPointer                              m_MovingPointSet;
  PointsContainerPointer                       m_MovingPointSetContainer;
  
  PointSetMetricPointer                        m_PointSetMetric;
  
  unsigned int                                 m_MaximumNumberOfIterationsRoundMainLoop;
  
  double                                       m_BlockSize;
  
  double                                       m_BlockHalfWidth;
  
  double                                       m_BlockSpacing;
  
  double                                       m_BlockSubSampling;
  
  double                                       m_Epsilon;
  
  double                                       m_ParameterReductionFactor;
  
  double                                       m_MinimumBlockSize;
  
  int                                          m_PercentageOfPointsToKeep;
  
  int                                          m_PercentageOfPointsInLeastTrimmedSquares;
  
  std::string                                  m_TransformedMovingImageFileName;
  
  std::string                                  m_TransformedMovingImageFileExt;
  
  bool                                         m_WriteTransformedMovingImage;

  std::string                                  m_PointSetFileNameWithoutExtension;
  
  bool                                         m_WritePointSet;
  
  bool                                         m_ScaleByMillimetres;
  
  bool                                         m_UseGradientMagnitudeVariance;
  
  bool                                         m_NoZero;
  
  ImagePixelType                               m_TransformedMovingImagePadValue;
  
  MinimumMaximumImageCalculatorPointer         m_MinMaxCalculator;
};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkBlockMatchingMethod.txx"
#endif

#endif