ITK  6.0.0
Insight Toolkit
SphinxExamples/src/Registration/Common/GlobalRegistrationOfTwoImages/Code.py
1 import itk
2 
3 Dimension = 2
4 PixelType = itk.UC
5 Double = itk.D
6 Float = itk.F
7 ImageType = itk.Image[PixelType, Dimension]
8 
9 
10 def main():
11  # The transform that will map the fixed image into the moving image.
12  TransformType = itk.TranslationTransform[Double, Dimension]
13 
14  # An optimizer is required to explore the parameter space of the transform
15  # in search of optimal values of the metric.
17 
18  # The metric will compare how well the two images match each other. Metric
19  # types are usually parameterized by the image types as it can be seen in
20  # the following type declaration.
21  MetricType = itk.MeanSquaresImageToImageMetric[ImageType, ImageType]
22 
23  # Finally, the type of the interpolator is declared. The interpolator will
24  # evaluate the intensities of the moving image at non-grid positions.
25  InterpolatorType = itk.LinearInterpolateImageFunction[ImageType, Double]
26 
27  # The registration method type is instantiated using the types of the
28  # fixed and moving images. This class is responsible for interconnecting
29  # all the components that we have described so far.
30  RegistrationType = itk.ImageRegistrationMethod[ImageType, ImageType]
31 
32  # Create components
33  metric = MetricType.New()
34  transform = TransformType.New()
35  optimizer = OptimizerType.New()
36  interpolator = InterpolatorType.New()
37  registration = RegistrationType.New()
38 
39  # Each component is now connected to the instance of the registration method.
40  registration.SetMetric(metric)
41  registration.SetOptimizer(optimizer)
42  registration.SetTransform(transform)
43  registration.SetInterpolator(interpolator)
44 
45  # Get the two images
46  fixedImage = ImageType.New()
47  movingImage = ImageType.New()
48 
49  CreateSphereImage(fixedImage)
50  CreateEllipseImage(movingImage)
51 
52  # Write the two synthetic inputs
53  itk.imwrite(fixedImage, "fixed.png")
54  itk.imwrite(movingImage, "moving.png")
55 
56  # Set the registration inputs
57  registration.SetFixedImage(fixedImage)
58  registration.SetMovingImage(movingImage)
59 
60  registration.SetFixedImageRegion(fixedImage.GetLargestPossibleRegion())
61 
62  # Initialize the transform
63  initialParameters = itk.OptimizerParameters[itk.D](
64  transform.GetNumberOfParameters()
65  )
66 
67  initialParameters[0] = 0.0 # Initial offset along X
68  initialParameters[1] = 0.0 # Initial offset along Y
69 
70  registration.SetInitialTransformParameters(initialParameters)
71 
72  optimizer.SetMaximumStepLength(4.00)
73  optimizer.SetMinimumStepLength(0.01)
74 
75  # Set a stopping criterion
76  optimizer.SetNumberOfIterations(200)
77 
78  # Connect an observer
79  # surface = dict()
80  # def print_iteration():
81  # surface[tuple(optimizer.GetCurrentPosition())] = optimizer.GetValue()
82  # print(surface)
83 
84  # optimizer.AddObserver(itk.IterationEvent(), print_iteration)
85 
86  try:
87  registration.Update()
88  except:
89  print("ExceptionObject caught !")
90  return
91 
92  # The result of the registration process is an array of parameters that
93  # defines the spatial transformation in an unique way. This final result is
94  # obtained using the \code{GetLastTransformParameters()} method.
95 
96  finalParameters = registration.GetLastTransformParameters()
97 
98  # In the case of the \doxygen{TranslationTransform}, there is a
99  # straightforward interpretation of the parameters. Each element of the
100  # array corresponds to a translation along one spatial dimension.
101 
102  TranslationAlongX = finalParameters[0]
103  TranslationAlongY = finalParameters[1]
104 
105  # The optimizer can be queried for the actual number of iterations
106  # performed to reach convergence. The \code{GetCurrentIteration()}
107  # method returns this value. A large number of iterations may be an
108  # indication that the maximum step length has been set too small, which
109  # is undesirable since it results in long computational times.
110 
111  numberOfIterations = optimizer.GetCurrentIteration()
112 
113  # The value of the image metric corresponding to the last set of parameters
114  # can be obtained with the \code{GetValue()} method of the optimizer.
115 
116  bestValue = optimizer.GetValue()
117 
118  # Print out results
119  print("Result = ")
120  print(" Translation X = {}".format(TranslationAlongX))
121  print(" Translation Y = {}".format(TranslationAlongY))
122  print(" Iterations = {}".format(numberOfIterations))
123  print(" Metric value = {}".format(bestValue))
124 
125  # It is common, as the last step of a registration task, to use the
126  # resulting transform to map the moving image into the fixed image space.
127  # This is easily done with the \doxygen{ResampleImageFilter}. Please
128  # refer to Section~\ref{sec:ResampleImageFilter} for details on the use
129  # of this filter. First, a ResampleImageFilter type is instantiated
130  # using the image types. It is convenient to use the fixed image type as
131  # the output type since it is likely that the transformed moving image
132  # will be compared with the fixed image.
133  ResampleFilterType = itk.ResampleImageFilter[ImageType, ImageType]
134 
135  # A resampling filter is created and the moving image is connected as its input.
136 
137  resampler = ResampleFilterType.New()
138  resampler.SetInput(movingImage)
139 
140  # The Transform that is produced as output of the Registration method is
141  # also passed as input to the resampling filter. Note the use of the
142  # methods \code{GetOutput()} and \code{Get()}. This combination is needed
143  # here because the registration method acts as a filter whose output is a
144  # transform decorated in the form of a \doxygen{DataObject}. For details in
145  # this construction you may want to read the documentation of the
146  # \doxygen{DataObjectDecorator}.
147 
148  resampler.SetTransform(registration.GetOutput().Get())
149 
150  # As described in Section \ref{sec:ResampleImageFilter}, the
151  # ResampleImageFilter requires additional parameters to be specified, in
152  # particular, the spacing, origin and size of the output image. The default
153  # pixel value is also set to a distinct gray level in order to highlight
154  # the regions that are mapped outside of the moving image.
155 
156  resampler.SetSize(itk.size(fixedImage))
157  resampler.SetOutputOrigin(fixedImage.GetOrigin())
158  resampler.SetOutputSpacing(fixedImage.GetSpacing())
159  resampler.SetOutputDirection(fixedImage.GetDirection())
160  resampler.SetDefaultPixelValue(100)
161 
162  # The output of the filter is passed to a writer that will store the
163  # image in a file. An \doxygen{CastImageFilter} is used to convert the
164  # pixel type of the resampled image to the final type used by the
165  # writer. The cast and writer filters are instantiated below.
166  CastFilterType = itk.CastImageFilter[ImageType, ImageType]
167 
168  caster = CastFilterType.New()
169  caster.SetInput(resampler.GetOutput())
170  itk.imwrite(caster.GetOutput(), "outputPython.png")
171 
172  # The fixed image and the transformed moving image can easily be compared
173  # using the \doxygen{SubtractImageFilter}. This pixel-wise filter computes
174  # the difference between homologous pixels of its two input images.
175 
176  DifferenceFilterType = itk.SubtractImageFilter[ImageType, ImageType, ImageType]
177 
178  difference = DifferenceFilterType.New()
179  difference.SetInput1(fixedImage)
180  difference.SetInput2(resampler.GetOutput())
181 
182  return
183 
184 
185 def CreateEllipseImage(image):
186  EllipseType = itk.EllipseSpatialObject[Dimension]
187 
188  SpatialObjectToImageFilterType = itk.SpatialObjectToImageFilter[
189  itk.SpatialObject[2], ImageType
190  ]
191 
192  imageFilter = SpatialObjectToImageFilterType.New()
193 
194  size = itk.Size[Dimension]()
195  size[0] = 100
196  size[1] = 100
197 
198  imageFilter.SetSize(size)
199 
200  spacing = [1, 1]
201  imageFilter.SetSpacing(spacing)
202 
203  ellipse = EllipseType.New()
204  radiusArray = itk.Array[Float]()
205  radiusArray.SetSize(Dimension)
206  radiusArray[0] = 10
207  radiusArray[1] = 20
208  ellipse.SetRadiusInObjectSpace(radiusArray)
209 
210  TransformType = itk.AffineTransform[Double, Dimension]
211  transform = TransformType.New()
212  transform.SetIdentity()
213 
214  translation = itk.Vector[Float, Dimension]()
215  translation[0] = 65
216  translation[1] = 45
217  transform.Translate(translation, False)
218 
219  ellipse.SetObjectToParentTransform(transform)
220 
221  imageFilter.SetInput(ellipse)
222 
223  ellipse.SetDefaultInsideValue(255)
224  ellipse.SetDefaultOutsideValue(0)
225  imageFilter.SetUseObjectValue(True)
226  imageFilter.SetOutsideValue(0)
227 
228  imageFilter.Update()
229  image.Graft(imageFilter.GetOutput())
230 
231 
232 def CreateSphereImage(image):
233  EllipseType = itk.EllipseSpatialObject[Dimension]
234 
235  SpatialObjectToImageFilterType = itk.SpatialObjectToImageFilter[
236  itk.SpatialObject[2], itk.Image[itk.UC, 2]
237  ]
238 
239  imageFilter = SpatialObjectToImageFilterType.New()
240 
241  size = itk.Size[Dimension]()
242  size[0] = 100
243  size[1] = 100
244  imageFilter.SetSize(size)
245 
246  spacing = [1, 1]
247  imageFilter.SetSpacing(spacing)
248 
249  ellipse = EllipseType.New()
250  radiusArray = itk.Array[Float]()
251  radiusArray.SetSize(Dimension)
252  radiusArray[0] = 10
253  radiusArray[1] = 10
254  ellipse.SetRadiusInObjectSpace(radiusArray)
255 
256  TransformType = itk.AffineTransform[Double, Dimension]
257  transform = TransformType.New()
258  transform.SetIdentity()
259 
260  translation = itk.Vector[PixelType, Dimension]()
261  translation[0] = 50
262  translation[1] = 50
263  transform.Translate(translation, False)
264 
265  ellipse.SetObjectToParentTransform(transform)
266 
267  imageFilter.SetInput(ellipse)
268 
269  ellipse.SetDefaultInsideValue(255)
270  ellipse.SetDefaultOutsideValue(0)
271  imageFilter.SetUseObjectValue(True)
272  imageFilter.SetOutsideValue(0)
273 
274  imageFilter.Update()
275  image.Graft(imageFilter.GetOutput())
276 
277 
278 if __name__ == "__main__":
279  main()
itk::CastImageFilter
Casts input pixels to output pixel type.
Definition: itkCastImageFilter.h:100
itk::OptimizerParameters
Class to hold and manage different parameter types used during optimization.
Definition: itkOptimizerParameters.h:36
itk::Size
Represent a n-dimensional size (bounds) of a n-dimensional image.
Definition: itkSize.h:69
itk::Vector
A templated class holding a n-Dimensional vector.
Definition: itkVector.h:62
itk::ImageRegistrationMethod
Base class for Image Registration Methods.
Definition: itkImageRegistrationMethod.h:70
itk::SpatialObjectToImageFilter
Base class for filters that take a SpatialObject as input and produce an image as output....
Definition: itkSpatialObjectToImageFilter.h:41
itk::AffineTransform
Definition: itkAffineTransform.h:101
itk::EllipseSpatialObject
Definition: itkEllipseSpatialObject.h:38
itk::RegularStepGradientDescentOptimizer
Implement a gradient descent optimizer.
Definition: itkRegularStepGradientDescentOptimizer.h:33
itk::LinearInterpolateImageFunction
Linearly interpolate an image at specified positions.
Definition: itkLinearInterpolateImageFunction.h:51
itk::TranslationTransform
Translation transformation of a vector space (e.g. space coordinates)
Definition: itkTranslationTransform.h:43
itk::SpatialObject
Implementation of the composite pattern.
Definition: itkSpatialObject.h:58
itk::SubtractImageFilter
Pixel-wise subtraction of two images.
Definition: itkSubtractImageFilter.h:68
itk::MeanSquaresImageToImageMetric
TODO.
Definition: itkMeanSquaresImageToImageMetric.h:41
itk::Array
Array class with size defined at construction time.
Definition: itkArray.h:47
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88