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