ITK  4.9.0
Insight Segmentation and Registration Toolkit
SphinxExamples/src/Registration/Common/Perform2DTranslationRegistrationWithMeanSquares/Code.py
1 #!/usr/bin/env python
2 
3 #==========================================================================
4 #
5 # Copyright Insight Software Consortium
6 #
7 # Licensed under the Apache License, Version 2.0 (the "License");
8 # you may not use this file except in compliance with the License.
9 # You may obtain a copy of the License at
10 #
11 # http://www.apache.org/licenses/LICENSE-2.0.txt
12 #
13 # Unless required by applicable law or agreed to in writing, software
14 # distributed under the License is distributed on an "AS IS" BASIS,
15 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 # See the License for the specific language governing permissions and
17 # limitations under the License.
18 #
19 #==========================================================================*/
20 
21 import sys
22 import itk
23 
24 from distutils.version import StrictVersion as VS
25 if VS(itk.Version.GetITKVersion()) < VS("4.9.0"):
26  print("ITK 4.9.0 is required.")
27  sys.exit(1)
28 
29 if len(sys.argv) != 6:
30  print("Usage: " + sys.argv[0] + " <fixedImageFile> <movingImageFile> "
31  "<outputImagefile> <differenceImageAfter> <differenceImageBefore>")
32  sys.exit(1)
33 
34 fixedImage = sys.argv[1]
35 movingImage = sys.argv[2]
36 outputImage = sys.argv[3]
37 differenceImageAfter = sys.argv[4]
38 differenceImageBefore = sys.argv[5]
39 
40 Dimension = 2
41 PixelType = itk.F
42 
43 FixedImageType = itk.Image[PixelType, Dimension]
44 MovingImageType = itk.Image[PixelType, Dimension]
45 
46 TransformType = itk.TranslationTransform[itk.D, Dimension]
47 OptimizerType = itk.RegularStepGradientDescentOptimizerv4[itk.D]
48 
50  FixedImageType, MovingImageType]
51 
52 RegistrationType = itk.ImageRegistrationMethodv4[
53  FixedImageType, MovingImageType]
54 
55 metric = MetricType.New()
56 optimizer = OptimizerType.New()
57 transform = TransformType.New()
58 registration = RegistrationType.New()
59 
60 registration.SetMetric(metric)
61 registration.SetOptimizer(optimizer)
62 registration.SetInitialTransform(transform)
63 
64 FixedLinearInterpolatorType = itk.LinearInterpolateImageFunction[
65  FixedImageType, itk.D]
66 MovingLinearInterpolatorType = itk.LinearInterpolateImageFunction[
67  MovingImageType, itk.D]
68 
69 fixedInterpolator = FixedLinearInterpolatorType.New()
70 movingInterpolator = MovingLinearInterpolatorType.New()
71 
72 metric.SetFixedInterpolator(fixedInterpolator)
73 metric.SetMovingInterpolator(movingInterpolator)
74 
75 FixedImageReaderType = itk.ImageFileReader[FixedImageType]
76 MovingImageReaderType = itk.ImageFileReader[MovingImageType]
77 fixedImageReader = FixedImageReaderType.New()
78 movingImageReader = MovingImageReaderType.New()
79 
80 fixedImageReader.SetFileName(fixedImage)
81 movingImageReader.SetFileName(movingImage)
82 
83 registration.SetFixedImage(fixedImageReader.GetOutput())
84 registration.SetMovingImage(movingImageReader.GetOutput())
85 
86 movingInitialTransform = TransformType.New()
87 initialParameters = movingInitialTransform.GetParameters()
88 initialParameters[0] = 0
89 initialParameters[1] = 0
90 movingInitialTransform.SetParameters(initialParameters)
91 
92 registration.SetMovingInitialTransform(movingInitialTransform)
93 
94 identityTransform = TransformType.New()
95 identityTransform.SetIdentity()
96 
97 registration.SetFixedInitialTransform(identityTransform)
98 
99 optimizer.SetLearningRate(4)
100 optimizer.SetMinimumStepLength(0.001)
101 optimizer.SetRelaxationFactor(0.5)
102 optimizer.SetNumberOfIterations(200)
103 
104 registration.SetNumberOfLevels(1)
105 registration.SetSmoothingSigmasPerLevel([0])
106 registration.SetShrinkFactorsPerLevel([1])
107 
108 registration.Update()
109 
110 transform = registration.GetTransform()
111 
112 finalParameters = transform.GetParameters()
113 TranslationAlongX = finalParameters.GetElement(0)
114 TranslationAlongY = finalParameters.GetElement(1)
115 
116 numberOfIterations = optimizer.GetCurrentIteration()
117 
118 bestValue = optimizer.GetValue()
119 
120 print("Result = ")
121 print(" Translation X = " + str(TranslationAlongX))
122 print(" Translation Y = " + str(TranslationAlongY))
123 print(" Iterations = " + str(numberOfIterations))
124 print(" Metric value = " + str(bestValue))
125 
126 CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
127 outputCompositeTransform = CompositeTransformType.New()
128 outputCompositeTransform.AddTransform(movingInitialTransform)
129 outputCompositeTransform.AddTransform(registration.GetModifiableTransform())
130 
131 ResampleFilterType = itk.ResampleImageFilter[MovingImageType, FixedImageType]
132 resampler = ResampleFilterType.New()
133 resampler.SetInput(movingImageReader.GetOutput())
134 resampler.SetTransform(outputCompositeTransform)
135 
136 fixedImage = fixedImageReader.GetOutput()
137 resampler.SetSize(fixedImage.GetLargestPossibleRegion().GetSize())
138 resampler.SetOutputOrigin(fixedImage.GetOrigin())
139 resampler.SetOutputSpacing(fixedImage.GetSpacing())
140 resampler.SetOutputDirection(fixedImage.GetDirection())
141 resampler.SetDefaultPixelValue(100)
142 
143 OutputPixelType = itk.UC
144 OutputImageType = itk.Image[OutputPixelType, Dimension]
145 
146 CastFilterType = itk.CastImageFilter[FixedImageType, OutputImageType]
147 WriterType = itk.ImageFileWriter[OutputImageType]
148 
149 writer = WriterType.New()
150 caster = CastFilterType.New()
151 
152 writer.SetFileName(outputImage)
153 
154 caster.SetInput(resampler.GetOutput())
155 writer.SetInput(caster.GetOutput())
156 writer.Update()
157 
158 DifferenceFilterType = itk.SubtractImageFilter[
159  FixedImageType, FixedImageType, FixedImageType]
160 
161 difference = DifferenceFilterType.New()
162 
163 difference.SetInput1(fixedImageReader.GetOutput())
164 difference.SetInput2(resampler.GetOutput())
165 
166 RescalerType = itk.RescaleIntensityImageFilter[FixedImageType, OutputImageType]
167 intensityRescaler = RescalerType.New()
168 
169 intensityRescaler.SetInput(difference.GetOutput())
170 intensityRescaler.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
171 intensityRescaler.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
172 
173 resampler.SetDefaultPixelValue(1)
174 
175 writer2 = WriterType.New()
176 writer2.SetInput(intensityRescaler.GetOutput())
177 
178 writer2.SetFileName(differenceImageAfter)
179 writer2.Update()
180 
181 resampler.SetTransform(identityTransform)
182 
183 writer2.SetFileName(differenceImageBefore)
184 writer2.Update()