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 #ifndef itkMeanSquareRegistrationFunction_h
19 #define itkMeanSquareRegistrationFunction_h
20 
21 #include "itkPDEDeformableRegistrationFunction.h"
22 #include "itkPoint.h"
23 #include "itkLinearInterpolateImageFunction.h"
24 #include "itkCentralDifferenceImageFunction.h"
25 
26 namespace itk
27 {
28 /**
29  * \class MeanSquareRegistrationFunction
30  *
31  * This class encapsulate the PDE which drives the demons registration
32  * algorithm. It is used by MeanSquareRegistrationFilter to compute the
33  * output displacement field which will map a moving image onto a
34  * a fixed image.
35  *
36  * Non-integer moving image values are obtained by using
37  * interpolation. The default interpolator is of type
38  * LinearInterpolateImageFunction. The user may set other
39  * interpolators via method SetMovingImageInterpolator. Note that the input
40  * interpolator must derive from baseclass InterpolateImageFunction.
41  *
42  * This class is templated over the fixed image type, moving image type,
43  * and the displacement field type.
44  *
45  * \warning This filter assumes that the fixed image type, moving image type
46  * and displacement field type all have the same number of dimensions.
47  *
48  * \sa MeanSquareRegistrationFilter
49  * \ingroup FiniteDifferenceFunctions
50  * \ingroup ITKRegistrationCommon
51  */
52 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
53 class ITK_TEMPLATE_EXPORT MeanSquareRegistrationFunction:
54   public PDEDeformableRegistrationFunction< TFixedImage,
55                                             TMovingImage, TDisplacementField >
56 {
57 public:
58   ITK_DISALLOW_COPY_AND_ASSIGN(MeanSquareRegistrationFunction);
59 
60   /** Standard class type aliases. */
61   using Self = MeanSquareRegistrationFunction;
62   using Superclass = PDEDeformableRegistrationFunction< TFixedImage,
63                                              TMovingImage, TDisplacementField >;
64   using Pointer = SmartPointer< Self >;
65   using ConstPointer = SmartPointer< const Self >;
66 
67   /** Method for creation through the object factory. */
68   itkNewMacro(Self);
69 
70   /** Run-time type information (and related methods). */
71   itkTypeMacro(MeanSquareRegistrationFunction,
72                PDEDeformableRegistrationFunction);
73 
74   /** MovingImage image type. */
75   using MovingImageType = typename Superclass::MovingImageType;
76   using MovingImagePointer = typename Superclass::MovingImagePointer;
77 
78   /** FixedImage image type. */
79   using FixedImageType = typename Superclass::FixedImageType;
80   using FixedImagePointer = typename Superclass::FixedImagePointer;
81   using IndexType = typename FixedImageType::IndexType;
82   using SizeType = typename FixedImageType::SizeType;
83   using SpacingType = typename FixedImageType::SpacingType;
84 
85   /** Displacement field type. */
86   using DisplacementFieldType = typename Superclass::DisplacementFieldType;
87   using DisplacementFieldPixelType = typename DisplacementFieldType::PixelType;
88   using DisplacementFieldTypePointer = typename Superclass::DisplacementFieldTypePointer;
89 
90   /** Inherit some enums from the superclass. */
91   static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
92 
93   /** Inherit some enums from the superclass. */
94   using PixelType = typename Superclass::PixelType;
95   using RadiusType = typename Superclass::RadiusType;
96   using NeighborhoodType = typename Superclass::NeighborhoodType;
97   using FloatOffsetType = typename Superclass::FloatOffsetType;
98   using TimeStepType = typename Superclass::TimeStepType;
99 
100   /** Interpolator type. */
101   using CoordRepType = double;
102   using InterpolatorType = InterpolateImageFunction< MovingImageType, CoordRepType >;
103   using InterpolatorPointer = typename InterpolatorType::Pointer;
104   using PointType = typename InterpolatorType::PointType;
105   using DefaultInterpolatorType = LinearInterpolateImageFunction< MovingImageType, CoordRepType >;
106 
107   /** Covariant vector type. */
108   using CovariantVectorType = CovariantVector< double, Self::ImageDimension >;
109 
110   /** Gradient calculator type. */
111   using GradientCalculatorType = CentralDifferenceImageFunction< FixedImageType >;
112   using GradientCalculatorPointer = typename GradientCalculatorType::Pointer;
113 
114   /** Set the moving image interpolator. */
SetMovingImageInterpolator(InterpolatorType * ptr)115   void SetMovingImageInterpolator(InterpolatorType *ptr)
116   { m_MovingImageInterpolator = ptr; }
117 
118   /** Get the moving image interpolator. */
GetMovingImageInterpolator()119   InterpolatorType * GetMovingImageInterpolator()
120   { return m_MovingImageInterpolator; }
121 
122   /** This class uses a constant timestep of 1. */
ComputeGlobalTimeStep(void * itkNotUsed (GlobalData))123   TimeStepType ComputeGlobalTimeStep( void *itkNotUsed(GlobalData) ) const override
124   { return m_TimeStep; }
125 
126   /** Return a pointer to a global data structure that is passed to
127    * this object from the solver at each calculation.  */
GetGlobalDataPointer()128   void * GetGlobalDataPointer() const override
129   {
130     auto * global = new GlobalDataStruct();
131 
132     return global;
133   }
134 
135   /** Release memory for global data structure. */
ReleaseGlobalDataPointer(void * GlobalData)136   void ReleaseGlobalDataPointer(void *GlobalData) const override
137   { delete (GlobalDataStruct *)GlobalData;  }
138 
139   /** Set the object's state before each iteration. */
140   void InitializeIteration() override;
141 
142   /** This method is called by a finite difference solver image filter at
143    * each pixel that does not lie on a data set boundary */
144   PixelType  ComputeUpdate( const NeighborhoodType & neighborhood,
145                                     void *globalData,
146                                     const FloatOffsetType & offset = FloatOffsetType(0.0) ) override;
147 
148 protected:
149   MeanSquareRegistrationFunction();
150   ~MeanSquareRegistrationFunction() override = default;
151   void PrintSelf(std::ostream & os, Indent indent) const override;
152 
153   /** FixedImage image neighborhood iterator type. */
154   using FixedImageNeighborhoodIteratorType = ConstNeighborhoodIterator< FixedImageType >;
155 
156   /** A global data type for this class of equation. Used to store
157    * iterators for the fixed image. */
158   struct GlobalDataStruct {
159     FixedImageNeighborhoodIteratorType m_FixedImageIterator;
160   };
161 
162 private:
163   /** Cache fixed image information. */
164   SpacingType m_FixedImageSpacing;
165 
166   /** Function to compute derivatives of the fixed image. */
167   GradientCalculatorPointer m_FixedImageGradientCalculator;
168 
169   /** Function to interpolate the moving image. */
170   InterpolatorPointer m_MovingImageInterpolator;
171 
172   /** The global timestep. */
173   TimeStepType m_TimeStep;
174 
175   /** Threshold below which the denominator term is considered zero. */
176   double m_DenominatorThreshold;
177 
178   /** Threshold below which two intensity value are assumed to match. */
179   double m_IntensityDifferenceThreshold;
180 };
181 } // end namespace itk
182 
183 #ifndef ITK_MANUAL_INSTANTIATION
184 #include "itkMeanSquareRegistrationFunction.hxx"
185 #endif
186 
187 #endif
188