ITK  4.9.0
Insight Segmentation and Registration Toolkit
Examples/RegistrationITKv4/ImageRegistration5.py
1 #==========================================================================
2 #
3 # Copyright Insight Software Consortium
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 'Usage: ImageRegistration3.py fixedImageFile movingImageFile outputImagefile'
31  exit()
32 
33 
34 #
35 # Define data types
36 #
37 FixedImageType = itk.Image[itk.F, 2]
38 MovingImageType = itk.Image[itk.F, 2]
39 TransformType = itk.CenteredRigid2DTransform[itk.D]
40 OptimizerType = itk.RegularStepGradientDescentOptimizerv4[itk.D]
41 RegistrationType = itk.ImageRegistrationMethodv4[FixedImageType,
42  MovingImageType]
43 MetricType = itk.MeanSquaresImageToImageMetricv4[FixedImageType,
44  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 = ( fixedOrigin.GetElement(0) + fixedSpacing.GetElement(0) * fixedSize.GetElement(0) / 2.0,
90  fixedOrigin.GetElement(1) + fixedSpacing.GetElement(1) * fixedSize.GetElement(1) / 2.0 )
91 
92 # center of the moving image
93 movingSpacing = movingImage.GetSpacing()
94 movingOrigin = movingImage.GetOrigin()
95 movingSize = movingImage.GetLargestPossibleRegion().GetSize()
96 
97 centerMoving = ( movingOrigin.GetElement(0) + movingSpacing.GetElement(0) * movingSize.GetElement(0) / 2.0,
98  movingOrigin.GetElement(1) + movingSpacing.GetElement(1) * movingSize.GetElement(1) / 2.0 )
99 
100 # transform center
101 center = transform.GetCenter()
102 center.SetElement( 0, centerFixed[0] )
103 center.SetElement( 1, centerFixed[1] )
104 
105 # transform translation
106 translation = transform.GetTranslation()
107 translation.SetElement( 0, centerMoving[0] - centerFixed[0] )
108 translation.SetElement( 1, centerMoving[1] - centerFixed[1] )
109 
110 registration.SetInitialTransform(transform)
111 
112 initialParameters = transform.GetParameters()
113 
114 print "Initial Parameters: "
115 print "Angle: %f" % (initialParameters.GetElement(0), )
116 print "Center: %f, %f" % ( initialParameters.GetElement(1), initialParameters.GetElement(2) )
117 print "Translation: %f, %f" % (initialParameters.GetElement(3), initialParameters.GetElement(4))
118 
119 
120 #
121 # Define optimizer parameters
122 #
123 
124 # optimizer scale
125 translationScale = 1.0 / 1000.0
126 
127 optimizerScales = itk.OptimizerParameters[itk.D](transform.GetNumberOfParameters())
128 optimizerScales.SetElement(0, 1.0)
129 optimizerScales.SetElement(1, translationScale)
130 optimizerScales.SetElement(2, translationScale)
131 optimizerScales.SetElement(3, translationScale)
132 optimizerScales.SetElement(4, translationScale)
133 
134 optimizer.SetScales( optimizerScales )
135 
136 optimizer.SetRelaxationFactor( 0.6 );
137 optimizer.SetLearningRate( 0.1 );
138 optimizer.SetMinimumStepLength( 0.001 );
139 optimizer.SetNumberOfIterations( 200 );
140 
141 
142 #
143 # One level registration process without shrinking and smoothing.
144 #
145 registration.SetNumberOfLevels(1)
146 registration.SetSmoothingSigmasPerLevel([0])
147 registration.SetShrinkFactorsPerLevel([1])
148 
149 
150 #
151 # Iteration Observer
152 #
153 def iterationUpdate():
154  currentParameter = transform.GetParameters()
155  print "M: %f P: %f %f %f %f %f " % ( optimizer.GetValue(),
156  currentParameter.GetElement(0),
157  currentParameter.GetElement(1),
158  currentParameter.GetElement(2),
159  currentParameter.GetElement(3),
160  currentParameter.GetElement(4) )
161 
162 iterationCommand = itk.PyCommand.New()
163 iterationCommand.SetCommandCallable( iterationUpdate )
164 optimizer.AddObserver( itk.IterationEvent(), iterationCommand )
165 
166 print "Starting registration"
167 
168 
169 #
170 # Start the registration process
171 #
172 registration.Update()
173 
174 
175 #
176 # Get the final parameters of the transformation
177 #
178 finalParameters = registration.GetOutput().Get().GetParameters()
179 
180 print "Final Registration Parameters "
181 print "Angle in radians = %f" % finalParameters.GetElement(0)
182 print "Rotation Center X = %f" % finalParameters.GetElement(1)
183 print "Rotation Center Y = %f" % finalParameters.GetElement(2)
184 print "Translation in X = %f" % finalParameters.GetElement(3)
185 print "Translation in Y = %f" % finalParameters.GetElement(4)
186 
187 
188 #
189 # Now, we use the final transform for resampling the
190 # moving image.
191 #
192 resampler = itk.ResampleImageFilter[MovingImageType,FixedImageType].New()
193 resampler.SetTransform(registration.GetTransform())
194 resampler.SetInput(movingImageReader.GetOutput())
195 
196 region = fixedImage.GetLargestPossibleRegion()
197 
198 resampler.SetSize(region.GetSize())
199 resampler.SetOutputOrigin(fixedImage.GetOrigin())
200 resampler.SetOutputSpacing(fixedImage.GetSpacing())
201 resampler.SetOutputDirection(fixedImage.GetDirection())
202 resampler.SetDefaultPixelValue(100)
203 
204 OutputImageType = itk.Image[itk.UC, 2]
205 outputCast = itk.CastImageFilter[FixedImageType, OutputImageType].New()
206 outputCast.SetInput(resampler.GetOutput())
207 
208 
209 #
210 # Write the resampled image
211 #
212 writer = itk.ImageFileWriter[OutputImageType].New()
213 writer.SetFileName( argv[3] )
214 writer.SetInput( outputCast.GetOutput() )
215 writer.Update()