ITK  4.6.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 <typename TMovingImage, typename TFixedImage, typename TFemObjectType>
118 class 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 
194 
201  {
202  return m_MovingImage;
203  }
204 
207  {
208  return m_OriginalMovingImage;
209  }
210 
212  void SetFixedImage(FixedImageType* T);
213 
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 
243  {
244  return m_WarpedImage;
245  }
246 
248  void ComputeJacobian( );
249 
252  {
253  return m_FloatImage;
254  }
255 
258  {
259  return m_Field;
260  }
261 
264  {
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  {
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);
357  {
358  SetUseLandmarks(false);
359  }
361 
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 
407  {
409  }
410 
413  {
415  }
416 
421  itkSetMacro(DoLineSearchOnImageEnergy, unsigned int);
422  itkGetMacro(DoLineSearchOnImageEnergy, unsigned int);
424 
425 
430  itkSetMacro(UseNormalizedGradient, bool);
431  itkGetMacro(UseNormalizedGradient, bool);
433  {
435  }
437 
439  {
441  }
442 
446  itkSetMacro(EmployRegridding, unsigned int);
447  itkGetMacro(EmployRegridding, unsigned int);
449 
454  itkSetMacro(LineSearchMaximumIterations, unsigned int);
455  itkGetMacro(LineSearchMaximumIterations, unsigned int);
457 
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 
494  {
495  m_Element = e;
496  }
497 
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);
520  {
522  }
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 = ITK_NULLPTR);
573 
576 
577  // FIXME - Not implemented
581 
582  // FIXME - Not implemented
585 
591 
592  // FIXME - Not implemented
596 
599 
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 
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
void SetElasticity(Float i, unsigned int which=0)
InterpolationGridType::SizeType InterpolationGridSizeType
Superclass::RegionType RegionType
Definition: itkImage.h:140
void FindBracketingTriplet(SolverType *mySolver, Float *a, Float *b, Float *c)
Image< float, itkGetStaticConstMacro(ImageDimension)> FloatImageType
Light weight base class for most itk classes.
FieldPointer ExpandVectorField(ExpandFactorsType *expandFactors, FieldType *f)
SmartPointer< const Self > ConstPointer
itk::Image< Element::ConstPointer, ImageDimension > InterpolationGridType
void GetLandmark(unsigned int i, PointType &source, PointType &target)
Represent the size (bounds) of a n-dimensional image.
Definition: itkSize.h:52
void operator=(const Self &)
void SetMeshPixelsPerElementAtEachResolution(unsigned int i, unsigned int which=0)
itk::Image< VectorType, itkGetStaticConstMacro(ImageDimension)> FieldType
void PrintVectorField(unsigned int modnum=1000)
virtual const RegionType & GetLargestPossibleRegion() const
Definition: itkImageBase.h:254
itk::VectorIndexSelectionCastImageFilter< FieldType, FloatImageType > IndexSelectCasterType
FEM Solver for time dependent problems; uses Crank-Nicolson implicit discretization scheme...
InterpolationGridType::PointType InterpolationGridPointType
FixedImageType::SizeType m_CurrentLevelImageSize
itk::ImageRegionIteratorWithIndex< FloatImageType > FloatImageIterator
void SetNumberOfIntegrationPoints(unsigned int i, unsigned int which=0)
Generate a rectilinar mesh from an image. The result is stored in a FEMObject.
void SetMovingImage(MovingImageType *R)
FEMObjectType * GetInputFEMObject(unsigned int level=0)
static const double e
The base of the natural logarithm or Euler&#39;s number
Definition: itkMath.h:45
void SetWidthOfMetricRegion(unsigned int i, unsigned int which=0)
FixedArray< double, ImageDimension > StandardDeviationsType
FixedImageType::RegionType m_Wregion
Float GetElasticity(unsigned int which=0)
itk::fem::ImageToRectilinearFEMObjectFilter< TMovingImage > ImageToMeshType
vnl_vector< unsigned int > m_MeshPixelsPerElementAtEachResolution
void ChooseMetric(unsigned int whichmetric)
itk::Vector< float, itkGetStaticConstMacro(ImageDimension)> VectorType
itk::VectorExpandImageFilter< FieldType, FieldType > ExpanderType
void SetRho(Float r, unsigned int which=0)
void AddLandmark(PointType source, PointType target)
Expand the size of a vector image by an integer factor in each dimension.
Linear elasticity material class.
void ApplyImageLoads(MovingImageType *i1, FixedImageType *i2)
vnl_vector< unsigned int > m_NumberOfIntegrationPoints
void SetMaxLevel(unsigned int level)
ExpanderType::ExpandFactorsType ExpandFactorsType
General image pair load that uses the itkFiniteDifferenceFunctions.
A templated class holding a n-Dimensional vector.
Definition: itkVector.h:62
void PrintSelf(std::ostream &os, Indent indent) const
FloatImageType * GetMetricImage(FieldType *F)
void SetMaterial(MaterialType::Pointer m)
A multi-dimensional iterator templated over image type that walks pixels within a region and is speci...
InterpolatorType::Pointer InterpolatorPointer
void CreateMesh(unsigned int ElementsPerSide, SolverType *solver)
ImageToImageFilter< TMovingImage, TFixedImage > Superclass
FiniteDifferenceFunctionLoad< MovingImageType, FixedImageType > ImageMetricLoadType
void InterpolateVectorField(SolverType *S)
const SizeType & GetSize() const
Warps an image using an input displacement field.
FEM Image registration filter. The image registration problem is modelled here with the finite elemen...
VectorInterpolateImageFunction< FieldType, CoordRepType > InterpolatorType
vnl_vector< unsigned int > m_Maxiters
Float EvaluateResidual(SolverType *mySolver, Float t)
void WarpImage(const MovingImageType *R)
itk::ImageRegionIteratorWithIndex< FixedImageType > ImageIterator
SolverCrankNicolson< ImageDimension > SolverType
unsigned int GetWidthOfMetricRegion(unsigned int which=0)
virtual void SetStandardDeviations(StandardDeviationsType _arg)
void IterativeSolve(SolverType *S)
Array for FEMP objects.
Definition: itkFEMPArray.h:40
MovingImageType::Pointer m_OriginalMovingImage
itk::WarpImageFilter< MovingImageType, FixedImageType, FieldType > WarperType
void InsertLandmark(unsigned int i, PointType source, PointType target)
VectorLinearInterpolateImageFunction< FieldType, CoordRepType > DefaultInterpolatorType
Float GoldenSection(SolverType *mySolver, Float tol=0.01, unsigned int MaxIters=25)
Base class for filters that take an image as input and produce an image as output.
static const unsigned int ImageDimension
itk::ImageRegionIteratorWithIndex< FieldType > FieldIterator
LinearSystemWrapperItpack LinearSystemSolverType
Control indentation during Print() invocation.
Definition: itkIndent.h:49
void SetFixedImage(FixedImageType *T)
void ApplyLoads(ImageSizeType Isz, double *spacing=ITK_NULLPTR)
Base class for all vector image interpolaters.
Extracts the selected index of the vector that is the input pixel type.
General abstract load base class.
ImageMetricLoadType::Pointer m_Load
void DeleteLandmark(unsigned int i)
LinearSystemWrapper class that uses Itpack numeric library functions to define and solve a sparse lin...
virtual void SetCreateMeshFromImage(bool _arg)
MetricBaseType::Pointer MetricBaseTypePointer
virtual void SetUseNormalizedGradient(bool _arg)
std::vector< typename LoadLandmark::Pointer > LandmarkArrayType
void SetInputFEMObject(FEMObjectType *F, unsigned int level=0)
void SetGamma(Float r, unsigned int which=0)
void SampleVectorFieldAtNodes(SolverType *S)
void SetMaximumIterations(unsigned int i, unsigned int which)
virtual void SetUseLandmarks(bool _arg)
Templated n-dimensional image class.
Definition: itkImage.h:75
PDEDeformableRegistrationFunction< FixedImageType, MovingImageType, FieldType > MetricBaseType
void EnforceDiffeomorphism(float thresh, SolverType *S, bool onlywriteimages)
Linearly interpolate a vector image at specified positions.
vnl_vector< unsigned int > m_MetricWidth