From f50c2e9fcfa164d3eb0b8ade90624401d09c6605 Mon Sep 17 00:00:00 2001
From: eDgar <edgar@openmail.cc>
Date: Wed, 2 Feb 2022 11:16:25 -0600
Subject: [RFC PATCH] Simple 3D registration example with Python

---
 .../hello_3d_registration.py                  | 206 ++++++++++++++++++
 1 file changed, 206 insertions(+)
 create mode 100644 Examples/RegistrationITKv4/hello_3d_registration.py

diff --git a/Examples/RegistrationITKv4/hello_3d_registration.py b/Examples/RegistrationITKv4/hello_3d_registration.py
new file mode 100644
index 0000000000..bb62e9e1c7
--- /dev/null
+++ b/Examples/RegistrationITKv4/hello_3d_registration.py
@@ -0,0 +1,206 @@
+#!/usr/bin/env python
+# This code is based on
+# https://itk.org/ITKExamples/src/Registration/Common/Perform2DTranslationRegistrationWithMeanSquares/Documentation.html
+# which holds the Apache 2.0 license (see below).
+#
+# Copyright 2022 edgar
+# License:
+#
+#    This program is free software: you can redistribute it
+#    and/or modify it under the terms of the GNU General
+#    Public License as published by the Free Software
+#    Foundation, version 3 of the License.
+#
+#    This program is distributed in the hope that it will be
+#    useful, but WITHOUT ANY WARRANTY; without even the
+#    implied warranty of MERCHANTABILITY or FITNESS FOR A
+#    PARTICULAR PURPOSE. See the GNU General Public License
+#    for more details.
+#
+#    You should have received a copy of the GNU General
+#    Public License along with this program. If not, see
+#    <https://www.gnu.org/licenses/>.
+#
+# Copyright Insight Software Consortium
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import sys
+import itk
+import numpy as np
+from distutils.version import StrictVersion as VS
+from skimage.morphology import disk
+from scipy.ndimage import shift
+
+
+if VS(itk.Version.GetITKVersion()) < VS("4.9.0"):
+    print("ITK 4.9.0 is required.")
+    sys.exit(1)
+
+# Set outputs
+output_image = "itk_hello_registration.out.vtk"
+first_frame_vtk = "itk_hello_registration.in01.vtk"
+last_frame_vtk = "itk_hello_registration.in04.vtk"
+diff_image_after = "itk_hello_registration.diff_after.vtk"
+diff_image_before = "itk_hello_registration.diff_before.vtk"
+
+# Draw a series with a 3D and moving disk
+series_length = 4
+side_size = 10
+obj2d = disk(side_size / 2)
+obj3d = np.asarray(
+    np.pad([disk(side_size / 2, dtype=np.uint8)
+            for idx in range(5)],
+           ((0, series_length + 1),
+            (0, series_length + 1),
+            (0, 0))))
+holder = np.zeros((series_length, *obj3d.shape),
+                  dtype=np.uint8, order="F")
+# Move the last frame
+holder[0] = obj3d
+holder[-1] = shift(obj3d,
+                   (series_length - 1, 0, 0))
+
+# # Uncomment to see the data
+# import matplotlib.pyplot as plt
+# from mpl_toolkits.mplot3d import Axes3D
+# ax = plt.figure().add_subplot(1, 1, 1,
+#                               projection=Axes3D.name)
+# _ = ax.voxels(holder[0])
+# _ = ax.voxels(holder[-1], alpha=0.5)
+# plt.show()
+
+# Set the right type
+pixel_type = itk.ctype("float")
+dimension = 3
+image_type = itk.Image[pixel_type, dimension]
+
+# Load first and last frames as ITK data (images);
+# swap axes ITK: i, j, k, Numpy: k, i, j
+# (https://discourse.itk.org/t/
+#   importing-image-from-array-and-axis-reorder/1192/)
+fixed_image_arr = holder[0].transpose((2, 1, 0))
+moving_image_arr = holder[-1].transpose((2, 1, 0))
+fixed_image = itk.image_from_array(
+    fixed_image_arr.astype(np.float32))
+moving_image = itk.image_from_array(
+    moving_image_arr.astype(np.float32))
+
+dimension = fixed_image.GetImageDimension()
+# fixed_image_type = itk.Image[pixel_type, dimension]
+# moving_image_type = itk.Image[pixel_type, dimension]
+fixed_image_type = type(fixed_image)
+moving_image_type = type(moving_image)
+
+transform_type = itk.TranslationTransform[itk.D, dimension]
+initial_transform = transform_type.New()
+
+optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
+    LearningRate=4,
+    MinimumStepLength=0.001,
+    RelaxationFactor=0.5,
+    NumberOfIterations=200)
+
+metric = itk.MeanSquaresImageToImageMetricv4[
+    fixed_image_type, moving_image_type].New()
+
+registration_method = itk.ImageRegistrationMethodv4
+registration = registration_method[
+    fixed_image_type, moving_image_type].New(
+        FixedImage=fixed_image,
+        MovingImage=moving_image,
+        Metric=metric,
+        Optimizer=optimizer,
+        InitialTransform=initial_transform)
+
+moving_initial_transform = transform_type.New()
+initial_parameters = moving_initial_transform.GetParameters()
+initial_parameters[0] = 0
+initial_parameters[1] = 0
+moving_initial_transform.SetParameters(initial_parameters)
+registration.SetMovingInitialTransform(moving_initial_transform)
+
+identity_transform = transform_type.New()
+identity_transform.SetIdentity()
+registration.SetFixedInitialTransform(identity_transform)
+
+registration.SetNumberOfLevels(1)
+registration.SetSmoothingSigmasPerLevel([0])
+registration.SetShrinkFactorsPerLevel([1])
+
+registration.Update()
+
+transform = registration.GetTransform()
+final_parameters = transform.GetParameters()
+translation_along_X = final_parameters.GetElement(0)
+translation_along_Y = final_parameters.GetElement(1)
+translation_along_Z = final_parameters.GetElement(2)
+
+number_of_iterations = optimizer.GetCurrentIteration()
+
+best_value = optimizer.GetValue()
+
+print("Result = ")
+print(" Translation X = " + str(translation_along_X))
+print(" Translation Y = " + str(translation_along_Y))
+print(" Translation Z = " + str(translation_along_Z))
+print(" Iterations    = " + str(number_of_iterations))
+print(" Metric value  = " + str(best_value))
+
+composite_transform_type = itk.CompositeTransform[
+    itk.D, dimension]
+output_composite_transform = composite_transform_type.New()
+output_composite_transform.AddTransform(
+    moving_initial_transform)
+output_composite_transform.AddTransform(
+    registration.GetModifiableTransform())
+
+resampler = itk.ResampleImageFilter.New(
+    Input=moving_image,
+    Transform=output_composite_transform,
+    UseReferenceImage=True,
+    ReferenceImage=fixed_image)
+resampler.SetDefaultPixelValue(100)
+
+output_pixel_type = itk.ctype('unsigned char')
+output_image_type = itk.Image[output_pixel_type, dimension]
+
+caster = itk.CastImageFilter[
+    fixed_image_type, output_image_type].New(Input=resampler)
+
+writer = itk.ImageFileWriter.New(
+    Input=caster, FileName=output_image)
+writer.SetFileName(output_image)
+writer.Update()
+
+difference = itk.SubtractImageFilter.New(
+    Input1=fixed_image, Input2=resampler)
+
+intensity_filter = itk.RescaleIntensityImageFilter
+intensity_rescaler = intensity_filter[
+    fixed_image_type, output_image_type].New(
+        Input=difference,
+        OutputMinimum=itk.NumericTraits[output_pixel_type].min(),
+        OutputMaximum=itk.NumericTraits[output_pixel_type].max())
+
+resampler.SetDefaultPixelValue(1)
+writer.SetInput(intensity_rescaler.GetOutput())
+writer.SetFileName(diff_image_after)
+writer.Update()
+
+resampler.SetTransform(identity_transform)
+writer.SetFileName(diff_image_before)
+writer.Update()
+
+itk.imwrite(fixed_image, first_frame_vtk)
+itk.imwrite(moving_image, last_frame_vtk)
-- 
2.34.1

