ITK  4.6.0
Insight Segmentation and Registration Toolkit
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 
20 from InsightToolkit import *
21 
22 from sys import argv
23 
24 #
25 # Read the fixed and moving images using filenames
26 # from the command line arguments
27 #
28 fixedImageReader = itkImageFileReaderF2_New()
29 movingImageReader = itkImageFileReaderF2_New()
30 
31 fixedImageReader.SetFileName( argv[1] )
32 movingImageReader.SetFileName( argv[2] )
33 
34 fixedImageReader.Update()
35 movingImageReader.Update()
36 
37 fixedImage = fixedImageReader.GetOutput()
38 movingImage = movingImageReader.GetOutput()
39 
40 #
41 # Instantiate the classes for the registration framework
42 #
43 registration = itkImageRegistrationMethodF2F2_New()
44 imageMetric = itkMeanSquaresImageToImageMetricF2F2_New()
45 transform = itkCenteredRigid2DTransform_New()
46 optimizer = itkRegularStepGradientDescentOptimizer_New()
47 interpolator = itkLinearInterpolateImageFunctionF2D_New()
48 
49 registration.SetOptimizer( optimizer.GetPointer() )
50 registration.SetTransform( transform.GetPointer() )
51 registration.SetInterpolator( interpolator.GetPointer() )
52 registration.SetMetric( imageMetric.GetPointer() )
53 registration.SetFixedImage( fixedImage )
54 registration.SetMovingImage( movingImage )
55 registration.SetFixedImageRegion( fixedImage.GetBufferedRegion() )
56 
57 
58 #
59 # Initial transform parameters
60 #
61 transform.SetAngle( 0.0 );
62 
63 # center of the fixed image
64 fixedSpacing = fixedImage.GetSpacing()
65 fixedOrigin = fixedImage.GetOrigin()
66 fixedSize = fixedImage.GetLargestPossibleRegion().GetSize()
67 
68 centerFixed = ( fixedOrigin.GetElement(0) + fixedSpacing.GetElement(0) * fixedSize.GetElement(0) / 2.0,
69  fixedOrigin.GetElement(1) + fixedSpacing.GetElement(1) * fixedSize.GetElement(1) / 2.0 )
70 
71 # center of the moving image
72 movingSpacing = movingImage.GetSpacing()
73 movingOrigin = movingImage.GetOrigin()
74 movingSize = movingImage.GetLargestPossibleRegion().GetSize()
75 
76 centerMoving = ( movingOrigin.GetElement(0) + movingSpacing.GetElement(0) * movingSize.GetElement(0) / 2.0,
77  movingOrigin.GetElement(1) + movingSpacing.GetElement(1) * movingSize.GetElement(1) / 2.0 )
78 
79 # transform center
80 center = transform.GetCenter()
81 center.SetElement( 0, centerFixed[0] )
82 center.SetElement( 1, centerFixed[1] )
83 
84 # transform translation
85 translation = transform.GetTranslation()
86 translation.SetElement( 0, centerMoving[0] - centerFixed[0] )
87 translation.SetElement( 1, centerMoving[1] - centerFixed[1] )
88 
89 initialParameters = transform.GetParameters()
90 
91 print "Initial Parameters: "
92 print "Angle: %f" % (initialParameters.GetElement(0), )
93 print "Center: %f, %f" % ( initialParameters.GetElement(1), initialParameters.GetElement(2) )
94 print "Translation: %f, %f" % (initialParameters.GetElement(3), initialParameters.GetElement(4))
95 
96 registration.SetInitialTransformParameters( initialParameters )
97 
98 #
99 # Define optimizer parameters
100 #
101 
102 # optimizer scale
103 translationScale = 1.0 / 1000.0
104 
105 optimizerScales = itkArrayD( transform.GetNumberOfParameters() )
106 optimizerScales.SetElement(0, 1.0)
107 optimizerScales.SetElement(1, translationScale)
108 optimizerScales.SetElement(2, translationScale)
109 optimizerScales.SetElement(3, translationScale)
110 optimizerScales.SetElement(4, translationScale)
111 
112 optimizer.SetScales( optimizerScales )
113 optimizer.SetMaximumStepLength( 0.1 )
114 optimizer.SetMinimumStepLength( 0.001 )
115 optimizer.SetNumberOfIterations( 200 )
116 
117 #
118 # Iteration Observer
119 #
120 def iterationUpdate():
121  currentParameter = transform.GetParameters()
122  print "M: %f P: %f %f %f %f %f " % ( optimizer.GetValue(),
123  currentParameter.GetElement(0),
124  currentParameter.GetElement(1),
125  currentParameter.GetElement(2),
126  currentParameter.GetElement(3),
127  currentParameter.GetElement(4) )
128 
129 iterationCommand = itkPyCommand_New()
130 iterationCommand.SetCommandCallable( iterationUpdate )
131 optimizer.AddObserver( itkIterationEvent(), iterationCommand.GetPointer() )
132 
133 print "Starting registration"
134 
135 #
136 # Start the registration process
137 #
138 
139 registration.Update()
140 
141 #
142 # Get the final parameters of the transformation
143 #
144 finalParameters = registration.GetLastTransformParameters()
145 
146 print "Final Registration Parameters "
147 print "Angle in radians = %f" % finalParameters.GetElement(0)
148 print "Rotation Center X = %f" % finalParameters.GetElement(1)
149 print "Rotation Center Y = %f" % finalParameters.GetElement(2)
150 print "Translation in X = %f" % finalParameters.GetElement(3)
151 print "Translation in Y = %f" % finalParameters.GetElement(4)
152 
153 # Now, we use the final transform for resampling the moving image.
154 resampler = itkResampleImageFilterF2F2_New()
155 
156 resampler.SetTransform( transform.GetPointer() )
157 resampler.SetInput( movingImage )
158 
159 region = fixedImage.GetLargestPossibleRegion()
160 
161 resampler.SetSize( region.GetSize() )
162 resampler.SetOutputSpacing( fixedImage.GetSpacing() )
163 resampler.SetOutputDirection( fixedImage.GetDirection() )
164 resampler.SetOutputOrigin( fixedImage.GetOrigin() )
165 resampler.SetDefaultPixelValue( 100 )
166 
167 #
168 # Cast for output
169 #
170 outputCast = itkRescaleIntensityImageFilterF2US2_New()
171 outputCast.SetInput( resampler.GetOutput() )
172 outputCast.SetOutputMinimum( 0 )
173 outputCast.SetOutputMaximum( 65535 )
174 
175 writer = itkImageFileWriterUS2_New()
176 
177 writer.SetFileName( argv[3] )
178 writer.SetInput( outputCast.GetOutput() )
179 writer.Update()