First reconstruction

RTK is a library, therefore it’s meant to be integrated into application. This tutorial shows how to create a simple FirstReconstruction project that links with RTK. The source code for this tutorial is located in RTK/examples/FirstReconstruction.

C++

CMakeLists

• First you need to create a CMakeLists.txt

cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR)

# This project is designed to be built outside the RTK source tree.
project(FirstReconstruction)

# Find ITK with RTK
find_package(ITK REQUIRED COMPONENTS RTK)
include(${ITK_USE_FILE})

# Executable(s)
add_executable(FirstReconstruction FirstReconstruction.cxx )
target_link_libraries(FirstReconstruction ${ITK_LIBRARIES})
if(${RTK_USE_CUDA})
  add_executable(FirstCudaReconstruction FirstCudaReconstruction.cxx )
  target_link_libraries(FirstCudaReconstruction ${ITK_LIBRARIES})
endif()

• Run CMake on the FirstReconstruction directory and create a FirstReconstruction-bin,

• Configure and build the project using your favorite compiler,

• Run FirstReconstruction image.mha geometry.xml. If everything runs correctly, you should see a few messages ending with Done! and two new files in the current directory, image.mha and geometry.xml. image.mha is the reconstruction of a sphere from 360 projections and geometry.xml is the geometry file in RTK’s geometry format.

Code

// RTK includes
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
#include <rtkRayEllipsoidIntersectionImageFilter.h>
#include <rtkFDKConeBeamReconstructionFilter.h>
#include <rtkFieldOfViewImageFilter.h>

// ITK includes
#include <itkImageFileWriter.h>

int
main(int argc, char ** argv)
{
  if (argc < 3)
  {
    std::cout << "Usage: FirstReconstruction <outputimage> <outputgeometry>" << std::endl;
    return EXIT_FAILURE;
  }

  // Defines the image type
  using ImageType = itk::Image<float, 3>;

  // Defines the RTK geometry object
  using GeometryType = rtk::ThreeDCircularProjectionGeometry;
  GeometryType::Pointer geometry = GeometryType::New();
  unsigned int          numberOfProjections = 360;
  double                firstAngle = 0;
  double                angularArc = 360;
  unsigned int          sid = 600;  // source to isocenter distance
  unsigned int          sdd = 1200; // source to detector distance
  for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
  {
    double angle = firstAngle + noProj * angularArc / numberOfProjections;
    geometry->AddProjection(sid, sdd, angle);
  }

  // Write the geometry to disk
  rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
  xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
  xmlWriter->SetFilename(argv[2]);
  xmlWriter->SetObject(geometry);
  xmlWriter->WriteFile();

  // Create a stack of empty projection images
  using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
  ConstantImageSourceType::Pointer     constantImageSource = ConstantImageSourceType::New();
  ConstantImageSourceType::PointType   origin;
  ConstantImageSourceType::SpacingType spacing;
  ConstantImageSourceType::SizeType    sizeOutput;

  origin[0] = -127;
  origin[1] = -127;
  origin[2] = 0.;

  sizeOutput[0] = 128;
  sizeOutput[1] = 128;
  sizeOutput[2] = numberOfProjections;

  spacing.Fill(2.);

  constantImageSource->SetOrigin(origin);
  constantImageSource->SetSpacing(spacing);
  constantImageSource->SetSize(sizeOutput);
  constantImageSource->SetConstant(0.);

  // Create projections of an ellipse
  using REIType = rtk::RayEllipsoidIntersectionImageFilter<ImageType, ImageType>;
  REIType::Pointer    rei = REIType::New();
  REIType::VectorType semiprincipalaxis, center;
  semiprincipalaxis.Fill(50.);
  center.Fill(0.);
  center[2] = 10.;
  rei->SetDensity(2.);
  rei->SetAngle(0.);
  rei->SetCenter(center);
  rei->SetAxis(semiprincipalaxis);
  rei->SetGeometry(geometry);
  rei->SetInput(constantImageSource->GetOutput());

  // Create reconstructed image
  ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();
  sizeOutput.Fill(128);
  origin.Fill(-63.5);
  spacing.Fill(1.);
  constantImageSource2->SetOrigin(origin);
  constantImageSource2->SetSpacing(spacing);
  constantImageSource2->SetSize(sizeOutput);
  constantImageSource2->SetConstant(0.);

  // FDK reconstruction
  std::cout << "Reconstructing..." << std::endl;
  using FDKCPUType = rtk::FDKConeBeamReconstructionFilter<ImageType>;
  FDKCPUType::Pointer feldkamp = FDKCPUType::New();
  feldkamp->SetInput(0, constantImageSource2->GetOutput());
  feldkamp->SetInput(1, rei->GetOutput());
  feldkamp->SetGeometry(geometry);
  feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
  feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);

  // Field-of-view masking
  using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
  FOVFilterType::Pointer fieldofview = FOVFilterType::New();
  fieldofview->SetInput(0, feldkamp->GetOutput());
  fieldofview->SetProjectionsStack(rei->GetOutput());
  fieldofview->SetGeometry(geometry);

  // Writer
  std::cout << "Writing output image..." << std::endl;
  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(argv[1]);
  writer->SetInput(fieldofview->GetOutput());
  writer->Update();

  std::cout << "Done!" << std::endl;
  return EXIT_SUCCESS;
}

C++ Cuda

CMakeLists

• First you need to create a CMakeLists.txt

cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR)

# This project is designed to be built outside the RTK source tree.
project(FirstReconstruction)

# Find ITK with RTK
find_package(ITK REQUIRED COMPONENTS RTK)
include(${ITK_USE_FILE})

# Executable(s)
add_executable(FirstReconstruction FirstReconstruction.cxx )
target_link_libraries(FirstReconstruction ${ITK_LIBRARIES})
if(${RTK_USE_CUDA})
  add_executable(FirstCudaReconstruction FirstCudaReconstruction.cxx )
  target_link_libraries(FirstCudaReconstruction ${ITK_LIBRARIES})
endif()

• Run CMake on the FirstReconstruction directory and create a FirstReconstruction-bin,

• Configure and build the project using your favorite compiler,

• Run FirstReconstruction image.mha geometry.xml. If everything runs correctly, you should see a few messages ending with Done! and two new files in the current directory, image.mha and geometry.xml. image.mha is the reconstruction of a sphere from 360 projections and geometry.xml is the geometry file in RTK’s geometry format.

Code

// RTK includes
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
#include <rtkRayEllipsoidIntersectionImageFilter.h>
#include <rtkCudaFDKConeBeamReconstructionFilter.h>
#include <rtkFieldOfViewImageFilter.h>

// ITK includes
#include <itkImageFileWriter.h>

int
main(int argc, char ** argv)
{
  if (argc < 3)
  {
    std::cout << "Usage: FirstReconstruction <outputimage> <outputgeometry>" << std::endl;
    return EXIT_FAILURE;
  }

  // Defines the image type
  using ImageType = itk::CudaImage<float, 3>;

  // Defines the RTK geometry object
  using GeometryType = rtk::ThreeDCircularProjectionGeometry;
  GeometryType::Pointer geometry = GeometryType::New();
  unsigned int          numberOfProjections = 360;
  double                firstAngle = 0;
  double                angularArc = 360;
  unsigned int          sid = 600;  // source to isocenter distance
  unsigned int          sdd = 1200; // source to detector distance
  for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
  {
    double angle = firstAngle + noProj * angularArc / numberOfProjections;
    geometry->AddProjection(sid, sdd, angle);
  }

  // Write the geometry to disk
  rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
  xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
  xmlWriter->SetFilename(argv[2]);
  xmlWriter->SetObject(geometry);
  xmlWriter->WriteFile();

  // Create a stack of empty projection images
  using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
  ConstantImageSourceType::Pointer     constantImageSource = ConstantImageSourceType::New();
  ConstantImageSourceType::PointType   origin;
  ConstantImageSourceType::SpacingType spacing;
  ConstantImageSourceType::SizeType    sizeOutput;

  origin[0] = -127;
  origin[1] = -127;
  origin[2] = 0.;

  sizeOutput[0] = 128;
  sizeOutput[1] = 128;
  sizeOutput[2] = numberOfProjections;

  spacing.Fill(2.);

  constantImageSource->SetOrigin(origin);
  constantImageSource->SetSpacing(spacing);
  constantImageSource->SetSize(sizeOutput);
  constantImageSource->SetConstant(0.);

  // Create projections of an ellipse
  using REIType = rtk::RayEllipsoidIntersectionImageFilter<ImageType, ImageType>;
  REIType::Pointer    rei = REIType::New();
  REIType::VectorType semiprincipalaxis, center;
  semiprincipalaxis.Fill(50.);
  center.Fill(0.);
  center[2] = 10.;
  rei->SetDensity(2.);
  rei->SetAngle(0.);
  rei->SetCenter(center);
  rei->SetAxis(semiprincipalaxis);
  rei->SetGeometry(geometry);
  rei->SetInput(constantImageSource->GetOutput());

  // Create reconstructed image
  ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();
  sizeOutput.Fill(128);
  origin.Fill(-63.5);
  spacing.Fill(1.);
  constantImageSource2->SetOrigin(origin);
  constantImageSource2->SetSpacing(spacing);
  constantImageSource2->SetSize(sizeOutput);
  constantImageSource2->SetConstant(0.);

  // FDK reconstruction
  std::cout << "Reconstructing..." << std::endl;
  using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
  FDKGPUType::Pointer feldkamp = FDKGPUType::New();
  feldkamp->SetInput(0, constantImageSource2->GetOutput());
  feldkamp->SetInput(1, rei->GetOutput());
  feldkamp->SetGeometry(geometry);
  feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
  feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);

  // Field-of-view masking
  using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
  FOVFilterType::Pointer fieldofview = FOVFilterType::New();
  fieldofview->SetInput(0, feldkamp->GetOutput());
  fieldofview->SetProjectionsStack(rei->GetOutput());
  fieldofview->SetGeometry(geometry);

  // Writer
  std::cout << "Writing output image..." << std::endl;
  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(argv[1]);
  writer->SetInput(fieldofview->GetOutput());
  writer->Update();

  std::cout << "Done!" << std::endl;
  return EXIT_SUCCESS;
}

Python

Code

#!/usr/bin/env python
import sys
import itk
from itk import RTK as rtk

if len ( sys.argv ) < 3:
  print( "Usage: FirstReconstruction <outputimage> <outputgeometry>" )
  sys.exit ( 1 )

# Defines the image type
ImageType = itk.Image[itk.F,3]

# Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New()
numberOfProjections = 360
firstAngle = 0.
angularArc = 360.
sid = 600 # source to isocenter distance
sdd = 1200 # source to detector distance
for x in range(0,numberOfProjections):
  angle = firstAngle + x * angularArc / numberOfProjections
  geometry.AddProjection(sid,sdd,angle)

# Writing the geometry to disk
rtk.write_geometry(geometry, sys.argv[2])

# Create a stack of empty projection images
ConstantImageSourceType = rtk.ConstantImageSource[ImageType]
constantImageSource = ConstantImageSourceType.New()
origin = [ -127, -127, 0. ]
sizeOutput = [ 128, 128,  numberOfProjections ]
spacing = [ 2.0, 2.0, 2.0 ]
constantImageSource.SetOrigin( origin )
constantImageSource.SetSpacing( spacing )
constantImageSource.SetSize( sizeOutput )
constantImageSource.SetConstant(0.)

REIType = rtk.RayEllipsoidIntersectionImageFilter[ImageType, ImageType]
rei = REIType.New()
semiprincipalaxis = [ 50, 50, 50]
center = [ 0, 0, 10]
# Set GrayScale value, axes, center...
rei.SetDensity(2)
rei.SetAngle(0)
rei.SetCenter(center)
rei.SetAxis(semiprincipalaxis)
rei.SetGeometry( geometry )
rei.SetInput(constantImageSource.GetOutput())

# Create reconstructed image
constantImageSource2 = ConstantImageSourceType.New()
sizeOutput = [ 128, 128, 128 ]
origin = [ -63.5, -63.5, -63.5 ]
spacing = [ 1.0, 1.0, 1.0 ]
constantImageSource2.SetOrigin( origin )
constantImageSource2.SetSpacing( spacing )
constantImageSource2.SetSize( sizeOutput )
constantImageSource2.SetConstant(0.)

# FDK reconstruction
print("Reconstructing...")
FDKCPUType = rtk.FDKConeBeamReconstructionFilter[ImageType]
feldkamp = FDKCPUType.New()
feldkamp.SetInput(0, constantImageSource2.GetOutput())
feldkamp.SetInput(1, rei.GetOutput())
feldkamp.SetGeometry(geometry)
feldkamp.GetRampFilter().SetTruncationCorrection(0.0)
feldkamp.GetRampFilter().SetHannCutFrequency(0.0)

# Field-of-view masking
FOVFilterType = rtk.FieldOfViewImageFilter[ImageType, ImageType]
fieldofview = FOVFilterType.New()
fieldofview.SetInput(0, feldkamp.GetOutput())
fieldofview.SetProjectionsStack(rei.GetOutput())
fieldofview.SetGeometry(geometry)

# Writer
print("Writing output image...")
WriterType = rtk.ImageFileWriter[ImageType]
writer = WriterType.New()
writer.SetFileName(sys.argv[1])
writer.SetInput(fieldofview.GetOutput())
writer.Update()

print("Done!")

Python Cuda

Code

#!/usr/bin/env python
import os
import sys
import itk
from itk import RTK as rtk

if len ( sys.argv ) < 3:
  print( "Usage: FirstReconstruction <outputimage> <outputgeometry>" )
  sys.exit ( 1 )

# Import Windows CUDA_PATH for dll (required for some Python versions,
# https://docs.python.org/3/whatsnew/3.8.html#bpo-36085-whatsnew)
if sys.platform == 'win32':
  os.add_dll_directory(os.path.join(os.environ['CUDA_PATH'], 'bin'))

# Defines the image type
GPUImageType = rtk.CudaImage[itk.F,3]
CPUImageType = rtk.Image[itk.F,3]

# Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New()
numberOfProjections = 360
firstAngle = 0.
angularArc = 360.
sid = 600 # source to isocenter distance
sdd = 1200 # source to detector distance
for x in range(0,numberOfProjections):
  angle = firstAngle + x * angularArc / numberOfProjections
  geometry.AddProjection(sid,sdd,angle)

# Writing the geometry to disk
rtk.write_geometry(geometry, sys.argv[2])

# Create a stack of empty projection images
ConstantImageSourceType = rtk.ConstantImageSource[GPUImageType]
constantImageSource = ConstantImageSourceType.New()
origin = [ -127, -127, 0. ]
sizeOutput = [ 128, 128,  numberOfProjections ]
spacing = [ 2.0, 2.0, 2.0 ]
constantImageSource.SetOrigin( origin )
constantImageSource.SetSpacing( spacing )
constantImageSource.SetSize( sizeOutput )
constantImageSource.SetConstant(0.)

REIType = rtk.RayEllipsoidIntersectionImageFilter[CPUImageType, CPUImageType]
rei = REIType.New()
semiprincipalaxis = [ 50, 50, 50]
center = [ 0, 0, 10]
# Set GrayScale value, axes, center...
rei.SetDensity(2)
rei.SetAngle(0)
rei.SetCenter(center)
rei.SetAxis(semiprincipalaxis)
rei.SetGeometry( geometry )
rei.SetInput(constantImageSource.GetOutput())

# Create reconstructed image
constantImageSource2 = ConstantImageSourceType.New()
sizeOutput = [ 128, 128, 128 ]
origin = [ -63.5, -63.5, -63.5 ]
spacing = [ 1.0, 1.0, 1.0 ]
constantImageSource2.SetOrigin( origin )
constantImageSource2.SetSpacing( spacing )
constantImageSource2.SetSize( sizeOutput )
constantImageSource2.SetConstant(0.)

# Graft the projections to an itk::CudaImage
projections = GPUImageType.New()
rei.Update()
projections.SetPixelContainer(rei.GetOutput().GetPixelContainer())
projections.CopyInformation(rei.GetOutput())
projections.SetBufferedRegion(rei.GetOutput().GetBufferedRegion())
projections.SetRequestedRegion(rei.GetOutput().GetRequestedRegion())

# FDK reconstruction
print("Reconstructing...")
FDKGPUType = rtk.CudaFDKConeBeamReconstructionFilter
feldkamp = FDKGPUType.New()
feldkamp.SetInput(0, constantImageSource2.GetOutput())
feldkamp.SetInput(1, projections)
feldkamp.SetGeometry(geometry)
feldkamp.GetRampFilter().SetTruncationCorrection(0.0)
feldkamp.GetRampFilter().SetHannCutFrequency(0.0)

# Field-of-view masking
FOVFilterType = rtk.FieldOfViewImageFilter[CPUImageType, CPUImageType]
fieldofview = FOVFilterType.New()
fieldofview.SetInput(0, feldkamp.GetOutput())
fieldofview.SetProjectionsStack(rei.GetOutput())
fieldofview.SetGeometry(geometry)

# Writer
print("Writing output image...")
WriterType = rtk.ImageFileWriter[CPUImageType]
writer = WriterType.New()
writer.SetFileName(sys.argv[1])
writer.SetInput(fieldofview.GetOutput())
writer.Update()

print("Done!")