ITK  4.3.0
Insight Segmentation and Registration Toolkit
itkFEMRegistrationFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #ifndef __itkFEMRegistrationFilter_h
20 #define __itkFEMRegistrationFilter_h
21 
24 #include "itkFEMObject.h"
27 #include "itkFEMImageMetricLoad.h"
29 
30 #include "itkImage.h"
31 #include "itkVector.h"
36 #include "itkWarpImageFilter.h"
37 #include "itkImageToImageMetric.h"
40 #include "itkFixedArray.h"
41 
43 #include "itkFEMLoadLandmark.h"
44 
45 #include "vnl/vnl_vector.h"
46 #include "vnl/vnl_math.h"
47 #include "vnl/vnl_vector_fixed.h"
48 
49 #include <iostream>
50 #include <string>
51 
52 namespace itk
53 {
54 namespace fem
55 {
56 
117 template <class TMovingImage, class TFixedImage, class TFemObjectType>
118 class ITK_EXPORT FEMRegistrationFilter : public ImageToImageFilter<TMovingImage, TFixedImage>
119 {
120 public:
125 
127  itkNewMacro(Self);
128 
131 
132  typedef TMovingImage MovingImageType;
133  typedef TFixedImage FixedImageType;
134  typedef TFemObjectType FEMObjectType;
135  typedef typename FixedImageType::PixelType PixelType;
136  typedef typename FixedImageType::SizeType ImageSizeType;
137  typedef typename FixedImageType::PointType PointType;
138 
140  itkStaticConstMacro(ImageDimension, unsigned int,
141  FixedImageType::ImageDimension);
142 
146 
147  enum Sign { positive = 1, negative = -1 };
148  typedef double Float;
150 
151  typedef std::vector<typename LoadLandmark::Pointer> LandmarkArrayType;
155 
160 
162 
164  typedef double CoordRepType;
167 
169 
173 
175 
178 
180 
186 
189 
190  /*---------------------- Main functions ----------------------*/
191 
193  void SetMovingImage(MovingImageType* R);
194 
200  MovingImageType * GetMovingImage()
201  {
202  return m_MovingImage;
203  }
204 
206  MovingImageType * GetOriginalMovingImage()
207  {
208  return m_OriginalMovingImage;
209  }
210 
212  void SetFixedImage(FixedImageType* T);
213 
214  FixedImageType * GetFixedImage()
215  {
216  return m_FixedImage;
217  }
218 
224  void SetInputFEMObject(FEMObjectType* F, unsigned int level = 0);
225 
226  FEMObjectType * GetInputFEMObject(unsigned int level = 0);
227 
229  void RunRegistration(void);
230 
232  void IterativeSolve(SolverType *S);
233 
235  void MultiResSolve();
236 
238  void WarpImage(const MovingImageType * R);
239 
242  FixedImageType * GetWarpedImage()
243  {
244  return m_WarpedImage;
245  }
246 
248  void ComputeJacobian( );
249 
251  FloatImageType * GetJacobianImage()
252  {
253  return m_FloatImage;
254  }
255 
257  FieldType * GetDisplacementField()
258  {
259  return m_Field;
260  }
261 
263  void SetDisplacementField(FieldType* F)
264  {
265  m_FieldSize = F->GetLargestPossibleRegion().GetSize();
266  m_Field = F;
267  }
269 
271  void AddLandmark(PointType source, PointType target);
272 
273  void InsertLandmark(unsigned int i, PointType source, PointType target);
274 
275  void DeleteLandmark(unsigned int i);
276 
277  void ClearLandmarks();
278 
279  void GetLandmark(unsigned int i, PointType& source, PointType& target);
280 
288  void EnforceDiffeomorphism(float thresh, SolverType *S, bool onlywriteimages);
289 
296  void SetMeshPixelsPerElementAtEachResolution(unsigned int i, unsigned int which = 0)
297  {
298  m_MeshPixelsPerElementAtEachResolution[which] = i;
299  }
301 
305  void SetNumberOfIntegrationPoints(unsigned int i, unsigned int which = 0)
306  {
307  m_NumberOfIntegrationPoints[which] = i;
308  }
309 
315  void SetWidthOfMetricRegion(unsigned int i, unsigned int which = 0)
316  {
317  m_MetricWidth[which] = i;
318  }
319 
320  unsigned int GetWidthOfMetricRegion(unsigned int which = 0)
321  {
322  return m_MetricWidth[which];
323  }
324 
329  void SetMaximumIterations(unsigned int i, unsigned int which)
330  {
331  m_Maxiters[which] = i;
332  }
333 
339  itkSetMacro(TimeStep, Float);
340  itkGetMacro(TimeStep, Float);
342 
347  itkSetMacro(Alpha, Float);
348  itkGetMacro(Alpha, Float);
350 
354  itkSetMacro(UseLandmarks, bool);
355  itkGetMacro(UseLandmarks, bool);
356  void SetUseLandmarksOff()
357  {
358  SetUseLandmarks(false);
359  }
361 
362  void SetUseLandmarksOn()
363  {
364  SetUseLandmarks(true);
365  }
366 
370  itkSetMacro(UseMassMatrix, bool);
371  itkGetMacro(UseMassMatrix, bool);
373 
377  itkSetMacro(EnergyReductionFactor, Float);
378  itkGetMacro(EnergyReductionFactor, Float);
380 
382  void SetElasticity(Float i, unsigned int which = 0)
383  {
384  m_E[which] = i;
385  }
386 
388  Float GetElasticity(unsigned int which = 0)
389  {
390  return m_E[which];
391  }
392 
394  void SetRho(Float r, unsigned int which = 0)
395  {
396  m_Rho[which] = r;
397  }
398 
400  void SetGamma(Float r, unsigned int which = 0)
401  {
402  m_Gamma[which] = r;
403  }
404 
406  void SetDescentDirectionMinimize()
407  {
408  m_DescentDirection = positive;
409  }
410 
412  void SetDescentDirectionMaximize()
413  {
414  m_DescentDirection = negative;
415  }
416 
421  itkSetMacro(DoLineSearchOnImageEnergy, unsigned int);
422  itkGetMacro(DoLineSearchOnImageEnergy, unsigned int);
424 
425 
430  itkSetMacro(UseNormalizedGradient, bool);
431  itkGetMacro(UseNormalizedGradient, bool);
432  void SetUseNormalizedGradientOff()
433  {
434  SetUseNormalizedGradient(false);
435  }
437 
438  void SetUseNormalizedGradientOn()
439  {
440  SetUseNormalizedGradient(true);
441  }
442 
446  itkSetMacro(EmployRegridding, unsigned int);
447  itkGetMacro(EmployRegridding, unsigned int);
449 
454  itkSetMacro(LineSearchMaximumIterations, unsigned int);
455  itkGetMacro(LineSearchMaximumIterations, unsigned int);
457 
461  ImageSizeType GetImageSize()
462  {
463  return m_FullImageSize;
464  }
465 
471  {
472  return m_Metric;
473  }
474 
475  void SetMetric(MetricBaseTypePointer MP)
476  {
477  m_Metric = MP;
478  }
479 
488  void ChooseMetric( unsigned int whichmetric);
489 
493  unsigned int GetMetricType()
494  {
495  return m_WhichMetric;
496  }
497 
499  void SetElement(Element::Pointer e)
500  {
501  m_Element = e;
502  }
503 
505  void SetMaterial(MaterialType::Pointer m)
506  {
507  m_Material = m;
508  }
509 
511  void PrintVectorField(unsigned int modnum = 1000);
512 
516  void SetMaxLevel(unsigned int level);
517  itkGetMacro(MaxLevel, unsigned int);
519 
524  itkSetMacro(CreateMeshFromImage, bool);
525  void SetCreateMeshFromImageOn()
526  {
527  SetCreateMeshFromImage(true);
528  }
529  void SetCreateMeshFromImageOff()
530  {
531  SetCreateMeshFromImage(false);
532  }
533  itkGetMacro(CreateMeshFromImage, bool);
535 
537  itkSetObjectMacro( Interpolator, InterpolatorType );
538 
540  itkGetObjectMacro( Interpolator, InterpolatorType );
541 
545  itkSetMacro(StandardDeviations, StandardDeviationsType);
546  virtual void SetStandardDeviations(double value);
548 
551  itkGetConstReferenceMacro(StandardDeviations, StandardDeviationsType);
552 
555  itkSetMacro(MaximumKernelWidth, unsigned int);
556  itkGetConstMacro(MaximumKernelWidth, unsigned int);
558 
561  itkSetMacro(MaximumError, double);
562  itkGetConstMacro(MaximumError, double);
564 
565 
566 // HELPER FUNCTIONS
567 
568 protected:
571 
572  void PrintSelf(std::ostream & os, Indent indent) const;
573 
575  void CreateMesh(unsigned int ElementsPerSide, SolverType *solver);
576 
578  void ApplyLoads(ImageSizeType Isz, double* spacing = NULL);
579 
581  void ApplyImageLoads(MovingImageType* i1, FixedImageType* i2);
582 
583  // FIXME - Not implemented
586  void CreateLinearSystemSolver();
587 
588  // FIXME - Not implemented
590  Float EvaluateEnergy();
591 
596  void InterpolateVectorField(SolverType *S);
597 
598  // FIXME - Not implemented
601  FloatImageType * GetMetricImage(FieldType* F);
602 
604  FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);
605 
607  void SampleVectorFieldAtNodes(SolverType *S);
608 
610  Float EvaluateResidual(SolverType *mySolver, Float t);
611 
612  // FIXME - Replace with BSD Code
613  /* Finds a triplet that brackets the energy minimum. From Numerical Recipes.*/
614  // void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
615  void FindBracketingTriplet(SolverType *mySolver, Float* a, Float* b, Float* c);
616 
621  Float GoldenSection(SolverType *mySolver, Float tol = 0.01, unsigned int MaxIters = 25);
622 
624  itkGetConstObjectMacro( Load, ImageMetricLoadType );
625  itkSetObjectMacro( Load, ImageMetricLoadType );
627 
629  void SmoothDisplacementField();
630 
631 private:
632 
633  void InitializeField();
634 
635  FEMRegistrationFilter(const Self &); // purposely not implemented
636  void operator=(const Self &); // purposely not implemented
637 
640 
641  /* Parameters used to define Multi-resolution Registration */
643  vnl_vector<unsigned int> m_MetricWidth; // number of iterations at each level
645  unsigned int m_TotalIterations; // total number of iterations that were run
646  unsigned int m_MaxLevel; // Number of Resolution Levels for registration.
647  unsigned int m_FileCount; // keeps track of number of files written
648  unsigned int m_CurrentLevel; // current resolution level
649 
650  typename FixedImageType::SizeType m_CurrentLevelImageSize;
651 
652  unsigned int m_WhichMetric; // Metric used for image registration
653  // 0 = Mean Squares
654  // 1 = Normalized Correlation
655  // 2 = Mutual Information
656  // 3 = Demons Metric
657 
661 
662  Float m_TimeStep; // time step
663  vnl_vector<Float> m_E; // elasticity
664  vnl_vector<Float> m_Rho; // mass matrix weight
665  vnl_vector<Float> m_Gamma; // image similarity weight
666  Float m_Energy; // current value of energy
667  Float m_MinE; // minimum recorded energy
668  Float m_MinJacobian; // minimum recorded energy
669  Float m_Alpha; // difference parameter
670 
671  bool m_UseLandmarks; // Use landmark points
672  bool m_UseMassMatrix; // Use Mass matrix in FEM solution
673  bool m_UseNormalizedGradient; // Use normalized gradient magnitude in the metric
674  bool m_CreateMeshFromImage; // Create the mesh based on the fixed image
675  unsigned int m_EmployRegridding; // Use regridding
676  Sign m_DescentDirection; // Metric minimizes or maximizes
677  Float m_EnergyReductionFactor; // Factor we want to reduce the energy - Helps with convergence
679  ImageSizeType m_ImageOrigin; // image origin
680 
687 
688  // only use TotalField if re-gridding is employed.
690  typename ImageMetricLoadType::Pointer m_Load; // Defines the load to use
691 
692  // define the warper
694 
695  // declare a new image to hold the warped reference
696  typename FixedImageType::Pointer m_WarpedImage;
698  typename FixedImageType::RegionType m_Wregion;
699  typename FixedImageType::IndexType m_Windex;
700 
701  // declare images for target and reference
702  typename MovingImageType::Pointer m_MovingImage;
703  typename MovingImageType::Pointer m_OriginalMovingImage;
704  typename FixedImageType::Pointer m_FixedImage;
705 
706  // element and metric pointers
710  typename FEMObjectType::Pointer m_FEMObject;
711 
714 
717 
719  unsigned int m_MaximumKernelWidth;
720 
722 
723 };
724 
725 }
726 } // end namespace itk::fem
727 
728 #ifndef ITK_MANUAL_INSTANTIATION
729 #include "itkFEMRegistrationFilter.hxx"
730 #endif
731 
732 #endif
733