ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkFEMRegistrationFilter.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 
00019 #ifndef __itkFEMRegistrationFilter_h
00020 #define __itkFEMRegistrationFilter_h
00021 
00022 #include "itkFEMLinearSystemWrapperItpack.h"
00023 #include "itkFEMLinearSystemWrapperDenseVNL.h"
00024 #include "itkFEMObject.h"
00025 #include "itkFEMSolverCrankNicolson.h"
00026 #include "itkFEMMaterialLinearElasticity.h"
00027 #include "itkFEMImageMetricLoad.h"
00028 #include "itkFEMFiniteDifferenceFunctionLoad.h"
00029 
00030 #include "itkImage.h"
00031 #include "itkVector.h"
00032 #include "itkImageRegionIteratorWithIndex.h"
00033 #include "itkImageToRectilinearFEMObjectFilter.h"
00034 #include "itkVectorCastImageFilter.h"
00035 #include "itkVectorIndexSelectionCastImageFilter.h"
00036 #include "itkWarpImageFilter.h"
00037 #include "itkImageToImageMetric.h"
00038 #include "itkTranslationTransform.h"
00039 #include "itkVectorExpandImageFilter.h"
00040 #include "itkFixedArray.h"
00041 
00042 #include "itkRecursiveMultiResolutionPyramidImageFilter.h"
00043 #include "itkFEMLoadLandmark.h"
00044 
00045 #include "vnl/vnl_vector.h"
00046 #include "vnl/vnl_math.h"
00047 #include "vnl/vnl_vector_fixed.h"
00048 
00049 #include <iostream>
00050 #include <string>
00051 
00052 namespace itk
00053 {
00054 namespace fem
00055 {
00056 
00117 template <class TMovingImage, class TFixedImage, class TFemObjectType>
00118 class ITK_EXPORT  FEMRegistrationFilter : public ImageToImageFilter<TMovingImage, TFixedImage>
00119 {
00120 public:
00121   typedef FEMRegistrationFilter                         Self;
00122   typedef ImageToImageFilter<TMovingImage, TFixedImage> Superclass;
00123   typedef SmartPointer<Self>                            Pointer;
00124   typedef SmartPointer<const Self>                      ConstPointer;
00125 
00127   itkNewMacro(Self);
00128 
00130   itkTypeMacro(FEMRegistrationFilter, ImageToImageFilter );
00131 
00132   typedef TMovingImage                       MovingImageType;
00133   typedef TFixedImage                        FixedImageType;
00134   typedef TFemObjectType                     FEMObjectType;
00135   typedef typename FixedImageType::PixelType PixelType;
00136   typedef typename FixedImageType::SizeType  ImageSizeType;
00137   typedef typename FixedImageType::PointType PointType;
00138 
00140   itkStaticConstMacro(ImageDimension, unsigned int,
00141                       FixedImageType::ImageDimension);
00142 
00143   typedef Image<float, itkGetStaticConstMacro(ImageDimension)> FloatImageType;
00144   typedef LinearSystemWrapperItpack                            LinearSystemSolverType;
00145   typedef SolverCrankNicolson<ImageDimension>                  SolverType;
00146 
00147   enum Sign { positive = 1, negative = -1 };
00148   typedef double          Float;
00149   typedef Load::ArrayType LoadArray;
00150 
00151   typedef std::vector<typename LoadLandmark::Pointer>                      LandmarkArrayType;
00152   typedef itk::Vector<float, itkGetStaticConstMacro(ImageDimension)>       VectorType;
00153   typedef itk::Image<VectorType, itkGetStaticConstMacro(ImageDimension)>   FieldType;
00154   typedef itk::WarpImageFilter<MovingImageType, FixedImageType, FieldType> WarperType;
00155 
00156   typedef MaterialLinearElasticity                          MaterialType;
00157   typedef itk::ImageRegionIteratorWithIndex<FixedImageType> ImageIterator;
00158   typedef itk::ImageRegionIteratorWithIndex<FloatImageType> FloatImageIterator;
00159   typedef itk::ImageRegionIteratorWithIndex<FieldType>      FieldIterator;
00160 
00161   typedef itk::VectorIndexSelectionCastImageFilter<FieldType, FloatImageType> IndexSelectCasterType;
00162 
00164   typedef double                                                  CoordRepType;
00165   typedef VectorInterpolateImageFunction<FieldType, CoordRepType> InterpolatorType;
00166   typedef typename InterpolatorType::Pointer                      InterpolatorPointer;
00167 
00168   typedef VectorLinearInterpolateImageFunction<FieldType, CoordRepType> DefaultInterpolatorType;
00169 
00170   typedef typename itk::Image<Element::ConstPointer, ImageDimension> InterpolationGridType;
00171   typedef typename InterpolationGridType::SizeType                   InterpolationGridSizeType;
00172   typedef typename InterpolationGridType::PointType                  InterpolationGridPointType;
00173 
00174   typedef itk::fem::ImageToRectilinearFEMObjectFilter<TMovingImage> ImageToMeshType;
00175 
00176   typedef itk::VectorExpandImageFilter<FieldType, FieldType> ExpanderType;
00177   typedef typename ExpanderType::ExpandFactorsType           ExpandFactorsType;
00178 
00179   typedef  typename FieldType::Pointer FieldPointer;
00180 
00182   typedef FiniteDifferenceFunctionLoad<MovingImageType, FixedImageType>
00183   ImageMetricLoadType;
00184   typedef PDEDeformableRegistrationFunction<FixedImageType, MovingImageType, FieldType>
00185   MetricBaseType;
00186 
00187   typedef typename MetricBaseType::Pointer     MetricBaseTypePointer;
00188   typedef FixedArray< double, ImageDimension > StandardDeviationsType;
00189 
00190   /*----------------------  Main functions ----------------------*/
00191 
00193   void SetMovingImage(MovingImageType* R);
00194 
00200   MovingImageType * GetMovingImage()
00201   {
00202     return m_MovingImage;
00203   }
00204 
00206   MovingImageType * GetOriginalMovingImage()
00207   {
00208     return m_OriginalMovingImage;
00209   }
00210 
00212   void SetFixedImage(FixedImageType* T);
00213 
00214   FixedImageType * GetFixedImage()
00215   {
00216     return m_FixedImage;
00217   }
00218 
00224   void SetInputFEMObject(FEMObjectType* F, unsigned int level = 0);
00225 
00226   FEMObjectType * GetInputFEMObject(unsigned int level = 0);
00227 
00229   void RunRegistration(void);
00230 
00232   void IterativeSolve(SolverType *S);
00233 
00235   void MultiResSolve();
00236 
00238   void WarpImage(const MovingImageType * R);
00239 
00242   FixedImageType * GetWarpedImage()
00243   {
00244     return m_WarpedImage;
00245   }
00246 
00248   void ComputeJacobian( );
00249 
00251   FloatImageType * GetJacobianImage()
00252   {
00253     return m_FloatImage;
00254   }
00255 
00257   FieldType * GetDisplacementField()
00258   {
00259     return m_Field;
00260   }
00261 
00263   void SetDisplacementField(FieldType* F)
00264   {
00265     m_FieldSize = F->GetLargestPossibleRegion().GetSize();
00266     m_Field = F;
00267   }
00269 
00271   void AddLandmark(PointType source, PointType target);
00272 
00273   void InsertLandmark(unsigned int i, PointType source, PointType target);
00274 
00275   void DeleteLandmark(unsigned int i);
00276 
00277   void ClearLandmarks();
00278 
00279   void GetLandmark(unsigned int i, PointType& source, PointType& target);
00280 
00288   void EnforceDiffeomorphism(float thresh, SolverType *S,  bool onlywriteimages);
00289 
00296   void SetMeshPixelsPerElementAtEachResolution(unsigned int i, unsigned int which = 0)
00297   {
00298     m_MeshPixelsPerElementAtEachResolution[which] = i;
00299   }
00301 
00305   void SetNumberOfIntegrationPoints(unsigned int i, unsigned int which = 0)
00306   {
00307     m_NumberOfIntegrationPoints[which] = i;
00308   }
00309 
00315   void SetWidthOfMetricRegion(unsigned int i, unsigned int which = 0)
00316   {
00317     m_MetricWidth[which] = i;
00318   }
00319 
00320   unsigned int GetWidthOfMetricRegion(unsigned int which = 0)
00321   {
00322     return m_MetricWidth[which];
00323   }
00324 
00329   void SetMaximumIterations(unsigned int i, unsigned int which)
00330   {
00331     m_Maxiters[which] = i;
00332   }
00333 
00339   itkSetMacro(TimeStep, Float);
00340   itkGetMacro(TimeStep, Float);
00342 
00347   itkSetMacro(Alpha, Float);
00348   itkGetMacro(Alpha, Float);
00350 
00354   itkSetMacro(UseLandmarks, bool);
00355   itkGetMacro(UseLandmarks, bool);
00356   void SetUseLandmarksOff()
00357   {
00358     SetUseLandmarks(false);
00359   }
00361 
00362   void SetUseLandmarksOn()
00363   {
00364     SetUseLandmarks(true);
00365   }
00366 
00370   itkSetMacro(UseMassMatrix, bool);
00371   itkGetMacro(UseMassMatrix, bool);
00373 
00377   itkSetMacro(EnergyReductionFactor, Float);
00378   itkGetMacro(EnergyReductionFactor, Float);
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 
00421   itkSetMacro(DoLineSearchOnImageEnergy, unsigned int);
00422   itkGetMacro(DoLineSearchOnImageEnergy, unsigned int);
00424 
00425 
00430   itkSetMacro(UseNormalizedGradient, bool);
00431   itkGetMacro(UseNormalizedGradient, bool);
00432   void SetUseNormalizedGradientOff()
00433   {
00434     SetUseNormalizedGradient(false);
00435   }
00437 
00438   void SetUseNormalizedGradientOn()
00439   {
00440     SetUseNormalizedGradient(true);
00441   }
00442 
00446   itkSetMacro(EmployRegridding, unsigned int);
00447   itkGetMacro(EmployRegridding, unsigned int);
00449 
00454   itkSetMacro(LineSearchMaximumIterations, unsigned int);
00455   itkGetMacro(LineSearchMaximumIterations, unsigned int);
00457 
00461   ImageSizeType GetImageSize()
00462   {
00463     return m_FullImageSize;
00464   }
00465 
00470   MetricBaseTypePointer GetMetric()
00471   {
00472     return m_Metric;
00473   }
00474 
00475   void SetMetric(MetricBaseTypePointer MP)
00476   {
00477     m_Metric = MP;
00478   }
00479 
00488   void ChooseMetric( unsigned int whichmetric);
00489 
00493   unsigned int GetMetricType()
00494   {
00495     return m_WhichMetric;
00496   }
00497 
00499   void SetElement(Element::Pointer e)
00500   {
00501     m_Element = e;
00502   }
00503 
00505   void SetMaterial(MaterialType::Pointer m)
00506   {
00507     m_Material = m;
00508   }
00509 
00511   void PrintVectorField(unsigned int modnum = 1000);
00512 
00516   void SetMaxLevel(unsigned int level);
00517   itkGetMacro(MaxLevel, unsigned int);
00519 
00524   itkSetMacro(CreateMeshFromImage, bool);
00525   void SetCreateMeshFromImageOn()
00526   {
00527     SetCreateMeshFromImage(true);
00528   }
00529   void SetCreateMeshFromImageOff()
00530   {
00531     SetCreateMeshFromImage(false);
00532   }
00533   itkGetMacro(CreateMeshFromImage, bool);
00535 
00537   itkSetObjectMacro( Interpolator, InterpolatorType );
00538 
00540   itkGetObjectMacro( Interpolator, InterpolatorType );
00541 
00545   itkSetMacro(StandardDeviations, StandardDeviationsType);
00546   virtual void SetStandardDeviations(double value);
00548 
00551   itkGetConstReferenceMacro(StandardDeviations, StandardDeviationsType);
00552 
00555   itkSetMacro(MaximumKernelWidth, unsigned int);
00556   itkGetConstMacro(MaximumKernelWidth, unsigned int);
00558 
00561   itkSetMacro(MaximumError, double);
00562   itkGetConstMacro(MaximumError, double);
00564 
00565 
00566 // HELPER FUNCTIONS
00567 protected:
00568   FEMRegistrationFilter();
00569   ~FEMRegistrationFilter();
00570 
00571   void PrintSelf(std::ostream & os, Indent indent) const;
00572 
00574   void CreateMesh(unsigned int ElementsPerSide, SolverType *solver);
00575 
00577   void ApplyLoads(ImageSizeType Isz, double* spacing = NULL);
00578 
00580   void ApplyImageLoads(MovingImageType* i1, FixedImageType* i2);
00581 
00582   // FIXME - Not implemented
00585   void  CreateLinearSystemSolver();
00586 
00587   // FIXME - Not implemented
00589   Float EvaluateEnergy();
00590 
00595   void InterpolateVectorField(SolverType *S);
00596 
00597   // FIXME - Not implemented
00600   FloatImageType * GetMetricImage(FieldType* F);
00601 
00603   FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);
00604 
00606   void SampleVectorFieldAtNodes(SolverType *S);
00607 
00609   Float EvaluateResidual(SolverType *mySolver, Float t);
00610 
00611   // FIXME - Replace with BSD Code
00612   /* Finds a triplet that brackets the energy minimum.  From Numerical Recipes.*/
00613   // void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
00614   void FindBracketingTriplet(SolverType *mySolver, Float* a, Float* b, Float* c);
00615 
00620   Float GoldenSection(SolverType *mySolver, Float tol = 0.01, unsigned int MaxIters = 25);
00621 
00623   itkGetConstObjectMacro( Load, ImageMetricLoadType );
00624   itkSetObjectMacro( Load, ImageMetricLoadType );
00626 
00628   void SmoothDisplacementField();
00629 
00630 private:
00631 
00632   void InitializeField();
00633 
00634   FEMRegistrationFilter(const Self &); // purposely not implemented
00635   void operator=(const Self &);        // purposely not implemented
00636 
00637   unsigned int m_DoLineSearchOnImageEnergy;
00638   unsigned int m_LineSearchMaximumIterations;
00639 
00640   /* Parameters used to define Multi-resolution Registration */
00641   vnl_vector<unsigned int> m_NumberOfIntegrationPoints; // resolution of integration
00642   vnl_vector<unsigned int> m_MetricWidth;               // number of iterations at each level
00643   vnl_vector<unsigned int> m_Maxiters;                  // max iterations
00644   unsigned int             m_TotalIterations;           // total number of iterations that were run
00645   unsigned int             m_MaxLevel;                  // Number of Resolution Levels for registration.
00646   unsigned int             m_FileCount;                 // keeps track of number of files written
00647   unsigned int             m_CurrentLevel;              // current resolution level
00648 
00649   typename FixedImageType::SizeType     m_CurrentLevelImageSize;
00650 
00651   unsigned int m_WhichMetric;                 // Metric used for image registration
00652                                               // 0 = Mean Squares
00653                                               // 1 = Normalized Correlation
00654                                               // 2 = Mutual Information
00655                                               // 3 = Demons Metric
00656 
00659   vnl_vector<unsigned int> m_MeshPixelsPerElementAtEachResolution;
00660 
00661   Float             m_TimeStep;        // time step
00662   vnl_vector<Float> m_E;               // elasticity
00663   vnl_vector<Float> m_Rho;             // mass matrix weight
00664   vnl_vector<Float> m_Gamma;           // image similarity weight
00665   Float             m_Energy;          // current value of energy
00666   Float             m_MinE;            // minimum recorded energy
00667   Float             m_MinJacobian;     // minimum recorded energy
00668   Float             m_Alpha;           // difference parameter
00669 
00670   bool          m_UseLandmarks;               // Use landmark points
00671   bool          m_UseMassMatrix;              // Use Mass matrix in FEM solution
00672   bool          m_UseNormalizedGradient;      // Use normalized gradient magnitude in the metric
00673   bool          m_CreateMeshFromImage;        // Create the mesh based on the fixed image
00674   unsigned int  m_EmployRegridding;           // Use regridding
00675   Sign          m_DescentDirection;           // Metric minimizes or maximizes
00676   Float         m_EnergyReductionFactor;      // Factor we want to reduce the energy - Helps with convergence
00677   ImageSizeType m_FullImageSize;              // image size
00678   ImageSizeType m_ImageOrigin;                // image origin
00679 
00681   ImageSizeType                  m_ImageScaling;
00682   ImageSizeType                  m_CurrentImageScaling;
00683   typename FieldType::RegionType m_FieldRegion;
00684   typename FieldType::SizeType   m_FieldSize;
00685   typename FieldType::Pointer    m_Field;
00686 
00687   // only use TotalField if re-gridding is employed.
00688   typename FieldType::Pointer               m_TotalField;
00689   typename ImageMetricLoadType::Pointer     m_Load; // Defines the load to use
00690 
00691   // define the warper
00692   typename WarperType::Pointer          m_Warper;
00693 
00694   // declare a new image to hold the warped  reference
00695   typename FixedImageType::Pointer      m_WarpedImage;
00696   typename FloatImageType::Pointer      m_FloatImage;
00697   typename FixedImageType::RegionType   m_Wregion;
00698   typename FixedImageType::IndexType    m_Windex;
00699 
00700   // declare images for target and reference
00701   typename MovingImageType::Pointer     m_MovingImage;
00702   typename MovingImageType::Pointer     m_OriginalMovingImage;
00703   typename FixedImageType::Pointer      m_FixedImage;
00704 
00705   // element and metric pointers
00706   typename Element::Pointer             m_Element;
00707   typename MaterialType::Pointer        m_Material;
00708   MetricBaseTypePointer                 m_Metric;
00709   typename FEMObjectType::Pointer       m_FEMObject;
00710 
00711   LandmarkArrayType   m_LandmarkArray;
00712   InterpolatorPointer m_Interpolator;
00713 
00715   double m_MaximumError;
00716 
00718   unsigned int m_MaximumKernelWidth;
00719 
00720   StandardDeviationsType m_StandardDeviations;
00721 
00722 };
00723 
00724 }
00725 }  // end namespace itk::fem
00726 
00727 #ifndef ITK_MANUAL_INSTANTIATION
00728 #include "itkFEMRegistrationFilter.hxx"
00729 #endif
00730 
00731 #endif
00732