Inline Reconstruction

InlineReconstruction demonstrates how to perform cone-beam CT reconstruction using RTK’s implementation of FDK while projections are being acquired.

The program simulates the acquisition of X-ray projections in a separate thread, gradually reconstructing the 3D volume as new projections become available.

C++

#include <atomic>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <thread>
#include <vector>

#include <itkExtractImageFilter.h>
#include <itkImageFileWriter.h>

#include <rtkConstantImageSource.h>
#include <rtkFDKConeBeamReconstructionFilter.h>
#include <rtkParkerShortScanImageFilter.h>
#include <rtkProjectionsReader.h>
#include <rtkSheppLoganPhantomFilter.h>

// CUDA headers (if available)
#ifdef RTK_USE_CUDA
#  include <itkCudaImage.h>
#  include <rtkCudaFDKConeBeamReconstructionFilter.h>
#  include <rtkCudaParkerShortScanImageFilter.h>
#endif

// Type definitions
using ImageType = itk::Image<float, 3>;

// Constants for acquisition and reconstruction
const int                       numProjections = 64;
const double                    sid = 1000.0;
const double                    sdd = 1500.0;
const double                    arc = 200.0;
const double                    spacing = 1.0;
const int                       imageSize = 256;
const double                    origin = -0.5 * (imageSize - 1);
const std::chrono::milliseconds acquisitionSleepDuration(100);

std::atomic<int> AcquiredProjectionsCount(0);

// Function for the thread simulating an acquisition
// Simulates and writes a projection every 100 ms.
void
Acquisition()
{
  for (int i = 0; i < numProjections; ++i)
  {
    // Create geometry for the acquisition
    auto geometryAcq = rtk::ThreeDCircularProjectionGeometry::New();
    geometryAcq->AddProjection(sid, sdd, i * arc / numProjections);

    // Create a constant image source for the projection
    auto projection = rtk::ConstantImageSource<ImageType>::New();
    projection->SetOrigin(itk::MakePoint(origin, origin, 0.0));
    projection->SetSize(itk::MakeSize(imageSize, imageSize, 1));
    projection->SetSpacing(itk::MakeVector(spacing, spacing, spacing));
    projection->Update();

    // Create Shepp-Logan phantom filter
    auto phantom = rtk::SheppLoganPhantomFilter<ImageType, ImageType>::New();
    phantom->SetInput(projection->GetOutput());
    phantom->SetGeometry(geometryAcq);
    phantom->SetPhantomScale(100);

    // Create filename for the output projection
    std::ostringstream filename;
    filename << "projection_" << std::setw(3) << std::setfill('0') << i << ".mha";
    phantom->Update();

    // Write the projection image to file
    itk::WriteImage(phantom->GetOutput(), filename.str());
    AcquiredProjectionsCount++;
    std::this_thread::sleep_for(acquisitionSleepDuration);
  }
}

int
main(int, char *[])
{
  // Launch the simulated acquisition
  std::thread acquisitionThread(Acquisition);

  // Create the geometry for reconstruction
  auto geometryRec = rtk::ThreeDCircularProjectionGeometry::New();
  for (int i = 0; i < numProjections; ++i)
  {
    geometryRec->AddProjection(sid, sdd, i * arc / numProjections);
  }

  // Create the reconstruction pipeline
  auto reader = rtk::ProjectionsReader<ImageType>::New();
  auto extractor = itk::ExtractImageFilter<ImageType, ImageType>::New();
  extractor->SetInput(reader->GetOutput());

  double originValue = origin * sid / sdd;
  double spacingValue = spacing * sid / sdd;

#ifdef RTK_USE_CUDA
  using CudaImageType = itk::CudaImage<float, 3>;
  // Use CUDA for Parker short scan image filter
  auto parker = rtk::CudaParkerShortScanImageFilter::New();
  parker->SetGeometry(geometryRec);

  auto reconstructionSource = rtk::ConstantImageSource<CudaImageType>::New();
  reconstructionSource->SetOrigin(itk::MakePoint(originValue, originValue, originValue));
  reconstructionSource->SetSpacing(itk::MakeVector(spacingValue, spacingValue, spacingValue));
  reconstructionSource->SetSize(itk::MakeSize(imageSize, imageSize, imageSize));
  auto fdk = rtk::CudaFDKConeBeamReconstructionFilter::New();
#else
  // Use CPU for Parker short scan image filter
  auto parker = rtk::ParkerShortScanImageFilter<ImageType>::New();
  parker->SetInput(extractor->GetOutput());
  parker->SetGeometry(geometryRec);

  auto reconstructionSource = rtk::ConstantImageSource<ImageType>::New();
  reconstructionSource->SetOrigin(itk::MakePoint(originValue, originValue, originValue));
  reconstructionSource->SetSpacing(itk::MakeVector(spacingValue, spacingValue, spacingValue));
  reconstructionSource->SetSize(itk::MakeSize(imageSize, imageSize, imageSize));
  auto fdk = rtk::FDKConeBeamReconstructionFilter<ImageType>::New();
#endif

  fdk->SetGeometry(geometryRec);
  fdk->SetInput(0, reconstructionSource->GetOutput());
  fdk->SetInput(1, parker->GetOutput());

  // Online reconstruction loop
  int                      reconstructedProjectionsCount = 0;
  std::vector<std::string> projectionFileNames(numProjections, "projection_000.mha");

  while (reconstructedProjectionsCount != numProjections)
  {
    if (reconstructedProjectionsCount < AcquiredProjectionsCount)
    {
      std::cout << "Processing projection #" << reconstructedProjectionsCount << "\r";

      if (reconstructedProjectionsCount == 0)
      {
        reader->SetFileNames(projectionFileNames);
        reader->UpdateOutputInformation();
      }
      else
      {
        // Create filename for the next projection
        std::ostringstream filename;
        filename << "projection_" << std::setw(3) << std::setfill('0') << reconstructedProjectionsCount << ".mha";
        projectionFileNames[reconstructedProjectionsCount] = filename.str();
        reader->SetFileNames(projectionFileNames);

        // Disconnect the pipeline to allow for proper reconstruction
#ifdef RTK_USE_CUDA
        CudaImageType::Pointer reconstructedImage = fdk->GetOutput();
#else
        ImageType::Pointer reconstructedImage = fdk->GetOutput();
#endif
        reconstructedImage->DisconnectPipeline();
        fdk->SetInput(reconstructedImage);
      }

      // Extract and read the new projection
      auto extractedRegion = reader->GetOutput()->GetLargestPossibleRegion();
      extractedRegion.SetIndex(2, reconstructedProjectionsCount);
      extractedRegion.SetSize(2, 1);
      extractor->SetExtractionRegion(extractedRegion);
      extractor->UpdateLargestPossibleRegion();

#ifdef RTK_USE_CUDA
      // Prepare projection for CUDA processing
      auto projection = CudaImageType::New();
      projection->SetPixelContainer(extractor->GetOutput()->GetPixelContainer());
      projection->CopyInformation(extractor->GetOutput());
      projection->SetBufferedRegion(extractor->GetOutput()->GetBufferedRegion());
      projection->SetRequestedRegion(extractor->GetOutput()->GetRequestedRegion());
      parker->SetInput(projection);
#endif

      fdk->Update();
      reconstructedProjectionsCount++;
    }
    else
    {
      // Sleep briefly to avoid busy waiting
      std::this_thread::sleep_for(std::chrono::milliseconds(10));
    }
  }

  acquisitionThread.join();
  itk::WriteImage(fdk->GetOutput(), "fdk.mha");

  return 0;
}

Python

import sys
import os
import time
import threading
import itk
from itk import RTK as rtk

# Parameters controlling the geometry of the simulated acquisition
image_type = itk.Image[itk.F, 3]
has_gpu_capability = hasattr(itk, "CudaImage")
if has_gpu_capability:
    cuda_image_type = itk.CudaImage[itk.F, 3]
nproj = 64
sid = 1000
sdd = 1500
arc = 200
spacing = 1
size = 256
origin = -0.5 * (size - 1)
number_of_acquired_projections = 0


# Function for the thread simulating an acquisition: simulates and writes a
# projection every 100 ms.
def acquisition():
    global number_of_acquired_projections
    for i in range(nproj):
        geometry_acq = rtk.ThreeDCircularProjectionGeometry.New()
        geometry_acq.AddProjection(sid, sdd, i * arc / nproj)
        projection = rtk.constant_image_source(
            origin=[origin, origin, 0.0],
            size=[size, size, 1],
            spacing=[spacing] * 3,
            ttype=image_type,
        )
        projection = rtk.shepp_logan_phantom_filter(
            projection,
            geometry=geometry_acq,
            phantom_scale=100,
        )
        itk.imwrite(projection, f"projection_{i:03d}.mha")

        # The code assumes that incrementing an integer is thread safe, which
        # is the case with CPython
        # https://docs.python.org/3/library/threading.html
        number_of_acquired_projections += 1
        time.sleep(0.1)


# Launches the acquisition thread, then waits for new projections and reconstructs inline / on-the-fly
# Launch the simulated acquisition
thread = threading.Thread(target=acquisition)
thread.start()

# Create the (expected) geometry and filenames
geometry_rec = rtk.ThreeDCircularProjectionGeometry.New()
for i in range(nproj):
    geometry_rec.AddProjection(sid, sdd, i * arc / nproj)

# Create the reconstruction pipeline
reader = rtk.ProjectionsReader[image_type].New()
extractor = itk.ExtractImageFilter[image_type, image_type].New(Input=reader.GetOutput())
if has_gpu_capability:
    parker = rtk.CudaParkerShortScanImageFilter.New(Geometry=geometry_rec)
    reconstruction_source = rtk.ConstantImageSource[cuda_image_type].New(
        Origin=[origin * sid / sdd] * 3,
        Spacing=[spacing * sid / sdd] * 3,
        Size=[size] * 3,
    )
    fdk = rtk.CudaFDKConeBeamReconstructionFilter.New(Geometry=geometry_rec)
else:
    parker = rtk.ParkerShortScanImageFilter[image_type, image_type].New(
        Input=extractor.GetOutput(), Geometry=geometry_rec
    )
    reconstruction_source = rtk.ConstantImageSource[image_type].New(
        Origin=[origin * sid / sdd] * 3,
        Spacing=[spacing * sid / sdd] * 3,
        Size=[size] * 3,
    )
    fdk = rtk.FDKConeBeamReconstructionFilter[image_type, image_type, itk.F].New(
        Geometry=geometry_rec
    )
fdk.SetInput(0, reconstruction_source.GetOutput())
fdk.SetInput(1, parker.GetOutput())

# Do the online / on-the-fly reconstruction: wait for acquired projections
# and use them as soon as a new one is available
number_of_reconstructed_projections = 0
while number_of_reconstructed_projections != nproj:
    new_projection_available = bool(
        number_of_reconstructed_projections < number_of_acquired_projections
    )
    if not new_projection_available:
        time.sleep(0.01)
    else:
        print(
            f"Processing projection #{number_of_reconstructed_projections}",
            end="\r",
        )
        if number_of_reconstructed_projections == 0:
            # First projection, mimick a stack from one file and prepare extracted region
            projection_file_names = ["projection_000.mha"] * nproj
            reader.SetFileNames(projection_file_names)
            reader.UpdateOutputInformation()
            extracted_region = reader.GetOutput().GetLargestPossibleRegion()
            extracted_region.SetSize(2, 1)
        else:
            # Update file name list with the new projection
            projection_file_names[number_of_reconstructed_projections] = (
                f"projection_{number_of_reconstructed_projections:03d}.mha"
            )
            reader.SetFileNames(projection_file_names)

            # Reconnect FDK output to
            reconstructed_image = fdk.GetOutput()
            reconstructed_image.DisconnectPipeline()
            fdk.SetInput(reconstructed_image)

        # Only extract and read the new projection
        extracted_region.SetIndex(2, number_of_reconstructed_projections)
        extractor.SetExtractionRegion(extracted_region)
        extractor.UpdateLargestPossibleRegion()
        if has_gpu_capability:
            projection = cuda_image_type.New()
            projection.SetPixelContainer(extractor.GetOutput().GetPixelContainer())
            projection.CopyInformation(extractor.GetOutput())
            projection.SetBufferedRegion(extractor.GetOutput().GetBufferedRegion())
            projection.SetRequestedRegion(extractor.GetOutput().GetRequestedRegion())
            parker.SetInput(projection)
        fdk.Update()
        number_of_reconstructed_projections += 1
thread.join()
writer = itk.ImageFileWriter[image_type].New(Input=fdk.GetOutput(), FileName="fdk.mha")
writer.Update()