4D FDK

This example shows how to perform filtered backprojection cone-beam CT reconstruction using either CPU or GPU resources depending on how RTK was compiled. It reads its data from disk. The data can be generated by the Generate 4D data example or downloaded from Girder.

C++

#include "rtkConfiguration.h"
#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h"
#include "rtkFDKConeBeamReconstructionFilter.h"
#ifdef RTK_USE_CUDA
#  include "rtkCudaFDKConeBeamReconstructionFilter.h"
#endif
#include "rtkSelectOneProjectionPerCycleImageFilter.h"
#include "rtkProjectionsReader.h"

#ifdef RTK_USE_CUDA
#  include <itkCudaImage.h>
#endif
#include <itkImageFileWriter.h>

int
main(int, char *[])
{
  using OutputPixelType = float;
  constexpr unsigned int Dimension = 3;

  using CPUOutputImageType = itk::Image<OutputPixelType, Dimension>;
#ifdef RTK_USE_CUDA
  using OutputImageType = itk::CudaImage<OutputPixelType, Dimension>;
#else
  using OutputImageType = CPUOutputImageType;
#endif

  // Read geometry, projections and signal
  rtk::ThreeDCircularProjectionGeometry::Pointer geometry;
  TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml"));

  using ReaderType = rtk::ProjectionsReader<OutputImageType>;
  auto                     projectionsReader = ReaderType::New();
  std::vector<std::string> fileNames = std::vector<std::string>();
  fileNames.push_back("four_d_projections.mha");
  projectionsReader->SetFileNames(fileNames);
  TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update());

  // Part specific to 4D
  auto selector = rtk::SelectOneProjectionPerCycleImageFilter<OutputImageType>::New();
  selector->SetInput(projectionsReader->GetOutput());
  selector->SetInputGeometry(geometry);
  selector->SetSignalFilename("four_d_signal.txt");

  // Create one frame of the reconstructed image
  using ConstantImageSourceType = rtk::ConstantImageSource<OutputImageType>;
  auto constantImageSource = ConstantImageSourceType::New();
  constantImageSource->SetOrigin(itk::MakePoint(-63., -31., -63.));
  constantImageSource->SetSpacing(itk::MakeVector(4., 4., 4.));
  constantImageSource->SetSize(itk::MakeSize(32, 16, 32));
  constantImageSource->SetConstant(0.);
  constantImageSource->Update();
  TRY_AND_EXIT_ON_ITK_EXCEPTION(constantImageSource->Update())

  // FDK reconstruction filtering
#ifdef RTK_USE_CUDA
  using FDKType = rtk::CudaFDKConeBeamReconstructionFilter;
#else
  using FDKType = rtk::FDKConeBeamReconstructionFilter<OutputImageType>;
#endif
  auto feldkamp = FDKType::New();
  feldkamp->SetInput(0, constantImageSource->GetOutput());
  feldkamp->SetInput(1, selector->GetOutput());
  feldkamp->SetGeometry(selector->GetOutputGeometry());
  TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->Update());

  // Create empty 4D image
  using FourDOutputImageType = itk::Image<OutputPixelType, Dimension + 1>;
  using ConstantFourDSourceType = rtk::ConstantImageSource<FourDOutputImageType>;
  auto fourDSource = ConstantFourDSourceType::New();
  fourDSource->SetOrigin(itk::MakePoint(-63., -31., -63., 0.));
  fourDSource->SetSpacing(itk::MakeVector(4., 4., 4., 1.));
  fourDSource->SetSize(itk::MakeSize(32, 16, 32, 8));
  fourDSource->SetConstant(0.);
  fourDSource->Update();

  // Go over each frame, reconstruct 3D frame and paste with iterators in 4D image
  for (int f = 0; f < 8; f++)
  {
    selector->SetPhase(f / (double)8);
    TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->UpdateLargestPossibleRegion())

    ConstantFourDSourceType::OutputImageRegionType region;
    region = fourDSource->GetOutput()->GetLargestPossibleRegion();
    region.SetIndex(3, f);
    region.SetSize(3, 1);

    itk::ImageRegionIterator<FourDOutputImageType> it4D(fourDSource->GetOutput(), region);
    itk::ImageRegionIterator<CPUOutputImageType>   it3D(feldkamp->GetOutput(),
                                                      feldkamp->GetOutput()->GetLargestPossibleRegion());
    while (!it3D.IsAtEnd())
    {
      it4D.Set(it3D.Get());
      ++it4D;
      ++it3D;
    }
  }

  // Write
  TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(fourDSource->GetOutput(), "fourdfdk.mha"));

  return EXIT_SUCCESS;
}

Python

import itk
from itk import RTK as rtk

OutputPixelType = itk.F
Dimension = 3

CPUOutputImageType = itk.Image[OutputPixelType, Dimension]

if hasattr(itk, "CudaImage"):
    OutputImageType = itk.CudaImage[OutputPixelType, Dimension]
else:
    OutputImageType = CPUOutputImageType

# Read geometry, projections and signal
geometry = rtk.read_geometry("four_d_geometry.xml")

projectionsReader = rtk.ProjectionsReader[CPUOutputImageType].New()
fileNames = ["four_d_projections.mha"]
projectionsReader.SetFileNames(fileNames)
projectionsReader.Update()

# Part specific to 4D
selector = rtk.SelectOneProjectionPerCycleImageFilter[CPUOutputImageType].New()
selector.SetInput(projectionsReader.GetOutput())
selector.SetInputGeometry(geometry)
selector.SetSignalFilename("four_d_signal.txt")

# Create one frame of the reconstructed image
constantImageSource = rtk.ConstantImageSource[OutputImageType].New()

constantImageSource.SetOrigin([-63.0, -31.0, -63.0])
constantImageSource.SetSpacing([4.0, 4.0, 4.0])
constantImageSource.SetSize([32, 16, 32])
constantImageSource.SetConstant(0.0)
constantImageSource.Update()

# FDK reconstruction filtering
if hasattr(itk, "CudaImage"):
    cuda_projections = itk.CudaImageFromImageFilter[CPUOutputImageType].New()
    cuda_projections.SetInput(selector.GetOutput())

    feldkamp = rtk.CudaFDKConeBeamReconstructionFilter.New()
    feldkamp.SetInput(1, cuda_projections.GetOutput())
else:
    feldkamp = rtk.FDKConeBeamReconstructionFilter[
        OutputImageType, OutputImageType, OutputPixelType
    ].New()
    feldkamp.SetInput(1, selector.GetOutput())

feldkamp.SetInput(0, constantImageSource.GetOutput())
feldkamp.SetGeometry(selector.GetOutputGeometry())

# Create empty 4D image
FourDOutputImageType = itk.Image[OutputPixelType, Dimension + 1]
fourDSource = rtk.ConstantImageSource[FourDOutputImageType].New()
fourDSource.SetOrigin([-63.0, -31.0, -63.0, 0.0])
fourDSource.SetSpacing([4.0, 4.0, 4.0, 1.0])
fourDSource.SetSize([32, 16, 32, 8])
fourDSource.SetConstant(0.0)
fourDSource.Update()

# Go over each frame, reconstruct 3D frame and paste with iterators in 4D image
it4D = itk.GetArrayViewFromImage(fourDSource.GetOutput())
for f in range(8):
    selector.SetPhase(f / 8.0)
    feldkamp.UpdateLargestPossibleRegion()

    it3D = itk.GetArrayFromImage(feldkamp.GetOutput(), ttype=CPUOutputImageType)
    it4D[f] = it3D

# Write
writer = itk.ImageFileWriter[FourDOutputImageType].New()
writer.SetFileName("fourdfdk.mha")
writer.SetInput(fourDSource.GetOutput())
writer.Update()