itkRegistrationBasedCTEFilter.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 itkRegistrationBasedCTEFilter_h
#define itkRegistrationBasedCTEFilter_h

#include <itkImageToImageFilter.h>
#include <itkVector.h>
#include <itkImage.h>
#include <itkImageFileWriter.h>
#include <itkLinearInterpolateImageFunction.h>
#include <itkGaussianSmoothVectorFieldFilter.h>
#include <itkDisplacementFieldJacobianDeterminantFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkSetOutputVectorToCurrentPositionFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkVectorMagnitudeImageFilter.h>
#include <itkVectorVPlusLambdaUImageFilter.h>
#include "itkVectorPhiPlusDeltaTTimesVFilter.h"
#include "itkDasGradientFilter.h"
#include <itkDiscreteGaussianImageFilter.h>

namespace itk {

/*
 * \class RegistrationBasedCTEFilter
 * \brief This class implements Das et al Neuroimage 45 (2009) 867-879.
 * 
 * The inputs to this filter should be exactly as follows.
 * <pre>
 * input1 = white matter pv image, set using SetWhiteMatterPVMap(image)
 * input2 = white + grey pv image, set using SetWhitePlusGreyMatterPVMap(image)
 * input3 = thickness prior image, set using SetThicknessPriorMap(image)
 * input4 = the grey white border, as a binary image where 1 = boundary and 0 is background, set using SetGWI(image)
 * </pre>
 * 
 * The output is the DiReCT thickness map.
 * 
 * After registration, we can additionally output the displacement field, 
 * using the method WriteDisplacementField(filename). 
 */
template< class TInputImage, typename TScalarType>
class ITK_EXPORT RegistrationBasedCTEFilter :
  public ImageToImageFilter<TInputImage, TInputImage>
{
public:

  typedef TScalarType                                                   VectorDataType;
  typedef RegistrationBasedCTEFilter                                    Self;
  typedef ImageToImageFilter<TInputImage, TInputImage>                  Superclass;
  typedef SmartPointer<Self>                                            Pointer;
  typedef SmartPointer<const Self>                                      ConstPointer;

  itkNewMacro(Self);
  
  itkTypeMacro(RegistrationBasedCTEFilter, ImageToImageFilter);

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

  typedef TInputImage                                                   ImageType;
  typedef typename ImageType::PixelType                                 PixelType;
  typedef typename ImageType::Pointer                                   ImagePointer;
  typedef typename ImageType::SizeType                                  SizeType;
  typedef typename ImageType::IndexType                                 IndexType;
  typedef typename ImageType::RegionType                                RegionType;
  typedef typename ImageType::SpacingType                               SpacingType;
  typedef typename ImageType::PointType                                 OriginType;
  typedef typename ImageType::DirectionType                             DirectionType;
  typedef Vector< VectorDataType, itkGetStaticConstMacro(Dimension) >   VectorPixelType;
  typedef Image< VectorPixelType, itkGetStaticConstMacro(Dimension) >   VectorImageType;
  typedef typename VectorImageType::Pointer                             VectorImagePointer;
  typedef typename VectorImageType::SizeType                            VectorImageSizeType;
  typedef typename VectorImageType::IndexType                           VectorImageIndexType;
  typedef typename VectorImageType::RegionType                          VectorImageRegionType;
  typedef typename VectorImageType::SpacingType                         VectorImageSpacingType;
  typedef typename VectorImageType::PointType                           VectorImagePointType;
  typedef typename VectorImageType::PointType                           VectorImageOriginType;
  typedef typename VectorImageType::DirectionType                       VectorImageDirectionType;
  typedef ImageFileReader< VectorImageType  >                           VectorImageReaderType;
  typedef ImageFileWriter< VectorImageType  >                           VectorImageWriterType;
  typedef ImageFileWriter< ImageType >                                  ScalarImageWriterType;
  
  typedef LinearInterpolateImageFunction< TInputImage, TScalarType >    LinearInterpolatorType;
  typedef GaussianSmoothVectorFieldFilter<TScalarType, 
                     itkGetStaticConstMacro(Dimension),
                     itkGetStaticConstMacro(Dimension)
                     >                                                  ConvolveFilterType;
  typedef DisplacementFieldJacobianDeterminantFilter<VectorImageType, 
                                                     TScalarType,
                                                     TInputImage>       JacobianFilterType;
  typedef MinimumMaximumImageCalculator<TInputImage>                    MinMaxJacobianType;
  typedef SetOutputVectorToCurrentPositionFilter<TScalarType, 
                                     itkGetStaticConstMacro(Dimension)> InitializePhiFilterType;
  typedef SubtractImageFilter<VectorImageType, VectorImageType>         SubtractImageFilterType;
  typedef VectorMagnitudeImageFilter<VectorImageType, TInputImage>      VectorMagnitudeFilterType; 
  typedef VectorVPlusLambdaUImageFilter<TScalarType,
                                     itkGetStaticConstMacro(Dimension)> VectorVPlusLambdaUFilterType;

  typedef VectorPhiPlusDeltaTTimesVFilter<TScalarType,
                                     itkGetStaticConstMacro(Dimension)> VectorPhiPlusDeltaTTimesVFilterType;
  typedef DasGradientFilter<TScalarType, 
                            itkGetStaticConstMacro(Dimension)>          DasGradientFilterType;
  typedef DasTransformImageFilter<TScalarType, 
                                  itkGetStaticConstMacro(Dimension)>    DasTransformImageFilterType;                               
  typedef DiscreteGaussianImageFilter<ImageType,ImageType>              GaussianSmoothImageFilterType;
  
  void SetWhiteMatterPVMap(ImagePointer image) { this->SetInput(0, image); }
  
  void SetWhitePlusGreyMatterPVMap(ImagePointer image) { this->SetInput(1, image); }
  
  void SetThicknessPriorMap(ImagePointer image) { this->SetInput(2, image); }

  void SetGWI(ImagePointer image) { this->SetInput(3, image); }
  
  itkSetMacro(MaxIterations, unsigned int);
  itkGetMacro(MaxIterations, unsigned int);

  itkSetMacro(M, unsigned int);
  itkGetMacro(M, unsigned int);

  itkSetMacro(N, unsigned int);
  itkGetMacro(N, unsigned int);
  
  itkSetMacro(Lambda, double);
  itkGetMacro(Lambda, double);
  
  itkSetMacro(Sigma, double);
  itkGetMacro(Sigma, double);

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

  itkSetMacro(Alpha, double);
  itkGetMacro(Alpha, double);

  void WriteDisplacementField(std::string filename);
  
  itkSetMacro(OutputAbsoluteLocation, bool);
  itkGetMacro(OutputAbsoluteLocation, bool);

  itkSetMacro(TrackJacobian, bool);
  itkGetMacro(TrackJacobian, bool);

  itkSetMacro(WriteMagnitudeOfDisplacementImage, bool);
  itkGetMacro(WriteMagnitudeOfDisplacementImage, bool);

  itkSetMacro(WriteTSurfImage, bool);
  itkGetMacro(WriteTSurfImage, bool);

  itkSetMacro(WriteGradientImage, bool);
  itkGetMacro(WriteGradientImage, bool);

  itkSetMacro(WriteVelocityImage, bool);
  itkGetMacro(WriteVelocityImage, bool);

  itkSetMacro(WriteTransformedMovingImage, bool);
  itkGetMacro(WriteTransformedMovingImage, bool);

  itkSetMacro(SmoothPVMaps, bool);
  itkGetMacro(SmoothPVMaps, bool);

  itkSetMacro(SmoothPVMapSigma, double);
  itkGetMacro(SmoothPVMapSigma, double);

  itkSetMacro(UseGradientMovingImage, bool);
  itkGetMacro(UseGradientMovingImage, bool);

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

  /* The main filter method. Note, single threaded. */
  virtual void GenerateData();

private:
  
  RegistrationBasedCTEFilter(const Self&); 
  void operator=(const Self&); 
  
  unsigned int       m_MaxIterations;
  unsigned int       m_M;
  unsigned int       m_N;
  double             m_Epsilon;
  double             m_Sigma;
  double             m_Lambda;
  double             m_Alpha;
  double             m_SmoothPVMapSigma;
  bool               m_OutputAbsoluteLocation;
  bool               m_TrackJacobian;
  bool               m_WriteMagnitudeOfDisplacementImage;
  bool               m_WriteTSurfImage;
  bool               m_WriteGradientImage;
  bool               m_WriteVelocityImage;
  bool               m_WriteTransformedMovingImage;
  bool               m_SmoothPVMaps;
  bool               m_UseGradientMovingImage;
  
  VectorImagePointer m_PhiZeroImageUninitialized;
  VectorImagePointer m_PhiZeroImage;
  VectorImagePointer m_FinalPhiImage;
  
  typename InitializePhiFilterType::Pointer m_InitializePhiZeroFilter;
  typename GaussianSmoothImageFilterType::Pointer m_SmoothWMPVMapFilter;
  typename GaussianSmoothImageFilterType::Pointer m_SmoothGMWMPVMapFilter;
  
  void InitializeScalarImage(ImageType *image,
                             RegionType& region,
                             SpacingType& spacing,
                             OriginType& origin,
                             DirectionType& direction);
  
  void InitializeVectorImage(VectorImageType* image, 
                             const VectorImageRegionType& region, 
                             const VectorImageSpacingType& spacing, 
                             const VectorImageOriginType& origin, 
                             const VectorImageDirectionType& direction);
  
  double EvaluateVelocityField(VectorImageType* velocityField, double dt);
  
  double EvaluateRegistrationSimilarity(VectorImageType* phi, ImageType* target, ImageType* source);
  
  double EvaluateCostFunction(double velocityFieldEnergy, double imageSimilarity);
  
  void CalculateMinAndMaxJacobian(VectorImageType* phi, double &min, double& max);
  
  double CalculateMaxDisplacement(VectorImageType* phi);

  double CalculateMaxMagnitude(VectorImageType* vec);
  
  void CopyVectorField(VectorImageType* a, VectorImageType* b);
  
  void RemoveFile(std::string filename);
  
  std::string GetFileName(int i, std::string, std::string);
  
}; // end class
  
} // end namespace

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

#endif