itkUCLBSplineTransform.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 itkUCLBSplineTransform_h
#define itkUCLBSplineTransform_h

#include <iostream>
#include "itkDeformableTransform.h"
#include <itkImage.h>
#include <itkImageRegion.h>
#include <itkVector.h>
#include <itkArray2D.h>

#include <itkContinuousIndex.h>
#include <itkImageRegionIterator.h>
#include <itkSingleValuedCostFunction.h>
#include <itkScalarImageToNormalizedGradientVectorImageFilter.h>

namespace itk
{

template <
    class TFixedImage,                   // Templated over the image type.
    class TScalarType,                   // Data type for scalars
    unsigned int NDimensions,            // Number of Dimensions i.e. 2D or 3D
    class TDeformationScalar             // Data type in the deformation field. 
    >            
class ITK_EXPORT UCLBSplineTransform : 
          public DeformableTransform< TFixedImage, TScalarType, NDimensions, TDeformationScalar >
{
public:
  
  typedef UCLBSplineTransform                                            Self;
  typedef DeformableTransform< TFixedImage, TScalarType, NDimensions, TDeformationScalar> Superclass;
  typedef SmartPointer<Self>                                          Pointer;
  typedef SmartPointer<const Self>                                    ConstPointer;
      
  itkNewMacro( Self );

  itkTypeMacro( UCLBSplineTransform, DeformableTransform );

  itkStaticConstMacro(SpaceDimension, unsigned int, NDimensions);

  typedef typename Superclass::ScalarType                         ScalarType;

  typedef typename Superclass::ParametersType                     ParametersType;
  typedef typename Superclass::DerivativeType                     DerivativeType;
  typedef typename Superclass::JacobianType                       JacobianType;

  typedef typename Superclass::InputPointType                     InputPointType;
  typedef typename Superclass::OutputPointType                    OutputPointType;

  typedef typename Superclass::DeformationFieldPixelType          DeformationFieldPixelType;
  typedef typename Superclass::DeformationFieldType               DeformationFieldType;
  typedef typename Superclass::DeformationFieldPointer            DeformationFieldPointer;
  typedef typename Superclass::DeformationFieldRegionType         DeformationFieldRegionType;
  typedef typename Superclass::DeformationFieldIndexType          DeformationFieldIndexType;
  typedef typename Superclass::DeformationFieldSizeType           DeformationFieldSizeType;
  typedef typename Superclass::DeformationFieldSpacingType        DeformationFieldSpacingType;
  typedef typename Superclass::DeformationFieldDirectionType      DeformationFieldDirectionType;
  typedef typename Superclass::DeformationFieldOriginType         DeformationFieldOriginType;
  
  typedef Vector< TDeformationScalar, NDimensions >               GridPixelType;
  typedef Image<GridPixelType, NDimensions>                       GridImageType;
  typedef typename GridImageType::Pointer                         GridImagePointer;
  typedef ImageRegion<NDimensions>                                GridRegionType;
  typedef typename GridRegionType::IndexType                      GridIndexType;
  typedef typename GridRegionType::SizeType                       GridSizeType;
  typedef typename GridImageType::SpacingType                     GridSpacingType;
  typedef typename GridImageType::DirectionType                   GridDirectionType;
  typedef typename GridImageType::PointType                       GridOriginType;
  
  typedef TDeformationScalar                                      BendingEnergyPixelType;
  typedef Image<TDeformationScalar, NDimensions>                  BendingEnergyImageType;
  typedef ImageRegion<NDimensions>                                BendingEnergyImageRegionType;
  typedef typename BendingEnergyImageRegionType::SizeType         BendingEnergyImageSizeType;
  typedef typename BendingEnergyImageType::Pointer                BendingEnergyImagePointer;
  typedef const BendingEnergyImageType*                           BendingEnergyImageConstPointer;  
  typedef ImageRegionIterator<BendingEnergyImageType>             BendingEnergyIteratorType;
  typedef ScalarImageToNormalizedGradientVectorImageFilter<BendingEnergyImageType, 
                                                     TDeformationScalar> BendingEnergyDerivativeFilterType;
  typedef typename BendingEnergyDerivativeFilterType::Pointer     BendingEnergyDerivativeFilterPointer;
  typedef typename BendingEnergyDerivativeFilterType::OutputPixelType BendingEnergyDerivativePixelType;
  typedef typename BendingEnergyDerivativeFilterType::OutputImageType BendingEnergyDerivativeImageType;
  typedef ImageRegionConstIterator<BendingEnergyDerivativeImageType>  BendingEnergyDerivativeIteratorType;
  
  typedef TFixedImage                                             FixedImageType;
  typedef typename FixedImageType::ConstPointer                   FixedImagePointer;
  
  typedef itk::ContinuousIndex<TDeformationScalar, NDimensions>   GridVoxelCoordinateType;
  typedef itk::ImageRegionIterator<DeformationFieldType>          DeformationFieldIteratorType;    
  typedef itk::ImageRegionConstIteratorWithIndex<GridImageType>   GridConstIteratorType; 
  typedef itk::ImageRegionIterator<GridImageType>                 GridIteratorType;
  
  typedef SingleValuedCostFunction::MeasureType                   MeasureType;

  virtual void SetIdentity();

  void Initialize(FixedImagePointer fixedImage, GridSpacingType finalGridSpacingInMillimetres, int numberOfLevels);
  
  void InterpolateNextGrid(FixedImagePointer fixedImage);
  
  void SetParameters(const ParametersType & parameters);

  virtual const ParametersType& GetFixedParameters(void) const;
  
  virtual void SetFixedParameters(const ParametersType& parameters);

  MeasureType GetBendingEnergy();

  itkGetMacro(BendingEnergyHasBeenUpdatedFlag, bool);
  
  void GetBendingEnergyDerivative(DerivativeType & derivative);
  
  void GetBendingEnergyDerivativeDaniel(DerivativeType & derivative);
  
  void GetBendingEnergyDerivativeMarc(DerivativeType & derivative);
  
  itkGetObjectMacro(Grid, GridImageType);

  virtual void WriteParameters(std::string filename) { this->WriteControlPointImage(filename); }

  void WriteControlPointImage(std::string filename);

  virtual OutputPointType  TransformPoint(const InputPointType  &point ) const;

protected:

  UCLBSplineTransform();
  virtual ~UCLBSplineTransform();

  void PrintSelf(std::ostream &os, Indent indent) const;

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

  const static unsigned int s_LookupTableRows = 1000;
  const static unsigned int s_LookupTableSize = s_LookupTableRows - 1;
  const static unsigned int s_LookupTableCols = 4;
  
  GridImagePointer          m_Grid;
  GridImagePointer          m_OldGrid;
  
  bool                                 m_BendingEnergyHasBeenUpdatedFlag;
  BendingEnergyImagePointer            m_BendingEnergyGrid;
  BendingEnergyDerivativeFilterPointer m_BendingEnergyDerivativeFilter;
  
  void InterpolateDeformationField2D();
  
  void InterpolateDeformationField3DDaniel();

  void InterpolateDeformationField3DMarc();

  MeasureType GetBendingEnergy2D(TScalarType divisor) const;
  
  MeasureType GetBendingEnergy3DDaniel(TScalarType divisor) const;

  MeasureType GetBendingEnergy3DMarc(TScalarType divisor) const;

  void InterpolateNextGrid2D(GridImagePointer& oldGrid, GridImagePointer &newGrid);
  
  void InterpolateNextGrid3D(GridImagePointer& oldGrid, GridImagePointer &newGrid);
  
  TScalarType GetNewOrigin(TScalarType oldSize, TScalarType oldSpacing, TScalarType oldOrigin, TScalarType newSize, TScalarType newSpacing);
  
  void InitializeGrid(FixedImagePointer fixedImage, 
      GridRegionType gridRegion,
      GridSpacingType gridSpacing,
      GridDirectionType gridDirection,
      GridOriginType gridOrigin);
  
  Array2D<TScalarType>    m_Lookup;
  Array2D<TScalarType>    m_Lookup1stDerivative;
  Array2D<TScalarType>    m_Lookup2ndDerivative;

  TScalarType B0(TScalarType u)  { return ((1-u)*(1-u)*(1-u))/6.0; } 
  TScalarType B1(TScalarType u)  { return (3*u*u*u - 6*u*u + 4)/6.0; }
  TScalarType B2(TScalarType u)  { return (-3*u*u*u + 3*u*u + 3*u + 1)/6.0; }
  TScalarType B3(TScalarType u)  { return (u*u*u)/6.0; }
  TScalarType B01(TScalarType u) { return -(1-u)*(1-u)/2.0; }
  TScalarType B11(TScalarType u) { return (9*u*u - 12*u)/6.0; }
  TScalarType B21(TScalarType u) { return (-9*u*u + 6*u + 3)/6.0; }
  TScalarType B31(TScalarType u) { return (u*u)/2.0; }
  TScalarType B02(TScalarType u) { return 1-u; }
  TScalarType B12(TScalarType u) { return 3*u - 2; }
  TScalarType B22(TScalarType u) { return -3*u + 1; }
  TScalarType B32(TScalarType u) { return u; }

  TScalarType B(int i, TScalarType u)
    {
      switch (i) {
        case 0: return B0(u);
        case 1: return B1(u);
        case 2: return B2(u);
        case 3: return B3(u);
      }
      return 0;
    }

  TScalarType B1(int i, TScalarType u)
    {
      switch (i) {
        case 0: return B01(u);
        case 1: return B11(u);
        case 2: return B21(u);
        case 3: return B31(u);
      }
      return 0;
    }

  TScalarType B2(int i, TScalarType u)
    {
      switch (i) {
        case 0: return B02(u);
        case 1: return B12(u);
        case 2: return B22(u);
        case 3: return B32(u);
      }
      return 0;
    }

  TScalarType B(int i, int derivative, TScalarType u)
    {
      switch(derivative) {
        case 0: return B(i, u);
        case 1: return B1(i, u);
        case 2: return B2(i, u);
      }
      return 0;
    }
  
}; // class UCLBSplineTransform


}  // namespace itk

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

#endif /* itkUCLBSplineTransform_h */