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.
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 withDone!
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;
}
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 withDone!
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;
}
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!")
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!")