ITK  4.4.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 
470  itkGetModifiableObjectMacro(Metric, MetricBaseType);
471  itkSetObjectMacro(Metric, MetricBaseType);
473 
482  void ChooseMetric( unsigned int whichmetric);
483 
487  unsigned int GetMetricType()
488  {
489  return m_WhichMetric;
490  }
491 
493  void SetElement(Element::Pointer e)
494  {
495  m_Element = e;
496  }
497 
499  void SetMaterial(MaterialType::Pointer m)
500  {
501  m_Material = m;
502  }
503 
505  void PrintVectorField(unsigned int modnum = 1000);
506 
510  void SetMaxLevel(unsigned int level);
511  itkGetMacro(MaxLevel, unsigned int);
513 
518  itkSetMacro(CreateMeshFromImage, bool);
519  void SetCreateMeshFromImageOn()
520  {
521  SetCreateMeshFromImage(true);
522  }
523  void SetCreateMeshFromImageOff()
524  {
525  SetCreateMeshFromImage(false);
526  }
527  itkGetMacro(CreateMeshFromImage, bool);
529 
531  itkSetObjectMacro( Interpolator, InterpolatorType );
532 
534  itkGetModifiableObjectMacro( Interpolator, InterpolatorType );
535 
539  itkSetMacro(StandardDeviations, StandardDeviationsType);
540  virtual void SetStandardDeviations(double value);
542 
545  itkGetConstReferenceMacro(StandardDeviations, StandardDeviationsType);
546 
549  itkSetMacro(MaximumKernelWidth, unsigned int);
550  itkGetConstMacro(MaximumKernelWidth, unsigned int);
552 
555  itkSetMacro(MaximumError, double);
556  itkGetConstMacro(MaximumError, double);
558 
559 
560 // HELPER FUNCTIONS
561 
562 protected:
565 
566  void PrintSelf(std::ostream & os, Indent indent) const;
567 
569  void CreateMesh(unsigned int ElementsPerSide, SolverType *solver);
570 
572  void ApplyLoads(ImageSizeType Isz, double* spacing = NULL);
573 
575  void ApplyImageLoads(MovingImageType* i1, FixedImageType* i2);
576 
577  // FIXME - Not implemented
580  void CreateLinearSystemSolver();
581 
582  // FIXME - Not implemented
584  Float EvaluateEnergy();
585 
590  void InterpolateVectorField(SolverType *S);
591 
592  // FIXME - Not implemented
595  FloatImageType * GetMetricImage(FieldType* F);
596 
598  FieldPointer ExpandVectorField(ExpandFactorsType* expandFactors, FieldType* f);
599 
601  void SampleVectorFieldAtNodes(SolverType *S);
602 
604  Float EvaluateResidual(SolverType *mySolver, Float t);
605 
606  // FIXME - Replace with BSD Code
607  /* Finds a triplet that brackets the energy minimum. From Numerical Recipes.*/
608  // void FindBracketingTriplet(SolverType& mySolver,Float* a, Float* b, Float* c);
609  void FindBracketingTriplet(SolverType *mySolver, Float* a, Float* b, Float* c);
610 
615  Float GoldenSection(SolverType *mySolver, Float tol = 0.01, unsigned int MaxIters = 25);
616 
618  itkSetObjectMacro( Load, ImageMetricLoadType );
619  itkGetModifiableObjectMacro(Load, ImageMetricLoadType );
621 
623  void SmoothDisplacementField();
624 
625 private:
626 
627  void InitializeField();
628 
629  FEMRegistrationFilter(const Self &); // purposely not implemented
630  void operator=(const Self &); // purposely not implemented
631 
634 
635  /* Parameters used to define Multi-resolution Registration */
636  vnl_vector<unsigned int> m_NumberOfIntegrationPoints; // resolution of integration
637  vnl_vector<unsigned int> m_MetricWidth; // number of iterations at each level
638  vnl_vector<unsigned int> m_Maxiters; // max iterations
639  unsigned int m_TotalIterations; // total number of iterations that were run
640  unsigned int m_MaxLevel; // Number of Resolution Levels for registration.
641  unsigned int m_FileCount; // keeps track of number of files written
642  unsigned int m_CurrentLevel; // current resolution level
643 
644  typename FixedImageType::SizeType m_CurrentLevelImageSize;
645 
646  unsigned int m_WhichMetric; // Metric used for image registration
647  // 0 = Mean Squares
648  // 1 = Normalized Correlation
649  // 2 = Mutual Information
650  // 3 = Demons Metric
651 
654  vnl_vector<unsigned int> m_MeshPixelsPerElementAtEachResolution;
655 
656  Float m_TimeStep; // time step
657  vnl_vector<Float> m_E; // elasticity
658  vnl_vector<Float> m_Rho; // mass matrix weight
659  vnl_vector<Float> m_Gamma; // image similarity weight
660  Float m_Energy; // current value of energy
661  Float m_MinE; // minimum recorded energy
662  Float m_MinJacobian; // minimum recorded energy
663  Float m_Alpha; // difference parameter
664 
665  bool m_UseLandmarks; // Use landmark points
666  bool m_UseMassMatrix; // Use Mass matrix in FEM solution
667  bool m_UseNormalizedGradient; // Use normalized gradient magnitude in the metric
668  bool m_CreateMeshFromImage; // Create the mesh based on the fixed image
669  unsigned int m_EmployRegridding; // Use regridding
670  Sign m_DescentDirection; // Metric minimizes or maximizes
671  Float m_EnergyReductionFactor; // Factor we want to reduce the energy - Helps with convergence
673  ImageSizeType m_ImageOrigin; // image origin
674 
681 
682  // only use TotalField if re-gridding is employed.
684  typename ImageMetricLoadType::Pointer m_Load; // Defines the load to use
685 
686  // define the warper
688 
689  // declare a new image to hold the warped reference
690  typename FixedImageType::Pointer m_WarpedImage;
692  typename FixedImageType::RegionType m_Wregion;
693  typename FixedImageType::IndexType m_Windex;
694 
695  // declare images for target and reference
696  typename MovingImageType::Pointer m_MovingImage;
697  typename MovingImageType::Pointer m_OriginalMovingImage;
698  typename FixedImageType::Pointer m_FixedImage;
699 
700  // element and metric pointers
704  typename FEMObjectType::Pointer m_FEMObject;
705 
708 
711 
713  unsigned int m_MaximumKernelWidth;
714 
716 
717 };
718 
719 }
720 } // end namespace itk::fem
721 
722 #ifndef ITK_MANUAL_INSTANTIATION
723 #include "itkFEMRegistrationFilter.hxx"
724 #endif
725 
726 #endif
727