ITK  4.8.0
Insight Segmentation and Registration Toolkit
Examples/RegistrationITKv4/ImageRegistration4.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 from InsightToolkit import *
20 
21 from sys import argv
22 
23 
24 
25 #
26 # Read the fixed and moving images using filenames
27 # from the command line arguments
28 #
29 fixedImageReader = itkImageFileReaderF2_New()
30 movingImageReader = itkImageFileReaderF2_New()
31 
32 fixedImageReader.SetFileName( argv[1] )
33 movingImageReader.SetFileName( argv[2] )
34 
35 fixedImageReader.Update()
36 movingImageReader.Update()
37 
38 fixedImage = fixedImageReader.GetOutput()
39 movingImage = movingImageReader.GetOutput()
40 
41 
42 
43 
44 #
45 # Instantiate the classes for the registration framework
46 #
47 registration = itkImageRegistrationMethodF2F2_New()
48 imageMetric = itkMattesMutualInformationImageToImageMetricF2F2_New()
49 transform = itkTranslationTransform2_New()
50 optimizer = itkRegularStepGradientDescentOptimizer_New()
51 interpolator = itkLinearInterpolateImageFunctionF2D_New()
52 
53 
54 imageMetric.SetNumberOfHistogramBins( 20 );
55 imageMetric.SetNumberOfSpatialSamples( 10000 );
56 
57 registration.SetOptimizer( optimizer.GetPointer() )
58 registration.SetTransform( transform.GetPointer() )
59 registration.SetInterpolator( interpolator.GetPointer() )
60 registration.SetMetric( imageMetric.GetPointer() )
61 registration.SetFixedImage( fixedImage )
62 registration.SetMovingImage( movingImage )
63 
64 registration.SetFixedImageRegion( fixedImage.GetBufferedRegion() )
65 
66 transform.SetIdentity()
67 initialParameters = transform.GetParameters()
68 
69 registration.SetInitialTransformParameters( initialParameters )
70 
71 #
72 # Iteration Observer
73 #
74 def iterationUpdate():
75  currentParameter = transform.GetParameters()
76  print "M: %f P: %f %f " % ( optimizer.GetValue(),
77  currentParameter.GetElement(0),
78  currentParameter.GetElement(1) )
79 
80 iterationCommand = itkPyCommand_New()
81 iterationCommand.SetCommandCallable( iterationUpdate )
82 optimizer.AddObserver( itkIterationEvent(), iterationCommand.GetPointer() )
83 
84 
85 
86 #
87 # Define optimizer parameters
88 #
89 optimizer.SetMaximumStepLength( 4.00 )
90 optimizer.SetMinimumStepLength( 0.01 )
91 optimizer.SetNumberOfIterations( 200 )
92 
93 
94 print "Starting registration"
95 
96 #
97 # Start the registration process
98 #
99 
100 registration.Update()
101 
102 
103 #
104 # Get the final parameters of the transformation
105 #
106 finalParameters = registration.GetLastTransformParameters()
107 
108 print "Final Registration Parameters "
109 print "Translation X = %f" % (finalParameters.GetElement(0),)
110 print "Translation Y = %f" % (finalParameters.GetElement(1),)
111 
112 
113 
114 
115 #
116 # Now, we use the final transform for resampling the
117 # moving image.
118 #
119 resampler = itkResampleImageFilterF2F2_New()
120 resampler.SetTransform( transform.GetPointer() )
121 resampler.SetInput( movingImage )
122 
123 region = fixedImage.GetLargestPossibleRegion()
124 
125 resampler.SetSize( region.GetSize() )
126 
127 resampler.SetOutputSpacing( fixedImage.GetSpacing() )
128 resampler.SetOutputDirection( fixedImage.GetDirection() )
129 resampler.SetOutputOrigin( fixedImage.GetOrigin() )
130 resampler.SetDefaultPixelValue( 100 )
131 
132 outputCast = itkRescaleIntensityImageFilterF2US2_New()
133 outputCast.SetOutputMinimum( 0 )
134 outputCast.SetOutputMaximum( 65535 )
135 outputCast.SetInput(resampler.GetOutput())
136 
137 #
138 # Write the resampled image
139 #
140 writer = itkImageFileWriterUS2_New()
141 
142 writer.SetFileName( argv[3] )
143 writer.SetInput( outputCast.GetOutput() )
144 writer.Update()