SimpleITK  
ImageRegistrationMethodDisplacement1/ImageRegistrationMethodDisplacement1.py
1 #!/usr/bin/env python
2 # =========================================================================
3 #
4 # Copyright NumFOCUS
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 import SimpleITK as sitk
21 import sys
22 import os
23 
24 
25 def command_iteration(method):
26  if method.GetOptimizerIteration() == 0:
27  print(f"\tLevel: {method.GetCurrentLevel()}")
28  print(f"\tScales: {method.GetOptimizerScales()}")
29  print(f"#{method.GetOptimizerIteration()}")
30  print(f"\tMetric Value: {method.GetMetricValue():10.5f}")
31  print(f"\tLearningRate: {method.GetOptimizerLearningRate():10.5f}")
32  if method.GetOptimizerConvergenceValue() != sys.float_info.max:
33  print(
34  "\tConvergence Value: "
35  + f"{method.GetOptimizerConvergenceValue():.5e}"
36  )
37 
38 
39 def command_multiresolution_iteration(method):
40  print(f"\tStop Condition: {method.GetOptimizerStopConditionDescription()}")
41  print("============= Resolution Change =============")
42 
43 
44 if len(sys.argv) < 4:
45  print(
46  "Usage:",
47  sys.argv[0],
48  "<fixedImageFilter> <movingImageFile>",
49  "<outputTransformFile>",
50  )
51  sys.exit(1)
52 
53 fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)
54 
55 moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)
56 
58  fixed, moving, sitk.AffineTransform(fixed.GetDimension())
59 )
60 
62 
63 R.SetShrinkFactorsPerLevel([3, 2, 1])
64 R.SetSmoothingSigmasPerLevel([2, 1, 1])
65 
66 R.SetMetricAsJointHistogramMutualInformation(20)
67 R.MetricUseFixedImageGradientFilterOff()
68 
69 R.SetOptimizerAsGradientDescent(
70  learningRate=1.0,
71  numberOfIterations=100,
72  estimateLearningRate=R.EachIteration,
73 )
74 R.SetOptimizerScalesFromPhysicalShift()
75 
76 R.SetInitialTransform(initialTx)
77 
78 R.SetInterpolator(sitk.sitkLinear)
79 
80 R.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(R))
81 R.AddCommand(
82  sitk.sitkMultiResolutionIterationEvent,
83  lambda: command_multiresolution_iteration(R),
84 )
85 
86 outTx1 = R.Execute(fixed, moving)
87 
88 print("-------")
89 print(outTx1)
90 print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
91 print(f" Iteration: {R.GetOptimizerIteration()}")
92 print(f" Metric value: {R.GetMetricValue()}")
93 
94 displacementField = sitk.Image(fixed.GetSize(), sitk.sitkVectorFloat64)
95 displacementField.CopyInformation(fixed)
96 displacementTx = sitk.DisplacementFieldTransform(displacementField)
97 del displacementField
98 displacementTx.SetSmoothingGaussianOnUpdate(
99  varianceForUpdateField=0.0, varianceForTotalField=1.5
100 )
101 
102 R.SetMovingInitialTransform(outTx1)
103 R.SetInitialTransform(displacementTx, inPlace=True)
104 
105 R.SetMetricAsANTSNeighborhoodCorrelation(4)
106 R.MetricUseFixedImageGradientFilterOff()
107 
108 R.SetShrinkFactorsPerLevel([3, 2, 1])
109 R.SetSmoothingSigmasPerLevel([2, 1, 1])
110 
111 R.SetOptimizerScalesFromPhysicalShift()
112 R.SetOptimizerAsGradientDescent(
113  learningRate=1,
114  numberOfIterations=300,
115  estimateLearningRate=R.EachIteration,
116 )
117 
118 R.Execute(fixed, moving)
119 
120 print("-------")
121 print(displacementTx)
122 print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
123 print(f" Iteration: {R.GetOptimizerIteration()}")
124 print(f" Metric value: {R.GetMetricValue()}")
125 
126 compositeTx = sitk.CompositeTransform([outTx1, displacementTx])
127 sitk.WriteTransform(compositeTx, sys.argv[3])
128 
129 if "SITK_NOSHOW" not in os.environ:
130  sitk.Show(displacementTx.GetDisplacementField(), "Displacement Field")
131 
132  resampler = sitk.ResampleImageFilter()
133  resampler.SetReferenceImage(fixed)
134  resampler.SetInterpolator(sitk.sitkLinear)
135  resampler.SetDefaultPixelValue(100)
136  resampler.SetTransform(compositeTx)
137 
138  out = resampler.Execute(moving)
139  simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8)
140  simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)
141  cimg = sitk.Compose(simg1, simg2, simg1 // 2.0 + simg2 // 2.0)
142  sitk.Show(cimg, "ImageRegistration1 Composition")
itk::simple::Image
The Image class for SimpleITK.
Definition: sitkImage.h:76
itk::simple::WriteTransform
SITKCommon_EXPORT void WriteTransform(const Transform &transform, const std::string &filename)
itk::simple::DisplacementFieldTransform
A dense deformable transform over a bounded spatial domain for 2D or 3D coordinates space.
Definition: sitkDisplacementFieldTransform.h:34
itk::simple::RescaleIntensity
Image RescaleIntensity(const Image &image1, double outputMinimum=0, double outputMaximum=255)
Applies a linear transformation to the intensity levels of the input Image .
itk::simple::Show
void SITKIO_EXPORT Show(const Image &image, const std::string &title="", const bool debugOn=ProcessObject::GetGlobalDefaultDebug())
itk::simple::ReadImage
SITKIO_EXPORT Image ReadImage(const std::vector< std::string > &fileNames, PixelIDValueEnum outputPixelType=sitkUnknown, const std::string &imageIO="")
ReadImage is a procedural interface to the ImageSeriesReader class which is convenient for most image...
itk::simple::Cast
Image Cast(const Image &image, PixelIDValueEnum pixelID)
itk::simple::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: sitkResampleImageFilter.h:52
itk::simple::CompositeTransform
This class contains a stack of transforms and concatenates them by composition.
Definition: sitkCompositeTransform.h:58
itk::simple::Compose
Image Compose(const Image &image1, const Image &image2, const Image &image3, const Image &image4, const Image &image5)
ComposeImageFilter combine several scalar images into a multicomponent image.
itk::simple::ImageRegistrationMethod
An interface method to the modular ITKv4 registration framework.
Definition: sitkImageRegistrationMethod.h:89
itk::simple::AffineTransform
An affine transformation about a fixed center with translation for a 2D or 3D coordinate.
Definition: sitkAffineTransform.h:33
itk::simple::CenteredTransformInitializer
Transform CenteredTransformInitializer(const Image &fixedImage, const Image &movingImage, const Transform &transform, CenteredTransformInitializerFilter::OperationModeType operationMode=itk::simple::CenteredTransformInitializerFilter::MOMENTS)
CenteredTransformInitializer is a helper class intended to initialize the center of rotation and the ...