itkFEMRegistrationFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMRegistrationFilter.h,v $
  Language:  C++
  Date:      $Date: 2002/11/14 21:02:30 $
  Version:   $Revision: 1.8 $

  Copyright (c) 2002 Insight 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 "itkImage.h"
#include "itkVector.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkVectorCastImageFilter.h"
#include "itkVectorIndexSelectionCastImageFilter.h"
#include "itkWarpImageFilter.h"
#include "itkImageToImageMetric.h"
#include "itkTranslationTransform.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 TReference,class TTarget> 
class  ITK_EXPORT  FEMRegistrationFilter : public ImageToImageFilter<TReference, TTarget>
{
public:
  typedef FEMRegistrationFilter                              Self;
  typedef ImageToImageFilter<TReference, TTarget> Superclass;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  itkNewMacro(Self);
  
  itkTypeMacro(FEMRegistrationFilter, ImageToImageFilter );
  
  typedef TReference                                ImageType;
  typedef TTarget                                   TargetImageType;
  typedef typename ImageType::PixelType             ImageDataType;
  typedef typename ImageType::PixelType             PixelType;
  typedef typename ImageType::SizeType              ImageSizeType;

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

  typedef Image< float, itkGetStaticConstMacro(ImageDimension) >            FloatImageType;
  typedef LinearSystemWrapperItpack                 LinearSystemSolverType;
  typedef SolverCrankNicolson                       SolverType;
  enum Sign { positive = 1, negative = -1 };
  typedef double                                    Float;
   

  typedef MaterialLinearElasticity                  MaterialType;
  typedef ImageToImageMetric<ImageType,TargetImageType>   MetricBaseType;
  typedef typename MetricBaseType::Pointer          MetricBaseTypePointer;
  typedef itk::Vector<float,itkGetStaticConstMacro(ImageDimension)>         VectorType;
  typedef itk::Image<VectorType,itkGetStaticConstMacro(ImageDimension)>     FieldType;
  typedef itk::WarpImageFilter<ImageType,ImageType, FieldType> WarperType;
  typedef itk::ImageRegionIteratorWithIndex<ImageType>         ImageIterator; 
  typedef itk::ImageRegionIteratorWithIndex<FieldType>         FieldIterator; 
  typedef itk::VectorIndexSelectionCastImageFilter<FieldType,FloatImageType> IndexSelectCasterType;

  typedef  ImageMetricLoad<ImageType,ImageType>     ImageMetricLoadType;
  
  class FEMOF : public FEMObjectFactory<FEMLightObject>{};
  
  /* Main functions */
 
  bool      ReadConfigFile(const char*);

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


  void      CreateMesh(double ElementsPerSide, Solver& S);

  void      ApplyLoads(SolverType& S,ImageSizeType Isz); 

  void      CreateLinearSystemSolver();

  void      IterativeSolve(SolverType& S);  

  void      MultiResSolve();

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

  FloatImageType*      GetMetricImage(FieldType* F); 
  
  void      SampleVectorFieldAtNodes(SolverType& S);

  void      WarpImage(ImageType* R);      

  int       WriteDisplacementField(unsigned int index);

  void      SetReferenceFile(const char* r) {m_ReferenceFileName=r;}

  std::string GetReferenceFile() {return m_ReferenceFileName;}
  
  void      SetTargetFile(const char* t) {m_TargetFileName=t;}
  
  std::string GetTargetFile() {return m_TargetFileName;}

  
  void SetReferenceImage(ImageType* R);
  
  void SetTargetImage(TargetImageType* T);
  
  ImageType* GetReferenceImage(){return m_RefImg;}
  
  TargetImageType* GetTargetImage(){return m_TarImg;}
  
  
  ImageType* GetWarpedImage(){return m_WarpedImage;}

  void ComputeJacobian();

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

  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();
  
  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* );

  FEMRegistrationFilter( ); 
  ~FEMRegistrationFilter(); 

protected :


  void PrintSelf(std::ostream& os, Indent indent) const 
  { 
      std::cout << "FEM registration filter "<< std::endl;
      Superclass::PrintSelf( os, indent );
  }

private :

  FEMRegistrationFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented
    
  std::string      m_ConfigFileName;
  std::string      m_ResultsFileName;
  std::string      m_ReferenceFileName;  
  std::string      m_TargetFileName;
  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;
  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_Alpha; // difference parameter 
  Float     m_EnergyReductionFactor; 
  Float     m_Temp;

  bool  m_WriteDisplacementField;
  bool  m_DoMultiRes;
  bool  m_UseLandmarks;
  bool  m_ReadMeshFile;
  bool  m_UseMassMatrix;
  Sign  m_DescentDirection;

  ImageSizeType     m_ImageSize; // image size
  ImageSizeType     m_ImageOrigin; // image size
  ImageSizeType     m_ImageScaling; 
  typename ImageType::RegionType   m_FieldRegion;
  typename ImageType::SizeType     m_FieldSize;
  typename FieldType::Pointer      m_Field;

  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 ImageType::Pointer      m_WarpedImage;
  typename FloatImageType::Pointer      m_FloatImage;
  typename ImageType::RegionType   m_Wregion; 
  typename ImageType::IndexType    m_Windex;
 
 // declare images for target and reference
  typename ImageType::Pointer      m_RefImg;
  typename ImageType::Pointer      m_TarImg;
  typename ImageType::RegionType   m_Rregion;
  typename ImageType::RegionType   m_Tregion;
  typename ImageType::IndexType    m_Rindex;
  typename ImageType::IndexType    m_Tindex;

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

};

}} // end namespace itk::fem

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

#endif

---