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: 2010-03-02 03:40:35 $
00007   Version:   $Revision: 1.29 $
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 
00216   void      SetMovingFile(const char* r)
00217     {
00218     m_MovingFileName=r;
00219     }
00220 
00221   std::string GetMovingFile()
00222     {
00223     return m_MovingFileName;
00224     }
00225   
00227   void      SetFixedFile(const char* t) {m_FixedFileName=t;}
00228 
00229   std::string GetFixedFile() {return m_FixedFileName;}
00230 
00231   
00235   void SetMovingImage(MovingImageType* R);
00236 
00238   void SetFixedImage(FixedImageType* T);
00239 
00240   MovingImageType* GetMovingImage(){return m_MovingImage;}
00241   MovingImageType* GetOriginalMovingImage(){return m_OriginalMovingImage;}
00242   
00243   FixedImageType* GetFixedImage(){return m_FixedImage;}
00244   
00245   
00248   FixedImageType* GetWarpedImage(){return m_WarpedImage;}
00249 
00251   void ComputeJacobian(float sign=1.0, FieldType* field=NULL , float smooth=0.0);
00252 
00254   FloatImageType* GetJacobianImage(){return m_FloatImage;}
00255 
00256 
00259   FieldType* GetDeformationField(){return m_Field;}
00260 
00262   void SetDeformationField(FieldType* F)
00263     {  
00264     m_FieldSize=F->GetLargestPossibleRegion().GetSize();
00265     m_Field=F;
00266     }
00268 
00271   void      SetLandmarkFile(const char* l)
00272     {
00273     m_LandmarkFileName=l;
00274     }
00275 
00277   void      UseLandmarks(bool b)
00278     {
00279     m_UseLandmarks=b;
00280     }
00281 
00289   void      EnforceDiffeomorphism(float thresh , SolverType& S ,  bool onlywriteimages);
00290 
00291 
00296   void      SetResultsFile(const char* r)
00297     {
00298     m_ResultsFileName=r;
00299     }
00300 
00301   void      SetResultsFileName (const char* f)
00302     {
00303     m_ResultsFileName=f;
00304     }
00305 
00306   std::string GetResultsFileName ()
00307     {
00308     return m_ResultsFileName;
00309     } 
00310 
00312   void      SetDisplacementsFile(const char* r)
00313     {
00314     m_DisplacementsFileName=r;
00315     }
00316 
00323   void      SetMeshPixelsPerElementAtEachResolution(unsigned int i,unsigned int which=0)
00324     {
00325     m_MeshPixelsPerElementAtEachResolution[which]=i;
00326     }
00328 
00332   void      SetNumberOfIntegrationPoints(unsigned int i,unsigned int which=0)
00333     {
00334     m_NumberOfIntegrationPoints[which]=i;
00335     }
00336 
00342   void      SetWidthOfMetricRegion(unsigned int i,unsigned int which=0)
00343     {
00344     m_MetricWidth[which]=i;
00345     }
00346   unsigned int      GetWidthOfMetricRegion(unsigned int which=0)
00347     {
00348     return m_MetricWidth[which];
00349     }
00351 
00356   void      SetMaximumIterations(unsigned int i,unsigned int which)
00357     {
00358     m_Maxiters[which]=i;
00359     }
00360 
00363   void      SetTimeStep(Float i)
00364     {
00365     m_Dt=i;
00366     }
00367 
00369   void      SetAlpha(Float a)
00370     {
00371     m_Alpha=a;
00372     }
00373 
00376   void      SetEnergyReductionFactor(Float i)
00377     {
00378     m_EnergyReductionFactor=i;
00379     }
00380 
00382   void      SetElasticity(Float i,unsigned int which=0)
00383     {
00384     m_E[which]=i;
00385     } 
00386 
00388   Float     GetElasticity(unsigned int which=0)
00389     {
00390      return m_E[which];
00391     }
00392 
00394   void      SetRho(Float r,unsigned int which=0)
00395     {
00396     m_Rho[which]=r;
00397     } 
00398 
00400   void      SetGamma(Float r,unsigned int which=0)
00401     {
00402     m_Gamma[which]=r;
00403     } 
00404 
00406   void      SetDescentDirectionMinimize()
00407     {
00408     m_DescentDirection=positive;
00409     }
00410 
00412   void      SetDescentDirectionMaximize()
00413     {
00414     m_DescentDirection=negative;
00415     } 
00416 
00419   void      DoLineSearch(unsigned int b)
00420     {
00421     m_DoLineSearchOnImageEnergy=b;
00422     } 
00423 
00425   void      DoMultiRes(bool b)
00426     {
00427     m_DoMultiRes=b;
00428     } 
00429 
00431   void      EmployRegridding(unsigned int b)
00432     {
00433     m_EmployRegridding=b;
00434     } 
00435 
00437   void      SetLineSearchMaximumIterations(unsigned int f)
00438     {
00439     m_LineSearchMaximumIterations=f;
00440     } 
00441 
00443   void      SetWriteDisplacements(bool b)
00444     {
00445     m_WriteDisplacementField=b;
00446     }
00447 
00449   bool      GetWriteDisplacements()
00450     {
00451     return m_WriteDisplacementField;
00452     }
00453 
00456   void      SetConfigFileName (const char* f){m_ConfigFileName=f;}
00457 
00458   std::string GetConfigFileName () {return m_ConfigFileName; }
00459   
00460   ImageSizeType GetImageSize(){ return m_FullImageSize; }
00461 
00463   MetricBaseTypePointer    GetMetric() { return m_Metric; }
00464   void      SetMetric(MetricBaseTypePointer MP) { m_Metric=MP; }
00466 
00469   void      ChooseMetric( float whichmetric); 
00470 
00472   void      SetElement(Element::Pointer e) {m_Element=e;}
00473 
00475   void      SetMaterial(MaterialType::Pointer m) {m_Material=m;}
00476 
00477   void      PrintVectorField(unsigned int modnum=1000);
00478 
00479   void      SetNumLevels(unsigned int i) { m_NumLevels=i; }
00480   void      SetMaxLevel(unsigned int i) { m_MaxLevel=i; }
00481 
00482   void      SetTemp(Float i) { m_Temp=i; }
00483 
00485   FEMRegistrationFilter( ); 
00486   ~FEMRegistrationFilter(); 
00488 
00489 // HELPER FUNCTIONS
00490 protected:
00491 
00492   
00499   class FEMOF : public FEMObjectFactory<FEMLightObject>{
00500   protected:
00501     FEMOF();
00502     ~FEMOF();
00503     };
00505 
00506  
00508   void      CreateMesh(double ElementsPerSide, Solver& S, ImageSizeType sz);
00509 
00511   void      ApplyLoads(SolverType& S,ImageSizeType Isz,double* spacing=NULL); 
00512 
00514   void      ApplyImageLoads(SolverType& S, MovingImageType* i1, FixedImageType* i2); 
00515 
00516   
00519   void      CreateLinearSystemSolver();
00520 
00521   
00523   Float     EvaluateEnergy();
00524 
00529   void      InterpolateVectorField(SolverType& S); 
00530 
00533   FloatImageType*      GetMetricImage(FieldType* F); 
00534 
00535  
00537   typedef  typename FieldType::Pointer FieldPointer;
00538   FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);
00540 
00541 
00543   void      SampleVectorFieldAtNodes(SolverType& S);
00544 
00545   Float EvaluateResidual(SolverType& mySolver,Float t);
00546 
00547   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00548   void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
00549   
00553   Float GoldenSection(SolverType& mySolver,Float tol=0.01,unsigned int MaxIters=25);
00554 
00556 //  itkSetMacro( Load, ImageMetricLoadType* );
00557   itkGetConstMacro( Load, ImageMetricLoadType* );
00559 
00560 
00561   void PrintSelf(std::ostream& os, Indent indent) const;
00562 
00563 private:
00564 
00565   void InitializeField();
00566 
00567   FEMRegistrationFilter(const Self&); //purposely not implemented
00568   void operator=(const Self&); //purposely not implemented
00569     
00570   std::string      m_ConfigFileName;
00571   std::string      m_ResultsFileName;
00572   std::string      m_MovingFileName;   // This variable is currently not being used.
00573   std::string      m_FixedFileName;    // This variable is currently not being used.
00574   std::string      m_LandmarkFileName;
00575   std::string      m_DisplacementsFileName;
00576   std::string      m_MeshFileName;
00577 
00578   unsigned int     m_DoLineSearchOnImageEnergy; 
00579   unsigned int     m_LineSearchMaximumIterations;
00580 
00581   vnl_vector<unsigned int>     m_NumberOfIntegrationPoints;// resolution of integration
00582   vnl_vector<unsigned int>     m_MetricWidth;
00583   vnl_vector<unsigned int>     m_Maxiters; // max iterations
00584   unsigned int                 m_TotalIterations;
00585   unsigned int                 m_NumLevels; // Number of Resolution Levels
00586   unsigned int                 m_MaxLevel;  // Maximum Level (NumLevels is original resolution).
00587   unsigned int                 m_MeshLevels;// Number of Mesh Resolutions ( should be >= 1)
00588   unsigned int                 m_MeshStep;  // Ratio Between Mesh Resolutions ( currently set to 2, should be >= 1)
00589   unsigned int                 m_FileCount; // keeps track of number of files written
00590   unsigned int                 m_CurrentLevel;
00591 
00592   typename FixedImageType::SizeType     m_CurrentLevelImageSize; 
00593 
00594   unsigned int     m_WhichMetric;
00595  
00598   vnl_vector<unsigned int> m_MeshPixelsPerElementAtEachResolution;
00599 
00600   Float                 m_Dt; // time step
00601   vnl_vector<Float>     m_E;  // elasticity 
00602   vnl_vector<Float>     m_Rho;   // mass matrix weight
00603   vnl_vector<Float>     m_Gamma;   // image similarity weight
00604   Float                 m_Energy; // current value of energy
00605   Float                 m_MinE;  // minimum recorded energy
00606   Float                 m_MinJacobian;  // minimum recorded energy
00607   Float                 m_Alpha; // difference parameter 
00609   Float     m_EnergyReductionFactor; 
00610   Float     m_Temp;
00611 
00612   bool         m_WriteDisplacementField;
00613   bool         m_DoMultiRes;
00614   bool         m_UseLandmarks;
00615   bool         m_ReadMeshFile;
00616   bool         m_UseMassMatrix;
00617   unsigned int m_EmployRegridding;
00618   Sign         m_DescentDirection;
00619 
00620   ImageSizeType     m_FullImageSize; // image size
00621   ImageSizeType     m_ImageOrigin; // image size
00623   ImageSizeType                    m_ImageScaling; 
00624   ImageSizeType                    m_CurrentImageScaling; 
00625   typename FieldType::RegionType   m_FieldRegion;
00626   typename FieldType::SizeType     m_FieldSize;
00627   typename FieldType::Pointer      m_Field;
00628   // only use TotalField if re-gridding is employed.
00629   typename FieldType::Pointer      m_TotalField;
00630   ImageMetricLoadType*             m_Load; // Defines the load to use
00631 
00632   // define the warper
00633   typename WarperType::Pointer m_Warper; 
00634 
00635   // declare a new image to hold the warped  reference
00636   typename FixedImageType::Pointer      m_WarpedImage;
00637   typename FloatImageType::Pointer      m_FloatImage;
00638   typename FixedImageType::RegionType   m_Wregion; 
00639   typename FixedImageType::IndexType    m_Windex;
00640  
00641   // declare images for target and reference
00642   typename MovingImageType::Pointer     m_MovingImage;
00643   typename MovingImageType::Pointer     m_OriginalMovingImage;
00644   typename FixedImageType::Pointer      m_FixedImage;
00645 
00646   // element and metric pointers
00647   typename Element::Pointer        m_Element;
00648   typename MaterialType::Pointer   m_Material;
00649   MetricBaseTypePointer            m_Metric;
00650 
00651   LandmarkArrayType      m_LandmarkArray;
00652   InterpolatorPointer    m_Interpolator;
00653 
00654  
00655 
00656 };
00657 
00658 }} // end namespace itk::fem
00659 
00660 #ifndef ITK_MANUAL_INSTANTIATION
00661 #include "itkFEMRegistrationFilter.txx"
00662 #endif
00663 
00664 #endif
00665 

Generated at Fri Apr 16 18:19:22 2010 for ITK by doxygen 1.6.1 written by Dimitri van Heesch, © 1997-2000