SimpleITK  
ImageRegistrationMethodBSpline3/ImageRegistrationMethodBSpline3.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 BSplineTransform and the JointHistogramMutualInformation metric. """
22
23import sys
24import os
25import SimpleITK as sitk
26
27
28def command_iteration(method, bspline_transform):
29 """ Callback invoked each iteration """
30 if method.GetOptimizerIteration() == 0:
31 # The BSpline is resized before the first optimizer
32 # iteration is completed per level. Print the transform object
33 # to show the adapted BSpline transform.
34 print(bspline_transform)
35
36 print(f"{method.GetOptimizerIteration():3} " + f"= {method.GetMetricValue():10.5f}")
37
38
39def command_multi_iteration(method):
40 """ Callback invoked before starting a multi-resolution level.
41 The sitkMultiResolutionIterationEvent occurs before the
42 resolution of the transform. This event is used here to print
43 the status of the optimizer from the previous registration level.
44 """
45 if method.GetCurrentLevel() > 0:
46 print(
47 "Optimizer stop condition: " + f"{R.GetOptimizerStopConditionDescription()}"
48 )
49 print(f" Iteration: {R.GetOptimizerIteration()}")
50 print(f" Metric value: {R.GetMetricValue()}")
51
52 print("--------- Resolution Changing ---------")
53
54
55if len(sys.argv) < 4:
56 print(
57 "Usage:",
58 sys.argv[0],
59 "<fixedImageFilter> <movingImageFile>",
60 "<outputTransformFile>",
61 )
62 sys.exit(1)
63
64fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)
65
66moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)
67
68transformDomainMeshSize = [2] * fixed.GetDimension()
69tx = sitk.BSplineTransformInitializer(fixed, transformDomainMeshSize)
70
71print(f"Initial Number of Parameters: {tx.GetNumberOfParameters()}")
72
74R.SetMetricAsJointHistogramMutualInformation()
75
76R.SetOptimizerAsGradientDescentLineSearch(
77 5.0, 100, convergenceMinimumValue=1e-4, convergenceWindowSize=5
78)
79
80R.SetInterpolator(sitk.sitkLinear)
81
82R.SetInitialTransformAsBSpline(tx, inPlace=True, scaleFactors=[1, 2, 5])
83R.SetShrinkFactorsPerLevel([4, 2, 1])
84R.SetSmoothingSigmasPerLevel([4, 2, 1])
85
86R.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(R, tx))
87R.AddCommand(sitk.sitkMultiResolutionIterationEvent, lambda: command_multi_iteration(R))
88
89outTx = R.Execute(fixed, moving)
90
91print("-------")
92print(tx)
93print(outTx)
94print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
95print(f" Iteration: {R.GetOptimizerIteration()}")
96print(f" Metric value: {R.GetMetricValue()}")
97
98sitk.WriteTransform(outTx, sys.argv[3])
99
100if "SITK_NOSHOW" not in os.environ:
101 resampler = sitk.ResampleImageFilter()
102 resampler.SetReferenceImage(fixed)
103 resampler.SetInterpolator(sitk.sitkLinear)
104 resampler.SetDefaultPixelValue(100)
105 resampler.SetTransform(outTx)
106
107 out = resampler.Execute(moving)
108 simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8)
109 simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)
110 cimg = sitk.Compose(simg1, simg2, simg1 // 2.0 + simg2 // 2.0)
111 sitk.Show(cimg, "Image Registration Composition")
An interface method to the modular ITKv4 registration framework.
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())
BSplineTransform BSplineTransformInitializer(const Image &image1, const std::vector< uint32_t > &transformDomainMeshSize=std::vector< uint32_t >(3, 1u), unsigned int order=3u)
BSplineTransformInitializerFilter is a helper class intended to initialize the control point grid suc...
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)