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: 2009-04-23 03:53:35 $
00007   Version:   $Revision: 1.26 $
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 namespace itk {
00051 namespace fem {
00052 
00105 template<class TMovingImage,class TFixedImage> 
00106 class  ITK_EXPORT  FEMRegistrationFilter : public ImageToImageFilter<TMovingImage, TFixedImage>
00107 {
00108 public:
00109   typedef FEMRegistrationFilter                         Self;
00110   typedef ImageToImageFilter<TMovingImage, TFixedImage> Superclass;
00111   typedef SmartPointer<Self>                            Pointer;
00112   typedef SmartPointer<const Self>                      ConstPointer;
00113 
00115   itkNewMacro(Self);
00116 
00118   itkTypeMacro(FEMRegistrationFilter, ImageToImageFilter );
00119 
00120   typedef TMovingImage                        MovingImageType;
00121   typedef TFixedImage                         FixedImageType;
00122   typedef typename FixedImageType::PixelType  PixelType;
00123   typedef typename FixedImageType::SizeType   ImageSizeType;
00124 
00126   itkStaticConstMacro(ImageDimension, unsigned int,
00127                       FixedImageType::ImageDimension);
00128 
00129   typedef Image< float, itkGetStaticConstMacro(ImageDimension) >            FloatImageType;
00130   typedef LinearSystemWrapperItpack            LinearSystemSolverType;
00131   typedef SolverCrankNicolson                  SolverType;
00132   enum Sign { positive = 1, negative = -1 };
00133   typedef double                               Float;
00134   typedef Load::ArrayType                      LoadArray;
00135     
00136 
00137   typedef std::vector<typename LoadLandmark::Pointer>       LandmarkArrayType;
00138   typedef itk::Vector<float,itkGetStaticConstMacro(ImageDimension)>
00139                                                             VectorType;
00140   typedef itk::Image<VectorType,itkGetStaticConstMacro(ImageDimension)>
00141                                                             FieldType;
00142   typedef itk::WarpImageFilter<MovingImageType,FixedImageType, FieldType>
00143                                                             WarperType;
00144   typedef MaterialLinearElasticity                          MaterialType;
00145   typedef itk::ImageRegionIteratorWithIndex<FixedImageType> ImageIterator; 
00146   typedef itk::ImageRegionIteratorWithIndex<FloatImageType> FloatImageIterator; 
00147   typedef itk::ImageRegionIteratorWithIndex<FieldType>      FieldIterator; 
00148   typedef itk::VectorIndexSelectionCastImageFilter<FieldType,FloatImageType>
00149                                                             IndexSelectCasterType;
00150 
00152   typedef double                                            CoordRepType;
00153   typedef VectorInterpolateImageFunction<FieldType,CoordRepType> 
00154                                                             InterpolatorType;
00155   typedef typename InterpolatorType::Pointer                InterpolatorPointer;
00156   typedef VectorLinearInterpolateImageFunction<FieldType,CoordRepType> 
00157                                                             DefaultInterpolatorType;
00158 
00160   itkSetObjectMacro( Interpolator, InterpolatorType );
00161 
00163   itkGetObjectMacro( Interpolator, InterpolatorType );
00164 
00165 
00166   typedef itk::VectorExpandImageFilter<FieldType,FieldType> ExpanderType;
00167   typedef typename ExpanderType::ExpandFactorsType          ExpandFactorsType;
00168 
00169   typedef itk::RecursiveMultiResolutionPyramidImageFilter<FixedImageType,FixedImageType>
00170     FixedPyramidType;
00171 
00173 //#define USEIMAGEMETRIC
00174 #ifdef  USEIMAGEMETRIC
00175   typedef ImageToImageMetric<ImageType,FixedImageType>   MetricBaseType;
00176   typedef  ImageMetricLoad<ImageType,ImageType>          ImageMetricLoadType;
00177 #else
00178   typedef  FiniteDifferenceFunctionLoad<MovingImageType,FixedImageType>
00179                                                          ImageMetricLoadType;
00180   typedef PDEDeformableRegistrationFunction<FixedImageType,MovingImageType,FieldType>
00181                                                          MetricBaseType;
00182 #endif
00183   typedef typename MetricBaseType::Pointer          MetricBaseTypePointer;
00184   /* Main functions */
00185 
00187   bool      ReadConfigFile(const char*);
00188 
00189 
00191   void      RunRegistration(void); 
00192 
00196   void      WriteWarpedImage(const char* fn);
00197 
00198 
00200   void      IterativeSolve(SolverType& S);  
00201 
00203   void      MultiResSolve();
00204 
00206   void      WarpImage(const MovingImageType * R);
00207 
00209   int       WriteDisplacementField(unsigned int index);
00210 
00212   int       WriteDisplacementFieldMultiComponent();
00213 
00215   void      SetMovingFile(const char* r)
00216     {
00217     m_MovingFileName=r;
00218     }
00219 
00220   std::string GetMovingFile()
00221     {
00222     return m_MovingFileName;
00223     }
00224   
00225   void      SetFixedFile(const char* t) {m_FixedFileName=t;}
00226   
00227   std::string GetFixedFile() {return m_FixedFileName;}
00228 
00229   
00233   void SetMovingImage(MovingImageType* R);
00234 
00236   void SetFixedImage(FixedImageType* T);
00237 
00238   MovingImageType* GetMovingImage(){return m_MovingImage;}
00239   MovingImageType* GetOriginalMovingImage(){return m_OriginalMovingImage;}
00240   
00241   FixedImageType* GetFixedImage(){return m_FixedImage;}
00242   
00243   
00246   FixedImageType* GetWarpedImage(){return m_WarpedImage;}
00247 
00249   void ComputeJacobian(float sign=1.0, FieldType* field=NULL , float smooth=0.0);
00250 
00252   FloatImageType* GetJacobianImage(){return m_FloatImage;}
00253 
00254 
00257   FieldType* GetDeformationField(){return m_Field;}
00258 
00260   void SetDeformationField(FieldType* F)
00261     {  
00262     m_FieldSize=F->GetLargestPossibleRegion().GetSize();
00263     m_Field=F;
00264     }
00266 
00269   void      SetLandmarkFile(const char* l)
00270     {
00271     m_LandmarkFileName=l;
00272     }
00273 
00275   void      UseLandmarks(bool b)
00276     {
00277     m_UseLandmarks=b;
00278     }
00279 
00287   void      EnforceDiffeomorphism(float thresh , SolverType& S ,  bool onlywriteimages);
00288 
00289 
00294   void      SetResultsFile(const char* r)
00295     {
00296     m_ResultsFileName=r;
00297     }
00298 
00299   void      SetResultsFileName (const char* f)
00300     {
00301     m_ResultsFileName=f;
00302     }
00303 
00304   std::string GetResultsFileName ()
00305     {
00306     return m_ResultsFileName;
00307     } 
00308 
00310   void      SetDisplacementsFile(const char* r)
00311     {
00312     m_DisplacementsFileName=r;
00313     }
00314 
00321   void      SetMeshPixelsPerElementAtEachResolution(unsigned int i,unsigned int which=0)
00322     {
00323     m_MeshPixelsPerElementAtEachResolution[which]=i;
00324     }
00326 
00330   void      SetNumberOfIntegrationPoints(unsigned int i,unsigned int which=0)
00331     {
00332     m_NumberOfIntegrationPoints[which]=i;
00333     }
00334 
00340   void      SetWidthOfMetricRegion(unsigned int i,unsigned int which=0)
00341     {
00342     m_MetricWidth[which]=i;
00343     }
00344   unsigned int      GetWidthOfMetricRegion(unsigned int which=0)
00345     {
00346     return m_MetricWidth[which];
00347     }
00349 
00354   void      SetMaximumIterations(unsigned int i,unsigned int which)
00355     {
00356     m_Maxiters[which]=i;
00357     }
00358 
00361   void      SetTimeStep(Float i)
00362     {
00363     m_Dt=i;
00364     }
00365 
00367   void      SetAlpha(Float a)
00368     {
00369     m_Alpha=a;
00370     }
00371 
00374   void      SetEnergyReductionFactor(Float i)
00375     {
00376     m_EnergyReductionFactor=i;
00377     }
00378 
00380   void      SetElasticity(Float i,unsigned int which=0)
00381     {
00382     m_E[which]=i;
00383     } 
00384 
00386   Float     GetElasticity(unsigned int which=0)
00387     {
00388      return m_E[which];
00389     }
00390 
00392   void      SetRho(Float r,unsigned int which=0)
00393     {
00394     m_Rho[which]=r;
00395     } 
00396 
00398   void      SetGamma(Float r,unsigned int which=0)
00399     {
00400     m_Gamma[which]=r;
00401     } 
00402 
00404   void      SetDescentDirectionMinimize()
00405     {
00406     m_DescentDirection=positive;
00407     }
00408 
00410   void      SetDescentDirectionMaximize()
00411     {
00412     m_DescentDirection=negative;
00413     } 
00414 
00417   void      DoLineSearch(unsigned int b)
00418     {
00419     m_DoLineSearchOnImageEnergy=b;
00420     } 
00421 
00423   void      DoMultiRes(bool b)
00424     {
00425     m_DoMultiRes=b;
00426     } 
00427 
00429   void      EmployRegridding(unsigned int b)
00430     {
00431     m_EmployRegridding=b;
00432     } 
00433 
00435   void      SetLineSearchMaximumIterations(unsigned int f)
00436     {
00437     m_LineSearchMaximumIterations=f;
00438     } 
00439 
00441   void      SetWriteDisplacements(bool b)
00442     {
00443     m_WriteDisplacementField=b;
00444     }
00445 
00447   bool      GetWriteDisplacements()
00448     {
00449     return m_WriteDisplacementField;
00450     }
00451 
00454   void      SetConfigFileName (const char* f){m_ConfigFileName=f;}
00455 
00456   std::string GetConfigFileName () {return m_ConfigFileName; }
00457   
00458   ImageSizeType GetImageSize(){ return m_FullImageSize; }
00459 
00461   MetricBaseTypePointer    GetMetric() { return m_Metric; }
00462   void      SetMetric(MetricBaseTypePointer MP) { m_Metric=MP; }
00464 
00467   void      ChooseMetric( float whichmetric); 
00468 
00470   void      SetElement(Element::Pointer e) {m_Element=e;}
00471 
00473   void      SetMaterial(MaterialType::Pointer m) {m_Material=m;}
00474 
00475   void      PrintVectorField(unsigned int modnum=1000);
00476 
00477   void      SetNumLevels(unsigned int i) { m_NumLevels=i; }
00478   void      SetMaxLevel(unsigned int i) { m_MaxLevel=i; }
00479 
00480   void      SetTemp(Float i) { m_Temp=i; }
00481 
00483   FEMRegistrationFilter( ); 
00484   ~FEMRegistrationFilter(); 
00486 
00487 // HELPER FUNCTIONS
00488 protected :
00489 
00490   
00495   class FEMOF : public FEMObjectFactory<FEMLightObject>{
00496   protected:
00497     FEMOF();
00498     ~FEMOF();
00499     };
00500 
00501  
00503   void      CreateMesh(double ElementsPerSide, Solver& S, ImageSizeType sz);
00504 
00506   void      ApplyLoads(SolverType& S,ImageSizeType Isz,double* spacing=NULL); 
00507 
00509   void      ApplyImageLoads(SolverType& S, MovingImageType* i1, FixedImageType* i2); 
00510 
00511   
00514   void      CreateLinearSystemSolver();
00515 
00516   
00518   Float     EvaluateEnergy();
00519 
00524   void      InterpolateVectorField(SolverType& S); 
00525 
00528   FloatImageType*      GetMetricImage(FieldType* F); 
00529 
00530  
00532   typedef  typename FieldType::Pointer FieldPointer;
00533   FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);
00535 
00536 
00538   void      SampleVectorFieldAtNodes(SolverType& S);
00539 
00540   Float EvaluateResidual(SolverType& mySolver,Float t);
00541 
00542   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00543   void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
00544   
00548   Float GoldenSection(SolverType& mySolver,Float tol=0.01,unsigned int MaxIters=25);
00549 
00551 //  itkSetMacro( Load, ImageMetricLoadType* );
00552   itkGetConstMacro( Load, ImageMetricLoadType* );
00554 
00555 
00556   void PrintSelf(std::ostream& os, Indent indent) const;
00557 
00558 private :
00559 
00560   void InitializeField();
00561 
00562   FEMRegistrationFilter(const Self&); //purposely not implemented
00563   void operator=(const Self&); //purposely not implemented
00564     
00565   std::string      m_ConfigFileName;
00566   std::string      m_ResultsFileName;
00567   std::string      m_MovingFileName;  
00568   std::string      m_FixedFileName;
00569   std::string      m_LandmarkFileName;
00570   std::string      m_DisplacementsFileName;
00571   std::string      m_MeshFileName;
00572 
00573   unsigned int     m_DoLineSearchOnImageEnergy; 
00574   unsigned int     m_LineSearchMaximumIterations;
00575 
00576   vnl_vector<unsigned int>     m_NumberOfIntegrationPoints;// resolution of integration
00577   vnl_vector<unsigned int>     m_MetricWidth;
00578   vnl_vector<unsigned int>     m_Maxiters; // max iterations
00579   unsigned int                 m_TotalIterations;
00580   unsigned int                 m_NumLevels; // Number of Resolution Levels
00581   unsigned int                 m_MaxLevel;  // Maximum Level (NumLevels is original resolution).
00582   unsigned int                 m_MeshLevels;// Number of Mesh Resolutions ( should be >= 1)
00583   unsigned int                 m_MeshStep;  // Ratio Between Mesh Resolutions ( currently set to 2, should be >= 1)
00584   unsigned int                 m_FileCount; // keeps track of number of files written
00585   unsigned int                 m_CurrentLevel;
00586 
00587   typename FixedImageType::SizeType     m_CurrentLevelImageSize; 
00588 
00589   unsigned int     m_WhichMetric;
00590  
00593   vnl_vector<unsigned int> m_MeshPixelsPerElementAtEachResolution;
00594 
00595   Float                 m_Dt; // time step
00596   vnl_vector<Float>     m_E;  // elasticity 
00597   vnl_vector<Float>     m_Rho;   // mass matrix weight
00598   vnl_vector<Float>     m_Gamma;   // image similarity weight
00599   Float                 m_Energy; // current value of energy
00600   Float                 m_MinE;  // minimum recorded energy
00601   Float                 m_MinJacobian;  // minimum recorded energy
00602   Float                 m_Alpha; // difference parameter 
00604   Float     m_EnergyReductionFactor; 
00605   Float     m_Temp;
00606 
00607   bool         m_WriteDisplacementField;
00608   bool         m_DoMultiRes;
00609   bool         m_UseLandmarks;
00610   bool         m_ReadMeshFile;
00611   bool         m_UseMassMatrix;
00612   unsigned int m_EmployRegridding;
00613   Sign         m_DescentDirection;
00614 
00615   ImageSizeType     m_FullImageSize; // image size
00616   ImageSizeType     m_ImageOrigin; // image size
00618   ImageSizeType                    m_ImageScaling; 
00619   ImageSizeType                    m_CurrentImageScaling; 
00620   typename FieldType::RegionType   m_FieldRegion;
00621   typename FieldType::SizeType     m_FieldSize;
00622   typename FieldType::Pointer      m_Field;
00623   // only use TotalField if re-gridding is employed.
00624   typename FieldType::Pointer      m_TotalField;
00625   ImageMetricLoadType*             m_Load; // Defines the load to use
00626 
00627   // define the warper
00628   typename WarperType::Pointer m_Warper; 
00629 
00630   // declare a new image to hold the warped  reference
00631   typename FixedImageType::Pointer      m_WarpedImage;
00632   typename FloatImageType::Pointer      m_FloatImage;
00633   typename FixedImageType::RegionType   m_Wregion; 
00634   typename FixedImageType::IndexType    m_Windex;
00635  
00636   // declare images for target and reference
00637   typename MovingImageType::Pointer     m_MovingImage;
00638   typename MovingImageType::Pointer     m_OriginalMovingImage;
00639   typename FixedImageType::Pointer      m_FixedImage;
00640 
00641   // element and metric pointers
00642   typename Element::Pointer        m_Element;
00643   typename MaterialType::Pointer   m_Material;
00644   MetricBaseTypePointer            m_Metric;
00645 
00646   LandmarkArrayType      m_LandmarkArray;
00647   InterpolatorPointer    m_Interpolator;
00648 
00649  
00650 
00651 };
00652 
00653 }} // end namespace itk::fem
00654 
00655 #ifndef ITK_MANUAL_INSTANTIATION
00656 #include "itkFEMRegistrationFilter.txx"
00657 #endif
00658 
00659 #endif
00660 

Generated at Tue Sep 15 02:55:35 2009 for ITK by doxygen 1.5.8 written by Dimitri van Heesch, © 1997-2000