SimpleITK  
DemonsRegistration2/DemonsRegistration2.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 fast symmetric forces Demons image
21 registation. """
22
23import sys
24import os
25import SimpleITK as sitk
26
27
28def command_iteration(sitk_filter):
29 """ Callback invoked when the filter progresses. """
30 print(f"{sitk_filter.GetElapsedIterations():3} = {sitk_filter.GetMetric():10.5f}")
31
32
33if len(sys.argv) < 4:
34 print(
35 "Usage:",
36 sys.argv[0],
37 "<fixedImageFilter> <movingImageFile>",
38 "[initialTransformFile] <outputTransformFile>",
39 )
40 sys.exit(1)
41
42fixed = sitk.ReadImage(sys.argv[1])
43
44moving = sitk.ReadImage(sys.argv[2])
45
47if fixed.GetPixelID() in (sitk.sitkUInt8, sitk.sitkInt8):
48 matcher.SetNumberOfHistogramLevels(128)
49else:
50 matcher.SetNumberOfHistogramLevels(1024)
51matcher.SetNumberOfMatchPoints(7)
52matcher.ThresholdAtMeanIntensityOn()
53moving = matcher.Execute(moving, fixed)
54
55# The fast symmetric forces Demons Registration Filter
56# Note there is a whole family of Demons Registration algorithms included in
57# SimpleITK
59demons.SetNumberOfIterations(200)
60# Standard deviation for Gaussian smoothing of displacement field
61demons.SetStandardDeviations(1.0)
62
63demons.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(demons))
64
65if len(sys.argv) > 4:
66 initialTransform = sitk.ReadTransform(sys.argv[3])
67 sys.argv[-1] = sys.argv.pop()
68
69 toDisplacementFilter = sitk.TransformToDisplacementFieldFilter()
70 toDisplacementFilter.SetReferenceImage(fixed)
71
72 displacementField = toDisplacementFilter.Execute(initialTransform)
73
74 displacementField = demons.Execute(fixed, moving, displacementField)
75
76else:
77
78 displacementField = demons.Execute(fixed, moving)
79
80print("-------")
81print(f"Number Of Iterations: {demons.GetElapsedIterations()}")
82print(f" RMS: {demons.GetRMSChange()}")
83
84outTx = sitk.DisplacementFieldTransform(displacementField)
85
86sitk.WriteTransform(outTx, sys.argv[3])
87
88if "SITK_NOSHOW" not in os.environ:
89 resampler = sitk.ResampleImageFilter()
90 resampler.SetReferenceImage(fixed)
91 resampler.SetInterpolator(sitk.sitkLinear)
92 resampler.SetDefaultPixelValue(100)
93 resampler.SetTransform(outTx)
94
95 out = resampler.Execute(moving)
96 simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8)
97 simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)
98 cimg = sitk.Compose(simg1, simg2, simg1 // 2.0 + simg2 // 2.0)
99 sitk.Show(cimg, "DeformableRegistration1 Composition")
A dense deformable transform over a bounded spatial domain for 2D or 3D coordinates space.
Deformably register two images using a symmetric forces demons algorithm.
Normalize the grayscale values for a source image by matching the shape of the source image histogram...
Resample an image via a coordinate transform.
Generate a displacement field from 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())
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 .
SITKCommon_EXPORT Transform ReadTransform(const PathType &filename)
Image Cast(const Image &image, PixelIDValueEnum pixelID)