SimpleITK  1.0.1
ImageRegistrationMethodExhaustive/ImageRegistrationMethodExhaustive.py
1 #!/usr/bin/env python
2 #=========================================================================
3 #
4 # Copyright Insight Software Consortium
5 #
6 # Licensed under the Apache License, Version 2.0 (the "License");
7 # you may not use this file except in compliance with the License.
8 # You may obtain a copy of the License at
9 #
10 # http://www.apache.org/licenses/LICENSE-2.0.txt
11 #
12 # Unless required by applicable law or agreed to in writing, software
13 # distributed under the License is distributed on an "AS IS" BASIS,
14 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 # See the License for the specific language governing permissions and
16 # limitations under the License.
17 #
18 #=========================================================================
19 
20 """
21 This script demonstrates the use of the Exhaustive optimizer in the
22 ImageRegistrationMethod to estimate a good initial rotation position.
23 
24 Because gradient descent base optimization can get stuck in local
25 minima, a good initial transform is critical for reasonable
26 results. Search a reasonable space on a grid with brute force may be a
27 reliable way to get a starting location for further optimization.
28 
29 The initial translation and center of rotation for the transform is
30 initialized based on the first principle moments of the intensities of
31 the image. Then in either 2D or 3D a Euler transform is used to
32 exhaustively search a grid of the rotation space at a certain step
33 size. The resulting transform is a reasonable guess where to start
34 further registration.
35 """
36 
37 from __future__ import print_function
38 from __future__ import division
39 
40 
41 import SimpleITK as sitk
42 import sys
43 import os
44 from math import pi
45 
46 
47 def command_iteration(method) :
48  if (method.GetOptimizerIteration()==0):
49  print("Scales: ", method.GetOptimizerScales())
50  print("{0:3} = {1:7.5f} : {2}".format(method.GetOptimizerIteration(),
51  method.GetMetricValue(),
52  method.GetOptimizerPosition()))
53 
54 
55 
56 if len ( sys.argv ) < 4:
57  print( "Usage: {0} <fixedImageFilter> <movingImageFile> <outputTransformFile>".format(sys.argv[0]))
58  sys.exit ( 1 )
59 
60 fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)
61 
62 moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)
63 
65 
66 R.SetMetricAsMattesMutualInformation(numberOfHistogramBins = 50)
67 
68 sample_per_axis=12
69 if fixed.GetDimension() == 2:
71  # Set the number of samples (radius) in each dimension, with a
72  # default step size of 1.0
73  R.SetOptimizerAsExhaustive([sample_per_axis//2,0,0])
74  # Utilize the scale to set the step size for each dimension
75  R.SetOptimizerScales([2.0*pi/sample_per_axis, 1.0,1.0])
76 elif fixed.GetDimension() == 3:
78  R.SetOptimizerAsExhaustive([sample_per_axis//2,sample_per_axis//2,sample_per_axis//4,0,0,0])
79  R.SetOptimizerScales([2.0*pi/sample_per_axis,2.0*pi/sample_per_axis,2.0*pi/sample_per_axis,1.0,1.0,1.0])
80 
81 # Initialize the transform with a translation and the center of
82 # rotation from the moments of intensity.
83 tx = sitk.CenteredTransformInitializer(fixed, moving, tx)
84 
85 R.SetInitialTransform(tx)
86 
87 R.SetInterpolator(sitk.sitkLinear)
88 
89 R.AddCommand( sitk.sitkIterationEvent, lambda: command_iteration(R) )
90 
91 outTx = R.Execute(fixed, moving)
92 
93 print("-------")
94 print(outTx)
95 print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
96 print(" Iteration: {0}".format(R.GetOptimizerIteration()))
97 print(" Metric value: {0}".format(R.GetMetricValue()))
98 
99 
100 sitk.WriteTransform(outTx, sys.argv[3])
101 
102 if ( not "SITK_NOSHOW" in os.environ ):
103 
104  resampler = sitk.ResampleImageFilter()
105  resampler.SetReferenceImage(fixed);
106  resampler.SetInterpolator(sitk.sitkLinear)
107  resampler.SetDefaultPixelValue(1)
108  resampler.SetTransform(outTx)
109 
110  out = resampler.Execute(moving)
111 
112  simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8)
113  simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)
114  cimg = sitk.Compose(simg1, simg2, simg1/2.+simg2/2.)
115  sitk.Show( cimg, "ImageRegistrationExhaustive Composition" )