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 itkESMDemonsRegistrationFunction_h
19 #define itkESMDemonsRegistrationFunction_h
20 
21 #include "itkPDEDeformableRegistrationFunction.h"
22 #include "itkCentralDifferenceImageFunction.h"
23 #include "itkWarpImageFilter.h"
24 #include <mutex>
25 
26 namespace itk
27 {
28 /**
29  * \class ESMDemonsRegistrationFunction
30  *
31  * \brief Fast implementation of the symmetric demons registration force
32  *
33  * This class provides a substantially faster implementation of the
34  * symmetric demons registration force. Speed is improved by keeping
35  * a deformed copy of the moving image for gradient evaluation.
36  *
37  * Symmetric forces simply means using the mean of the gradient
38  * of the fixed image and the gradient of the warped moving
39  * image.
40  *
41  * Note that this class also enables the use of fixed, mapped moving
42  * and warped moving images forces by using a call to SetUseGradientType
43  *
44  * The moving image should not be saturated. We indeed use
45  * NumericTraits<MovingPixelType>::Max() as a special value.
46  *
47  * \author Tom Vercauteren, INRIA & Mauna Kea Technologies
48  *
49  * This implementation was taken from the Insight Journal paper:
50  * https://hdl.handle.net/1926/510
51  *
52  * \sa SymmetricForcesDemonsRegistrationFunction
53  * \sa SymmetricForcesDemonsRegistrationFilter
54  * \sa DemonsRegistrationFilter
55  * \sa DemonsRegistrationFunction
56  * \ingroup FiniteDifferenceFunctions
57  *
58  * \ingroup ITKPDEDeformableRegistration
59  */
60 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
61 class ITK_TEMPLATE_EXPORT ESMDemonsRegistrationFunction:
62   public PDEDeformableRegistrationFunction< TFixedImage,
63                                             TMovingImage, TDisplacementField >
64 {
65 public:
66   ITK_DISALLOW_COPY_AND_ASSIGN(ESMDemonsRegistrationFunction);
67 
68   /** Standard class type aliases. */
69   using Self = ESMDemonsRegistrationFunction;
70   using Superclass = PDEDeformableRegistrationFunction<
71     TFixedImage, TMovingImage, TDisplacementField >;
72 
73   using Pointer = SmartPointer< Self >;
74   using ConstPointer = SmartPointer< const Self >;
75 
76   /** Method for creation through the object factory. */
77   itkNewMacro(Self);
78 
79   /** Run-time type information (and related methods). */
80   itkTypeMacro(ESMDemonsRegistrationFunction,
81                PDEDeformableRegistrationFunction);
82 
83   /** MovingImage image type. */
84   using MovingImageType = typename Superclass::MovingImageType;
85   using MovingImagePointer = typename Superclass::MovingImagePointer;
86   using MovingPixelType = typename MovingImageType::PixelType;
87 
88   /** FixedImage image type. */
89   using FixedImageType = typename Superclass::FixedImageType;
90   using FixedImagePointer = typename Superclass::FixedImagePointer;
91   using IndexType = typename FixedImageType::IndexType;
92   using SizeType = typename FixedImageType::SizeType;
93   using SpacingType = typename FixedImageType::SpacingType;
94   using DirectionType = typename FixedImageType::DirectionType;
95 
96   /** Deformation field type. */
97   using DisplacementFieldType = typename Superclass::DisplacementFieldType;
98   using DisplacementFieldTypePointer = typename Superclass::DisplacementFieldTypePointer;
99 
100   /** Inherit some enums from the superclass. */
101   static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
102 
103   /** Inherit some enums from the superclass. */
104   using PixelType = typename Superclass::PixelType;
105   using RadiusType = typename Superclass::RadiusType;
106   using NeighborhoodType = typename Superclass::NeighborhoodType;
107   using FloatOffsetType = typename Superclass::FloatOffsetType;
108   using TimeStepType = typename Superclass::TimeStepType;
109 
110   /** Interpolator type. */
111   using CoordRepType = double;
112   using InterpolatorType = InterpolateImageFunction<
113     MovingImageType, CoordRepType >;
114   using InterpolatorPointer = typename InterpolatorType::Pointer;
115   using PointType = typename InterpolatorType::PointType;
116   using DefaultInterpolatorType = LinearInterpolateImageFunction<
117     MovingImageType, CoordRepType >;
118 
119   /** Warper type */
120   using WarperType = WarpImageFilter<
121     MovingImageType,
122     MovingImageType, DisplacementFieldType >;
123 
124   using WarperPointer = typename WarperType::Pointer;
125 
126   /** Covariant vector type. */
127   using CovariantVectorType = CovariantVector< double, Self::ImageDimension >;
128 
129   /** Fixed image gradient calculator type. */
130   using GradientCalculatorType = CentralDifferenceImageFunction< FixedImageType >;
131   using GradientCalculatorPointer = typename GradientCalculatorType::Pointer;
132 
133   /** Moving image gradient (unwarped) calculator type. */
134   using MovingImageGradientCalculatorType =
135       CentralDifferenceImageFunction< MovingImageType, CoordRepType >;
136   using MovingImageGradientCalculatorPointer = typename MovingImageGradientCalculatorType::Pointer;
137 
138   /** Set the moving image interpolator. */
SetMovingImageInterpolator(InterpolatorType * ptr)139   void SetMovingImageInterpolator(InterpolatorType *ptr)
140   { m_MovingImageInterpolator = ptr; m_MovingImageWarper->SetInterpolator(ptr); }
141 
142   /** Get the moving image interpolator. */
GetMovingImageInterpolator()143   InterpolatorType * GetMovingImageInterpolator()
144   { return m_MovingImageInterpolator; }
145 
146   /** This class uses a constant timestep of 1. */
ComputeGlobalTimeStep(void * itkNotUsed (GlobalData))147   TimeStepType ComputeGlobalTimeStep( void *itkNotUsed(GlobalData) ) const override
148   { return m_TimeStep; }
149 
150   /** Return a pointer to a global data structure that is passed to
151    * this object from the solver at each calculation.  */
GetGlobalDataPointer()152   void * GetGlobalDataPointer() const override
153   {
154     auto * global = new GlobalDataStruct();
155 
156     global->m_SumOfSquaredDifference  = 0.0;
157     global->m_NumberOfPixelsProcessed = 0L;
158     global->m_SumOfSquaredChange      = 0;
159     return global;
160   }
161 
162   /** Release memory for global data structure. */
163   void ReleaseGlobalDataPointer(void *GlobalData) const override;
164 
165   /** Set the object's state before each iteration. */
166   void InitializeIteration() override;
167 
168   /** This method is called by a finite difference solver image filter at
169    * each pixel that does not lie on a data set boundary */
170   PixelType  ComputeUpdate( const NeighborhoodType & neighborhood,
171                                     void *globalData,
172                                     const FloatOffsetType & offset = FloatOffsetType(0.0) ) override;
173 
174   /** Get the metric value. The metric value is the mean square difference
175    * in intensity between the fixed image and transforming moving image
176    * computed over the the overlapping region between the two images. */
GetMetric()177   virtual double GetMetric() const
178   { return m_Metric; }
179 
180   /** Get the rms change in deformation field. */
GetRMSChange()181   virtual const double & GetRMSChange() const
182   { return m_RMSChange; }
183 
184   /** Set/Get the threshold below which the absolute difference of
185    * intensity yields a match. When the intensities match between a
186    * moving and fixed image pixel, the update vector (for that
187    * iteration) will be the zero vector. Default is 0.001. */
188   virtual void SetIntensityDifferenceThreshold(double);
189 
190   virtual double GetIntensityDifferenceThreshold() const;
191 
192   /** Set/Get the maximum update step length. In Thirion this is 0.5.
193    *  Setting it to 0 implies no restriction (beware of numerical
194    *  instability in this case. */
SetMaximumUpdateStepLength(double sm)195   virtual void SetMaximumUpdateStepLength(double sm)
196   {
197     this->m_MaximumUpdateStepLength = sm;
198   }
199 
GetMaximumUpdateStepLength()200   virtual double GetMaximumUpdateStepLength() const
201   {
202     return this->m_MaximumUpdateStepLength;
203   }
204 
205   /** Type of available image forces */
206   enum GradientType {
207     Symmetric = 0,
208     Fixed = 1,
209     WarpedMoving = 2,
210     MappedMoving = 3
211     };
212 
213   /** Set/Get the type of used image forces */
SetUseGradientType(GradientType gtype)214   virtual void SetUseGradientType(GradientType gtype)
215   { m_UseGradientType = gtype; }
GetUseGradientType()216   virtual GradientType GetUseGradientType() const
217   { return m_UseGradientType; }
218 
219 protected:
220   ESMDemonsRegistrationFunction();
221   ~ESMDemonsRegistrationFunction() override = default;
222   void PrintSelf(std::ostream & os, Indent indent) const override;
223 
224   /** FixedImage image neighborhood iterator type. */
225   using FixedImageNeighborhoodIteratorType = ConstNeighborhoodIterator< FixedImageType >;
226 
227   /** A global data type for this class of equation. Used to store
228    * iterators for the fixed image. */
229   struct GlobalDataStruct {
230     double m_SumOfSquaredDifference;
231     SizeValueType m_NumberOfPixelsProcessed;
232     double m_SumOfSquaredChange;
233   };
234 
235 private:
236   /** Cache fixed image information. */
237   PointType     m_FixedImageOrigin;
238   SpacingType   m_FixedImageSpacing;
239   DirectionType m_FixedImageDirection;
240   double        m_Normalizer;
241 
242   /** Function to compute derivatives of the fixed image. */
243   GradientCalculatorPointer m_FixedImageGradientCalculator;
244 
245   /** Function to compute derivatives of the moving image (unwarped). */
246   MovingImageGradientCalculatorPointer m_MappedMovingImageGradientCalculator;
247 
248   GradientType m_UseGradientType;
249 
250   /** Function to interpolate the moving image. */
251   InterpolatorPointer m_MovingImageInterpolator;
252 
253   /** Filter to warp moving image for fast gradient computation. */
254   WarperPointer m_MovingImageWarper;
255 
256   MovingImageType *m_MovingImageWarperOutput;
257 
258   /** The global timestep. */
259   TimeStepType m_TimeStep;
260 
261   /** Threshold below which the denominator term is considered zero. */
262   double m_DenominatorThreshold;
263 
264   /** Threshold below which two intensity value are assumed to match. */
265   double m_IntensityDifferenceThreshold;
266 
267   /** Maximum update step length in pixels (default is 0.5 as in Thirion). */
268   double m_MaximumUpdateStepLength;
269 
270   /** The metric value is the mean square difference in intensity between
271    * the fixed image and transforming moving image computed over the
272    * the overlapping region between the two images. */
273   mutable double        m_Metric;
274   mutable double        m_SumOfSquaredDifference;
275   mutable SizeValueType m_NumberOfPixelsProcessed;
276   mutable double        m_RMSChange;
277   mutable double        m_SumOfSquaredChange;
278 
279   /** Mutex lock to protect modification to metric. */
280   mutable std::mutex m_MetricCalculationLock;
281 };
282 } // end namespace itk
283 
284 #ifndef ITK_MANUAL_INSTANTIATION
285 #include "itkESMDemonsRegistrationFunction.hxx"
286 #endif
287 
288 #endif
289