Project, draw and reconstruct geometric phantom

GeometricPhantom demonstrates how to create 2D projections of a Forbild phantom, draw a reference 3D image and reconstruct from the projections with FDK. You can find how to create your own phantom file in the phantom definition page.

This example uses the Thorax Forbild phantom, available in the dedicated GitHub repository for Forbild phantom files compatible with RTK.

#include "rtkDrawGeometricPhantomImageFilter.h"
#include "rtkConstantImageSource.h"
#include "rtkProjectGeometricPhantomImageFilter.h"
#include "rtkFDKConeBeamReconstructionFilter.h"
#include "rtkThreeDCircularProjectionGeometryXMLFileWriter.h"

int
main()
{
  constexpr unsigned int Dimension = 3;
  using OutputPixelType = float;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;

  constexpr unsigned int numberOfProjections = 180;
  constexpr double       angularArc = 360.;
  constexpr unsigned int sid = 600;  // source to isocenter distance
  constexpr unsigned int sdd = 1200; // source to detector distance
  constexpr double       scale = 2.;
  const std::string      configFileName = "Thorax";

  itk::Matrix<double, Dimension, Dimension> rotation;
  rotation[0][0] = 1.;
  rotation[1][2] = 1.;
  rotation[2][1] = 1.;

  // Set up the geometry for the projections
  auto geometry = rtk::ThreeDCircularProjectionGeometry::New();
  for (unsigned int x = 0; x < numberOfProjections; x++)
  {
    geometry->AddProjection(sid, sdd, x * angularArc / numberOfProjections);
  }
  rtk::WriteGeometry(geometry, "geometry.xml");

  // Create a constant image source for the tomography
  auto tomographySource = rtk::ConstantImageSource<OutputImageType>::New();
  tomographySource->SetOrigin(itk::MakePoint(-63.5, -63.5, -63.5));
  tomographySource->SetSize(itk::MakeSize(128, 128, 128));

  // Create a constant image source for the projections
  auto projectionsSource = rtk::ConstantImageSource<OutputImageType>::New();
  projectionsSource->SetOrigin(itk::MakePoint(-127., -127., -127.));
  projectionsSource->SetSpacing(itk::MakeVector(2., 2., 2.));
  projectionsSource->SetSize(itk::MakeSize(128, 128, numberOfProjections));

  // Project the geometric phantom image
  auto pgp = rtk::ProjectGeometricPhantomImageFilter<OutputImageType, OutputImageType>::New();
  pgp->SetInput(projectionsSource->GetOutput());
  pgp->SetGeometry(geometry);
  pgp->SetPhantomScale(scale);
  pgp->SetRotationMatrix(rotation);
  pgp->SetConfigFile(configFileName);
  pgp->SetIsForbildConfigFile(true);
  TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(pgp->GetOutput(), "projections.mha"));

  // Draw the geometric phantom image
  auto dgp = rtk::DrawGeometricPhantomImageFilter<OutputImageType, OutputImageType>::New();
  dgp->SetInput(tomographySource->GetOutput());
  dgp->SetPhantomScale(scale);
  dgp->SetRotationMatrix(rotation);
  dgp->SetConfigFile(configFileName);
  dgp->SetIsForbildConfigFile(true);
  TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(dgp->GetOutput(), "ref.mha"));

  // Perform FDK reconstruction filtering
  auto feldkamp = rtk::FDKConeBeamReconstructionFilter<OutputImageType>::New();
  feldkamp->SetInput(0, tomographySource->GetOutput());
  feldkamp->SetInput(1, pgp->GetOutput());
  feldkamp->SetGeometry(geometry);
  TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(feldkamp->GetOutput(), "fdk.mha"));

  return EXIT_SUCCESS;
}

import itk
from itk import RTK as rtk

Dimension = 3
OutputPixelType = itk.F
OutputImageType = itk.Image[OutputPixelType, Dimension]

numberOfProjections = 180
angularArc = 360.0
sid = 600
sdd = 1200
scale = 2.0
config_file_name = "Thorax"
rotation = itk.matrix_from_array([[1.0, 0.0, 0.0], [0.0, 0.0, -1.0], [0.0, 1.0, 0.0]])

# Set up the geometry for the projections
geometry = rtk.ThreeDCircularProjectionGeometry.New()
for x in range(0, numberOfProjections):
    angle = x * angularArc / numberOfProjections
    geometry.AddProjection(sid, sdd, angle)
rtk.write_geometry(geometry, "geometry.xml")

# Project the geometric phantom image
projections_source = rtk.constant_image_source(
    ttype=[OutputImageType],
    origin=[-127.0, -127.0, 0.0],
    size=[128, 128, numberOfProjections],
    spacing=[2.0] * 3,
)
pgp = rtk.project_geometric_phantom_image_filter(
    projections_source,
    geometry=geometry,
    phantom_scale=scale,
    rotation_matrix=rotation,
    config_file=config_file_name,
    is_forbild_config_file=True,
)
itk.imwrite(pgp, "projections.mha")

# Draw the geometric phantom image
ref_source = rtk.constant_image_source(
    ttype=[OutputImageType], origin=[-63.5] * 3, size=[128] * 3
)
dgp = rtk.draw_geometric_phantom_image_filter(
    ref_source,
    phantom_scale=scale,
    rotation_matrix=rotation,
    config_file=config_file_name,
    is_forbild_config_file=True,
)
itk.imwrite(dgp, "ref.mha")

# Perform FDK reconstruction filtering
fdk_source = rtk.constant_image_source(
    ttype=[OutputImageType],
    origin=[-63.5] * 3,
    size=[128] * 3,
    spacing=[1.0] * 3,
    constant=0.0,
)
feldkamp = rtk.fdk_cone_beam_reconstruction_filter(fdk_source, pgp, geometry=geometry)
itk.imwrite(feldkamp, "fdk.mha")

The results displayed with rtkshowgeometry are:

visu