ITK  4.9.0
Insight Segmentation and Registration Toolkit
Examples/RegistrationITKv3/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 import itk
20 from sys import argv
21 
22 
23 #
24 # Check input parameters
25 # INPUTS(fixedImage): {BrainProtonDensitySliceBorder20.png}
26 # INPUTS(movingImage): {BrainProtonDensitySliceShifted13x17y.png}
27 #
28 if len(argv) < 4:
29  print 'Missing Parameters'
30  print 'Usage: ImageRegistration4.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.TranslationTransform[itk.D, 2]
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.MattesMutualInformationImageToImageMetric[FixedImageType, MovingImageType].New()
65 transform = TransformType.New()
67 interpolator = itk.LinearInterpolateImageFunction[FixedImageType, itk.D].New()
68 
69 imageMetric.SetNumberOfHistogramBins( 20 );
70 imageMetric.SetNumberOfSpatialSamples( 10000 );
71 
72 registration.SetOptimizer( optimizer )
73 registration.SetTransform( transform )
74 registration.SetInterpolator( interpolator )
75 registration.SetMetric( imageMetric )
76 
77 registration.SetFixedImage( fixedImage )
78 registration.SetMovingImage( movingImage )
79 
80 registration.SetFixedImageRegion( fixedImage.GetBufferedRegion() )
81 
82 transform.SetIdentity()
83 initialParameters = transform.GetParameters()
84 
85 registration.SetInitialTransformParameters( initialParameters )
86 
87 
88 #
89 # Iteration Observer
90 #
91 def iterationUpdate():
92  currentParameter = transform.GetParameters()
93  print "M: %f P: %f %f " % ( optimizer.GetValue(),
94  currentParameter.GetElement(0),
95  currentParameter.GetElement(1) )
96 
97 iterationCommand = itk.PyCommand.New()
98 iterationCommand.SetCommandCallable( iterationUpdate )
99 optimizer.AddObserver( itk.IterationEvent(), iterationCommand )
100 
101 
102 #
103 # Define optimizer parameters
104 #
105 optimizer.SetMaximumStepLength( 4.00 )
106 optimizer.SetMinimumStepLength( 0.01 )
107 optimizer.SetNumberOfIterations( 200 )
108 
109 print "Starting registration"
110 
111 
112 #
113 # Start the registration process
114 #
115 registration.Update()
116 
117 
118 #
119 # Get the final parameters of the transformation
120 #
121 finalParameters = registration.GetLastTransformParameters()
122 
123 print "Final Registration Parameters "
124 print "Translation X = %f" % (finalParameters.GetElement(0),)
125 print "Translation Y = %f" % (finalParameters.GetElement(1),)
126 
127 
128 #
129 # Now, we use the final transform for resampling the
130 # moving image.
131 #
132 resampler = itk.ResampleImageFilter[MovingImageType, FixedImageType].New()
133 resampler.SetTransform( transform )
134 resampler.SetInput( movingImage )
135 
136 region = fixedImage.GetLargestPossibleRegion()
137 
138 resampler.SetSize( region.GetSize() )
139 
140 resampler.SetOutputSpacing( fixedImage.GetSpacing() )
141 resampler.SetOutputDirection( fixedImage.GetDirection() )
142 resampler.SetOutputOrigin( fixedImage.GetOrigin() )
143 resampler.SetDefaultPixelValue( 100 )
144 
145 outputCast = itk.RescaleIntensityImageFilter[FixedImageType, OutputImageType].New()
146 outputCast.SetInput(resampler.GetOutput())
147 
148 
149 #
150 # Write the resampled image
151 #
152 writer = itk.ImageFileWriter[OutputImageType].New()
153 
154 writer.SetFileName( argv[3] )
155 writer.SetInput( outputCast.GetOutput() )
156 writer.Update()