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
00028 #include "itkImage.h"
00029 #include "itkVector.h"
00030 #include "itkImageRegionIteratorWithIndex.h"
00031 #include "itkVectorCastImageFilter.h"
00032 #include "itkVectorIndexSelectionCastImageFilter.h"
00033 #include "itkWarpImageFilter.h"
00034 #include "itkImageToImageMetric.h"
00035 #include "itkTranslationTransform.h"
00036
00037
00038 #include "vnl/vnl_vector.h"
00039 #include "vnl/vnl_math.h"
00040 #include "vnl/vnl_vector_fixed.h"
00041
00042
00043 #include <iostream>
00044 #include <string>
00045
00046
00047
00048 namespace itk {
00049 namespace fem {
00050
00097 template<class TReference,class TTarget>
00098 class ITK_EXPORT FEMRegistrationFilter : public ImageToImageFilter<TReference, TTarget>
00099 {
00100 public:
00101 typedef FEMRegistrationFilter Self;
00102 typedef ImageToImageFilter<TReference, TTarget> Superclass;
00103 typedef SmartPointer<Self> Pointer;
00104 typedef SmartPointer<const Self> ConstPointer;
00105
00107 itkNewMacro(Self);
00108
00110 itkTypeMacro(FEMRegistrationFilter, ImageToImageFilter );
00111
00112 typedef TReference ImageType;
00113 typedef TTarget TargetImageType;
00114 typedef typename ImageType::PixelType ImageDataType;
00115 typedef typename ImageType::PixelType PixelType;
00116 typedef typename ImageType::SizeType ImageSizeType;
00117
00119 itkStaticConstMacro(ImageDimension, unsigned int,
00120 ImageType::ImageDimension);
00121
00122 typedef Image< float, itkGetStaticConstMacro(ImageDimension) > FloatImageType;
00123 typedef LinearSystemWrapperItpack LinearSystemSolverType;
00124 typedef SolverCrankNicolson SolverType;
00125 enum Sign { positive = 1, negative = -1 };
00126 typedef double Float;
00127
00128
00129 typedef MaterialLinearElasticity MaterialType;
00130 typedef ImageToImageMetric<ImageType,TargetImageType> MetricBaseType;
00131 typedef typename MetricBaseType::Pointer MetricBaseTypePointer;
00132 typedef itk::Vector<float,itkGetStaticConstMacro(ImageDimension)> VectorType;
00133 typedef itk::Image<VectorType,itkGetStaticConstMacro(ImageDimension)> FieldType;
00134 typedef itk::WarpImageFilter<ImageType,ImageType, FieldType> WarperType;
00135 typedef itk::ImageRegionIteratorWithIndex<ImageType> ImageIterator;
00136 typedef itk::ImageRegionIteratorWithIndex<FieldType> FieldIterator;
00137 typedef itk::VectorIndexSelectionCastImageFilter<FieldType,FloatImageType> IndexSelectCasterType;
00138
00140 typedef ImageMetricLoad<ImageType,ImageType> ImageMetricLoadType;
00141
00146 class FEMOF : public FEMObjectFactory<FEMLightObject>{};
00147
00148
00149
00151 bool ReadConfigFile(const char*);
00152
00154 void RunRegistration();
00155
00159 void WriteWarpedImage(const char* fn);
00160
00161
00165 void CreateMesh(double ElementsPerSide, Solver& S);
00166
00168 void ApplyLoads(SolverType& S,ImageSizeType Isz);
00169
00172 void CreateLinearSystemSolver();
00173
00175 void IterativeSolve(SolverType& S);
00176
00178 void MultiResSolve();
00179
00181 Float EvaluateEnergy();
00182
00187 void GetVectorField(SolverType& S);
00188
00190 void SampleVectorFieldAtNodes(SolverType& S);
00191
00193 void WarpImage(ImageType* R);
00194
00196 int WriteDisplacementField(unsigned int index);
00197
00199 void SetReferenceFile(const char* r) {m_ReferenceFileName=r;}
00200
00201 std::string GetReferenceFile() {return m_ReferenceFileName;}
00202
00203 void SetTargetFile(const char* t) {m_TargetFileName=t;}
00204
00205 std::string GetTargetFile() {return m_TargetFileName;}
00206
00207
00211 void SetReferenceImage(ImageType* R);
00212
00214 void SetTargetImage(TargetImageType* T);
00215
00216 ImageType* GetReferenceImage(){return m_RefImg;}
00217
00218 TargetImageType* GetTargetImage(){return m_TarImg;}
00219
00220
00223 ImageType* GetWarpedImage(){return m_WarpedImage;}
00224
00226 FieldType* GetDeformationField(){return m_Field;}
00227
00230 void SetLandmarkFile(const char* l) {m_LandmarkFileName=l; }
00231
00233 void UseLandmarks(bool b) {m_UseLandmarks=b;}
00234
00239 void SetResultsFile(const char* r) {m_ResultsFileName=r;}
00240
00241 void SetResultsFileName (const char* f){m_ResultsFileName=f;}
00242
00243 std::string GetResultsFileName () {return m_ResultsFileName;}
00244
00246 void SetDisplacementsFile(const char* r) {m_DisplacementsFileName=r;}
00247
00254 void SetMeshElementsPerDimensionAtEachResolution(unsigned int i,unsigned int which=0){ m_MeshElementsPerDimensionAtEachResolution[which]=i;}
00255
00259 void SetNumberOfIntegrationPoints(unsigned int i,unsigned int which=0){ m_NumberOfIntegrationPoints[which]=i;}
00260
00266 void SetWidthOfMetricRegion(unsigned int i,unsigned int which=0) { m_MetricWidth[which]=i;}
00267
00272 void SetMaximumIterations(unsigned int i,unsigned int which) { m_Maxiters[which]=i;}
00273
00276 void SetTimeStep(Float i) { m_dT=i;}
00277
00279 void SetAlpha(Float a) { m_Alpha=a;}
00280
00283 void SetEnergyReductionFactor(Float i) { m_EnergyReductionFactor=i;}
00284
00286 void SetElasticity(Float i,unsigned int which=0) { m_E[which]=i;}
00287
00289 Float GetElasticity(unsigned int which=0) { return m_E[which];}
00290
00292 void SetRho(Float r,unsigned int which=0) { m_Rho[which]=r;}
00293
00295 void SetGamma(Float r,unsigned int which=0) { m_Gamma[which]=r;}
00296
00298 void SetDescentDirectionMinimize() { m_DescentDirection=positive;}
00299
00301 void SetDescentDirectionMaximize() { m_DescentDirection=negative;}
00302
00304 void DoLineSearch(unsigned int b) { m_DoLineSearchOnImageEnergy=b; }
00305
00306
00308 void DoMultiRes(bool b) { m_DoMultiRes=b; }
00309
00311 void SetLineSearchFrequency(unsigned int f) { m_LineSearchFrequency=f; }
00312
00314 void SetLineSearchMaximumIterations(unsigned int f) { m_LineSearchMaximumIterations=f; }
00315
00317 void SetWriteDisplacements(bool b) {m_WriteDisplacementField=b;}
00318
00320 bool GetWriteDisplacements() {return m_WriteDisplacementField;}
00321
00324 void SetConfigFileName (const char* f){m_ConfigFileName=f;}
00325
00326 std::string GetConfigFileName () {return m_ConfigFileName; }
00327
00328 ImageSizeType GetImageSize(){ return m_ImageSize; }
00329
00331 void SetMetric(MetricBaseTypePointer MP) { m_Metric=MP; }
00332
00335 void ChooseMetric( float whichmetric);
00336
00338 void SetElement(Element::Pointer e) {m_Element=e;}
00339
00341 void SetMaterial(MaterialType::Pointer m) {m_Material=m;}
00342
00343 void PrintVectorField();
00344
00345 Float EvaluateResidual(SolverType& mySolver,Float t);
00346
00347
00348 void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
00349
00353 Float GoldenSection(SolverType& mySolver,Float tol=0.01,unsigned int MaxIters=25);
00354
00356 FEMRegistrationFilter( );
00357 ~FEMRegistrationFilter();
00358
00359 protected :
00360
00361
00362 void PrintSelf(std::ostream& os, Indent indent) const
00363 {
00364 std::cout << "FEM registration filter "<< std::endl;
00365 Superclass::PrintSelf( os, indent );
00366 }
00367
00368 private :
00369
00370 FEMRegistrationFilter(const Self&);
00371 void operator=(const Self&);
00372
00373 std::string m_ConfigFileName;
00374 std::string m_ResultsFileName;
00375 std::string m_ReferenceFileName;
00376 std::string m_TargetFileName;
00377 std::string m_LandmarkFileName;
00378 std::string m_DisplacementsFileName;
00379 std::string m_MeshFileName;
00380
00381 unsigned int m_DoLineSearchOnImageEnergy;
00382 unsigned int m_LineSearchFrequency;
00383 unsigned int m_LineSearchMaximumIterations;
00384
00385 vnl_vector<unsigned int> m_NumberOfIntegrationPoints;
00386 vnl_vector<unsigned int> m_MetricWidth;
00387 vnl_vector<unsigned int> m_Maxiters;
00388 unsigned int m_TotalIterations;
00389 unsigned int m_NumLevels;
00390 unsigned int m_MaxLevel;
00391 unsigned int m_MeshLevels;
00392 unsigned int m_MeshStep;
00393 unsigned int m_FileCount;
00394 unsigned int m_CurrentLevel;
00395
00398 vnl_vector<unsigned int> m_MeshElementsPerDimensionAtEachResolution;
00399
00400 Float m_dT;
00401 vnl_vector<Float> m_E;
00402 vnl_vector<Float> m_Rho;
00403 vnl_vector<Float> m_Gamma;
00404 Float m_Energy;
00405 Float m_MinE;
00406 Float m_Alpha;
00408 Float m_EnergyReductionFactor;
00409 Float m_Temp;
00410
00411 bool m_WriteDisplacementField;
00412 bool m_DoMultiRes;
00413 bool m_UseLandmarks;
00414 bool m_ReadMeshFile;
00415 Sign m_DescentDirection;
00416
00417 ImageSizeType m_ImageSize;
00418 ImageSizeType m_ImageOrigin;
00420 ImageSizeType m_ImageScaling;
00421 typename ImageType::RegionType m_FieldRegion;
00422 typename ImageType::SizeType m_FieldSize;
00423 typename FieldType::Pointer m_Field;
00424
00425 ImageMetricLoad<ImageType,ImageType>* m_Load;
00426
00427
00428 typename WarperType::Pointer m_Warper;
00429
00430
00431 typename ImageType::Pointer m_WarpedImage;
00432 typename ImageType::RegionType m_Wregion;
00433 typename ImageType::IndexType m_Windex;
00434
00435
00436 typename ImageType::Pointer m_RefImg;
00437 typename ImageType::Pointer m_TarImg;
00438 typename ImageType::RegionType m_Rregion;
00439 typename ImageType::RegionType m_Tregion;
00440 typename ImageType::IndexType m_Rindex;
00441 typename ImageType::IndexType m_Tindex;
00442
00443
00444 typename Element::Pointer m_Element;
00445 typename MaterialType::Pointer m_Material;
00446 MetricBaseTypePointer m_Metric;
00447
00448
00449 };
00450
00451 }}
00452
00453 #ifndef ITK_MANUAL_INSTANTIATION
00454 #include "itkFEMRegistrationFilter.txx"
00455 #endif
00456
00457 #endif