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 <iostream>
#include <thread>
#include <chrono>
#include <vector>
#include <iomanip>
#include <sstream>
#include <atomic>

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

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

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

// Type definitions
using ImageType = itk::Image<float, 3>;
using CudaImageType = itk::CudaImage<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);

    // Set the origin for the projection
    itk::Point<double, 3> point;
    point[0] = origin;
    point[1] = origin;
    point[2] = 0.0;

    // Create a constant image source for the projection
    auto projection = rtk::ConstantImageSource<ImageType>::New();
    projection->SetOrigin(point);
    projection->SetSize({ imageSize, imageSize, 1 });

    // Set spacing for the projection
    itk::Vector<double, 3> spacingVec;
    spacingVec.Fill(spacing);
    projection->SetSpacing(spacingVec);
    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()
{
  // 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());

  // Set the origin for the reconstruction
  itk::Point<double, 3> originPoint;
  originPoint.Fill(origin * sid / sdd);

  // Set spacing for the reconstruction
  itk::Vector<double, 3> spacingVec;
  spacingVec.Fill(spacing * sid / sdd);

#ifdef RTK_USE_CUDA
  // Use CUDA for Parker short scan image filter
  auto parker = rtk::CudaParkerShortScanImageFilter::New();
  parker->SetGeometry(geometryRec);

  auto reconstructionSource = rtk::ConstantImageSource<CudaImageType>::New();
  reconstructionSource->SetOrigin(originPoint);
  reconstructionSource->SetSpacing(spacingVec);
  reconstructionSource->SetSize({ 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(originPoint);
  reconstructionSource->Setspacing(spacingVec);
  reconstructionSource->SetSize({ 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
        ImageType::Pointer reconstructedImage = fdk->GetOutput();
        reconstructedImage->DisconnectPipeline();

        fdk->SetInput(0, reconstructionSource->GetOutput());
      }

      // 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.], 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(.1)

# Main program: launches the acquition thread, then waits for new projections an reconstructs inline / on-the-fly
if __name__ == "__main__":
    # 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].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].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(.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()