00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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
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
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
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
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
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&);
00568 void operator=(const Self&);
00569
00570 std::string m_ConfigFileName;
00571 std::string m_ResultsFileName;
00572 std::string m_MovingFileName;
00573 std::string m_FixedFileName;
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;
00582 vnl_vector<unsigned int> m_MetricWidth;
00583 vnl_vector<unsigned int> m_Maxiters;
00584 unsigned int m_TotalIterations;
00585 unsigned int m_NumLevels;
00586 unsigned int m_MaxLevel;
00587 unsigned int m_MeshLevels;
00588 unsigned int m_MeshStep;
00589 unsigned int m_FileCount;
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;
00601 vnl_vector<Float> m_E;
00602 vnl_vector<Float> m_Rho;
00603 vnl_vector<Float> m_Gamma;
00604 Float m_Energy;
00605 Float m_MinE;
00606 Float m_MinJacobian;
00607 Float m_Alpha;
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;
00621 ImageSizeType m_ImageOrigin;
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
00629 typename FieldType::Pointer m_TotalField;
00630 ImageMetricLoadType* m_Load;
00631
00632
00633 typename WarperType::Pointer m_Warper;
00634
00635
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
00642 typename MovingImageType::Pointer m_MovingImage;
00643 typename MovingImageType::Pointer m_OriginalMovingImage;
00644 typename FixedImageType::Pointer m_FixedImage;
00645
00646
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 }}
00659
00660 #ifndef ITK_MANUAL_INSTANTIATION
00661 #include "itkFEMRegistrationFilter.txx"
00662 #endif
00663
00664 #endif
00665