ITK  5.4.0
Insight Toolkit
Examples/RegistrationITKv4/ImageRegistration5.py
1# ==========================================================================
2#
3# Copyright NumFOCUS
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# https://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
19import itk
20from sys import argv
21
22
23#
24# Check input parameters
25# INPUTS(fixedImage): {BrainProtonDensitySliceBorder20.png}
26# INPUTS(movingImage): {BrainProtonDensitySliceRotated10.png}
27#
28if len(argv) < 4:
29 print("Missing Parameters")
30 print(
31 "Usage: ImageRegistration5.py fixedImageFile movingImageFile outputImagefile"
32 )
33 exit()
34
35
36#
37# Define data types
38#
39FixedImageType = itk.Image[itk.F, 2]
40MovingImageType = itk.Image[itk.F, 2]
41TransformType = itk.CenteredRigid2DTransform[itk.D]
43RegistrationType = itk.ImageRegistrationMethodv4[FixedImageType, MovingImageType]
44MetricType = itk.MeanSquaresImageToImageMetricv4[FixedImageType, MovingImageType]
45
46
47#
48# Read the fixed and moving images using filenames
49# from the command line arguments
50#
51fixedImageReader = itk.ImageFileReader[FixedImageType].New()
52movingImageReader = itk.ImageFileReader[MovingImageType].New()
53
54fixedImageReader.SetFileName(argv[1])
55movingImageReader.SetFileName(argv[2])
56
57fixedImageReader.Update()
58movingImageReader.Update()
59
60fixedImage = fixedImageReader.GetOutput()
61movingImage = movingImageReader.GetOutput()
62
63
64#
65# Instantiate the classes for the registration framework
66#
67registration = RegistrationType.New()
68imageMetric = MetricType.New()
69transform = TransformType.New()
70optimizer = OptimizerType.New()
71
72registration.SetOptimizer(optimizer)
73registration.SetMetric(imageMetric)
74
75registration.SetFixedImage(fixedImage)
76registration.SetMovingImage(movingImage)
77
78
79#
80# Initial transform parameters
81#
82transform.SetAngle(0.0)
83
84# center of the fixed image
85fixedSpacing = fixedImage.GetSpacing()
86fixedOrigin = fixedImage.GetOrigin()
87fixedSize = fixedImage.GetLargestPossibleRegion().GetSize()
88
89centerFixed = (
90 fixedOrigin.GetElement(0)
91 + fixedSpacing.GetElement(0) * fixedSize.GetElement(0) / 2.0,
92 fixedOrigin.GetElement(1)
93 + fixedSpacing.GetElement(1) * fixedSize.GetElement(1) / 2.0,
94)
95
96# center of the moving image
97movingSpacing = movingImage.GetSpacing()
98movingOrigin = movingImage.GetOrigin()
99movingSize = movingImage.GetLargestPossibleRegion().GetSize()
100
101centerMoving = (
102 movingOrigin.GetElement(0)
103 + movingSpacing.GetElement(0) * movingSize.GetElement(0) / 2.0,
104 movingOrigin.GetElement(1)
105 + movingSpacing.GetElement(1) * movingSize.GetElement(1) / 2.0,
106)
107
108# transform center
109center = transform.GetCenter()
110center.SetElement(0, centerFixed[0])
111center.SetElement(1, centerFixed[1])
112
113# transform translation
114translation = transform.GetTranslation()
115translation.SetElement(0, centerMoving[0] - centerFixed[0])
116translation.SetElement(1, centerMoving[1] - centerFixed[1])
117
118registration.SetInitialTransform(transform)
119
120initialParameters = transform.GetParameters()
121
122print("Initial Parameters: ")
123print(f"Angle: {initialParameters.GetElement(0):f}")
124print(
125 f"Center: {initialParameters.GetElement(1):f}, {initialParameters.GetElement(2):f}"
126)
127print(
128 "Translation: %f, %f"
129 % (initialParameters.GetElement(3), initialParameters.GetElement(4))
130)
131
132
133#
134# Define optimizer parameters
135#
136
137# optimizer scale
138translationScale = 1.0 / 1000.0
139
140optimizerScales = itk.OptimizerParameters[itk.D](transform.GetNumberOfParameters())
141optimizerScales.SetElement(0, 1.0)
142optimizerScales.SetElement(1, translationScale)
143optimizerScales.SetElement(2, translationScale)
144optimizerScales.SetElement(3, translationScale)
145optimizerScales.SetElement(4, translationScale)
146
147optimizer.SetScales(optimizerScales)
148
149optimizer.SetRelaxationFactor(0.6)
150optimizer.SetLearningRate(0.1)
151optimizer.SetMinimumStepLength(0.001)
152optimizer.SetNumberOfIterations(200)
153
154
155#
156# One level registration process without shrinking and smoothing.
157#
158registration.SetNumberOfLevels(1)
159registration.SetSmoothingSigmasPerLevel([0])
160registration.SetShrinkFactorsPerLevel([1])
161
162
163#
164# Iteration Observer
165#
166def iterationUpdate():
167 currentParameter = transform.GetParameters()
168 print(
169 "M: %f P: %f %f %f %f %f "
170 % (
171 optimizer.GetValue(),
172 currentParameter.GetElement(0),
173 currentParameter.GetElement(1),
174 currentParameter.GetElement(2),
175 currentParameter.GetElement(3),
176 currentParameter.GetElement(4),
177 )
178 )
179
180
181iterationCommand = itk.PyCommand.New()
182iterationCommand.SetCommandCallable(iterationUpdate)
183optimizer.AddObserver(itk.IterationEvent(), iterationCommand)
184
185print("Starting registration")
186
187
188#
189# Start the registration process
190#
191registration.Update()
192
193
194#
195# Get the final parameters of the transformation
196#
197finalParameters = registration.GetOutput().Get().GetParameters()
198
199print("Final Registration Parameters ")
200print(f"Angle in radians = {finalParameters.GetElement(0):f}")
201print(f"Rotation Center X = {finalParameters.GetElement(1):f}")
202print(f"Rotation Center Y = {finalParameters.GetElement(2):f}")
203print(f"Translation in X = {finalParameters.GetElement(3):f}")
204print(f"Translation in Y = {finalParameters.GetElement(4):f}")
205
206
207#
208# Now, we use the final transform for resampling the
209# moving image.
210#
211resampler = itk.ResampleImageFilter[MovingImageType, FixedImageType].New()
212resampler.SetTransform(registration.GetTransform())
213resampler.SetInput(movingImageReader.GetOutput())
214
215region = fixedImage.GetLargestPossibleRegion()
216
217resampler.SetSize(region.GetSize())
218resampler.SetOutputOrigin(fixedImage.GetOrigin())
219resampler.SetOutputSpacing(fixedImage.GetSpacing())
220resampler.SetOutputDirection(fixedImage.GetDirection())
221resampler.SetDefaultPixelValue(100)
222
223OutputImageType = itk.Image[itk.UC, 2]
224outputCast = itk.CastImageFilter[FixedImageType, OutputImageType].New()
225outputCast.SetInput(resampler.GetOutput())
226
227
228#
229# Write the resampled image
230#
231writer = itk.ImageFileWriter[OutputImageType].New()
232writer.SetFileName(argv[3])
233writer.SetInput(outputCast.GetOutput())
234writer.Update()
Casts input pixels to output pixel type.
CenteredRigid2DTransform of a vector space (e.g. space coordinates)
Data source that reads image data from a single file.
Writes image data to a single file.
Interface method for the current registration framework.
Templated n-dimensional image class.
Definition: itkImage.h:89
Class implementing a mean squares metric.
Class to hold and manage different parameter types used during optimization.
Resample an image via a coordinate transform.
static Pointer New()