itkFEMRegistrationFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMRegistrationFilter.h,v $
  Language:  C++
  Date:      $Date: 2010-03-02 03:40:35 $
  Version:   $Revision: 1.29 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

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

#ifndef __itkFEMRegistrationFilter_h
#define __itkFEMRegistrationFilter_h

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

#include "itkImage.h"
#include "itkVector.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkVectorCastImageFilter.h"
#include "itkVectorIndexSelectionCastImageFilter.h"
#include "itkWarpImageFilter.h"
#include "itkImageToImageMetric.h"
#include "itkTranslationTransform.h"
#include "itkVectorExpandImageFilter.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<class TMovingImage,class TFixedImage> 
class  ITK_EXPORT  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 typename FixedImageType::PixelType  PixelType;
  typedef typename FixedImageType::SizeType   ImageSizeType;

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

  typedef Image< float, itkGetStaticConstMacro(ImageDimension) >            FloatImageType;
  typedef LinearSystemWrapperItpack            LinearSystemSolverType;
  typedef SolverCrankNicolson                  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;

  itkSetObjectMacro( Interpolator, InterpolatorType );

  itkGetObjectMacro( Interpolator, InterpolatorType );


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

  typedef itk::RecursiveMultiResolutionPyramidImageFilter<FixedImageType,FixedImageType>
    FixedPyramidType;

//#define USEIMAGEMETRIC
#ifdef  USEIMAGEMETRIC
  typedef ImageToImageMetric<ImageType,FixedImageType>   MetricBaseType;
  typedef  ImageMetricLoad<ImageType,ImageType>          ImageMetricLoadType;
#else
  typedef  FiniteDifferenceFunctionLoad<MovingImageType,FixedImageType>
                                                         ImageMetricLoadType;
  typedef PDEDeformableRegistrationFunction<FixedImageType,MovingImageType,FieldType>
                                                         MetricBaseType;
#endif
  typedef typename MetricBaseType::Pointer          MetricBaseTypePointer;
  /* Main functions */
 
  bool      ReadConfigFile(const char*);
  

  void      RunRegistration(void); 
  
  void      WriteWarpedImage(const char* fn);


  void      IterativeSolve(SolverType& S);  

  void      MultiResSolve();

  void      WarpImage(const MovingImageType * R);

  int       WriteDisplacementField(unsigned int index);

  int       WriteDisplacementFieldMultiComponent();

  void      SetMovingFile(const char* r)
    {
    m_MovingFileName=r;
    }

  std::string GetMovingFile()
    {
    return m_MovingFileName;
    }
  
  void      SetFixedFile(const char* t) {m_FixedFileName=t;}
  
  std::string GetFixedFile() {return m_FixedFileName;}

  
  void SetMovingImage(MovingImageType* R);
  
  void SetFixedImage(FixedImageType* T);
  
  MovingImageType* GetMovingImage(){return m_MovingImage;}
  MovingImageType* GetOriginalMovingImage(){return m_OriginalMovingImage;}
  
  FixedImageType* GetFixedImage(){return m_FixedImage;}
  
  
  FixedImageType* GetWarpedImage(){return m_WarpedImage;}

  void ComputeJacobian(float sign=1.0, FieldType* field=NULL , float smooth=0.0);

  FloatImageType* GetJacobianImage(){return m_FloatImage;}
  

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

  void      SetLandmarkFile(const char* l)
    {
    m_LandmarkFileName=l;
    }

  void      UseLandmarks(bool b)
    {
    m_UseLandmarks=b;
    }

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


  void      SetResultsFile(const char* r)
    {
    m_ResultsFileName=r;
    }

  void      SetResultsFileName (const char* f)
    {
    m_ResultsFileName=f;
    }

  std::string GetResultsFileName ()
    {
    return m_ResultsFileName;
    } 

  void      SetDisplacementsFile(const char* r)
    {
    m_DisplacementsFileName=r;
    }
  
  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;
    }
  
  void      SetTimeStep(Float i)
    {
    m_Dt=i;
    }
  
  void      SetAlpha(Float a)
    {
    m_Alpha=a;
    }
  
  void      SetEnergyReductionFactor(Float i)
    {
    m_EnergyReductionFactor=i;
    }

  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;
    } 

  void      DoLineSearch(unsigned int b)
    {
    m_DoLineSearchOnImageEnergy=b;
    } 

  void      DoMultiRes(bool b)
    {
    m_DoMultiRes=b;
    } 

  void      EmployRegridding(unsigned int b)
    {
    m_EmployRegridding=b;
    } 

  void      SetLineSearchMaximumIterations(unsigned int f)
    {
    m_LineSearchMaximumIterations=f;
    } 
  
  void      SetWriteDisplacements(bool b)
    {
    m_WriteDisplacementField=b;
    }
  
  bool      GetWriteDisplacements()
    {
    return m_WriteDisplacementField;
    }
    
  void      SetConfigFileName (const char* f){m_ConfigFileName=f;}
  
  std::string GetConfigFileName () {return m_ConfigFileName; }
  
  ImageSizeType GetImageSize(){ return m_FullImageSize; }

  MetricBaseTypePointer    GetMetric() { return m_Metric; }
  void      SetMetric(MetricBaseTypePointer MP) { m_Metric=MP; }
  
  void      ChooseMetric( float whichmetric); 
  
  void      SetElement(Element::Pointer e) {m_Element=e;}

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

  void      PrintVectorField(unsigned int modnum=1000);

  void      SetNumLevels(unsigned int i) { m_NumLevels=i; }
  void      SetMaxLevel(unsigned int i) { m_MaxLevel=i; }

  void      SetTemp(Float i) { m_Temp=i; }

  FEMRegistrationFilter( ); 
  ~FEMRegistrationFilter(); 

// HELPER FUNCTIONS
protected:

  
  class FEMOF : public FEMObjectFactory<FEMLightObject>{
  protected:
    FEMOF();
    ~FEMOF();
    };

 
  void      CreateMesh(double ElementsPerSide, Solver& S, ImageSizeType sz);

  void      ApplyLoads(SolverType& S,ImageSizeType Isz,double* spacing=NULL); 

  void      ApplyImageLoads(SolverType& S, MovingImageType* i1, FixedImageType* i2); 

  
  void      CreateLinearSystemSolver();

  
  Float     EvaluateEnergy();
 
  void      InterpolateVectorField(SolverType& S); 

  FloatImageType*      GetMetricImage(FieldType* F); 

 
  typedef  typename FieldType::Pointer FieldPointer;
  FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);
  

  void      SampleVectorFieldAtNodes(SolverType& S);
  
  Float EvaluateResidual(SolverType& mySolver,Float t);

  /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
  void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
  
  Float GoldenSection(SolverType& mySolver,Float tol=0.01,unsigned int MaxIters=25);

//  itkSetMacro( Load, ImageMetricLoadType* );
  itkGetConstMacro( Load, ImageMetricLoadType* );


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

private:

  void InitializeField();

  FEMRegistrationFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented
    
  std::string      m_ConfigFileName;
  std::string      m_ResultsFileName;
  std::string      m_MovingFileName;   // This variable is currently not being used.
  std::string      m_FixedFileName;    // This variable is currently not being used.
  std::string      m_LandmarkFileName;
  std::string      m_DisplacementsFileName;
  std::string      m_MeshFileName;

  unsigned int     m_DoLineSearchOnImageEnergy; 
  unsigned int     m_LineSearchMaximumIterations;

  vnl_vector<unsigned int>     m_NumberOfIntegrationPoints;// resolution of integration
  vnl_vector<unsigned int>     m_MetricWidth;
  vnl_vector<unsigned int>     m_Maxiters; // max iterations
  unsigned int                 m_TotalIterations;
  unsigned int                 m_NumLevels; // Number of Resolution Levels
  unsigned int                 m_MaxLevel;  // Maximum Level (NumLevels is original resolution).
  unsigned int                 m_MeshLevels;// Number of Mesh Resolutions ( should be >= 1)
  unsigned int                 m_MeshStep;  // Ratio Between Mesh Resolutions ( currently set to 2, should be >= 1)
  unsigned int                 m_FileCount; // keeps track of number of files written
  unsigned int                 m_CurrentLevel;

  typename FixedImageType::SizeType     m_CurrentLevelImageSize; 

  unsigned int     m_WhichMetric;
 
  vnl_vector<unsigned int> m_MeshPixelsPerElementAtEachResolution;

  Float                 m_Dt; // 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 
  Float     m_EnergyReductionFactor; 
  Float     m_Temp;

  bool         m_WriteDisplacementField;
  bool         m_DoMultiRes;
  bool         m_UseLandmarks;
  bool         m_ReadMeshFile;
  bool         m_UseMassMatrix;
  unsigned int m_EmployRegridding;
  Sign         m_DescentDirection;

  ImageSizeType     m_FullImageSize; // image size
  ImageSizeType     m_ImageOrigin; // image size
  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;
  ImageMetricLoadType*             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;

  LandmarkArrayType      m_LandmarkArray;
  InterpolatorPointer    m_Interpolator;

 

};

}} // end namespace itk::fem

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

#endif