itkFEMRegistrationFilter.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 itkFEMRegistrationFilter_h
#define itkFEMRegistrationFilter_h

#include "itkFEMLinearSystemWrapperItpack.h"
#include "itkFEMLinearSystemWrapperDenseVNL.h"
#include "itkFEMObject.h"
#include "itkFEMSolverCrankNicolson.h"
#include "itkFEMMaterialLinearElasticity.h"
#include "itkFEMImageMetricLoad.h"
#include "itkFEMFiniteDifferenceFunctionLoad.h"

#include "itkImage.h"
#include "itkVector.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkImageToRectilinearFEMObjectFilter.h"
#include "itkVectorCastImageFilter.h"
#include "itkVectorIndexSelectionCastImageFilter.h"
#include "itkWarpImageFilter.h"
#include "itkImageToImageMetric.h"
#include "itkTranslationTransform.h"
#include "itkVectorExpandImageFilter.h"
#include "itkFixedArray.h"

#include "itkRecursiveMultiResolutionPyramidImageFilter.h"
#include "itkFEMLoadLandmark.h"

#include "vnl/vnl_vector.h"
#include "vnl/vnl_math.h"
#include "vnl/vnl_vector_fixed.h"

#include <iostream>
#include <string>

namespace itk
{
namespace fem
{

template <typename TMovingImage, typename TFixedImage, typename TFemObjectType>
class FEMRegistrationFilter : public ImageToImageFilter<TMovingImage, TFixedImage>
{
public:
  typedef FEMRegistrationFilter                         Self;
  typedef ImageToImageFilter<TMovingImage, TFixedImage> Superclass;
  typedef SmartPointer<Self>                            Pointer;
  typedef SmartPointer<const Self>                      ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(FEMRegistrationFilter, ImageToImageFilter );

  typedef TMovingImage                       MovingImageType;
  typedef TFixedImage                        FixedImageType;
  typedef TFemObjectType                     FEMObjectType;
  typedef typename FixedImageType::PixelType PixelType;
  typedef typename FixedImageType::SizeType  ImageSizeType;
  typedef typename FixedImageType::PointType PointType;

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

  typedef Image<float, itkGetStaticConstMacro(ImageDimension)> FloatImageType;
  typedef LinearSystemWrapperItpack                            LinearSystemSolverType;
  typedef SolverCrankNicolson<ImageDimension>                  SolverType;

  enum Sign { positive = 1, negative = -1 };
  typedef double          Float;
  typedef Load::ArrayType LoadArray;

  typedef std::vector<typename LoadLandmark::Pointer>                      LandmarkArrayType;
  typedef itk::Vector<float, itkGetStaticConstMacro(ImageDimension)>       VectorType;
  typedef itk::Image<VectorType, itkGetStaticConstMacro(ImageDimension)>   FieldType;
  typedef itk::WarpImageFilter<MovingImageType, FixedImageType, FieldType> WarperType;

  typedef MaterialLinearElasticity                          MaterialType;
  typedef itk::ImageRegionIteratorWithIndex<FixedImageType> ImageIterator;
  typedef itk::ImageRegionIteratorWithIndex<FloatImageType> FloatImageIterator;
  typedef itk::ImageRegionIteratorWithIndex<FieldType>      FieldIterator;

  typedef itk::VectorIndexSelectionCastImageFilter<FieldType, FloatImageType> IndexSelectCasterType;

  typedef double                                                  CoordRepType;
  typedef VectorInterpolateImageFunction<FieldType, CoordRepType> InterpolatorType;
  typedef typename InterpolatorType::Pointer                      InterpolatorPointer;

  typedef VectorLinearInterpolateImageFunction<FieldType, CoordRepType> DefaultInterpolatorType;

  typedef typename itk::Image<Element::ConstPointer, ImageDimension> InterpolationGridType;
  typedef typename InterpolationGridType::SizeType                   InterpolationGridSizeType;
  typedef typename InterpolationGridType::PointType                  InterpolationGridPointType;

  typedef itk::fem::ImageToRectilinearFEMObjectFilter<TMovingImage> ImageToMeshType;

  typedef itk::VectorExpandImageFilter<FieldType, FieldType> ExpanderType;
  typedef typename ExpanderType::ExpandFactorsType           ExpandFactorsType;

  typedef  typename FieldType::Pointer FieldPointer;

  typedef FiniteDifferenceFunctionLoad<MovingImageType, FixedImageType>
  ImageMetricLoadType;
  typedef PDEDeformableRegistrationFunction<FixedImageType, MovingImageType, FieldType>
  MetricBaseType;

  typedef typename MetricBaseType::Pointer     MetricBaseTypePointer;
  typedef FixedArray< double, ImageDimension > StandardDeviationsType;

  /*----------------------  Main functions ----------------------*/

  void SetMovingImage(MovingImageType* R);

  MovingImageType * GetMovingImage()
  {
    return m_MovingImage;
  }

  MovingImageType * GetOriginalMovingImage()
  {
    return m_OriginalMovingImage;
  }

  void SetFixedImage(FixedImageType* T);

  FixedImageType * GetFixedImage()
  {
    return m_FixedImage;
  }

  void SetInputFEMObject(FEMObjectType* F, unsigned int level = 0);

  FEMObjectType * GetInputFEMObject(unsigned int level = 0);

  void RunRegistration();

  void IterativeSolve(SolverType *S);

  void MultiResSolve();

  void WarpImage(const MovingImageType * R);

  FixedImageType * GetWarpedImage()
  {
    return m_WarpedImage;
  }

  void ComputeJacobian( );

  FloatImageType * GetJacobianImage()
  {
    return m_FloatImage;
  }

  FieldType * GetDisplacementField()
  {
    return m_Field;
  }

  void SetDisplacementField(FieldType* F)
  {
    m_FieldSize = F->GetLargestPossibleRegion().GetSize();
    m_Field = F;
  }

  void AddLandmark(PointType source, PointType target);

  void InsertLandmark(unsigned int i, PointType source, PointType target);

  void DeleteLandmark(unsigned int i);

  void ClearLandmarks();

  void GetLandmark(unsigned int i, PointType& source, PointType& target);

  void EnforceDiffeomorphism(float thresh, SolverType *S,  bool onlywriteimages);

  void SetMeshPixelsPerElementAtEachResolution(unsigned int i, unsigned int which = 0)
  {
    m_MeshPixelsPerElementAtEachResolution[which] = i;
  }

  void SetNumberOfIntegrationPoints(unsigned int i, unsigned int which = 0)
  {
    m_NumberOfIntegrationPoints[which] = i;
  }

  void SetWidthOfMetricRegion(unsigned int i, unsigned int which = 0)
  {
    m_MetricWidth[which] = i;
  }

  unsigned int GetWidthOfMetricRegion(unsigned int which = 0)
  {
    return m_MetricWidth[which];
  }

  void SetMaximumIterations(unsigned int i, unsigned int which)
  {
    m_Maxiters[which] = i;
  }

  itkSetMacro(TimeStep, Float);
  itkGetMacro(TimeStep, Float);

  itkSetMacro(Alpha, Float);
  itkGetMacro(Alpha, Float);

  itkSetMacro(UseLandmarks, bool);
  itkGetMacro(UseLandmarks, bool);
  void SetUseLandmarksOff()
  {
    SetUseLandmarks(false);
  }

  void SetUseLandmarksOn()
  {
    SetUseLandmarks(true);
  }

  itkSetMacro(UseMassMatrix, bool);
  itkGetMacro(UseMassMatrix, bool);

  itkSetMacro(EnergyReductionFactor, Float);
  itkGetMacro(EnergyReductionFactor, Float);

  void SetElasticity(Float i, unsigned int which = 0)
  {
    m_E[which] = i;
  }

  Float GetElasticity(unsigned int which = 0)
  {
    return m_E[which];
  }

  void  SetRho(Float r, unsigned int which = 0)
  {
    m_Rho[which] = r;
  }

  void SetGamma(Float r, unsigned int which = 0)
  {
    m_Gamma[which] = r;
  }

  void  SetDescentDirectionMinimize()
  {
    m_DescentDirection = positive;
  }

  void SetDescentDirectionMaximize()
  {
    m_DescentDirection = negative;
  }

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


  itkSetMacro(UseNormalizedGradient, bool);
  itkGetMacro(UseNormalizedGradient, bool);
  void SetUseNormalizedGradientOff()
  {
    SetUseNormalizedGradient(false);
  }

  void SetUseNormalizedGradientOn()
  {
    SetUseNormalizedGradient(true);
  }

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

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

  ImageSizeType GetImageSize()
  {
    return m_FullImageSize;
  }

  itkGetModifiableObjectMacro(Metric, MetricBaseType);
  itkSetObjectMacro(Metric, MetricBaseType);

  void ChooseMetric( unsigned int whichmetric);

  unsigned int GetMetricType()
  {
    return m_WhichMetric;
  }

  void SetElement(Element::Pointer e)
  {
    m_Element = e;
  }

  void SetMaterial(MaterialType::Pointer m)
  {
    m_Material = m;
  }

  void PrintVectorField(unsigned int modnum = 1000);

  void SetMaxLevel(unsigned int level);
  itkGetMacro(MaxLevel, unsigned int);

  itkSetMacro(CreateMeshFromImage, bool);
  void SetCreateMeshFromImageOn()
  {
    SetCreateMeshFromImage(true);
  }
  void SetCreateMeshFromImageOff()
  {
    SetCreateMeshFromImage(false);
  }
  itkGetMacro(CreateMeshFromImage, bool);

  itkSetObjectMacro( Interpolator, InterpolatorType );

  itkGetModifiableObjectMacro( Interpolator, InterpolatorType );

  itkSetMacro(StandardDeviations, StandardDeviationsType);
  virtual void SetStandardDeviations(double value);

  itkGetConstReferenceMacro(StandardDeviations, StandardDeviationsType);

  itkSetMacro(MaximumKernelWidth, unsigned int);
  itkGetConstMacro(MaximumKernelWidth, unsigned int);

  itkSetMacro(MaximumError, double);
  itkGetConstMacro(MaximumError, double);


// HELPER FUNCTIONS

protected:
  FEMRegistrationFilter();
  ~FEMRegistrationFilter();

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

  void CreateMesh(unsigned int ElementsPerSide, SolverType *solver);

  void ApplyLoads(ImageSizeType Isz, double* spacing = ITK_NULLPTR);

  void ApplyImageLoads(MovingImageType* i1, FixedImageType* i2);

  // FIXME - Not implemented
  void  CreateLinearSystemSolver();

  // FIXME - Not implemented
  Float EvaluateEnergy();

  void InterpolateVectorField(SolverType *S);

  // FIXME - Not implemented
  FloatImageType * GetMetricImage(FieldType* F);

  FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);

  void SampleVectorFieldAtNodes(SolverType *S);

  Float EvaluateResidual(SolverType *mySolver, Float t);

  // FIXME - Replace with BSD Code
  /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
  // void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
  void FindBracketingTriplet(SolverType *mySolver, Float* a, Float* b, Float* c);

  Float GoldenSection(SolverType *mySolver, Float tol = 0.01, unsigned int MaxIters = 25);

  itkSetObjectMacro( Load, ImageMetricLoadType );
  itkGetModifiableObjectMacro(Load, ImageMetricLoadType );

  void SmoothDisplacementField();

private:

  void InitializeField();

  FEMRegistrationFilter(const Self &) ITK_DELETE_FUNCTION;
  void operator=(const Self &) ITK_DELETE_FUNCTION;

  unsigned int m_DoLineSearchOnImageEnergy;
  unsigned int m_LineSearchMaximumIterations;

  /* Parameters used to define Multi-resolution Registration */
  vnl_vector<unsigned int> m_NumberOfIntegrationPoints; // resolution of integration
  vnl_vector<unsigned int> m_MetricWidth;               // number of iterations at each level
  vnl_vector<unsigned int> m_Maxiters;                  // max iterations
  unsigned int             m_TotalIterations;           // total number of iterations that were run
  unsigned int             m_MaxLevel;                  // Number of Resolution Levels for registration.
  unsigned int             m_FileCount;                 // keeps track of number of files written
  unsigned int             m_CurrentLevel;              // current resolution level

  typename FixedImageType::SizeType     m_CurrentLevelImageSize;

  unsigned int m_WhichMetric;                 // Metric used for image registration
                                              // 0 = Mean Squares
                                              // 1 = Normalized Correlation
                                              // 2 = Mutual Information
                                              // 3 = Demons Metric

  vnl_vector<unsigned int> m_MeshPixelsPerElementAtEachResolution;

  Float             m_TimeStep;        // time step
  vnl_vector<Float> m_E;               // elasticity
  vnl_vector<Float> m_Rho;             // mass matrix weight
  vnl_vector<Float> m_Gamma;           // image similarity weight
  Float             m_Energy;          // current value of energy
  Float             m_MinE;            // minimum recorded energy
  Float             m_MinJacobian;     // minimum recorded energy
  Float             m_Alpha;           // difference parameter

  bool          m_UseLandmarks;               // Use landmark points
  bool          m_UseMassMatrix;              // Use Mass matrix in FEM solution
  bool          m_UseNormalizedGradient;      // Use normalized gradient magnitude in the metric
  bool          m_CreateMeshFromImage;        // Create the mesh based on the fixed image
  unsigned int  m_EmployRegridding;           // Use regridding
  Sign          m_DescentDirection;           // Metric minimizes or maximizes
  Float         m_EnergyReductionFactor;      // Factor we want to reduce the energy - Helps with convergence
  ImageSizeType m_FullImageSize;              // image size
  ImageSizeType m_ImageOrigin;                // image origin

  ImageSizeType                  m_ImageScaling;
  ImageSizeType                  m_CurrentImageScaling;
  typename FieldType::RegionType m_FieldRegion;
  typename FieldType::SizeType   m_FieldSize;
  typename FieldType::Pointer    m_Field;

  // only use TotalField if re-gridding is employed.
  typename FieldType::Pointer               m_TotalField;
  typename ImageMetricLoadType::Pointer     m_Load; // Defines the load to use

  // define the warper
  typename WarperType::Pointer          m_Warper;

  // declare a new image to hold the warped  reference
  typename FixedImageType::Pointer      m_WarpedImage;
  typename FloatImageType::Pointer      m_FloatImage;
  typename FixedImageType::RegionType   m_Wregion;
  typename FixedImageType::IndexType    m_Windex;

  // declare images for target and reference
  typename MovingImageType::Pointer     m_MovingImage;
  typename MovingImageType::Pointer     m_OriginalMovingImage;
  typename FixedImageType::Pointer      m_FixedImage;

  // element and metric pointers
  typename Element::Pointer             m_Element;
  typename MaterialType::Pointer        m_Material;
  MetricBaseTypePointer                 m_Metric;
  typename FEMObjectType::Pointer       m_FEMObject;

  LandmarkArrayType   m_LandmarkArray;
  InterpolatorPointer m_Interpolator;

  double m_MaximumError;

  unsigned int m_MaximumKernelWidth;

  StandardDeviationsType m_StandardDeviations;

};

}
}  // end namespace itk::fem

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFEMRegistrationFilter.hxx"
#endif

#endif