1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
|
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkFastSymmetricForcesDemonsRegistrationFunction.h,v $
Language: C++
Date: $Date: 2005-05-03 23:05:49 $
Version: $Revision: 1.1 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itkFastSymmetricForcesDemonsRegistrationFunction_h_
#define _itkFastSymmetricForcesDemonsRegistrationFunction_h_
#include "itkPDEDeformableRegistrationFunction.h"
#include "itkPoint.h"
#include "itkCovariantVector.h"
#include "itkInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkCentralDifferenceImageFunction.h"
#include "itkWarpImageFilter.h"
namespace itk {
/**
* \class FastSymmetricForcesDemonsRegistrationFunction
*
* This class provides a substantially faster implementation of the
* symmetric demons registration force. Speed is improved by keeping
* a deformed copy of the moving image for gradient evaluation.
*
* \sa SymmetricForcesDemonsRegistrationFunction
* \sa SymmetricForcesDemonsRegistrationFilter
* \sa DemonsRegistrationFilter
* \sa DemonsRegistrationFunction
* \ingroup FiniteDifferenceFunctions
*
* \author Torsten Rohlfing, Neuroscience Program, SRI International
*/
template<class TFixedImage, class TMovingImage, class TDeformationField>
class ITK_EXPORT FastSymmetricForcesDemonsRegistrationFunction :
public PDEDeformableRegistrationFunction< TFixedImage,
TMovingImage, TDeformationField>
{
public:
/** Standard class typedefs. */
typedef FastSymmetricForcesDemonsRegistrationFunction Self;
typedef PDEDeformableRegistrationFunction< TFixedImage,
TMovingImage, TDeformationField > Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro( FastSymmetricForcesDemonsRegistrationFunction,
PDEDeformableRegistrationFunction );
/** MovingImage image type. */
typedef typename Superclass::MovingImageType MovingImageType;
typedef typename Superclass::MovingImagePointer MovingImagePointer;
/** FixedImage image type. */
typedef typename Superclass::FixedImageType FixedImageType;
typedef typename Superclass::FixedImagePointer FixedImagePointer;
typedef typename FixedImageType::IndexType IndexType;
typedef typename FixedImageType::SizeType SizeType;
typedef typename FixedImageType::SpacingType SpacingType;
/** Deformation field type. */
typedef typename Superclass::DeformationFieldType DeformationFieldType;
typedef typename Superclass::DeformationFieldTypePointer
DeformationFieldTypePointer;
/** Inherit some enums from the superclass. */
itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension);
/** Inherit some enums from the superclass. */
typedef typename Superclass::PixelType PixelType;
typedef typename Superclass::RadiusType RadiusType;
typedef typename Superclass::NeighborhoodType NeighborhoodType;
typedef typename Superclass::FloatOffsetType FloatOffsetType;
typedef typename Superclass::TimeStepType TimeStepType;
/** Interpolator type. */
typedef double CoordRepType;
typedef InterpolateImageFunction<MovingImageType,CoordRepType> InterpolatorType;
typedef typename InterpolatorType::Pointer InterpolatorPointer;
typedef typename InterpolatorType::PointType PointType;
typedef LinearInterpolateImageFunction<MovingImageType,CoordRepType>
DefaultInterpolatorType;
/** Warper type */
typedef WarpImageFilter<MovingImageType,MovingImageType,DeformationFieldType> WarperType;
typedef typename WarperType::Pointer WarperPointer;
/** Covariant vector type. */
typedef CovariantVector<double,itkGetStaticConstMacro(ImageDimension)> CovariantVectorType;
/** Gradient calculator types. */
typedef CentralDifferenceImageFunction<FixedImageType> GradientCalculatorType;
typedef typename GradientCalculatorType::Pointer GradientCalculatorPointer;
typedef CentralDifferenceImageFunction<MovingImageType> MovingGradientCalculatorType;
typedef typename MovingGradientCalculatorType::Pointer MovingGradientCalculatorPointer;
/** Set the moving image interpolator. */
void SetMovingImageInterpolator( InterpolatorType * ptr )
{ m_MovingImageInterpolator = ptr; m_MovingImageWarper->SetInterpolator( ptr ); }
/** Get the moving image interpolator. */
InterpolatorType * GetMovingImageInterpolator(void)
{ return m_MovingImageInterpolator; }
/** This class uses a constant timestep of 1. */
virtual TimeStepType ComputeGlobalTimeStep(void * itkNotUsed(GlobalData)) const
{ return m_TimeStep; }
/** Return a pointer to a global data structure that is passed to
* this object from the solver at each calculation. */
virtual void *GetGlobalDataPointer() const
{
GlobalDataStruct *global = new GlobalDataStruct();
global->m_SumOfSquaredDifference = 0.0;
global->m_NumberOfPixelsProcessed = 0L;
global->m_SumOfSquaredChange = 0;
return global;
}
/** Release memory for global data structure. */
virtual void ReleaseGlobalDataPointer( void *GlobalData ) const;
/** Set the object's state before each iteration. */
virtual void InitializeIteration();
/** This method is called by a finite difference solver image filter at
* each pixel that does not lie on a data set boundary */
virtual PixelType ComputeUpdate(const NeighborhoodType &neighborhood,
void *globalData,
const FloatOffsetType &offset = FloatOffsetType(0.0));
/** Get the metric value. The metric value is the mean square difference
* in intensity between the fixed image and transforming moving image
* computed over the the overlapping region between the two images. */
virtual double GetMetric() const
{ return m_Metric; }
/** Get the rms change in deformation field. */
virtual const double &GetRMSChange() const
{ return m_RMSChange; }
/** Set/Get the threshold below which the absolute difference of
* intensity yields a match. When the intensities match between a
* moving and fixed image pixel, the update vector (for that
* iteration) will be the zero vector. Default is 0.001. */
virtual void SetIntensityDifferenceThreshold(double);
virtual double GetIntensityDifferenceThreshold() const;
protected:
FastSymmetricForcesDemonsRegistrationFunction();
~FastSymmetricForcesDemonsRegistrationFunction() {}
void PrintSelf(std::ostream& os, Indent indent) const;
/** FixedImage image neighborhood iterator type. */
typedef ConstNeighborhoodIterator<FixedImageType> FixedImageNeighborhoodIteratorType;
/** A global data type for this class of equation. Used to store
* iterators for the fixed image. */
struct GlobalDataStruct
{
double m_SumOfSquaredDifference;
unsigned long m_NumberOfPixelsProcessed;
double m_SumOfSquaredChange;
};
private:
FastSymmetricForcesDemonsRegistrationFunction(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
/** Cache fixed image information. */
SpacingType m_FixedImageSpacing;
PointType m_FixedImageOrigin;
double m_Normalizer;
/** Function to compute derivatives of the fixed image. */
GradientCalculatorPointer m_FixedImageGradientCalculator;
/** Function to compute derivatives of the fixed image. */
MovingGradientCalculatorPointer m_WarpedMovingImageGradientCalculator;
/** Function to interpolate the moving image. */
InterpolatorPointer m_MovingImageInterpolator;
/** Filter to warp moving image for fast gradient computation. */
WarperPointer m_MovingImageWarper;
/** The global timestep. */
TimeStepType m_TimeStep;
/** Threshold below which the denominator term is considered zero. */
double m_DenominatorThreshold;
/** Threshold below which two intensity value are assumed to match. */
double m_IntensityDifferenceThreshold;
/** The metric value is the mean square difference in intensity between
* the fixed image and transforming moving image computed over the
* the overlapping region between the two images. */
mutable double m_Metric;
mutable double m_SumOfSquaredDifference;
mutable unsigned long m_NumberOfPixelsProcessed;
mutable double m_RMSChange;
mutable double m_SumOfSquaredChange;
/** Mutex lock to protect modification to metric. */
mutable SimpleFastMutexLock m_MetricCalculationLock;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFastSymmetricForcesDemonsRegistrationFunction.txx"
#endif
#endif
|
