itkFEMRegistrationFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMRegistrationFilter.h,v $
  Language:  C++
  Date:      $Date: 2004/02/06 16:55:17 $
  Version:   $Revision: 1.22 $

  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* );
  itkGetMacro( 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;  
  std::string      m_FixedFileName;
  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;


  // multi-resolution stuff
//  typename FixedPyramidType::Pointer   m_FixedPyramid;
//  typename FixedPyramidType::Pointer   m_MovingPyramid;
  
  LandmarkArrayType    m_LandmarkArray;

  InterpolatorPointer    m_Interpolator;

 

};

}} // end namespace itk::fem

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

#endif

---