ENH: Add DeformAVolumeWithAThinPlateSpline

Adapted from:

  https://itk.org/Doxygen/html/Examples_2RegistrationITKv4_2ThinPlateSplineWarp_8cxx-example.html

However, we use itk mesh IO to read the point sets, and use the
reference image interface to ResampleImageFilter. Both were added since
that example was created.

ITKRepositoryImport.todo

@@ -249,7 +249,6 @@
 * ../ITK/Examples/Registration/DeformationFieldJacobian.cxx
 * ../ITK/Examples/Registration/ImageRegistration5.tcl
 * ../ITK/Examples/Registration/ImageRegistration5.py
-* ../ITK/Examples/Registration/ThinPlateSplineWarp.cxx
 * ../ITK/Examples/Registration/ImageRegistration4.tcl
 * ../ITK/Examples/Registration/ImageRegistration16.cxx
 * ../ITK/Examples/Registration/ImageRegistration4.py

src/Core/Transform/CMakeLists.txt

@@ -39,3 +39,8 @@ add_example(GlobalRegistrationTwoImagesBSpline)
 compare_to_baseline(EXAMPLE_NAME GlobalRegistrationTwoImagesBSpline
   BASELINE_PREFIX GlobalRegistrationTwoImagesBSpline
   )
+
+add_example(DeformAVolumeWithAThinPlateSpline)
+compare_to_baseline(EXAMPLE_NAME DeformAVolumeWithAThinPlateSpline
+  BASELINE_PREFIX DeformedImageBaseline
+  )

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/CMakeLists.txt

@@ -0,0 +1,49 @@
+cmake_minimum_required(VERSION 3.16.3)
+
+project(DeformAVolumeWithAThinPlateSpline)
+
+find_package(ITK REQUIRED)
+include(${ITK_USE_FILE})
+
+add_executable(DeformAVolumeWithAThinPlateSpline Code.cxx)
+target_link_libraries(DeformAVolumeWithAThinPlateSpline ${ITK_LIBRARIES})
+
+install(TARGETS DeformAVolumeWithAThinPlateSpline
+  DESTINATION bin/ITKSphinxExamples/Core/Transform
+  COMPONENT Runtime
+)
+
+install(FILES Code.cxx CMakeLists.txt Code.py
+  DESTINATION share/ITKSphinxExamples/Code/Core/Transform/DeformAVolumeWithAThinPlateSpline
+  COMPONENT Code
+)
+
+enable_testing()
+
+set(source_landmarks ${CMAKE_CURRENT_BINARY_DIR}/SourceLandmarks.vtk)
+set(target_landmarks ${CMAKE_CURRENT_BINARY_DIR}/TargetLandmarks.vtk)
+set(input_image ${CMAKE_CURRENT_BINARY_DIR}/brainweb165a10f17.mha)
+set(deformed_image DeformedImage.mha)
+set(checkerboard_image CheckerBoardImage.mha)
+
+add_test(NAME DeformAVolumeWithAThinPlateSplineTest
+  COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/DeformAVolumeWithAThinPlateSpline
+  ${source_landmarks}
+  ${target_landmarks}
+  ${input_image}
+  ${deformed_image}
+  ${checkerboard_image}
+)
+
+if(ITK_WRAP_PYTHON)
+  find_package(PythonInterp REQUIRED)
+  string(REPLACE . "Python." deformed_image "${deformed_image}")
+  add_test(NAME DeformAVolumeWithAThinPlateSplineTestPython
+    COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Code.py
+      ${source_landmarks}
+      ${target_landmarks}
+      ${input_image}
+      ${deformed_image}
+      ${checkerboard_image}
+    )
+endif()

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/CheckerBoard.png.cid

@@ -0,0 +1 @@
+bafybeido67mh6ekbf4ksqnvohvhjeae3kps6lynmucjiz4goyjafujm6ki

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/Code.cxx

@@ -0,0 +1,119 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  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
+ *
+ *         https://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.
+ *
+ *=========================================================================*/
+
+#include "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include "itkImage.h"
+#include "itkResampleImageFilter.h"
+#include "itkThinPlateSplineKernelTransform.h"
+#include "itkPointSet.h"
+#include "itkMesh.h"
+#include "itkMeshFileReader.h"
+#include "itkCheckerBoardImageFilter.h"
+
+int
+main(int argc, char * argv[])
+{
+  if (argc != 6)
+  {
+    std::cerr << "Usage: " << argv[0];
+    std::cerr << " <SourceLandmarks>";
+    std::cerr << " <TargetLandmarks>";
+    std::cerr << " <InputImage>";
+    std::cerr << " <DeformedImage>";
+    std::cerr << " <CheckerBoardImage>";
+    std::cerr << std::endl;
+    return EXIT_FAILURE;
+  }
+  const char * sourceLandmarksFile = argv[1];
+  const char * targetLandmarksFile = argv[2];
+  const char * inputImageFile = argv[3];
+  const char * deformedImageFile = argv[4];
+  const char * checkerboardImageFile = argv[5];
+
+  constexpr unsigned int Dimension = 3;
+  using CoordinateRepType = double;
+
+  using TransformType = itk::ThinPlateSplineKernelTransform<CoordinateRepType, Dimension>;
+  using PointSetType = TransformType::PointSetType;
+
+  using PixelType = unsigned char;
+  using MeshType = itk::Mesh<PointSetType::PixelType, Dimension, PointSetType::MeshTraits>;
+  using ImageType = itk::Image<PixelType, Dimension>;
+
+  auto               sourceLandmarks = PointSetType::New();
+  auto               targetLandmarks = PointSetType::New();
+  ImageType::Pointer inputImage;
+  try
+  {
+    auto sourceLandmarksMesh = itk::ReadMesh<MeshType>(sourceLandmarksFile);
+    sourceLandmarks->SetPoints(sourceLandmarksMesh->GetPoints());
+
+    auto targetLandmarksMesh = itk::ReadMesh<MeshType>(targetLandmarksFile);
+    targetLandmarks->SetPoints(targetLandmarksMesh->GetPoints());
+
+    inputImage = itk::ReadImage<ImageType>(inputImageFile);
+  }
+  catch (itk::ExceptionObject & error)
+  {
+    std::cerr << "Error: " << error << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  using TransformType = itk::ThinPlateSplineKernelTransform<CoordinateRepType, Dimension>;
+
+  auto thinPlateSpline = TransformType::New();
+  thinPlateSpline->SetSourceLandmarks(sourceLandmarks);
+  thinPlateSpline->SetTargetLandmarks(targetLandmarks);
+  thinPlateSpline->ComputeWMatrix();
+
+  using ResamplerType = itk::ResampleImageFilter<ImageType, ImageType>;
+  auto resampler = ResamplerType::New();
+
+  resampler->SetInput(inputImage);
+
+  using InterpolatorType = itk::LinearInterpolateImageFunction<ImageType, CoordinateRepType>;
+  auto interpolator = InterpolatorType::New();
+  resampler->SetInterpolator(interpolator);
+
+  resampler->SetUseReferenceImage(true);
+  resampler->SetReferenceImage(inputImage);
+
+  resampler->SetTransform(thinPlateSpline);
+
+  using CheckerBoardFilterType = itk::CheckerBoardImageFilter<ImageType>;
+  auto checkerboardFilter = CheckerBoardFilterType::New();
+  checkerboardFilter->SetInput1(inputImage);
+  checkerboardFilter->SetInput2(resampler->GetOutput());
+
+  try
+  {
+    resampler->Update();
+    checkerboardFilter->Update();
+
+    itk::WriteImage(resampler->GetOutput(), deformedImageFile);
+    itk::WriteImage(checkerboardFilter->GetOutput(), checkerboardImageFile);
+  }
+  catch (itk::ExceptionObject & error)
+  {
+    std::cerr << "Error: " << error << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  return EXIT_SUCCESS;
+}

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/Code.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+
+# Copyright NumFOCUS
+#
+# 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
+#
+#        https://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 argparse
+import itk
+import numpy as np
+
+parser = argparse.ArgumentParser(
+    description="Deform a Volume With a Thin Plate Spline."
+)
+parser.add_argument("source_landmarks")
+parser.add_argument("target_landmarks")
+parser.add_argument("input_image")
+parser.add_argument("deformed_image")
+parser.add_argument("checker_board_image")
+args = parser.parse_args()
+
+Dimension = 3
+thin_plate_spline = itk.ThinPlateSplineKernelTransform[itk.D, Dimension].New()
+
+source_landmarks_mesh = itk.meshread(args.source_landmarks)
+# Cast points from float32 to float64
+points = itk.array_from_vector_container(source_landmarks_mesh.GetPoints())
+points = points.astype(np.float64)
+source_landmarks = thin_plate_spline.GetSourceLandmarks()
+source_landmarks.SetPoints(itk.vector_container_from_array(points.flatten()))
+
+target_landmarks_mesh = itk.meshread(args.target_landmarks)
+# Cast points from float32 to float64
+points = itk.array_from_vector_container(target_landmarks_mesh.GetPoints())
+points = points.astype(np.float64)
+target_landmarks = thin_plate_spline.GetTargetLandmarks()
+target_landmarks.SetPoints(itk.vector_container_from_array(points.flatten()))
+
+thin_plate_spline.ComputeWMatrix()
+
+input_image = itk.imread(args.input_image)
+
+deformed_image = itk.resample_image_filter(
+    input_image,
+    use_reference_image=True,
+    reference_image=input_image,
+    transform=thin_plate_spline,
+)
+itk.imwrite(deformed_image, args.deformed_image)
+
+checker_board = itk.checker_board_image_filter(input_image, deformed_image)
+itk.imwrite(checker_board, args.checker_board_image)

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/DeformedImage.png.cid

@@ -0,0 +1 @@
+bafybeif7peygrvroniic5yqavvxzuctdxqyuzz5tfmba5f63ew6ozpyl4u

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/DeformedImageBaseline.mha.cid

@@ -0,0 +1 @@
+bafybeih4kal3lobtzwsnsoc2vszskc46on3p5ojwa6tqriikca32c6ixum

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/Documentation.rst

@@ -0,0 +1,59 @@
+:name: DeformAVolumeWithAThinPlateSpline
+
+Deform a Volume With a Thin Plate Spline
+========================================
+
+.. index::
+   single: ThinPlateSplineKernelTransform
+
+Synopsis
+--------
+
+
+This example deforms a 3D volume with the thin plate spline.
+
+
+Results
+-------
+
+.. figure:: InputImage.png
+  :width: 640px
+  :alt: Input image
+
+  Input image
+
+.. figure:: DeformedImage.png
+  :width: 640px
+  :alt: Deformed image
+
+  Deformed image
+
+.. figure:: CheckerBoard.png
+  :width: 640px
+  :alt: CheckerBoard image
+
+  CheckerBoard image
+
+Code
+----
+
+Python
+......
+
+.. literalinclude:: Code.py
+   :language: python
+   :lines: 1,16-
+
+C++
+...
+
+.. literalinclude:: Code.cxx
+   :language: c++
+   :lines: 18-
+
+
+
+Classes demonstrated
+--------------------
+
+.. breathelink:: itk::ThinPlateSplineKernelTransform

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/InputImage.png.cid

@@ -0,0 +1 @@
+bafybeibzkxiedmtayg43pkunxie6hf6k3oqubepsh4luxjbx4onjeo5ep4

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/SourceLandmarks.vtk.cid

@@ -0,0 +1 @@
+bafkreicmsmsjebk2lvcwwqbz7u56tklm6qseu5q2dkn2rdidtegatbsc64

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/TargetLandmarks.vtk.cid

@@ -0,0 +1 @@
+bafkreiejdhhrghlfm2xwqzwdrwern3ssbguhp4t7g3sfen54gxwa3vlb3i

src/Core/Transform/DeformAVolumeWithAThinPlateSpline/brainweb165a10f17.mha.cid

@@ -0,0 +1 @@
+bafybeifjjxuuxkmfs3v2u763w76gklc5kbr6ecwxce4drafurjf44qys7a

src/Core/Transform/index.rst

@@ -14,3 +14,4 @@ Transform
   ScaleAnImage/Documentation.rst
   TranslateAVectorImage/Documentation.rst
   TranslateImage/Documentation.rst
+  DeformAVolumeWithAThinPlateSpline/Documentation.rst