ITK  4.13.0
Insight Segmentation and Registration Toolkit
Examples/RegistrationITKv3/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: ImageRegistration5.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 OutputImageType = itk.Image[itk.UC, 2]
40 TransformType = itk.CenteredRigid2DTransform
41 
42 
43 #
44 # Read the fixed and moving images using filenames
45 # from the command line arguments
46 #
47 fixedImageReader = itk.ImageFileReader[FixedImageType].New()
48 movingImageReader = itk.ImageFileReader[MovingImageType].New()
49 
50 fixedImageReader.SetFileName( argv[1] )
51 movingImageReader.SetFileName( argv[2] )
52 
53 fixedImageReader.Update()
54 movingImageReader.Update()
55 
56 fixedImage = fixedImageReader.GetOutput()
57 movingImage = movingImageReader.GetOutput()
58 
59 
60 #
61 # Instantiate the classes for the registration framework
62 #
63 registration = itk.ImageRegistrationMethod[FixedImageType, MovingImageType].New()
64 imageMetric = itk.MeanSquaresImageToImageMetric[FixedImageType, MovingImageType].New()
65 transform = TransformType.New()
67 interpolator = itk.LinearInterpolateImageFunction[FixedImageType, itk.D].New()
68 
69 registration.SetOptimizer( optimizer )
70 registration.SetTransform( transform )
71 registration.SetInterpolator( interpolator)
72 registration.SetMetric( imageMetric )
73 
74 registration.SetFixedImage( fixedImage )
75 registration.SetMovingImage( movingImage )
76 
77 registration.SetFixedImageRegion( fixedImage.GetBufferedRegion() )
78 
79 
80 #
81 # Initial transform parameters
82 #
83 transform.SetAngle( 0.0 );
84 
85 # center of the fixed image
86 fixedSpacing = fixedImage.GetSpacing()
87 fixedOrigin = fixedImage.GetOrigin()
88 fixedSize = fixedImage.GetLargestPossibleRegion().GetSize()
89 
90 centerFixed = ( fixedOrigin.GetElement(0) + fixedSpacing.GetElement(0) * fixedSize.GetElement(0) / 2.0,
91  fixedOrigin.GetElement(1) + fixedSpacing.GetElement(1) * fixedSize.GetElement(1) / 2.0 )
92 
93 # center of the moving image
94 movingSpacing= movingImage.GetSpacing()
95 movingOrigin = movingImage.GetOrigin()
96 movingSize = movingImage.GetLargestPossibleRegion().GetSize()
97 
98 centerMoving = ( movingOrigin.GetElement(0) + movingSpacing.GetElement(0) * movingSize.GetElement(0) / 2.0,
99  movingOrigin.GetElement(1) + movingSpacing.GetElement(1) * movingSize.GetElement(1) / 2.0 )
100 
101 # transform center
102 center = transform.GetCenter()
103 center.SetElement( 0, centerFixed[0] )
104 center.SetElement( 1, centerFixed[1] )
105 
106 # transform translation
107 translation = transform.GetTranslation()
108 translation.SetElement( 0, centerMoving[0] - centerFixed[0] )
109 translation.SetElement( 1, centerMoving[1] - centerFixed[1] )
110 
111 initialParameters = transform.GetParameters()
112 
113 print "Initial Parameters: "
114 print "Angle: %f" % (initialParameters.GetElement(0), )
115 print "Center: %f, %f" % ( initialParameters.GetElement(1), initialParameters.GetElement(2) )
116 print "Translation: %f, %f" % (initialParameters.GetElement(3), initialParameters.GetElement(4))
117 
118 registration.SetInitialTransformParameters( initialParameters )
119 
120 
121 #
122 # Define optimizer parameters
123 #
124 
125 # optimizer scale
126 translationScale = 1.0 / 1000.0
127 
128 optimizerScales = itk.Array[itk.D](transform.GetNumberOfParameters())
129 optimizerScales.SetElement(0, 1.0)
130 optimizerScales.SetElement(1, translationScale)
131 optimizerScales.SetElement(2, translationScale)
132 optimizerScales.SetElement(3, translationScale)
133 optimizerScales.SetElement(4, translationScale)
134 
135 optimizer.SetScales( optimizerScales )
136 optimizer.SetMaximumStepLength( 0.1 )
137 optimizer.SetMinimumStepLength( 0.001 )
138 optimizer.SetNumberOfIterations( 200 )
139 
140 
141 #
142 # Iteration Observer
143 #
144 def iterationUpdate():
145  currentParameter = transform.GetParameters()
146  print "M: %f P: %f %f %f %f %f " % ( optimizer.GetValue(),
147  currentParameter.GetElement(0),
148  currentParameter.GetElement(1),
149  currentParameter.GetElement(2),
150  currentParameter.GetElement(3),
151  currentParameter.GetElement(4) )
152 
153 iterationCommand = itk.PyCommand.New()
154 iterationCommand.SetCommandCallable( iterationUpdate )
155 optimizer.AddObserver( itk.IterationEvent(), iterationCommand )
156 
157 print "Starting registration"
158 
159 
160 #
161 # Start the registration process
162 #
163 registration.Update()
164 
165 
166 #
167 # Get the final parameters of the transformation
168 #
169 finalParameters = registration.GetLastTransformParameters()
170 
171 print "Final Registration Parameters "
172 print "Angle in radians = %f" % finalParameters.GetElement(0)
173 print "Rotation Center X = %f" % finalParameters.GetElement(1)
174 print "Rotation Center Y = %f" % finalParameters.GetElement(2)
175 print "Translation in X = %f" % finalParameters.GetElement(3)
176 print "Translation in Y = %f" % finalParameters.GetElement(4)
177 
178 # Now, we use the final transform for resampling the moving image.
179 resampler = itk.ResampleImageFilter[MovingImageType, FixedImageType].New()
180 resampler.SetTransform( transform )
181 resampler.SetInput( movingImage )
182 
183 region = fixedImage.GetLargestPossibleRegion()
184 
185 resampler.SetSize( region.GetSize() )
186 
187 resampler.SetOutputSpacing( fixedImage.GetSpacing() )
188 resampler.SetOutputOrigin( fixedImage.GetOrigin() )
189 resampler.SetOutputDirection( fixedImage.GetDirection() )
190 resampler.SetDefaultPixelValue( 100 )
191 
192 
193 #
194 # Cast for output
195 #
196 outputCast = itk.RescaleIntensityImageFilter[FixedImageType, OutputImageType].New()
197 outputCast.SetInput(resampler.GetOutput())
198 
199 
200 #
201 # Write the resampled image
202 #
203 writer = itk.ImageFileWriter[OutputImageType].New()
204 
205 writer.SetFileName( argv[3] )
206 writer.SetInput( outputCast.GetOutput() )
207 writer.Update()