Skip to content

Examples 3D+time #883

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
cyrilmory wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Examples 3D+time #883

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
537a9f2
ENH: Move 3D+time test data generation code into a helper file
cyrilmory Jan 8, 2026
a92760f
ENH: Add FourDConjugateGradient example and corresponding doc page
cyrilmory Jan 8, 2026
52d617b
ENH: Move 4D example data generation to a separate example
cyrilmory Jan 19, 2026
5fd3eba
DOC: Fix description of input data in FourDConjugateGradient example
cyrilmory Jan 30, 2026
ac45656
ENH: Add example FourDSART
cyrilmory Jan 30, 2026
7519030
WIP: Add FourDFDK example. Runs, but the reconstruction is terrible
cyrilmory Jan 31, 2026
dad9771
ENH: Add example FourDROOSTER
cyrilmory Jan 31, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@

#ifdef RTK_USE_CUDA
# include "itkCudaImage.h"
# include "rtkCudaConstantVolumeSeriesSource.h"
#endif

#include <itkImageFileWriter.h>
Expand Down
12 changes: 12 additions & 0 deletions examples/FourDConjugateGradient/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR)

# This project is designed to be built outside the RTK source tree.
project(FourDConjugateGradient)

# Find ITK with RTK
find_package(ITK REQUIRED COMPONENTS RTK)
include(${ITK_USE_FILE})

# Executable(s)
add_executable(FourDConjugateGradient FourDConjugateGradient.cxx)
target_link_libraries(FourDConjugateGradient ${ITK_LIBRARIES})
102 changes: 102 additions & 0 deletions examples/FourDConjugateGradient/FourDConjugateGradient.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
#include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h"
#include "rtkIterationCommands.h"
#include "rtkSignalToInterpolationWeights.h"
#include "rtkReorderProjectionsImageFilter.h"
#include "rtkProjectionsReader.h"
#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h"

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

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

#ifdef RTK_USE_CUDA
using VolumeSeriesType = itk::CudaImage<OutputPixelType, 4>;
using ProjectionStackType = itk::CudaImage<OutputPixelType, 3>;
using VolumeType = itk::CudaImage<OutputPixelType, 3>;
#else
using VolumeSeriesType = itk::Image<OutputPixelType, 4>;
using ProjectionStackType = itk::Image<OutputPixelType, 3>;
using VolumeType = itk::Image<OutputPixelType, 3>;
#endif

// Generate the input volume series, used as initial estimate by 4D conjugate gradient
auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.);
auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.);
auto fourDSize = itk::MakeSize(32, 16, 32, 8);

using ConstantFourDSourceType = rtk::ConstantImageSource<VolumeSeriesType>;
auto fourDSource = ConstantFourDSourceType::New();

fourDSource->SetOrigin(fourDOrigin);
fourDSource->SetSpacing(fourDSpacing);
fourDSource->SetSize(fourDSize);
fourDSource->SetConstant(0.);
fourDSource->Update();

// 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<ProjectionStackType>;
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());

// Re-order geometry and projections
// In the new order, projections with identical phases are packed together
std::vector<double> signal = rtk::ReadSignalFile("four_d_signal.txt");
auto reorder = rtk::ReorderProjectionsImageFilter<ProjectionStackType>::New();
reorder->SetInput(projectionsReader->GetOutput());
reorder->SetInputGeometry(geometry);
reorder->SetInputSignal(signal);
TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update())

// Release the memory holding the stack of original projections
projectionsReader->GetOutput()->ReleaseData();

// Compute the interpolation weights
auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New();
signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal());
signalToInterpolationWeights->SetNumberOfReconstructedFrames(fourDSize[3]);
TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update())

// Set the forward and back projection filters to be used
using ConjugateGradientFilterType =
rtk::FourDConjugateGradientConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>;
auto fourdconjugategradient = ConjugateGradientFilterType::New();

fourdconjugategradient->SetInputVolumeSeries(fourDSource->GetOutput());
fourdconjugategradient->SetNumberOfIterations(2);
#ifdef RTK_USE_CUDA
fourdconjugategradient->SetCudaConjugateGradient(true);
fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_CUDARAYCAST);
fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_CUDAVOXELBASED);
#else
fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_JOSEPH);
fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_VOXELBASED);
#endif

// Set the newly ordered arguments
fourdconjugategradient->SetInputProjectionStack(reorder->GetOutput());
fourdconjugategradient->SetGeometry(reorder->GetOutputGeometry());
fourdconjugategradient->SetWeights(signalToInterpolationWeights->GetOutput());
fourdconjugategradient->SetSignal(reorder->GetOutputSignal());

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

TRY_AND_EXIT_ON_ITK_EXCEPTION(fourdconjugategradient->Update())

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

return EXIT_SUCCESS;
}
18 changes: 18 additions & 0 deletions examples/FourDConjugateGradient/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
# 4D Conjugate Gradient

FourDConjugateGradient shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources.
You can refer to the [projectors documentation](../../documentation/docs/Projectors.md) to see all options available for the back and forwardprojections.

This example reads its data from disk in `your_rtk_build_folder/test/`. The data is generated by the example [GenerateFourDData](../GenerateFourDData/README.md). Therefore, you must run GenerateFourDData at least once before you run this example.
Alternatively, you can also [download the data on Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531).

`````{tab-set}

````{tab-item} C++

```{literalinclude} ./FourDConjugateGradient.cxx
:language: c++
```
````

`````
12 changes: 12 additions & 0 deletions examples/FourDFDK/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR)

# This project is designed to be built outside the RTK source tree.
project(FourDFDK)

# Find ITK with RTK
find_package(ITK REQUIRED COMPONENTS RTK)
include(${ITK_USE_FILE})

# Executable(s)
add_executable(FourDFDK FourDFDK.cxx)
target_link_libraries(FourDFDK ${ITK_LIBRARIES})
159 changes: 159 additions & 0 deletions examples/FourDFDK/FourDFDK.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,159 @@
#include "rtkConfiguration.h"
#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h"
#include "rtkDisplacedDetectorForOffsetFieldOfViewImageFilter.h"
#include "rtkParkerShortScanImageFilter.h"
#include "rtkFDKConeBeamReconstructionFilter.h"
#ifdef RTK_USE_CUDA
# include "rtkCudaDisplacedDetectorImageFilter.h"
// TODO # include "rtkCudaDisplacedDetectorForOffsetFieldOfViewImageFilter.h"
# include "rtkCudaParkerShortScanImageFilter.h"
# include "rtkCudaFDKConeBeamReconstructionFilter.h"
#endif
#include "rtkSelectOneProjectionPerCycleImageFilter.h"
#include "rtkProjectionsReader.h"

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

int
main(int argc, char * argv[])
{
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");

// Check on hardware parameter
#ifndef RTK_USE_CUDA
if (!strcmp(args_info.hardware_arg, "cuda"))
{
std::cerr << "The program has not been compiled with cuda option" << std::endl;
return EXIT_FAILURE;
}
#endif

// Displaced detector weighting
#ifdef RTK_USE_CUDA
using DDFType = rtk::CudaDisplacedDetectorImageFilter;
#else
using DDFType = rtk::DisplacedDetectorForOffsetFieldOfViewImageFilter<OutputImageType>;
#endif

rtk::DisplacedDetectorImageFilter<OutputImageType>::Pointer ddf;
#ifdef RTK_USE_CUDA
ddf = DDFType::New();
#else
ddf = rtk::DisplacedDetectorForOffsetFieldOfViewImageFilter<OutputImageType>::New();
#endif
ddf->SetInput(selector->GetOutput());
ddf->SetGeometry(selector->GetOutputGeometry());

// Short scan image filter
using PSSFCPUType = rtk::ParkerShortScanImageFilter<OutputImageType>;
#ifdef RTK_USE_CUDA
using PSSFType = rtk::CudaParkerShortScanImageFilter;
#else
using PSSFType = rtk::ParkerShortScanImageFilter<OutputImageType>;
#endif
PSSFCPUType::Pointer pssf;

#ifdef RTK_USE_CUDA
pssf = PSSFType::New();
#else
pssf = PSSFCPUType::New();
#endif
pssf->SetInput(ddf->GetOutput());
pssf->SetGeometry(selector->GetOutputGeometry());
pssf->InPlaceOff();

// 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, pssf->GetOutput());
feldkamp->SetGeometry(selector->GetOutputGeometry());
TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->Update());

// Streaming depending on streaming capability of writer
auto streamerBP = itk::StreamingImageFilter<CPUOutputImageType, CPUOutputImageType>::New();
streamerBP->SetInput(feldkamp->GetOutput());
streamerBP->SetNumberOfStreamDivisions(2);
TRY_AND_EXIT_ON_ITK_EXCEPTION(streamerBP->UpdateOutputInformation())

// 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, 4));
fourDSource->SetConstant(0.);
fourDSource->Update();

// Go over each frame, reconstruct 3D frame and paste with iterators in 4D image
for (int f = 0; f < 4; f++)
{
selector->SetPhase(f / (double)4);
TRY_AND_EXIT_ON_ITK_EXCEPTION(streamerBP->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(streamerBP->GetOutput(),
streamerBP->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;
}
18 changes: 18 additions & 0 deletions examples/FourDFDK/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
# 4D FDK

FourDFDK shows how to perform analytical cone-beam CT reconstruction using either CPU or GPU resources.
Only the [FDKBackProjection](https://www.openrtk.org/Doxygen/classrtk_1_1FDKBackProjectionImageFilter.html) and [CudaFDKBackProjection](https://www.openrtk.org/Doxygen/classrtk_1_1CudaFDKBackProjectionImageFilter.html) projectors can be used for back projection, depending on the `RTK_USE_CUDA` CMake variable.

This example reads its data from disk in `your_rtk_build_folder/test/`. The data is generated by the example [GenerateFourDData](../GenerateFourDData/README.md). Therefore, you must run GenerateFourDData at least once before you run this example.
Alternatively, you can also [download the data on Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531).

`````{tab-set}

````{tab-item} C++

```{literalinclude} ./FourDFDK.cxx
:language: c++
```
````

`````
12 changes: 12 additions & 0 deletions examples/FourDROOSTER/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR)

# This project is designed to be built outside the RTK source tree.
project(FourDROOSTER)

# Find ITK with RTK
find_package(ITK REQUIRED COMPONENTS RTK)
include(${ITK_USE_FILE})

# Executable(s)
add_executable(FourDROOSTER FourDROOSTER.cxx)
target_link_libraries(FourDROOSTER ${ITK_LIBRARIES})
Loading
Loading