Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkFEMRegistrationFilter.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkFEMRegistrationFilter.h,v $
00005   Language:  C++
00006   Date:      $Date: 2006/11/07 23:23:15 $
00007   Version:   $Revision: 1.24 $
00008 
00009   Copyright (c) Insight Software Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012      This software is distributed WITHOUT ANY WARRANTY; without even 
00013      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00014      PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #ifndef _itkFEMRegistrationFilter_h_
00019 #define _itkFEMRegistrationFilter_h_
00020 
00021 #include "itkFEMLinearSystemWrapperItpack.h"
00022 #include "itkFEMLinearSystemWrapperDenseVNL.h"
00023 #include "itkFEMGenerateMesh.h"
00024 #include "itkFEMSolverCrankNicolson.h"
00025 #include "itkFEMMaterialLinearElasticity.h"
00026 #include "itkFEMImageMetricLoad.h"
00027 #include "itkFEMFiniteDifferenceFunctionLoad.h"
00028 
00029 #include "itkImage.h"
00030 #include "itkVector.h"
00031 #include "itkImageRegionIteratorWithIndex.h"
00032 #include "itkVectorCastImageFilter.h"
00033 #include "itkVectorIndexSelectionCastImageFilter.h"
00034 #include "itkWarpImageFilter.h"
00035 #include "itkImageToImageMetric.h"
00036 #include "itkTranslationTransform.h"
00037 #include "itkVectorExpandImageFilter.h"
00038 
00039 #include "itkRecursiveMultiResolutionPyramidImageFilter.h"
00040 #include "itkFEMLoadLandmark.h"
00041 
00042 #include "vnl/vnl_vector.h"
00043 #include "vnl/vnl_math.h"
00044 #include "vnl/vnl_vector_fixed.h"
00045 
00046 
00047 #include <iostream>
00048 #include <string>
00049 
00050 
00051 
00052 namespace itk {
00053 namespace fem {
00054 
00101 template<class TMovingImage,class TFixedImage> 
00102 class  ITK_EXPORT  FEMRegistrationFilter : public ImageToImageFilter<TMovingImage, TFixedImage>
00103 {
00104 public:
00105   typedef FEMRegistrationFilter                              Self;
00106   typedef ImageToImageFilter<TMovingImage, TFixedImage> Superclass;
00107   typedef SmartPointer<Self> Pointer;
00108   typedef SmartPointer<const Self> ConstPointer;
00109 
00111   itkNewMacro(Self);
00112 
00114   itkTypeMacro(FEMRegistrationFilter, ImageToImageFilter );
00115 
00116   typedef TMovingImage                              MovingImageType;
00117   typedef TFixedImage                               FixedImageType;
00118   typedef typename FixedImageType::PixelType             PixelType;
00119   typedef typename FixedImageType::SizeType              ImageSizeType;
00120 
00122   itkStaticConstMacro(ImageDimension, unsigned int,
00123                       FixedImageType::ImageDimension);
00124 
00125   typedef Image< float, itkGetStaticConstMacro(ImageDimension) >            FloatImageType;
00126   typedef LinearSystemWrapperItpack                 LinearSystemSolverType;
00127   typedef SolverCrankNicolson                       SolverType;
00128   enum Sign { positive = 1, negative = -1 };
00129   typedef double                                    Float;
00130   typedef Load::ArrayType LoadArray;
00131     
00132 
00133   typedef std::vector<typename LoadLandmark::Pointer> LandmarkArrayType;
00134   typedef itk::Vector<float,itkGetStaticConstMacro(ImageDimension)>         VectorType;
00135   typedef itk::Image<VectorType,itkGetStaticConstMacro(ImageDimension)>     FieldType;
00136   typedef itk::WarpImageFilter<MovingImageType,FixedImageType, FieldType> WarperType;
00137   typedef MaterialLinearElasticity                  MaterialType;
00138   typedef itk::ImageRegionIteratorWithIndex<FixedImageType>         ImageIterator; 
00139   typedef itk::ImageRegionIteratorWithIndex<FloatImageType>         FloatImageIterator; 
00140   typedef itk::ImageRegionIteratorWithIndex<FieldType>         FieldIterator; 
00141   typedef itk::VectorIndexSelectionCastImageFilter<FieldType,FloatImageType> IndexSelectCasterType;
00142 
00143 
00145   typedef double CoordRepType;
00146   typedef VectorInterpolateImageFunction<FieldType,CoordRepType> 
00147     InterpolatorType;
00148   typedef typename InterpolatorType::Pointer InterpolatorPointer;
00149   typedef VectorLinearInterpolateImageFunction<FieldType,CoordRepType> 
00150     DefaultInterpolatorType;
00151 
00153   itkSetObjectMacro( Interpolator, InterpolatorType );
00154 
00156   itkGetObjectMacro( Interpolator, InterpolatorType );
00157 
00158 
00159   typedef itk::VectorExpandImageFilter<FieldType,FieldType> ExpanderType;
00160   typedef typename ExpanderType::ExpandFactorsType ExpandFactorsType;
00161 
00162   typedef itk::RecursiveMultiResolutionPyramidImageFilter<FixedImageType,FixedImageType>
00163     FixedPyramidType;
00164 
00166 //#define USEIMAGEMETRIC
00167 #ifdef  USEIMAGEMETRIC
00168   typedef ImageToImageMetric<ImageType,FixedImageType>   MetricBaseType;
00169   typedef  ImageMetricLoad<ImageType,ImageType>  ImageMetricLoadType;
00170 #else
00171   typedef  FiniteDifferenceFunctionLoad<MovingImageType,FixedImageType>  ImageMetricLoadType;
00172   typedef PDEDeformableRegistrationFunction<FixedImageType,MovingImageType,FieldType>   MetricBaseType;
00173 #endif
00174   typedef typename MetricBaseType::Pointer          MetricBaseTypePointer;
00175   /* Main functions */
00176 
00178   bool      ReadConfigFile(const char*);
00179 
00180 
00182   void      RunRegistration(void); 
00183 
00187   void      WriteWarpedImage(const char* fn);
00188 
00189 
00191   void      IterativeSolve(SolverType& S);  
00192 
00194   void      MultiResSolve();
00195 
00197   void      WarpImage(const MovingImageType * R);      
00198 
00200   int       WriteDisplacementField(unsigned int index);
00201 
00203   int       WriteDisplacementFieldMultiComponent();
00204 
00206   void      SetMovingFile(const char* r) {m_MovingFileName=r;}
00207 
00208   std::string GetMovingFile() {return m_MovingFileName;}
00209   
00210   void      SetFixedFile(const char* t) {m_FixedFileName=t;}
00211   
00212   std::string GetFixedFile() {return m_FixedFileName;}
00213 
00214   
00218   void SetMovingImage(MovingImageType* R);
00219 
00221   void SetFixedImage(FixedImageType* T);
00222 
00223   MovingImageType* GetMovingImage(){return m_MovingImage;}
00224   MovingImageType* GetOriginalMovingImage(){return m_OriginalMovingImage;}
00225   
00226   FixedImageType* GetFixedImage(){return m_FixedImage;}
00227   
00228   
00231   FixedImageType* GetWarpedImage(){return m_WarpedImage;}
00232 
00234   void ComputeJacobian(float sign=1.0, FieldType* field=NULL , float smooth=0.0);
00235 
00237   FloatImageType* GetJacobianImage(){return m_FloatImage;}
00238 
00239 
00241   FieldType* GetDeformationField(){return m_Field;}
00242 
00244   void SetDeformationField(FieldType* F)
00245   {  
00246     m_FieldSize=F->GetLargestPossibleRegion().GetSize();
00247     m_Field=F;
00248   }
00250 
00253   void      SetLandmarkFile(const char* l) {m_LandmarkFileName=l; }
00254 
00256   void      UseLandmarks(bool b) {m_UseLandmarks=b;}
00257 
00258 
00266   void      EnforceDiffeomorphism(float thresh , SolverType& S ,  bool onlywriteimages);
00267 
00268 
00273   void      SetResultsFile(const char* r) {m_ResultsFileName=r;}
00274 
00275   void      SetResultsFileName (const char* f){m_ResultsFileName=f;}
00276 
00277   std::string GetResultsFileName () {return m_ResultsFileName;} 
00278 
00280   void      SetDisplacementsFile(const char* r) {m_DisplacementsFileName=r;}
00281 
00288   void      SetMeshPixelsPerElementAtEachResolution(unsigned int i,unsigned int which=0){ m_MeshPixelsPerElementAtEachResolution[which]=i;}
00290 
00294   void      SetNumberOfIntegrationPoints(unsigned int i,unsigned int which=0){ m_NumberOfIntegrationPoints[which]=i;}
00295 
00301   void      SetWidthOfMetricRegion(unsigned int i,unsigned int which=0) { m_MetricWidth[which]=i;}
00302   unsigned int      GetWidthOfMetricRegion(unsigned int which=0) { return m_MetricWidth[which];}
00304 
00309   void      SetMaximumIterations(unsigned int i,unsigned int which) { m_Maxiters[which]=i;}
00310 
00313   void      SetTimeStep(Float i) { m_dT=i;}
00314 
00316   void      SetAlpha(Float a) { m_Alpha=a;}
00317 
00320   void      SetEnergyReductionFactor(Float i) { m_EnergyReductionFactor=i;}
00321 
00323   void      SetElasticity(Float i,unsigned int which=0) { m_E[which]=i;} 
00324 
00326   Float     GetElasticity(unsigned int which=0) { return m_E[which];}
00327 
00329   void      SetRho(Float r,unsigned int which=0) { m_Rho[which]=r;} 
00330 
00332   void      SetGamma(Float r,unsigned int which=0) { m_Gamma[which]=r;} 
00333 
00335   void      SetDescentDirectionMinimize() { m_DescentDirection=positive;}
00336 
00338   void      SetDescentDirectionMaximize() { m_DescentDirection=negative;} 
00339 
00341   void      DoLineSearch(unsigned int b) { m_DoLineSearchOnImageEnergy=b; } 
00342 
00343 
00345   void      DoMultiRes(bool b) { m_DoMultiRes=b; } 
00346 
00348   void      EmployRegridding(unsigned int b) { m_EmployRegridding=b; } 
00349 
00351   void      SetLineSearchMaximumIterations(unsigned int f) { m_LineSearchMaximumIterations=f; } 
00352 
00354   void      SetWriteDisplacements(bool b) {m_WriteDisplacementField=b;}
00355 
00357   bool      GetWriteDisplacements() {return m_WriteDisplacementField;}
00358 
00361   void      SetConfigFileName (const char* f){m_ConfigFileName=f;}
00362 
00363   std::string GetConfigFileName () {return m_ConfigFileName; }
00364   
00365   ImageSizeType GetImageSize(){ return m_FullImageSize; }
00366 
00368   MetricBaseTypePointer    GetMetric() { return m_Metric; }
00369   void      SetMetric(MetricBaseTypePointer MP) { m_Metric=MP; }
00371 
00374   void      ChooseMetric( float whichmetric); 
00375 
00377   void      SetElement(Element::Pointer e) {m_Element=e;}
00378 
00380   void      SetMaterial(MaterialType::Pointer m) {m_Material=m;}
00381 
00382   void      PrintVectorField(unsigned int modnum=1000);
00383 
00384   void      SetNumLevels(unsigned int i) { m_NumLevels=i; }
00385   void      SetMaxLevel(unsigned int i) { m_MaxLevel=i; }
00386 
00387   void      SetTemp(Float i) { m_Temp=i; }
00388 
00390   FEMRegistrationFilter( ); 
00391   ~FEMRegistrationFilter(); 
00393 
00394 // HELPER FUNCTIONS
00395 protected :
00396 
00397   
00402   class FEMOF : public FEMObjectFactory<FEMLightObject>{
00403   protected:
00404     FEMOF();
00405     ~FEMOF();
00406     };
00408 
00409  
00411   void      CreateMesh(double ElementsPerSide, Solver& S, ImageSizeType sz);
00412 
00414   void      ApplyLoads(SolverType& S,ImageSizeType Isz,double* spacing=NULL); 
00415 
00417   void      ApplyImageLoads(SolverType& S, MovingImageType* i1, FixedImageType* i2); 
00418 
00419   
00422   void      CreateLinearSystemSolver();
00423 
00424   
00426   Float     EvaluateEnergy();
00427 
00432   void      InterpolateVectorField(SolverType& S); 
00433 
00436   FloatImageType*      GetMetricImage(FieldType* F); 
00437 
00438  
00440   typedef  typename FieldType::Pointer FieldPointer;
00441   FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);
00443 
00444 
00446   void      SampleVectorFieldAtNodes(SolverType& S);
00447 
00448   Float EvaluateResidual(SolverType& mySolver,Float t);
00449 
00450   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00451   void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
00452   
00456   Float GoldenSection(SolverType& mySolver,Float tol=0.01,unsigned int MaxIters=25);
00457 
00459 //  itkSetMacro( Load, ImageMetricLoadType* );
00460   itkGetMacro( Load, ImageMetricLoadType* );
00462 
00463 
00464   void PrintSelf(std::ostream& os, Indent indent) const;
00465 
00466 private :
00467 
00468   void InitializeField();
00469 
00470   FEMRegistrationFilter(const Self&); //purposely not implemented
00471   void operator=(const Self&); //purposely not implemented
00472     
00473   std::string      m_ConfigFileName;
00474   std::string      m_ResultsFileName;
00475   std::string      m_MovingFileName;  
00476   std::string      m_FixedFileName;
00477   std::string      m_LandmarkFileName;
00478   std::string      m_DisplacementsFileName;
00479   std::string      m_MeshFileName;
00480 
00481   unsigned int     m_DoLineSearchOnImageEnergy; 
00482   unsigned int     m_LineSearchMaximumIterations;
00483 
00484   vnl_vector<unsigned int>     m_NumberOfIntegrationPoints;// resolution of integration
00485   vnl_vector<unsigned int>     m_MetricWidth;
00486   vnl_vector<unsigned int>     m_Maxiters; // max iterations
00487   unsigned int     m_TotalIterations;
00488   unsigned int     m_NumLevels; // Number of Resolution Levels
00489   unsigned int     m_MaxLevel;  // Maximum Level (NumLevels is original resolution).
00490   unsigned int     m_MeshLevels;// Number of Mesh Resolutions ( should be >= 1)
00491   unsigned int     m_MeshStep;  // Ratio Between Mesh Resolutions ( currently set to 2, should be >= 1)
00492   unsigned int     m_FileCount; // keeps track of number of files written
00493   unsigned int     m_CurrentLevel;
00494   typename FixedImageType::SizeType     m_CurrentLevelImageSize; 
00495   unsigned int     m_WhichMetric;
00496  
00499   vnl_vector<unsigned int> m_MeshPixelsPerElementAtEachResolution;
00500 
00501   Float     m_dT; // time step
00502   vnl_vector<Float>     m_E;  // elasticity 
00503   vnl_vector<Float>     m_Rho;   // mass matrix weight
00504   vnl_vector<Float>     m_Gamma;   // image similarity weight
00505   Float     m_Energy; // current value of energy
00506   Float     m_MinE;  // minimum recorded energy
00507   Float     m_MinJacobian;  // minimum recorded energy
00508   Float     m_Alpha; // difference parameter 
00510   Float     m_EnergyReductionFactor; 
00511   Float     m_Temp;
00512 
00513   bool  m_WriteDisplacementField;
00514   bool  m_DoMultiRes;
00515   bool  m_UseLandmarks;
00516   bool  m_ReadMeshFile;
00517   bool  m_UseMassMatrix;
00518   unsigned int m_EmployRegridding;
00519   Sign  m_DescentDirection;
00520 
00521   ImageSizeType     m_FullImageSize; // image size
00522   ImageSizeType     m_ImageOrigin; // image size
00524   ImageSizeType     m_ImageScaling; 
00525   ImageSizeType     m_CurrentImageScaling; 
00526   typename FieldType::RegionType   m_FieldRegion;
00527   typename FieldType::SizeType     m_FieldSize;
00528   typename FieldType::Pointer      m_Field;
00529   // only use TotalField if re-gridding is employed.
00530   typename FieldType::Pointer      m_TotalField;
00531   ImageMetricLoadType* m_Load; // Defines the load to use
00532 
00533   // define the warper
00534   typename WarperType::Pointer m_Warper; 
00535 
00536  // declare a new image to hold the warped  reference
00537   typename FixedImageType::Pointer      m_WarpedImage;
00538   typename FloatImageType::Pointer      m_FloatImage;
00539   typename FixedImageType::RegionType   m_Wregion; 
00540   typename FixedImageType::IndexType    m_Windex;
00541  
00542  // declare images for target and reference
00543   typename MovingImageType::Pointer      m_MovingImage;
00544   typename MovingImageType::Pointer      m_OriginalMovingImage;
00545   typename FixedImageType::Pointer      m_FixedImage;
00546 
00547   // element and metric pointers
00548   typename Element::Pointer        m_Element;
00549   typename MaterialType::Pointer   m_Material;
00550   MetricBaseTypePointer            m_Metric;
00551 
00552 
00553   // multi-resolution stuff
00554 //  typename FixedPyramidType::Pointer   m_FixedPyramid;
00555 //  typename FixedPyramidType::Pointer   m_MovingPyramid;
00556   
00557   LandmarkArrayType    m_LandmarkArray;
00558 
00559   InterpolatorPointer    m_Interpolator;
00560 
00561  
00562 
00563 };
00564 
00565 }} // end namespace itk::fem
00566 
00567 #ifndef ITK_MANUAL_INSTANTIATION
00568 #include "itkFEMRegistrationFilter.txx"
00569 #endif
00570 
00571 #endif
00572 
00573 

Generated at Tue Jul 29 20:01:03 2008 for ITK by doxygen 1.5.1 written by Dimitri van Heesch, © 1997-2000