Deform a Volume With a Thin Plate Spline#

Synopsis#

This example deforms a 3D volume with the thin plate spline.

Results#

Code#

Python#

#!/usr/bin/env python

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)

C++#

#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;
}

Classes demonstrated#

template<typename TParametersValueType, unsigned int NDimensions = 3> class ThinPlateSplineKernelTransform : public itk::KernelTransform<TParametersValueType, NDimensions>: This class defines the thin plate spline (TPS) transformation. It is implemented in as straightforward a manner as possible from the IEEE TMI paper by Davis, Khotanzad, Flamig, and Harms, Vol. 16 No. 3 June 1997

See itk::ThinPlateSplineKernelTransform for additional documentation.