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