forked from InsightSoftwareConsortium/ITK
-
Notifications
You must be signed in to change notification settings - Fork 0
/
itkSyNImageRegistrationMethod.h
269 lines (227 loc) · 11.2 KB
/
itkSyNImageRegistrationMethod.h
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
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkSyNImageRegistrationMethod_h
#define itkSyNImageRegistrationMethod_h
#include "itkImageRegistrationMethodv4.h"
#include "itkImageMaskSpatialObject.h"
#include "itkDisplacementFieldTransform.h"
namespace itk
{
/** \class SyNImageRegistrationMethod
* \brief Interface method for the performing greedy SyN image registration.
*
* For greedy SyN we use \c m_Transform to map the time-parameterized middle
* image to the fixed image (and vice versa using
* \c m_Transform->GetInverseDisplacementField() ). We employ another ivar,
* \c m_InverseTransform, to map the time-parameterized middle image to the
* moving image.
*
* Output: The output is the updated transform which has been added to the
* composite transform.
*
* This implementation is based on the source code in Advanced Normalization Tools (ANTs)
*
* Avants, B. B.; Tustison, N. J.; Song, G.; Cook, P. A.; Klein, A. & Gee, J. C.
* A reproducible evaluation of ANTs similarity metric performance in brain image registration.
* Neuroimage, Penn Image Computing and Science Laboratory, University of Pennsylvania,
* 2011, 54, 2033-2044
*
* The original paper discussing the method is here:
*
* Avants, B. B.; Epstein, C. L.; Grossman, M. & Gee, J. C.
* Symmetric diffeomorphic image registration with cross-correlation:
* evaluating automated labeling of elderly and neurodegenerative brain.
* Med Image Anal, Department of Radiology, University of Pennsylvania,
* 2008, 12, 26-41
*
* The method evolved since that time with crucial contributions from Gang Song and
* Nick Tustison. Though similar in spirit, this implementation is not identical.
*
* \todo Need to allow the fixed image to have a composite transform.
*
* \author Nick Tustison
* \author Brian Avants
*
* \ingroup ITKRegistrationMethodsv4
*/
template <typename TFixedImage,
typename TMovingImage,
typename TOutputTransform = DisplacementFieldTransform<double, TFixedImage::ImageDimension>,
typename TVirtualImage = TFixedImage,
typename TPointSet = PointSet<unsigned int, TFixedImage::ImageDimension>>
class ITK_TEMPLATE_EXPORT SyNImageRegistrationMethod
: public ImageRegistrationMethodv4<TFixedImage, TMovingImage, TOutputTransform, TVirtualImage, TPointSet>
{
public:
ITK_DISALLOW_COPY_AND_MOVE(SyNImageRegistrationMethod);
/** Standard class type aliases. */
using Self = SyNImageRegistrationMethod;
using Superclass = ImageRegistrationMethodv4<TFixedImage, TMovingImage, TOutputTransform, TVirtualImage, TPointSet>;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** ImageDimension constants */
static constexpr unsigned int ImageDimension = TFixedImage::ImageDimension;
/** \see LightObject::GetNameOfClass() */
itkOverrideGetNameOfClassMacro(SyNImageRegistrationMethod);
/** Input type alias for the images. */
using FixedImageType = TFixedImage;
using FixedImagePointer = typename FixedImageType::Pointer;
using typename Superclass::FixedImagesContainerType;
using MovingImageType = TMovingImage;
using MovingImagePointer = typename MovingImageType::Pointer;
using typename Superclass::MovingImagesContainerType;
using typename Superclass::PointSetType;
using PointSetPointer = typename PointSetType::Pointer;
using typename Superclass::PointSetsContainerType;
/** Metric and transform type alias */
using typename Superclass::ImageMetricType;
using ImageMetricPointer = typename ImageMetricType::Pointer;
using MeasureType = typename ImageMetricType::MeasureType;
using ImageMaskSpatialObjectType = ImageMaskSpatialObject<ImageDimension>;
using typename Superclass::FixedImageMaskType;
using FixedMaskImageType = typename ImageMaskSpatialObjectType::ImageType;
using typename Superclass::FixedImageMasksContainerType;
using typename Superclass::MovingImageMaskType;
using MovingMaskImageType = typename ImageMaskSpatialObjectType::ImageType;
using typename Superclass::MovingImageMasksContainerType;
using VirtualImageType = typename Superclass::VirtualImageType;
using typename Superclass::VirtualImageBaseType;
using typename Superclass::VirtualImageBaseConstPointer;
using typename Superclass::MultiMetricType;
using typename Superclass::MetricType;
using MetricPointer = typename MetricType::Pointer;
using typename Superclass::PointSetMetricType;
using typename Superclass::InitialTransformType;
using OutputTransformType = TOutputTransform;
using OutputTransformPointer = typename OutputTransformType::Pointer;
using RealType = typename OutputTransformType::ScalarType;
using DerivativeType = typename OutputTransformType::DerivativeType;
using DerivativeValueType = typename DerivativeType::ValueType;
using DisplacementFieldType = typename OutputTransformType::DisplacementFieldType;
using DisplacementFieldPointer = typename DisplacementFieldType::Pointer;
using DisplacementVectorType = typename DisplacementFieldType::PixelType;
using typename Superclass::CompositeTransformType;
using TransformBaseType = typename CompositeTransformType::TransformType;
using typename Superclass::DecoratedOutputTransformType;
using DecoratedOutputTransformPointer = typename DecoratedOutputTransformType::Pointer;
using DisplacementFieldTransformType = DisplacementFieldTransform<RealType, ImageDimension>;
using DisplacementFieldTransformPointer = typename DisplacementFieldTransformType::Pointer;
using NumberOfIterationsArrayType = Array<SizeValueType>;
/** Set/Get the learning rate. */
itkSetMacro(LearningRate, RealType);
itkGetConstMacro(LearningRate, RealType);
/** Set/Get the number of iterations per level. */
itkSetMacro(NumberOfIterationsPerLevel, NumberOfIterationsArrayType);
itkGetConstMacro(NumberOfIterationsPerLevel, NumberOfIterationsArrayType);
/** Set/Get the convergence threshold */
itkSetMacro(ConvergenceThreshold, RealType);
itkGetConstMacro(ConvergenceThreshold, RealType);
/** Set/Get the convergence window size */
itkSetMacro(ConvergenceWindowSize, unsigned int);
itkGetConstMacro(ConvergenceWindowSize, unsigned int);
/** Let the user control whether we compute metric derivatives in the downsampled or full-res space.
* The default is 'true' --- classic SyN --- but there may be advantages to the other approach.
* Classic SyN did not have this possibility. This implementation will let us explore the question.
*/
itkSetMacro(DownsampleImagesForMetricDerivatives, bool);
itkGetConstMacro(DownsampleImagesForMetricDerivatives, bool);
/** Allow the user to average the gradients in the mid-point domain. Default false.
* One might choose to do this to further reduce bias.
*/
itkSetMacro(AverageMidPointGradients, bool);
itkGetConstMacro(AverageMidPointGradients, bool);
/**
* Get/Set the Gaussian smoothing variance for the update field.
*/
itkSetMacro(GaussianSmoothingVarianceForTheUpdateField, RealType);
itkGetConstReferenceMacro(GaussianSmoothingVarianceForTheUpdateField, RealType);
/**
* Get/Set the Gaussian smoothing variance for the total field.
*/
itkSetMacro(GaussianSmoothingVarianceForTheTotalField, RealType);
itkGetConstReferenceMacro(GaussianSmoothingVarianceForTheTotalField, RealType);
/** Get modifiable FixedToMiddle and MovingToMiddle transforms to save the current state of the registration. */
itkGetModifiableObjectMacro(FixedToMiddleTransform, OutputTransformType);
itkGetModifiableObjectMacro(MovingToMiddleTransform, OutputTransformType);
/** Set FixedToMiddle and MovingToMiddle transforms to restore the registration from a saved state. */
itkSetObjectMacro(FixedToMiddleTransform, OutputTransformType);
itkSetObjectMacro(MovingToMiddleTransform, OutputTransformType);
protected:
SyNImageRegistrationMethod();
~SyNImageRegistrationMethod() override = default;
void
PrintSelf(std::ostream & os, Indent indent) const override;
/** Perform the registration. */
void
GenerateData() override;
/** Handle optimization internally.
* Starts the optimization at each level. Performs a basic gradient descent operation.
*/
virtual void
StartOptimization();
/**
* Initialize by setting the interconnects between the components. Need to override
* in the SyN class since we need to "adapt" the \c m_InverseTransform
*/
void
InitializeRegistrationAtEachLevel(const SizeValueType) override;
virtual DisplacementFieldPointer
ComputeUpdateField(const FixedImagesContainerType,
const PointSetsContainerType,
const TransformBaseType *,
const MovingImagesContainerType,
const PointSetsContainerType,
const TransformBaseType *,
const FixedImageMasksContainerType,
const MovingImageMasksContainerType,
MeasureType &);
virtual DisplacementFieldPointer
ComputeMetricGradientField(const FixedImagesContainerType,
const PointSetsContainerType,
const TransformBaseType *,
const MovingImagesContainerType,
const PointSetsContainerType,
const TransformBaseType *,
const FixedImageMasksContainerType,
const MovingImageMasksContainerType,
MeasureType &);
virtual DisplacementFieldPointer
ScaleUpdateField(const DisplacementFieldType *);
virtual DisplacementFieldPointer
GaussianSmoothDisplacementField(const DisplacementFieldType *, const RealType);
virtual DisplacementFieldPointer
InvertDisplacementField(const DisplacementFieldType *, const DisplacementFieldType * = nullptr);
RealType m_LearningRate{ 0.25 };
OutputTransformPointer m_MovingToMiddleTransform{ nullptr };
OutputTransformPointer m_FixedToMiddleTransform{ nullptr };
RealType m_ConvergenceThreshold{ static_cast<RealType>(1.0e-6) };
unsigned int m_ConvergenceWindowSize{ 10 };
NumberOfIterationsArrayType m_NumberOfIterationsPerLevel{};
bool m_DownsampleImagesForMetricDerivatives{ true };
bool m_AverageMidPointGradients{ false };
private:
RealType m_GaussianSmoothingVarianceForTheUpdateField{ 3.0 };
RealType m_GaussianSmoothingVarianceForTheTotalField{ 0.5 };
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
# include "itkSyNImageRegistrationMethod.hxx"
#endif
#endif