From 97f311b90be4e21e3df7a0a8eabd9dbbf78dd689 Mon Sep 17 00:00:00 2001 From: Cyril Mory Date: Thu, 8 Jan 2026 13:44:40 +0100 Subject: [PATCH 1/3] ENH: Add 4D C++ reconstruction examples and corresponding doc pages --- .../rtkfourdconjugategradient.cxx | 1 - .../FourDConjugateGradient/CMakeLists.txt | 30 +++ .../FourDConjugateGradient.cxx | 101 +++++++++ examples/FourDConjugateGradient/README.md | 15 ++ examples/FourDFDK/CMakeLists.txt | 30 +++ examples/FourDFDK/FourDFDK.cxx | 103 +++++++++ examples/FourDFDK/README.md | 15 ++ examples/FourDROOSTER/CMakeLists.txt | 30 +++ examples/FourDROOSTER/FourDROOSTER.cxx | 143 ++++++++++++ examples/FourDROOSTER/README.md | 15 ++ examples/GenerateFourDData/CMakeLists.txt | 12 + .../GenerateFourDData/GenerateFourDData.cxx | 214 ++++++++++++++++++ examples/GenerateFourDData/README.md | 21 ++ examples/README.md | 4 + include/rtkCudaExternTemplates.h | 4 + src/rtkCudaExternTemplates.cxx | 2 + test/CMakeLists.txt | 26 +++ 17 files changed, 765 insertions(+), 1 deletion(-) create mode 100644 examples/FourDConjugateGradient/CMakeLists.txt create mode 100644 examples/FourDConjugateGradient/FourDConjugateGradient.cxx create mode 100644 examples/FourDConjugateGradient/README.md create mode 100644 examples/FourDFDK/CMakeLists.txt create mode 100644 examples/FourDFDK/FourDFDK.cxx create mode 100644 examples/FourDFDK/README.md create mode 100644 examples/FourDROOSTER/CMakeLists.txt create mode 100644 examples/FourDROOSTER/FourDROOSTER.cxx create mode 100644 examples/FourDROOSTER/README.md create mode 100644 examples/GenerateFourDData/CMakeLists.txt create mode 100644 examples/GenerateFourDData/GenerateFourDData.cxx create mode 100644 examples/GenerateFourDData/README.md diff --git a/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx b/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx index fc57e78a0..f66820fbb 100644 --- a/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx +++ b/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx @@ -26,7 +26,6 @@ #ifdef RTK_USE_CUDA # include "itkCudaImage.h" -# include "rtkCudaConstantVolumeSeriesSource.h" #endif #include diff --git a/examples/FourDConjugateGradient/CMakeLists.txt b/examples/FourDConjugateGradient/CMakeLists.txt new file mode 100644 index 000000000..b3f694b08 --- /dev/null +++ b/examples/FourDConjugateGradient/CMakeLists.txt @@ -0,0 +1,30 @@ +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 + ITKImageIO +) +if("${ITK_VERSION_MAJOR}" VERSION_LESS "6") + include(${ITK_USE_FILE}) + set(RTK_EXAMPLE_TARGETS ${ITK_LIBRARIES}) +else() + itk_generate_factory_registration() + set( + RTK_EXAMPLE_TARGETS + ${ITK_INTERFACE_LIBRARIES} + ITK::ITKImageIO + ITK::ITKFFTImageFilterInit + ) + message(ITK_INTERFACE_LIBRARIES=${ITK_INTERFACE_LIBRARIES}) +endif() + +# Executable(s) +add_executable(FourDConjugateGradient FourDConjugateGradient.cxx) +target_link_libraries(FourDConjugateGradient ${RTK_EXAMPLE_TARGETS}) diff --git a/examples/FourDConjugateGradient/FourDConjugateGradient.cxx b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx new file mode 100644 index 000000000..95c2c77af --- /dev/null +++ b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx @@ -0,0 +1,101 @@ +#include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h" +#include "rtkIterationCommands.h" +#include "rtkSignalToInterpolationWeights.h" +#include "rtkReorderProjectionsImageFilter.h" +#include "rtkProjectionsReader.h" +#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h" + +#ifdef RTK_USE_CUDA +# include +#endif +#include + +int +main(int argc, char * argv[]) +{ + using OutputPixelType = float; + +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; +#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; + auto fourDSource = ConstantFourDSourceType::New(); + + fourDSource->SetOrigin(fourDOrigin); + fourDSource->SetSpacing(fourDSpacing); + fourDSource->SetSize(fourDSize); + 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; + auto projectionsReader = ReaderType::New(); + std::vector fileNames = std::vector(); + 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 signal = rtk::ReadSignalFile("four_d_signal.txt"); + auto reorder = rtk::ReorderProjectionsImageFilter::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; + 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::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; +} diff --git a/examples/FourDConjugateGradient/README.md b/examples/FourDConjugateGradient/README.md new file mode 100644 index 000000000..77c1a644d --- /dev/null +++ b/examples/FourDConjugateGradient/README.md @@ -0,0 +1,15 @@ +# 4D Conjugate Gradient + +This example shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources. +It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./FourDConjugateGradient.cxx +:language: c++ +``` +```` + +````` diff --git a/examples/FourDFDK/CMakeLists.txt b/examples/FourDFDK/CMakeLists.txt new file mode 100644 index 000000000..0b4cd4703 --- /dev/null +++ b/examples/FourDFDK/CMakeLists.txt @@ -0,0 +1,30 @@ +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 + ITKImageIO +) +if("${ITK_VERSION_MAJOR}" VERSION_LESS "6") + include(${ITK_USE_FILE}) + set(RTK_EXAMPLE_TARGETS ${ITK_LIBRARIES}) +else() + itk_generate_factory_registration() + set( + RTK_EXAMPLE_TARGETS + ${ITK_INTERFACE_LIBRARIES} + ITK::ITKImageIO + ITK::ITKFFTImageFilterInit + ) + message(ITK_INTERFACE_LIBRARIES=${ITK_INTERFACE_LIBRARIES}) +endif() + +# Executable(s) +add_executable(FourDFDK FourDFDK.cxx) +target_link_libraries(FourDFDK ${RTK_EXAMPLE_TARGETS}) diff --git a/examples/FourDFDK/FourDFDK.cxx b/examples/FourDFDK/FourDFDK.cxx new file mode 100644 index 000000000..7211622a4 --- /dev/null +++ b/examples/FourDFDK/FourDFDK.cxx @@ -0,0 +1,103 @@ +#include "rtkConfiguration.h" +#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h" +#include "rtkFDKConeBeamReconstructionFilter.h" +#ifdef RTK_USE_CUDA +# include "rtkCudaFDKConeBeamReconstructionFilter.h" +#endif +#include "rtkSelectOneProjectionPerCycleImageFilter.h" +#include "rtkProjectionsReader.h" + +#ifdef RTK_USE_CUDA +# include +#endif +#include + +int +main(int argc, char * argv[]) +{ + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#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; + auto projectionsReader = ReaderType::New(); + std::vector fileNames = std::vector(); + 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::New(); + selector->SetInput(projectionsReader->GetOutput()); + selector->SetInputGeometry(geometry); + selector->SetSignalFilename("four_d_signal.txt"); + + // Create one frame of the reconstructed image + using ConstantImageSourceType = rtk::ConstantImageSource; + 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; +#endif + auto feldkamp = FDKType::New(); + feldkamp->SetInput(0, constantImageSource->GetOutput()); + feldkamp->SetInput(1, selector->GetOutput()); + feldkamp->SetGeometry(selector->GetOutputGeometry()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->Update()); + + // Create empty 4D image + using FourDOutputImageType = itk::Image; + using ConstantFourDSourceType = rtk::ConstantImageSource; + 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(feldkamp->UpdateLargestPossibleRegion()) + + ConstantFourDSourceType::OutputImageRegionType region; + region = fourDSource->GetOutput()->GetLargestPossibleRegion(); + region.SetIndex(3, f); + region.SetSize(3, 1); + + itk::ImageRegionIterator it4D(fourDSource->GetOutput(), region); + itk::ImageRegionIterator it3D(feldkamp->GetOutput(), + feldkamp->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; +} diff --git a/examples/FourDFDK/README.md b/examples/FourDFDK/README.md new file mode 100644 index 000000000..0b06ab133 --- /dev/null +++ b/examples/FourDFDK/README.md @@ -0,0 +1,15 @@ +# 4D FDK + +This example shows how to perform filtered backprojection cone-beam CT reconstruction using either CPU or GPU resources depending on how RTK was compiled. +It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./FourDFDK.cxx +:language: c++ +``` +```` + +````` diff --git a/examples/FourDROOSTER/CMakeLists.txt b/examples/FourDROOSTER/CMakeLists.txt new file mode 100644 index 000000000..eaa95e17c --- /dev/null +++ b/examples/FourDROOSTER/CMakeLists.txt @@ -0,0 +1,30 @@ +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 + ITKImageIO +) +if("${ITK_VERSION_MAJOR}" VERSION_LESS "6") + include(${ITK_USE_FILE}) + set(RTK_EXAMPLE_TARGETS ${ITK_LIBRARIES}) +else() + itk_generate_factory_registration() + set( + RTK_EXAMPLE_TARGETS + ${ITK_INTERFACE_LIBRARIES} + ITK::ITKImageIO + ITK::ITKFFTImageFilterInit + ) + message(ITK_INTERFACE_LIBRARIES=${ITK_INTERFACE_LIBRARIES}) +endif() + +# Executable(s) +add_executable(FourDROOSTER FourDROOSTER.cxx) +target_link_libraries(FourDROOSTER ${RTK_EXAMPLE_TARGETS}) diff --git a/examples/FourDROOSTER/FourDROOSTER.cxx b/examples/FourDROOSTER/FourDROOSTER.cxx new file mode 100644 index 000000000..825641211 --- /dev/null +++ b/examples/FourDROOSTER/FourDROOSTER.cxx @@ -0,0 +1,143 @@ +#include "rtkFourDROOSTERConeBeamReconstructionFilter.h" +#include "rtkIterationCommands.h" +#include "rtkSignalToInterpolationWeights.h" +#include "rtkReorderProjectionsImageFilter.h" +#include "rtkProjectionsReader.h" +#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h" + +#ifdef RTK_USE_CUDA +# include +#endif +#include + +int +main(int argc, char * argv[]) +{ + using OutputPixelType = float; + using DVFVectorType = itk::CovariantVector; +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; + using DVFSequenceImageType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; + using DVFSequenceImageType = itk::Image; +#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; + auto fourDSource = ConstantFourDSourceType::New(); + + fourDSource->SetOrigin(fourDOrigin); + fourDSource->SetSpacing(fourDSpacing); + fourDSource->SetSize(fourDSize); + 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; + auto projectionsReader = ReaderType::New(); + std::vector fileNames = std::vector(); + 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 signal = rtk::ReadSignalFile("four_d_signal.txt"); + auto reorder = rtk::ReorderProjectionsImageFilter::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 ROOSTERFilterType = rtk::FourDROOSTERConeBeamReconstructionFilter; + auto rooster = ROOSTERFilterType::New(); + + rooster->SetInputVolumeSeries(fourDSource->GetOutput()); + rooster->SetCG_iterations(2); + rooster->SetMainLoop_iterations(2); +#ifdef RTK_USE_CUDA + rooster->SetCudaConjugateGradient(true); + rooster->SetUseCudaCyclicDeformation(true); + rooster->SetForwardProjectionFilter(ROOSTERFilterType::FP_CUDARAYCAST); + rooster->SetBackProjectionFilter(ROOSTERFilterType::BP_CUDAVOXELBASED); +#else + rooster->SetForwardProjectionFilter(ROOSTERFilterType::FP_JOSEPH); + rooster->SetBackProjectionFilter(ROOSTERFilterType::BP_VOXELBASED); +#endif + + // Set the newly ordered arguments + rooster->SetInputProjectionStack(reorder->GetOutput()); + rooster->SetGeometry(reorder->GetOutputGeometry()); + rooster->SetWeights(signalToInterpolationWeights->GetOutput()); + rooster->SetSignal(reorder->GetOutputSignal()); + + // For each optional regularization step, set whether or not + // it should be performed, and provide the necessary inputs + + // Positivity + rooster->SetPerformPositivity(true); + + // Motion mask + rooster->SetPerformMotionMask(false); + + // Spatial TV + rooster->SetGammaTVSpace(0.1); + rooster->SetTV_iterations(4); + rooster->SetPerformTVSpatialDenoising(true); + + // Spatial wavelets + rooster->SetPerformWaveletsSpatialDenoising(false); + + // Temporal TV + rooster->SetGammaTVTime(0.1); + rooster->SetTV_iterations(4); + rooster->SetPerformTVTemporalDenoising(true); + + // Temporal L0 + rooster->SetPerformL0TemporalDenoising(false); + + // Total nuclear variation + rooster->SetPerformTNVDenoising(false); + + // Warping + rooster->SetPerformWarping(true); + rooster->SetUseNearestNeighborInterpolationInWarping(false); + DVFSequenceImageType::Pointer dvf; + TRY_AND_EXIT_ON_ITK_EXCEPTION(dvf = itk::ReadImage("four_d_dvf.mha")) + rooster->SetDisplacementField(dvf); + rooster->SetComputeInverseWarpingByConjugateGradient(false); + DVFSequenceImageType::Pointer idvf; + TRY_AND_EXIT_ON_ITK_EXCEPTION(idvf = itk::ReadImage("four_d_idvf.mha")) + rooster->SetInverseDisplacementField(idvf); + + auto verboseIterationCommand = rtk::VerboseIterationCommand::New(); + rooster->AddObserver(itk::AnyEvent(), verboseIterationCommand); + + TRY_AND_EXIT_ON_ITK_EXCEPTION(rooster->Update()) + + // Write + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(rooster->GetOutput(), "fourdrooster.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/FourDROOSTER/README.md b/examples/FourDROOSTER/README.md new file mode 100644 index 000000000..477d17f1e --- /dev/null +++ b/examples/FourDROOSTER/README.md @@ -0,0 +1,15 @@ +# 4D ROOSTER + +This example shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources. +It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./FourDROOSTER.cxx +:language: c++ +``` +```` + +````` diff --git a/examples/GenerateFourDData/CMakeLists.txt b/examples/GenerateFourDData/CMakeLists.txt new file mode 100644 index 000000000..04b65011a --- /dev/null +++ b/examples/GenerateFourDData/CMakeLists.txt @@ -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(GenerateFourDData) + +# Find ITK with RTK +find_package(ITK REQUIRED COMPONENTS RTK) +include(${ITK_USE_FILE}) + +# Executable(s) +add_executable(GenerateFourDData GenerateFourDData.cxx) +target_link_libraries(GenerateFourDData ${ITK_LIBRARIES}) diff --git a/examples/GenerateFourDData/GenerateFourDData.cxx b/examples/GenerateFourDData/GenerateFourDData.cxx new file mode 100644 index 000000000..a4cda7af1 --- /dev/null +++ b/examples/GenerateFourDData/GenerateFourDData.cxx @@ -0,0 +1,214 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "rtkConstantImageSource.h" +#include "rtkRayEllipsoidIntersectionImageFilter.h" +#include "rtkDrawEllipsoidImageFilter.h" +#include "rtkThreeDCircularProjectionGeometry.h" + +int +main(int argc, char * argv[]) +{ + using PixelType = float; + using DVFVectorType = itk::CovariantVector; + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; + using DVFSequenceImageType = itk::Image; + using GeometryType = rtk::ThreeDCircularProjectionGeometry; + constexpr unsigned int NumberOfProjectionImages = 64; + + /* ========================= + * Constant volume source + * ========================= */ + using ConstantVolumeSourceType = rtk::ConstantImageSource; + auto volumeSource = ConstantVolumeSourceType::New(); + + auto origin = itk::MakePoint(-63., -31., -63.); + auto spacing = itk::MakeVector(4., 4., 4.); + auto size = itk::MakeSize(32, 16, 32); + + volumeSource->SetOrigin(origin); + volumeSource->SetSpacing(spacing); + volumeSource->SetSize(size); + volumeSource->Update(); + + /* ========================= + * Projection accumulation + * ========================= */ + using ConstantProjectionSourceType = rtk::ConstantImageSource; + auto projectionsSource = ConstantProjectionSourceType::New(); + + auto projOrigin = itk::MakePoint(-254., -254., -254.); + auto projSpacing = itk::MakeVector(8., 8., 1.); + auto projSize = itk::MakeSize(64, 64, NumberOfProjectionImages); + + projectionsSource->SetOrigin(projOrigin); + projectionsSource->SetSpacing(projSpacing); + projectionsSource->SetSize(projSize); + projectionsSource->Update(); + + auto oneProjectionSource = ConstantProjectionSourceType::New(); + projSize[2] = 1; + oneProjectionSource->SetOrigin(projOrigin); + oneProjectionSource->SetSpacing(projSpacing); + oneProjectionSource->SetSize(projSize); + + using REIType = rtk::RayEllipsoidIntersectionImageFilter; + using PasteType = itk::PasteImageFilter; + + /* Allocate explicit accumulator image */ + auto accumulated = ProjectionStackType::New(); + accumulated->SetOrigin(projectionsSource->GetOutput()->GetOrigin()); + accumulated->SetSpacing(projectionsSource->GetOutput()->GetSpacing()); + accumulated->SetDirection(projectionsSource->GetOutput()->GetDirection()); + accumulated->SetRegions(projectionsSource->GetOutput()->GetLargestPossibleRegion()); + accumulated->Allocate(); + accumulated->FillBuffer(0.); + + auto paste = PasteType::New(); + paste->SetDestinationImage(accumulated); + + itk::Index<3> destIndex; + destIndex.Fill(0); + + // Create the signal file and the geometry, to be filled in the for loop + std::string signalFileName = "four_d_signal.txt"; + std::ofstream signalFile(signalFileName.c_str()); + auto geometry = GeometryType::New(); + + for (unsigned int i = 0; i < NumberOfProjectionImages; ++i) + { + geometry->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); + + auto geom = GeometryType::New(); + geom->AddProjection(600., 1200., i * 360. / NumberOfProjectionImages, 0, 0, 0, 0, 20, 15); + + auto e1 = REIType::New(); + e1->SetInput(oneProjectionSource->GetOutput()); + e1->SetGeometry(geom); + e1->SetDensity(2.); + e1->SetAxis(itk::MakeVector(60., 30., 60.)); + e1->SetCenter(itk::MakePoint(0., 0., 0.)); + e1->InPlaceOff(); + e1->Update(); + + auto e2 = REIType::New(); + e2->SetInput(e1->GetOutput()); + e2->SetGeometry(geom); + e2->SetDensity(-1.); + e2->SetAxis(itk::MakeVector(8., 8., 8.)); + auto center = itk::MakePoint(4 * (itk::Math::abs((4 + i) % 8 - 4.) - 2.), 0., 0.); + e2->SetCenter(center); + e2->InPlaceOff(); + e2->Update(); + + paste->SetSourceImage(e2->GetOutput()); + paste->SetSourceRegion(e2->GetOutput()->GetLargestPossibleRegion()); + paste->SetDestinationIndex(destIndex); + paste->Update(); + + accumulated = paste->GetOutput(); + paste->SetDestinationImage(accumulated); + + destIndex[2]++; + + signalFile << (i % 8) / 8. << std::endl; + } + + signalFile.close(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(accumulated, "four_d_projections.mha")); + TRY_AND_EXIT_ON_ITK_EXCEPTION(rtk::WriteGeometry(geometry, "four_d_geometry.xml")) + + /* ========================= + * DVF & inverse DVF + * ========================= */ + auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.); + auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.); + auto fourDSize = itk::MakeSize(32, 16, 32, 8); + + auto dvf = DVFSequenceImageType::New(); + auto idvf = DVFSequenceImageType::New(); + + typename DVFSequenceImageType::RegionType region; + auto dvfSize = itk::MakeSize(fourDSize[0], fourDSize[1], fourDSize[2], 2); + region.SetSize(dvfSize); + + dvf->SetRegions(region); + dvf->SetOrigin(fourDOrigin); + dvf->SetSpacing(fourDSpacing); + dvf->Allocate(); + + idvf->SetRegions(region); + idvf->SetOrigin(fourDOrigin); + idvf->SetSpacing(fourDSpacing); + idvf->Allocate(); + + itk::ImageRegionIteratorWithIndex it(dvf, region); + itk::ImageRegionIteratorWithIndex iit(idvf, region); + + DVFVectorType v; + typename DVFSequenceImageType::IndexType centerIndex; + centerIndex.Fill(0); + centerIndex[0] = dvfSize[0] / 2; + centerIndex[1] = dvfSize[1] / 2; + centerIndex[2] = dvfSize[2] / 2; + + for (; !it.IsAtEnd(); ++it, ++iit) + { + v.Fill(0.); + auto d = it.GetIndex() - centerIndex; + if (0.3 * d[0] * d[0] + d[1] * d[1] + d[2] * d[2] < 40) + v[0] = (it.GetIndex()[3] == 0 ? -8. : 8.); + it.Set(v); + iit.Set(-v); + } + + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(dvf, "four_d_dvf.mha")); + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(idvf, "four_d_idvf.mha")); + + /* ========================= + * Ground truth + * ========================= */ + auto join = itk::JoinSeriesImageFilter::New(); + + for (unsigned int t = 0; t < fourDSize[3]; ++t) + { + using DEType = rtk::DrawEllipsoidImageFilter; + + auto de1 = DEType::New(); + de1->SetInput(volumeSource->GetOutput()); + de1->SetDensity(2.); + de1->SetAxis(itk::MakeVector(60., 30., 60.)); + de1->SetCenter(itk::MakePoint(0., 0., 0.)); + de1->InPlaceOff(); + de1->Update(); + + auto de2 = DEType::New(); + de2->SetInput(de1->GetOutput()); + de2->SetDensity(-1.); + de2->SetAxis(itk::MakeVector(8., 8., 8.)); + de2->SetCenter(itk::MakePoint(4 * (itk::Math::abs((4 + t) % 8 - 4.) - 2.), 0., 0.)); + de2->InPlaceOff(); + de2->Update(); + + using DuplicatorType = itk::ImageDuplicator; + auto duplicator = DuplicatorType::New(); + duplicator->SetInputImage(de2->GetOutput()); + duplicator->Update(); + + join->SetInput(t, duplicator->GetOutput()); + } + + join->Update(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(join->GetOutput(), "four_d_ground_truth.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/GenerateFourDData/README.md b/examples/GenerateFourDData/README.md new file mode 100644 index 000000000..1cc19be41 --- /dev/null +++ b/examples/GenerateFourDData/README.md @@ -0,0 +1,21 @@ +# Generate 4D data + +This example shows how to generate a full set of input data to run 4D examples. It contains: +- a moving phantom, made of two ellipsoids, one of which moves along the first dimension +- a geometry +- a phase signal +- the set of projections of the moving phantom +- a deformation vector field describing the motion of the moving ellipsoid +- the corresponding inverse deformation vector field + +You can also skip this part and [download the data](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./GenerateFourDData.cxx +:language: c++ +``` +```` +````` diff --git a/examples/README.md b/examples/README.md index d6d775f18..f67e085b3 100644 --- a/examples/README.md +++ b/examples/README.md @@ -22,4 +22,8 @@ providing efficient implementations for high-performance applications. ./InlineReconstruction/README.md ./AddNoise/README.md ./GeometricPhantom/README.md +./GenerateFourDData/README.md +./FourDFDK/README.md +./FourDConjugateGradient/README.md +./FourDROOSTER/README.md ``` diff --git a/include/rtkCudaExternTemplates.h b/include/rtkCudaExternTemplates.h index 1d7704ba1..1638acb50 100644 --- a/include/rtkCudaExternTemplates.h +++ b/include/rtkCudaExternTemplates.h @@ -22,7 +22,9 @@ #include "rtkConfiguration.h" #ifdef RTK_USE_CUDA +# include # include +# include # include # include # include @@ -49,6 +51,8 @@ extern template class RTK_EXPORT_EXPLICIT itk::ImageSource, itk::CudaImage>; extern template class RTK_EXPORT_EXPLICIT itk::InPlaceImageFilter, itk::CudaImage>; extern template class RTK_EXPORT_EXPLICIT itk::ImageSource, 3>>; +extern template class RTK_EXPORT_EXPLICIT itk::ExtractImageFilter, itk::CudaImage>; +extern template class RTK_EXPORT_EXPLICIT itk::CropImageFilter, itk::CudaImage>; ITK_GCC_PRAGMA_DIAG_POP() # endif diff --git a/src/rtkCudaExternTemplates.cxx b/src/rtkCudaExternTemplates.cxx index 594c9b396..9d6d8fd6f 100644 --- a/src/rtkCudaExternTemplates.cxx +++ b/src/rtkCudaExternTemplates.cxx @@ -25,6 +25,8 @@ template class itk::ImageSource>; template class itk::ImageToImageFilter, itk::CudaImage>; template class itk::InPlaceImageFilter, itk::CudaImage>; template class itk::ImageSource, 3>>; +template class itk::ExtractImageFilter, itk::CudaImage>; +template class itk::CropImageFilter, itk::CudaImage>; // RTK base class explicit instantiation definitions # include "rtkBackProjectionImageFilter.h" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fb744596c..c066b14ad 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -409,6 +409,32 @@ rtk_add_test(rtkConjugateGradientExampleTest ../examples/ConjugateGradient/ConjugateGradient.cxx DATA{Input/Forbild/Thorax} ) +rtk_add_test(rtkGenerateFourDDataExampleTest + ../examples/GenerateFourDData/GenerateFourDData.cxx +) +set_tests_properties( + rtkGenerateFourDDataExampleTest + PROPERTIES + FIXTURES_SETUP + FourDData +) +rtk_add_test(rtkFourDConjugateGradientExampleTest + ../examples/FourDConjugateGradient/FourDConjugateGradient.cxx +) +rtk_add_test(rtkFourDFDKExampleTest + ../examples/FourDFDK/FourDFDK.cxx +) +rtk_add_test(rtkFourDROOSTERExampleTest + ../examples/FourDROOSTER/FourDROOSTER.cxx +) +set_tests_properties( + rtkFourDConjugateGradientExampleTest + rtkFourDFDKExampleTest + rtkFourDROOSTERExampleTest + PROPERTIES + FIXTURES_REQUIRED + FourDData +) if(ITK_WRAP_PYTHON) itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py) From a3ee27437c3aced89570ee783bc6d217f0b8e7f1 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 10 Jun 2026 11:42:48 +0200 Subject: [PATCH 2/3] ENH: Add 4D python reconstruction examples --- .../FourDConjugateGradient.py | 114 +++++++++ examples/FourDFDK/FourDFDK.py | 87 +++++++ examples/FourDROOSTER/FourDROOSTER.py | 163 +++++++++++++ examples/FourDROOSTER/README.md | 7 + .../GenerateFourDData/GenerateFourDData.py | 227 ++++++++++++++++++ examples/GenerateFourDData/README.md | 7 + wrapping/CMakeLists.txt | 1 + wrapping/itkImageRegionRTK.wrap | 7 + wrapping/itkImageSourceRTK.wrap | 1 + wrapping/rtkConstantImageSource.wrap | 1 + 10 files changed, 615 insertions(+) create mode 100644 examples/FourDConjugateGradient/FourDConjugateGradient.py create mode 100644 examples/FourDFDK/FourDFDK.py create mode 100644 examples/FourDROOSTER/FourDROOSTER.py create mode 100644 examples/GenerateFourDData/GenerateFourDData.py diff --git a/examples/FourDConjugateGradient/FourDConjugateGradient.py b/examples/FourDConjugateGradient/FourDConjugateGradient.py new file mode 100644 index 000000000..c2147a3e9 --- /dev/null +++ b/examples/FourDConjugateGradient/FourDConjugateGradient.py @@ -0,0 +1,114 @@ +import itk +from itk import RTK as rtk + +OutputPixelType = itk.F +CpuProjectionStackType = itk.Image[OutputPixelType, 3] +CpuVolumeSeriesType = itk.Image[OutputPixelType, 4] + +if hasattr(itk, "CudaImage"): + ProjectionStackType = itk.CudaImage[OutputPixelType, 3] + VolumeSeriesType = itk.CudaImage[OutputPixelType, 4] +else: + ProjectionStackType = CpuProjectionStackType + VolumeSeriesType = CpuVolumeSeriesType + + +# Generate the input volume series, used as initial estimate by 4D conjugate gradient +four_d_source = rtk.constant_image_source( + ttype=[VolumeSeriesType], + origin=[-63.0, -31.0, -63.0, 0.0], + spacing=[4.0, 4.0, 4.0, 1.0], + size=[32, 16, 32, 8], +) + +# Read geometry, projections and signal +geometry = rtk.read_geometry("four_d_geometry.xml") + +projections_reader = rtk.ProjectionsReader[CpuProjectionStackType].New() +file_names = ["four_d_projections.mha"] +projections_reader.SetFileNames(file_names) +projections_reader.Update() + +signal = rtk.read_signal_file("four_d_signal.txt") + +# Re-order geometry and projections +# In the new order, projections with identical phases are packed together +reorder = rtk.ReorderProjectionsImageFilter[ + ProjectionStackType, ProjectionStackType +].New() +if hasattr(itk, "CudaImage"): + reorder.SetInput(itk.cuda_image_from_image(projections_reader.GetOutput())) +else: + reorder.SetInput(projections_reader.GetOutput()) +reorder.SetInputGeometry(geometry) +reorder.SetInputSignal(signal) +reorder.Update() + +# Release the memory holding the stack of original projections +projections_reader.GetOutput().ReleaseData() + +# Compute the interpolation weights +signal_to_interpolation_weights = rtk.SignalToInterpolationWeights.New() +signal_to_interpolation_weights.SetSignal(reorder.GetOutputSignal()) +signal_to_interpolation_weights.SetNumberOfReconstructedFrames( + four_d_source.GetLargestPossibleRegion().GetSize()[3] +) +signal_to_interpolation_weights.Update() + +# Set the forward and back projection filters to be used +ConjugateGradientFilterType = rtk.FourDConjugateGradientConeBeamReconstructionFilter[ + VolumeSeriesType, ProjectionStackType +] +fourdconjugategradient = ConjugateGradientFilterType.New() +if hasattr(itk, "CudaImage"): + fourdconjugategradient.SetCudaConjugateGradient(True) + fourdconjugategradient.SetForwardProjectionFilter( + ConjugateGradientFilterType.ForwardProjectionType_FP_CUDARAYCAST + ) + fourdconjugategradient.SetBackProjectionFilter( + ConjugateGradientFilterType.BackProjectionType_BP_CUDAVOXELBASED + ) +else: + fourdconjugategradient.SetForwardProjectionFilter( + ConjugateGradientFilterType.ForwardProjectionType_FP_JOSEPH + ) + fourdconjugategradient.SetBackProjectionFilter( + ConjugateGradientFilterType.BackProjectionType_BP_VOXELBASED + ) + +fourdconjugategradient.SetInputVolumeSeries(four_d_source) +fourdconjugategradient.SetNumberOfIterations(2) + +# Set the newly ordered arguments +fourdconjugategradient.SetInputProjectionStack(reorder.GetOutput()) +fourdconjugategradient.SetGeometry(reorder.GetOutputGeometry()) +fourdconjugategradient.SetWeights(signal_to_interpolation_weights.GetOutput()) +fourdconjugategradient.SetSignal(reorder.GetOutputSignal()) + + +class VerboseIterationCommand: + def __init__(self): + self.count = 0 + + def callback(self): + self.count += 1 + print(f"Iteration {self.count}", end="\r") + + +class VerboseEndCommand: + def callback(self): + print("") + + +cmd = VerboseIterationCommand() +fourdconjugategradient.AddObserver(itk.IterationEvent(), cmd.callback) +cmd_end = VerboseEndCommand() +fourdconjugategradient.AddObserver(itk.EndEvent(), cmd_end.callback) + +fourdconjugategradient.Update() + +# Write +writer = itk.ImageFileWriter[CpuVolumeSeriesType].New() +writer.SetFileName("fourdconjugategradient.mha") +writer.SetInput(fourdconjugategradient.GetOutput()) +writer.Update() diff --git a/examples/FourDFDK/FourDFDK.py b/examples/FourDFDK/FourDFDK.py new file mode 100644 index 000000000..2382dec5c --- /dev/null +++ b/examples/FourDFDK/FourDFDK.py @@ -0,0 +1,87 @@ +import itk +from itk import RTK as rtk + +OutputPixelType = itk.F +Dimension = 3 + +CPUOutputImageType = itk.Image[OutputPixelType, Dimension] + +if hasattr(itk, "CudaImage"): + OutputImageType = itk.CudaImage[OutputPixelType, Dimension] +else: + OutputImageType = CPUOutputImageType + +# Read geometry, projections and signal +geometry = rtk.read_geometry("four_d_geometry.xml") + +projections_reader = rtk.ProjectionsReader[CPUOutputImageType].New() +file_names = ["four_d_projections.mha"] +projections_reader.SetFileNames(file_names) +projections_reader.Update() + +# Part specific to 4D +selector = rtk.SelectOneProjectionPerCycleImageFilter[CPUOutputImageType].New() +selector.SetInput(projections_reader.GetOutput()) +selector.SetInputGeometry(geometry) +selector.SetSignalFilename("four_d_signal.txt") +selector.Update() + + +# Create one frame of the reconstructed image +constant_image_source = rtk.constant_image_source( + ttype=[OutputImageType], + origin=[-63.0, -31.0, -63.0], + spacing=[4.0, 4.0, 4.0], + size=[32, 16, 32], + constant=0.0, +) + +# FDK reconstruction filtering +if hasattr(itk, "CudaImage"): + cuda_projections = itk.CudaImageFromImageFilter[CPUOutputImageType].New() + cuda_projections.SetInput(selector.GetOutput()) + + feldkamp = rtk.CudaFDKConeBeamReconstructionFilter.New() + feldkamp.SetInput(0, constant_image_source) + feldkamp.SetInput(1, cuda_projections.GetOutput()) +else: + feldkamp = rtk.FDKConeBeamReconstructionFilter[ + OutputImageType, OutputImageType, OutputPixelType + ].New() + feldkamp.SetInput(0, constant_image_source) + feldkamp.SetInput(1, selector.GetOutput()) +feldkamp.SetGeometry(selector.GetOutputGeometry()) +feldkamp.InPlaceOff() + +# Create empty 4D image +FourDOutputImageType = itk.Image[OutputPixelType, Dimension + 1] +four_d_source = rtk.constant_image_source( + ttype=[FourDOutputImageType], + origin=[-63.0, -31.0, -63.0, 0.0], + spacing=[4.0, 4.0, 4.0, 1.0], + size=[32, 16, 32, 4], + constant=0.0, +) + +# Go over each frame, reconstruct 3D frame and paste with iterators in 4D image +for f in range(4): + selector.SetPhase(f / 4.0) + feldkamp.UpdateLargestPossibleRegion() + output_3d = feldkamp.GetOutput() + + region_4d = four_d_source.GetLargestPossibleRegion() + region_4d.SetIndex(3, f) + region_4d.SetSize(3, 1) + + size_3d = output_3d.GetLargestPossibleRegion().GetSize() + for z in range(size_3d[2]): + for y in range(size_3d[1]): + for x in range(size_3d[0]): + pixel = output_3d.GetPixel([x, y, z]) + four_d_source.SetPixel([x, y, z, f], pixel) + +# Write +writer = itk.ImageFileWriter[FourDOutputImageType].New() +writer.SetFileName("fourdfdk.mha") +writer.SetInput(four_d_source) +writer.Update() diff --git a/examples/FourDROOSTER/FourDROOSTER.py b/examples/FourDROOSTER/FourDROOSTER.py new file mode 100644 index 000000000..a6b6995d1 --- /dev/null +++ b/examples/FourDROOSTER/FourDROOSTER.py @@ -0,0 +1,163 @@ +import itk +from itk import RTK as rtk + +OutputPixelType = itk.F +CpuProjectionStackType = itk.Image[OutputPixelType, 3] +CpuVolumeSeriesType = itk.Image[OutputPixelType, 4] +DVFVectorType = itk.CovariantVector[OutputPixelType, 3] +CpuDVFSequenceImageType = itk.Image[DVFVectorType, 4] + +if hasattr(itk, "CudaImage"): + ProjectionStackType = itk.CudaImage[OutputPixelType, 3] + VolumeSeriesType = itk.CudaImage[OutputPixelType, 4] +else: + ProjectionStackType = CpuProjectionStackType + VolumeSeriesType = CpuVolumeSeriesType + + +# Generate the input volume series, used as initial estimate by 4D conjugate gradient +four_d_source = rtk.constant_image_source( + ttype=[VolumeSeriesType], + origin=[-63.0, -31.0, -63.0, 0.0], + spacing=[4.0, 4.0, 4.0, 1.0], + size=[32, 16, 32, 8], +) + +# Read geometry, projections and signal +geometry = rtk.read_geometry("four_d_geometry.xml") + +projections_reader = rtk.ProjectionsReader[CpuProjectionStackType].New() +file_names = ["four_d_projections.mha"] +projections_reader.SetFileNames(file_names) +projections_reader.Update() + +signal = rtk.read_signal_file("four_d_signal.txt") + +# Re-order geometry and projections +# In the new order, projections with identical phases are packed together +reorder = rtk.ReorderProjectionsImageFilter[ + ProjectionStackType, ProjectionStackType +].New() +if hasattr(itk, "CudaImage"): + reorder.SetInput(itk.cuda_image_from_image(projections_reader.GetOutput())) +else: + reorder.SetInput(projections_reader.GetOutput()) +reorder.SetInputGeometry(geometry) +reorder.SetInputSignal(signal) +reorder.Update() + +# Release the memory holding the stack of original projections +projections_reader.GetOutput().ReleaseData() + +# Compute the interpolation weights +signal_to_interpolation_weights = rtk.SignalToInterpolationWeights.New() +signal_to_interpolation_weights.SetSignal(reorder.GetOutputSignal()) +signal_to_interpolation_weights.SetNumberOfReconstructedFrames( + four_d_source.GetLargestPossibleRegion().GetSize()[3] +) +signal_to_interpolation_weights.Update() + +# Set the forward and back projection filters to be used +ROOSTERFilterType = rtk.FourDROOSTERConeBeamReconstructionFilter[ + VolumeSeriesType, ProjectionStackType +] +rooster = ROOSTERFilterType.New() + +rooster.SetInputVolumeSeries(four_d_source) +rooster.SetCG_iterations(2) +rooster.SetMainLoop_iterations(2) +if hasattr(itk, "CudaImage"): + rooster.SetCudaConjugateGradient(True) + rooster.SetUseCudaCyclicDeformation(True) + rooster.SetForwardProjectionFilter( + ROOSTERFilterType.ForwardProjectionType_FP_CUDARAYCAST + ) + rooster.SetBackProjectionFilter( + ROOSTERFilterType.BackProjectionType_BP_CUDAVOXELBASED + ) +else: + rooster.SetForwardProjectionFilter( + ROOSTERFilterType.ForwardProjectionType_FP_JOSEPH + ) + rooster.SetBackProjectionFilter(ROOSTERFilterType.BackProjectionType_BP_VOXELBASED) + +# Set the newly ordered arguments +rooster.SetInputProjectionStack(reorder.GetOutput()) +rooster.SetGeometry(reorder.GetOutputGeometry()) +rooster.SetWeights(signal_to_interpolation_weights.GetOutput()) +rooster.SetSignal(reorder.GetOutputSignal()) + +# For each optional regularization step, set whether or not +# it should be performed, and provide the necessary inputs + +# Positivity +rooster.SetPerformPositivity(True) + +# Motion mask +rooster.SetPerformMotionMask(False) + +# Spatial TV +rooster.SetGammaTVSpace(0.1) +rooster.SetTV_iterations(4) +rooster.SetPerformTVSpatialDenoising(True) + +# Spatial wavelets +rooster.SetPerformWaveletsSpatialDenoising(False) + +# Temporal TV +rooster.SetGammaTVTime(0.1) +rooster.SetTV_iterations(4) +rooster.SetPerformTVTemporalDenoising(True) + +# Temporal L0 +rooster.SetPerformL0TemporalDenoising(False) + +# Total nuclear variation +rooster.SetPerformTNVDenoising(False) + +# Warping +rooster.SetPerformWarping(True) +rooster.SetUseNearestNeighborInterpolationInWarping(False) + +dvf_reader = itk.ImageFileReader[CpuDVFSequenceImageType].New() +dvf_reader.SetFileName("four_d_dvf.mha") +dvf = dvf_reader.GetOutput() +if hasattr(itk, "CudaImage"): + dvf = itk.cuda_image_from_image(dvf) +rooster.SetDisplacementField(dvf) + +rooster.SetComputeInverseWarpingByConjugateGradient(False) +idvf_reader = itk.ImageFileReader[CpuDVFSequenceImageType].New() +idvf_reader.SetFileName("four_d_idvf.mha") +idvf = idvf_reader.GetOutput() +if hasattr(itk, "CudaImage"): + idvf = itk.cuda_image_from_image(idvf) +rooster.SetInverseDisplacementField(idvf) + + +class VerboseIterationCommand: + def __init__(self): + self.count = 0 + + def callback(self): + self.count += 1 + print(f"Iteration {self.count}", end="\r") + + +class VerboseEndCommand: + def callback(self): + print("") + + +cmd = VerboseIterationCommand() +rooster.AddObserver(itk.IterationEvent(), cmd.callback) +cmd_end = VerboseEndCommand() +rooster.AddObserver(itk.EndEvent(), cmd_end.callback) + +rooster.Update() + +# Write +writer = itk.ImageFileWriter[CpuVolumeSeriesType].New() +writer.SetFileName("fourdrooster.mha") +writer.SetInput(rooster.GetOutput()) +writer.Update() diff --git a/examples/FourDROOSTER/README.md b/examples/FourDROOSTER/README.md index 477d17f1e..6f5044d65 100644 --- a/examples/FourDROOSTER/README.md +++ b/examples/FourDROOSTER/README.md @@ -12,4 +12,11 @@ It reads its data from disk. The data can be generated by the [*Generate 4D data ``` ```` +````{tab-item} Python + +```{literalinclude} ./FourDROOSTER.py +:language: python +``` +```` + ````` diff --git a/examples/GenerateFourDData/GenerateFourDData.py b/examples/GenerateFourDData/GenerateFourDData.py new file mode 100644 index 000000000..c92682d51 --- /dev/null +++ b/examples/GenerateFourDData/GenerateFourDData.py @@ -0,0 +1,227 @@ +import itk +from itk import RTK as rtk + +PixelType = itk.F +DVFVectorType = itk.CovariantVector[PixelType, 3] +VolumeSeriesType = itk.Image[PixelType, 4] +ProjectionStackType = itk.Image[PixelType, 3] +VolumeType = itk.Image[PixelType, 3] +DVFSequenceImageType = itk.Image[DVFVectorType, 4] + +NumberOfProjectionImages = 64 +NumberOfFrames = 8 + +# /* ========================= +# * Constant volume source +# * ========================= */ +ConstantVolumeSourceType = rtk.ConstantImageSource[VolumeType] +volumeSource = ConstantVolumeSourceType.New() + +origin = [-63.0, -31.0, -63.0] +spacing = [4.0, 4.0, 4.0] +size = [32, 16, 32] + +volumeSource.SetOrigin(origin) +volumeSource.SetSpacing(spacing) +volumeSource.SetSize(size) +volumeSource.Update() + +# /* ========================= +# * Projection accumulation +# * ========================= */ +ConstantProjectionSourceType = rtk.ConstantImageSource[ProjectionStackType] +projectionsSource = ConstantProjectionSourceType.New() + +projOrigin = [-254.0, -254.0, -254.0] +projSpacing = [8.0, 8.0, 1.0] +projSize = [64, 64, NumberOfProjectionImages] + +projectionsSource.SetOrigin(projOrigin) +projectionsSource.SetSpacing(projSpacing) +projectionsSource.SetSize(projSize) +projectionsSource.Update() + +oneProjectionSource = ConstantProjectionSourceType.New() +oneProjSize = [projSize[0], projSize[1], 1] +oneProjectionSource.SetOrigin(projOrigin) +oneProjectionSource.SetSpacing(projSpacing) +oneProjectionSource.SetSize(oneProjSize) + +REIType = rtk.RayEllipsoidIntersectionImageFilter[VolumeType, ProjectionStackType] +PasteType = itk.PasteImageFilter[ProjectionStackType] + +# /* Allocate explicit accumulator image */ +accumulated = ProjectionStackType.New() +accumulated.SetOrigin(projectionsSource.GetOutput().GetOrigin()) +accumulated.SetSpacing(projectionsSource.GetOutput().GetSpacing()) +accumulated.SetDirection(projectionsSource.GetOutput().GetDirection()) +accumulated.SetRegions(projectionsSource.GetOutput().GetLargestPossibleRegion()) +accumulated.Allocate() +accumulated.FillBuffer(0.0) + +paste = PasteType.New() +paste.SetDestinationImage(accumulated) + +destIndex = [0, 0, 0] + +# Create the signal file and the geometry, to be filled in the for loop +signalFileName = "four_d_signal.txt" +geometry = rtk.ThreeDCircularProjectionGeometry.New() + +print("Generating projections, geometry, and signal...") +with open(signalFileName, "w", encoding="utf-8") as signalFile: + for i in range(NumberOfProjectionImages): + if i % 25 == 0: + print(f" projection {i + 1}/{NumberOfProjectionImages}") + + geometry.AddProjection( + 600.0, + 1200.0, + i * 360.0 / NumberOfProjectionImages, + 0, + 0, + 0, + 0, + 20, + 15, + ) + + geom = rtk.ThreeDCircularProjectionGeometry.New() + geom.AddProjection( + 600.0, + 1200.0, + i * 360.0 / NumberOfProjectionImages, + 0, + 0, + 0, + 0, + 20, + 15, + ) + + e1 = REIType.New() + e1.SetInput(oneProjectionSource.GetOutput()) + e1.SetGeometry(geom) + e1.SetDensity(2.0) + e1.SetAxis([60.0, 30.0, 60.0]) + e1.SetCenter([0.0, 0.0, 0.0]) + e1.InPlaceOff() + e1.Update() + + e2 = REIType.New() + e2.SetInput(e1.GetOutput()) + e2.SetGeometry(geom) + e2.SetDensity(-1.0) + e2.SetAxis([8.0, 8.0, 8.0]) + center = [4 * (abs((4 + i) % 8 - 4.0) - 2.0), 0.0, 0.0] + e2.SetCenter(center) + e2.InPlaceOff() + e2.Update() + + paste.SetSourceImage(e2.GetOutput()) + paste.SetSourceRegion(e2.GetOutput().GetLargestPossibleRegion()) + paste.SetDestinationIndex(destIndex) + paste.Update() + + accumulated = paste.GetOutput() + accumulated.DisconnectPipeline() + paste.SetDestinationImage(accumulated) + + destIndex[2] += 1 + + signalFile.write(f"{(i % 8) / 8.0}\n") + +itk.imwrite(accumulated, "four_d_projections.mha") +rtk.write_geometry(geometry, "four_d_geometry.xml") + +# /* ========================= +# * DVF & inverse DVF +# * ========================= */ +print("Generating DVF and inverse DVF...") + +fourDOrigin = [-63.0, -31.0, -63.0, 0.0] +fourDSpacing = [4.0, 4.0, 4.0, 1.0] +fourDSize = [32, 16, 32, NumberOfFrames] + +dvf = DVFSequenceImageType.New() +idvf = DVFSequenceImageType.New() + +region = itk.ImageRegion[4]() +dvfSize = [fourDSize[0], fourDSize[1], fourDSize[2], 2] +region.SetSize(dvfSize) + +dvf.SetRegions(region) +dvf.SetOrigin(fourDOrigin) +dvf.SetSpacing(fourDSpacing) +dvf.Allocate() + +idvf.SetRegions(region) +idvf.SetOrigin(fourDOrigin) +idvf.SetSpacing(fourDSpacing) +idvf.Allocate() + +centerIndex = [0, 0, 0, 0] +centerIndex[0] = dvfSize[0] // 2 +centerIndex[1] = dvfSize[1] // 2 +centerIndex[2] = dvfSize[2] // 2 + +for t in range(dvfSize[3]): + for z in range(dvfSize[2]): + for y in range(dvfSize[1]): + for x in range(dvfSize[0]): + v = DVFVectorType() + v.Fill(0.0) + + d0 = x - centerIndex[0] + d1 = y - centerIndex[1] + d2 = z - centerIndex[2] + if 0.3 * d0 * d0 + d1 * d1 + d2 * d2 < 40: + v[0] = -8.0 if t == 0 else 8.0 + + iv = DVFVectorType() + iv.Fill(0.0) + iv[0] = -v[0] + + index = [x, y, z, t] + dvf.SetPixel(index, v) + idvf.SetPixel(index, iv) + +itk.imwrite(dvf, "four_d_dvf.mha") +itk.imwrite(idvf, "four_d_idvf.mha") + +# /* ========================= +# * Ground truth +# * ========================= */ +print("Generating ground truth...") + +join = itk.JoinSeriesImageFilter[VolumeType, VolumeSeriesType].New() + +for t in range(fourDSize[3]): + DEType = rtk.DrawEllipsoidImageFilter[VolumeType, VolumeType] + + de1 = DEType.New() + de1.SetInput(volumeSource.GetOutput()) + de1.SetDensity(2.0) + de1.SetAxis([60.0, 30.0, 60.0]) + de1.SetCenter([0.0, 0.0, 0.0]) + de1.InPlaceOff() + de1.Update() + + de2 = DEType.New() + de2.SetInput(de1.GetOutput()) + de2.SetDensity(-1.0) + de2.SetAxis([8.0, 8.0, 8.0]) + de2.SetCenter([4 * (abs((4 + t) % 8 - 4.0) - 2.0), 0.0, 0.0]) + de2.InPlaceOff() + de2.Update() + + duplicator = itk.ImageDuplicator[VolumeType].New() + duplicator.SetInputImage(de2.GetOutput()) + duplicator.Update() + + join.SetInput(t, duplicator.GetOutput()) + +join.Update() +itk.imwrite(join.GetOutput(), "four_d_ground_truth.mha") + +print("Done.") diff --git a/examples/GenerateFourDData/README.md b/examples/GenerateFourDData/README.md index 1cc19be41..eeed217fa 100644 --- a/examples/GenerateFourDData/README.md +++ b/examples/GenerateFourDData/README.md @@ -18,4 +18,11 @@ You can also skip this part and [download the data](https://data.kitware.com/#co :language: c++ ``` ```` + +````{tab-item} Python + +```{literalinclude} ./GenerateFourDData.py +:language: python +``` +```` ````` diff --git a/wrapping/CMakeLists.txt b/wrapping/CMakeLists.txt index 332b24b4a..9a25c21c5 100644 --- a/wrapping/CMakeLists.txt +++ b/wrapping/CMakeLists.txt @@ -13,6 +13,7 @@ set( itkImageRTK itkCudaImageRTK itkVectorImageRTK + itkImageDuplicatorRTK itkImageSourceRTK itkImageToImageFilterRTK itkInPlaceImageFilterRTK diff --git a/wrapping/itkImageRegionRTK.wrap b/wrapping/itkImageRegionRTK.wrap index c424fa00e..b035b0e6c 100644 --- a/wrapping/itkImageRegionRTK.wrap +++ b/wrapping/itkImageRegionRTK.wrap @@ -4,3 +4,10 @@ if (${_index} EQUAL -1) itk_wrap_template(1 1) itk_end_wrap_class() endif() + +list(FIND ITK_WRAP_IMAGE_DIMS "4" _index) +if (${_index} EQUAL -1) + itk_wrap_class("itk::ImageRegion") + itk_wrap_template(4 4) + itk_end_wrap_class() +endif() diff --git a/wrapping/itkImageSourceRTK.wrap b/wrapping/itkImageSourceRTK.wrap index 8f40b0f5c..8749aef4b 100644 --- a/wrapping/itkImageSourceRTK.wrap +++ b/wrapping/itkImageSourceRTK.wrap @@ -38,6 +38,7 @@ itk_wrap_class("itk::ImageSource" POINTER) list(FIND ITK_WRAP_IMAGE_DIMS "4" _index) if (${_index} EQUAL -1) itk_wrap_template("I${ITKM_F}4" "itk::Image<${ITKT_F}, 4>") + itk_wrap_template("I${ITKM_CVF3}4" "itk::Image<${ITKT_CVF3}, 4>") if (ITK_WRAP_double) itk_wrap_template("I${ITKM_UD}4" "itk::Image<${ITKT_D}, 4>") endif() diff --git a/wrapping/rtkConstantImageSource.wrap b/wrapping/rtkConstantImageSource.wrap index f59c10894..e2e7080a0 100644 --- a/wrapping/rtkConstantImageSource.wrap +++ b/wrapping/rtkConstantImageSource.wrap @@ -6,6 +6,7 @@ itk_wrap_class("rtk::ConstantImageSource" POINTER) itk_wrap_image_filter("${WRAP_ITK_REAL}" 1) itk_wrap_image_filter("${WRAP_ITK_VECTOR_REAL}" 1) + itk_wrap_template("I${ITKM_CVF3}4" "itk::Image<${ITKT_CVF3}, 4>") if(RTK_USE_CUDA) foreach(d ${ITK_WRAP_IMAGE_DIMS}) From 5f8c92856c67a9c83ea013e38c631a08a09f5012 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 23 Jun 2026 15:27:02 +0200 Subject: [PATCH 3/3] Doc: Generate 4D example illustration data The illustration data were generated with 128 projections on a 512 x 512 x 128 projection stack with 1 mm detector spacing. The illustration reconstructions were generated on a 256 x 128 x 256 x 8 grid with 0.5 mm spatial spacing. --- .../docs/ExternalData/FourDConjugateGradient.gif.sha512 | 1 + documentation/docs/ExternalData/FourDFDK.gif.sha512 | 1 + documentation/docs/ExternalData/FourDRooster.gif.sha512 | 1 + examples/FourDConjugateGradient/README.md | 8 ++++++++ examples/FourDFDK/README.md | 9 +++++++++ examples/FourDROOSTER/README.md | 2 ++ 6 files changed, 22 insertions(+) create mode 100644 documentation/docs/ExternalData/FourDConjugateGradient.gif.sha512 create mode 100644 documentation/docs/ExternalData/FourDFDK.gif.sha512 create mode 100644 documentation/docs/ExternalData/FourDRooster.gif.sha512 diff --git a/documentation/docs/ExternalData/FourDConjugateGradient.gif.sha512 b/documentation/docs/ExternalData/FourDConjugateGradient.gif.sha512 new file mode 100644 index 000000000..b48691500 --- /dev/null +++ b/documentation/docs/ExternalData/FourDConjugateGradient.gif.sha512 @@ -0,0 +1 @@ +430dd919f1ad889d78a04cc70c32f1851b790f0423d6ee63ae3a9c626786a3fe9617a38b5655cecfb1890f56e6ffe298a01fe2d39210a7f7a16e539074b3f259 diff --git a/documentation/docs/ExternalData/FourDFDK.gif.sha512 b/documentation/docs/ExternalData/FourDFDK.gif.sha512 new file mode 100644 index 000000000..991cfb6ee --- /dev/null +++ b/documentation/docs/ExternalData/FourDFDK.gif.sha512 @@ -0,0 +1 @@ +6db1332af81b0ac92d9fbbf3ad42ce8d3f7c839923e979cd4aa7bb1285e7b1f482221d77c7983547020317ab652b1b05d998916f0952f96c741b58d37a0b2c80 diff --git a/documentation/docs/ExternalData/FourDRooster.gif.sha512 b/documentation/docs/ExternalData/FourDRooster.gif.sha512 new file mode 100644 index 000000000..991cfb6ee --- /dev/null +++ b/documentation/docs/ExternalData/FourDRooster.gif.sha512 @@ -0,0 +1 @@ +6db1332af81b0ac92d9fbbf3ad42ce8d3f7c839923e979cd4aa7bb1285e7b1f482221d77c7983547020317ab652b1b05d998916f0952f96c741b58d37a0b2c80 diff --git a/examples/FourDConjugateGradient/README.md b/examples/FourDConjugateGradient/README.md index 77c1a644d..279a3a8c6 100644 --- a/examples/FourDConjugateGradient/README.md +++ b/examples/FourDConjugateGradient/README.md @@ -3,6 +3,8 @@ This example shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources. It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). +![img_4D](../../documentation/docs/ExternalData/FourDConjugateGradient.gif){w=200px alt="conjugate gradient 4D reconstruction"} + `````{tab-set} ````{tab-item} C++ @@ -11,5 +13,11 @@ It reads its data from disk. The data can be generated by the [*Generate 4D data :language: c++ ``` ```` +````{tab-item} Python + +```{literalinclude} ./FourDConjugateGradient.py +:language: python +``` +```` ````` diff --git a/examples/FourDFDK/README.md b/examples/FourDFDK/README.md index 0b06ab133..5112aa8a7 100644 --- a/examples/FourDFDK/README.md +++ b/examples/FourDFDK/README.md @@ -3,6 +3,8 @@ This example shows how to perform filtered backprojection cone-beam CT reconstruction using either CPU or GPU resources depending on how RTK was compiled. It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). +![img_4D](../../documentation/docs/ExternalData/FourDFDK.gif){w=200px alt="Fdk 4D reconstruction"} + `````{tab-set} ````{tab-item} C++ @@ -12,4 +14,11 @@ It reads its data from disk. The data can be generated by the [*Generate 4D data ``` ```` +````{tab-item} Python + +```{literalinclude} ./FourDFDK.py +:language: python +``` +```` + ````` diff --git a/examples/FourDROOSTER/README.md b/examples/FourDROOSTER/README.md index 6f5044d65..7b4e0554e 100644 --- a/examples/FourDROOSTER/README.md +++ b/examples/FourDROOSTER/README.md @@ -3,6 +3,8 @@ This example shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources. It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). +![img_3D](../../documentation/docs/ExternalData/FourDRooster.gif){w=200px alt="Rooster 4D reconstruction"} + `````{tab-set} ````{tab-item} C++