Conjugate Gradient

ConjugateGradient shows how to perform iterative cone-beam CT reconstruction using either CPU and GPU resources. You can refer to the projectors documentation to see all options available for the back and forwardprojections.

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

C++

#include "rtkThreeDCircularProjectionGeometryXMLFile.h"
#include "rtkConjugateGradientConeBeamReconstructionFilter.h"
#include "rtkIterationCommands.h"
#include "rtkProjectGeometricPhantomImageFilter.h"

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

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

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

  constexpr unsigned int numberOfProjections = 60;
  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.;
  std::string            configFileName = "Thorax";
  if (argc > 1 && argv[1] && argv[1][0] != '\0')
  {
    configFileName = argv[1];
  }
  constexpr unsigned int niterations = 3;

  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 projections
  auto projectionsSource = rtk::ConstantImageSource<OutputImageType>::New();
  projectionsSource->SetOrigin(itk::MakePoint(-63., -63., -63.));
  projectionsSource->SetSpacing(itk::MakeVector(2., 2., 2.));
  projectionsSource->SetSize(itk::MakeSize(64, 64, numberOfProjections));

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

  // Create new empty volume
  auto conjugateGradientSource = rtk::ConstantImageSource<OutputImageType>::New();
  conjugateGradientSource->SetOrigin(itk::MakePoint(-63., -63., -63.));
  conjugateGradientSource->SetSpacing(itk::MakeVector(2., 2., 2.));
  conjugateGradientSource->SetSize(itk::MakeSize(64, 64, 64));

  // Set the forward and back projection filters to be used
  using ConjugateGradientFilterType = rtk::ConjugateGradientConeBeamReconstructionFilter<OutputImageType>;
  auto conjugategradient = ConjugateGradientFilterType::New();

  conjugategradient->SetInputVolume(conjugateGradientSource->GetOutput());
  conjugategradient->SetInputProjectionStack(pgp->GetOutput());
  conjugategradient->SetGeometry(geometry);
  conjugategradient->SetNumberOfIterations(niterations);
#ifdef RTK_USE_CUDA
  conjugategradient->SetCudaConjugateGradient(true);
  conjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_CUDARAYCAST);
  conjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_CUDAVOXELBASED);
#else
  conjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_JOSEPH);
  conjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_VOXELBASED);
#endif

  auto verboseIterationCommand = rtk::VerboseIterationCommand<ConjugateGradientFilterType>::New();
  conjugategradient->AddObserver(itk::AnyEvent(), verboseIterationCommand);

  TRY_AND_EXIT_ON_ITK_EXCEPTION(conjugategradient->Update())

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

  return EXIT_SUCCESS;
}

Python

import itk
from itk import RTK as rtk

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

niterations = 10
numberOfProjections = 180
angularArc = 360.0
sid = 600
sdd = 1200
scale = 2.0

# Using RTK applications python APIs to create the geometry and projections
rtk.rtksimulatedgeometry(
    nproj=numberOfProjections,
    arc=angularArc,
    sid=sid,
    sdd=sdd,
    output="geometry.xml",
)
rtk.rtkprojectgeometricphantom(
    geometry="geometry.xml",
    output="projections.mha",
    phantomfile="Thorax",
    phantomscale=scale,
    rotation="1,0,0,0,0,-1,0,1,0",
    spacing=2,
    size=128,
)

projections = itk.imread("projections.mha")
geometry = rtk.read_geometry("geometry.xml")

# Create a constant image source to initialize the reconstruction
conjugate_gradient_source = rtk.constant_image_source(
    ttype=[OutputImageType],
    origin=[-127.0, -127.0, -127.0],
    size=[128, 128, 128],
    spacing=[2.0] * 3,
)

# Set the forward and back projection filters to be used
if hasattr(itk, "CudaImage"):
    OutputCudaImageType = itk.CudaImage[itk.F, Dimension]
    ConjugateGradientFilterType = rtk.ConjugateGradientConeBeamReconstructionFilter[
        OutputCudaImageType, OutputCudaImageType, OutputCudaImageType
    ]
    conjugategradient = ConjugateGradientFilterType.New()
    conjugategradient.SetCudaConjugateGradient(True)
    conjugategradient.SetInput(itk.cuda_image_from_image(conjugate_gradient_source))
    conjugategradient.SetInput(1, itk.cuda_image_from_image(projections))
    conjugategradient.SetBackProjectionFilter(
        ConjugateGradientFilterType.BackProjectionType_BP_CUDAVOXELBASED
    )
    conjugategradient.SetForwardProjectionFilter(
        ConjugateGradientFilterType.ForwardProjectionType_FP_CUDARAYCAST
    )

else:
    ConjugateGradientFilterType = rtk.ConjugateGradientConeBeamReconstructionFilter[
        OutputImageType, OutputImageType, OutputImageType
    ]
    conjugategradient = ConjugateGradientFilterType.New()
    conjugategradient.SetInput(conjugate_gradient_source)
    conjugategradient.SetInput(1, projections)
    conjugategradient.SetBackProjectionFilter(
        ConjugateGradientFilterType.BackProjectionType_BP_VOXELBASED
    )
    conjugategradient.SetForwardProjectionFilter(
        ConjugateGradientFilterType.ForwardProjectionType_FP_JOSEPH
    )

conjugategradient.SetGeometry(geometry)
conjugategradient.SetNumberOfIterations(niterations)


# Create a command to observe the iterations
class VerboseIterationCommand:
    def __init__(self):
        self.count = 0

    def callback(self):
        self.count = self.count + 1
        print(f"Iteration {self.count}", end="\r")


class VerboseEndCommand:
    def callback(self):
        print("")


cmd = VerboseIterationCommand()
conjugategradient.AddObserver(itk.IterationEvent(), cmd.callback)
cmd = VerboseEndCommand()
conjugategradient.AddObserver(itk.EndEvent(), cmd.callback)
conjugategradient.Update()

# Write
writer = itk.ImageFileWriter[OutputImageType].New()
writer.SetFileName("conjugategradient.mha")
writer.SetInput(conjugategradient.GetOutput())
writer.Update()

The results displayed with rtkshowgeometry are: