From 8ae01f96258c90bb43991709121a0cf9a762d605 Mon Sep 17 00:00:00 2001 From: Cyril Mory Date: Thu, 8 Jan 2026 13:44:40 +0100 Subject: [PATCH 1/7] 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 | 12 + 17 files changed, 751 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..6bac805f9 --- /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, 8)); + fourDSource->SetConstant(0.); + fourDSource->Update(); + + // Go over each frame, reconstruct 3D frame and paste with iterators in 4D image + for (int f = 0; f < 8; f++) + { + selector->SetPhase(f / (double)8); + 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 18da9b27c..4fb498491 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -409,6 +409,18 @@ rtk_add_test(rtkConjugateGradientExampleTest ../examples/ConjugateGradient/ConjugateGradient.cxx DATA{Input/Forbild/Thorax} ) +rtk_add_test(rtkGenerateFourDDataExampleTest + ../examples/GenerateFourDData/GenerateFourDData.cxx +) +rtk_add_test(rtkFourDConjugateGradientExampleTest + ../examples/FourDConjugateGradient/FourDConjugateGradient.cxx +) +rtk_add_test(rtkFourDFDKExampleTest + ../examples/FourDFDK/FourDFDK.cxx +) +rtk_add_test(rtkFourDROOSTERExampleTest + ../examples/FourDROOSTER/FourDROOSTER.cxx +) #----------------------------------------------------------------------------- # Test drivers From 041a95288985b0a3118df5f808d10a5d3379d349 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 10 Jun 2026 11:42:48 +0200 Subject: [PATCH 2/7] ENH: Add 4D python reconstruction examples --- .../FourDConjugateGradient.py | 114 +++++++++ examples/FourDFDK/FourDFDK.py | 75 ++++++ 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, 603 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..c3a924480 --- /dev/null +++ b/examples/FourDFDK/FourDFDK.py @@ -0,0 +1,75 @@ +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") + +projectionsReader = rtk.ProjectionsReader[CPUOutputImageType].New() +fileNames = ["four_d_projections.mha"] +projectionsReader.SetFileNames(fileNames) +projectionsReader.Update() + +# Part specific to 4D +selector = rtk.SelectOneProjectionPerCycleImageFilter[CPUOutputImageType].New() +selector.SetInput(projectionsReader.GetOutput()) +selector.SetInputGeometry(geometry) +selector.SetSignalFilename("four_d_signal.txt") + +# Create one frame of the reconstructed image +constantImageSource = rtk.ConstantImageSource[OutputImageType].New() + +constantImageSource.SetOrigin([-63.0, -31.0, -63.0]) +constantImageSource.SetSpacing([4.0, 4.0, 4.0]) +constantImageSource.SetSize([32, 16, 32]) +constantImageSource.SetConstant(0.0) +constantImageSource.Update() + +# FDK reconstruction filtering +if hasattr(itk, "CudaImage"): + cuda_projections = itk.CudaImageFromImageFilter[CPUOutputImageType].New() + cuda_projections.SetInput(selector.GetOutput()) + + feldkamp = rtk.CudaFDKConeBeamReconstructionFilter.New() + feldkamp.SetInput(1, cuda_projections.GetOutput()) +else: + feldkamp = rtk.FDKConeBeamReconstructionFilter[ + OutputImageType, OutputImageType, OutputPixelType + ].New() + feldkamp.SetInput(1, selector.GetOutput()) + +feldkamp.SetInput(0, constantImageSource.GetOutput()) +feldkamp.SetGeometry(selector.GetOutputGeometry()) + +# Create empty 4D image +FourDOutputImageType = itk.Image[OutputPixelType, Dimension + 1] +fourDSource = rtk.ConstantImageSource[FourDOutputImageType].New() +fourDSource.SetOrigin([-63.0, -31.0, -63.0, 0.0]) +fourDSource.SetSpacing([4.0, 4.0, 4.0, 1.0]) +fourDSource.SetSize([32, 16, 32, 8]) +fourDSource.SetConstant(0.0) +fourDSource.Update() + +# Go over each frame, reconstruct 3D frame and paste with iterators in 4D image +it4D = itk.GetArrayViewFromImage(fourDSource.GetOutput()) +for f in range(8): + selector.SetPhase(f / 8.0) + feldkamp.UpdateLargestPossibleRegion() + + it3D = itk.GetArrayFromImage(feldkamp.GetOutput()) + it4D[f] = it3D + +# Write +writer = itk.ImageFileWriter[FourDOutputImageType].New() +writer.SetFileName("fourdfdk.mha") +writer.SetInput(fourDSource.GetOutput()) +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 ccf216c01c805119018dd09a5482b5c245088ad1 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 23 Jun 2026 15:27:02 +0200 Subject: [PATCH 3/7] 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..6d034c88c --- /dev/null +++ b/documentation/docs/ExternalData/FourDRooster.gif.sha512 @@ -0,0 +1 @@ +efb1e489c9eb1851438ccb3899e1b70fe5632104bd765fe57e974ecdb343a6dc00151f813f0a4d7aff02f2152ee62f0d24e07791aa10af7efa3e21610d8ad334 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++ From e30cbacffde4ae0fc75bf4306d5ec25e1b080e02 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 9 Jul 2026 11:31:29 +0200 Subject: [PATCH 4/7] ENH: Improve testing time Reduce the workload of long-running tests while keeping their existing coverage targets and quality checks. This lowers projection counts, sampled geometries, reconstruction sizes, iteration counts, and streaming splits in selected expensive tests. It also shortens the inline reconstruction example acquisition delay so the example test no longer spends most of its time sleeping. --- .../InlineReconstruction/InlineReconstruction.cxx | 8 ++++---- .../InlineReconstruction/InlineReconstruction.py | 4 ++-- test/rtkTestReg23ProjectionGeometry.cxx | 2 +- test/rtkadjointoperatorstest.cxx | 14 +++++++------- test/rtkadmmwaveletstest.cxx | 2 +- test/rtkcylindricaldetectorreconstructiontest.cxx | 4 ++-- test/rtkfdktest.cxx | 4 ++-- test/rtkforwardattenuatedprojectiontest.cxx | 4 ++-- test/rtkfourdroostertest.cxx | 14 +++++++------- test/rtkfourdsarttest.cxx | 10 +++++----- test/rtkgeometryclonetest.cxx | 14 +++++++------- test/rtkgeometryfrommatrixtest.cxx | 14 +++++++------- test/rtkosemtest.cxx | 8 ++++---- test/rtkrampfiltertest.cxx | 2 +- test/rtkregularizedconjugategradienttest.cxx | 10 +++++----- test/rtksarttest.cxx | 2 +- test/rtkspectralonesteptest.cxx | 6 +++--- test/rtkvariancereconstructiontest.cxx | 4 ++-- test/rtkzengforwardprojectiontest.cxx | 10 +++++++--- 19 files changed, 70 insertions(+), 66 deletions(-) diff --git a/examples/InlineReconstruction/InlineReconstruction.cxx b/examples/InlineReconstruction/InlineReconstruction.cxx index 2f879d37c..7aeb00f5a 100644 --- a/examples/InlineReconstruction/InlineReconstruction.cxx +++ b/examples/InlineReconstruction/InlineReconstruction.cxx @@ -26,19 +26,19 @@ using ImageType = itk::Image; // Constants for acquisition and reconstruction -const int numProjections = 64; +const int numProjections = 32; const double sid = 1000.0; const double sdd = 1500.0; const double arc = 200.0; const double spacing = 1.0; -const int imageSize = 256; +const int imageSize = 128; const double origin = -0.5 * (imageSize - 1); -const std::chrono::milliseconds acquisitionSleepDuration(100); +const std::chrono::milliseconds acquisitionSleepDuration(20); std::atomic AcquiredProjectionsCount(0); // Function for the thread simulating an acquisition -// Simulates and writes a projection every 100 ms. +// Simulates and writes a projection every 20 ms. void Acquisition() { diff --git a/examples/InlineReconstruction/InlineReconstruction.py b/examples/InlineReconstruction/InlineReconstruction.py index 3f1e8dcd6..8320a4f4c 100644 --- a/examples/InlineReconstruction/InlineReconstruction.py +++ b/examples/InlineReconstruction/InlineReconstruction.py @@ -21,7 +21,7 @@ # Function for the thread simulating an acquisition: simulates and writes a -# projection every 100 ms. +# projection every 20 ms. def acquisition(): global number_of_acquired_projections for i in range(nproj): @@ -44,7 +44,7 @@ def acquisition(): # is the case with CPython # https://docs.python.org/3/library/threading.html number_of_acquired_projections += 1 - time.sleep(0.1) + time.sleep(0.02) # Launches the acquisition thread, then waits for new projections and reconstructs inline / on-the-fly diff --git a/test/rtkTestReg23ProjectionGeometry.cxx b/test/rtkTestReg23ProjectionGeometry.cxx index 96badf637..bf5182c29 100644 --- a/test/rtkTestReg23ProjectionGeometry.cxx +++ b/test/rtkTestReg23ProjectionGeometry.cxx @@ -203,7 +203,7 @@ rtkTestReg23ProjectionGeometry(int argc, char * argv[]) gantryAngles.insert(gantryAngles.end(), criticalAngles.begin(), criticalAngles.end()); outOfPlaneAngles.insert(outOfPlaneAngles.end(), criticalAngles.begin(), criticalAngles.end()); // randomly sample further angles in [-180;+180] deg range: - constexpr int NUM_RANDOM_ANGLES = 100; + constexpr int NUM_RANDOM_ANGLES = 10; for (int u = 0; u < NUM_RANDOM_ANGLES; u++) { inPlaneAngles.push_back(generator->GetUniformVariate(-180, 180)); diff --git a/test/rtkadjointoperatorstest.cxx b/test/rtkadjointoperatorstest.cxx index ce1af3eb0..1bc4be60b 100644 --- a/test/rtkadjointoperatorstest.cxx +++ b/test/rtkadjointoperatorstest.cxx @@ -43,7 +43,7 @@ rtkadjointoperatorstest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 60; + constexpr unsigned int NumberOfProjectionImages = 14; #endif // Random image sources @@ -58,13 +58,13 @@ rtkadjointoperatorstest(int, char *[]) auto constantAttenuationSource = ConstantImageSourceType::New(); // Volume metadata - auto origin = itk::MakePoint(-128., -128., -128.); + auto origin = itk::MakePoint(-64., -64., -64.); #if FAST_TESTS_NO_CHECKS auto spacing = itk::MakeVector(252., 252., 252.); auto size = itk::MakeSize(2, 2, 2); #else auto spacing = itk::MakeVector(4., 4., 4.); - auto size = itk::MakeSize(65, 65, 65); + auto size = itk::MakeSize(32, 32, 32); #endif randomVolumeSource->SetOrigin(origin); randomVolumeSource->SetSpacing(spacing); @@ -81,16 +81,16 @@ rtkadjointoperatorstest(int, char *[]) constantAttenuationSource->SetOrigin(origin); constantAttenuationSource->SetSpacing(spacing); constantAttenuationSource->SetSize(size); - constantAttenuationSource->SetConstant(0.0154); + constantAttenuationSource->SetConstant(0.002); // Projections metadata - origin = itk::MakePoint(-128., -128., -128.); + origin = itk::MakePoint(-64., -64., -64.); #if FAST_TESTS_NO_CHECKS - spacing = itk::MakeVector(504., 504., 504.); + spacing = itk::MakeVector(126., 126., 126.); size = itk::MakeSize(2, 2, NumberOfProjectionImages); #else spacing = itk::MakeVector(4., 4., 1.); - size = itk::MakeSize(65, 65, NumberOfProjectionImages); + size = itk::MakeSize(32, 32, NumberOfProjectionImages); #endif randomProjectionsSource->SetOrigin(origin); randomProjectionsSource->SetSpacing(spacing); diff --git a/test/rtkadmmwaveletstest.cxx b/test/rtkadmmwaveletstest.cxx index 435a87151..7a1732e72 100644 --- a/test/rtkadmmwaveletstest.cxx +++ b/test/rtkadmmwaveletstest.cxx @@ -94,7 +94,7 @@ rtkadmmwaveletstest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 180; + constexpr unsigned int NumberOfProjectionImages = 60; #endif diff --git a/test/rtkcylindricaldetectorreconstructiontest.cxx b/test/rtkcylindricaldetectorreconstructiontest.cxx index 4ee752ffb..9e1ab8d34 100644 --- a/test/rtkcylindricaldetectorreconstructiontest.cxx +++ b/test/rtkcylindricaldetectorreconstructiontest.cxx @@ -37,7 +37,7 @@ rtkcylindricaldetectorreconstructiontest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 180; + constexpr unsigned int NumberOfProjectionImages = 40; #endif // Constant image sources @@ -110,7 +110,7 @@ rtkcylindricaldetectorreconstructiontest(int, char *[]) conjugategradient->SetInput(1, rei->GetOutput()); conjugategradient->SetInputWeights(uniformWeightsSource->GetOutput()); conjugategradient->SetGeometry(geometry); - conjugategradient->SetNumberOfIterations(5); + conjugategradient->SetNumberOfIterations(3); conjugategradient->SetDisableDisplacedDetectorFilter(true); std::cout << "\n\n****** Case 1: Joseph forward and back projectors ******" << std::endl; diff --git a/test/rtkfdktest.cxx b/test/rtkfdktest.cxx index ca7214b06..b1d8197fd 100644 --- a/test/rtkfdktest.cxx +++ b/test/rtkfdktest.cxx @@ -41,7 +41,7 @@ rtkfdktest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 180; + constexpr unsigned int NumberOfProjectionImages = 120; #endif // Constant image sources @@ -164,7 +164,7 @@ rtkfdktest(int, char *[]) auto streamer = itk::StreamingImageFilter::New(); streamer->SetInput(0, fov->GetOutput()); - streamer->SetNumberOfStreamDivisions(8); + streamer->SetNumberOfStreamDivisions(2); auto splitter = itk::ImageRegionSplitterDirection::New(); splitter->SetDirection(2); // Prevent splitting along z axis. As a result, splitting will be performed along y axis streamer->SetRegionSplitter(splitter); diff --git a/test/rtkforwardattenuatedprojectiontest.cxx b/test/rtkforwardattenuatedprojectiontest.cxx index 6a88338bd..bfe771420 100644 --- a/test/rtkforwardattenuatedprojectiontest.cxx +++ b/test/rtkforwardattenuatedprojectiontest.cxx @@ -45,7 +45,7 @@ rtkforwardattenuatedprojectiontest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 45; + constexpr unsigned int NumberOfProjectionImages = 24; #endif // Constant image sources @@ -191,7 +191,7 @@ rtkforwardattenuatedprojectiontest(int, char *[]) auto stream = itk::StreamingImageFilter::New(); stream->SetInput(jfp->GetOutput()); - stream->SetNumberOfStreamDivisions(9); + stream->SetNumberOfStreamDivisions(4); auto splitter = itk::ImageRegionSplitterDirection::New(); splitter->SetDirection(2); // Splitting along direction 1, NOT 2 stream->SetRegionSplitter(splitter); diff --git a/test/rtkfourdroostertest.cxx b/test/rtkfourdroostertest.cxx index fb7f57208..e3070553c 100644 --- a/test/rtkfourdroostertest.cxx +++ b/test/rtkfourdroostertest.cxx @@ -47,7 +47,7 @@ rtkfourdroostertest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 5; #else - constexpr unsigned int NumberOfProjectionImages = 64; + constexpr unsigned int NumberOfProjectionImages = 24; #endif // Constant image sources @@ -72,7 +72,7 @@ rtkfourdroostertest(int, char *[]) auto fourDSize = itk::MakeSize(8, 8, 8, 2); auto fourDSpacing = itk::MakeVector(16., 8., 16., 1.); #else - auto fourDSize = itk::MakeSize(32, 16, 32, 8); + auto fourDSize = itk::MakeSize(32, 16, 32, 4); auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.); #endif fourdSource->SetOrigin(fourDOrigin); @@ -142,7 +142,7 @@ rtkfourdroostertest(int, char *[]) // Ellipse 2 auto e2 = REIType::New(); semiprincipalaxis.Fill(8.); - center[0] = 4 * (std::abs((4 + noProj) % 8 - 4.) - 2.); + center[0] = 4 * (std::abs((fourDSize[3] / 2 + noProj) % fourDSize[3] - fourDSize[3] / 2.) - 2.); center[1] = 0.; center[2] = 0.; e2->SetInput(e1->GetOutput()); @@ -167,7 +167,7 @@ rtkfourdroostertest(int, char *[]) destinationIndex[2]++; // Signal - signalFile << (noProj % 8) / 8. << std::endl; + signalFile << (noProj % fourDSize[3]) / static_cast(fourDSize[3]) << std::endl; } signalFile.close(); @@ -251,7 +251,7 @@ rtkfourdroostertest(int, char *[]) axis2.Fill(8.); de2->SetAxis(axis2); DEType::PointType center2; - center2[0] = 4 * (std::abs((4 + n) % 8 - 4.) - 2.); + center2[0] = 4 * (std::abs((fourDSize[3] / 2 + n) % fourDSize[3] - fourDSize[3] / 2.) - 2.); center2[1] = 0.; center2[2] = 0.; de2->SetCenter(center2); @@ -298,7 +298,7 @@ rtkfourdroostertest(int, char *[]) rooster->SetCG_iterations(2); rooster->SetMainLoop_iterations(2); - rooster->SetTV_iterations(3); + rooster->SetTV_iterations(2); rooster->SetGammaTVSpace(1); rooster->SetGammaTVTime(0.1); @@ -307,7 +307,7 @@ rtkfourdroostertest(int, char *[]) rooster->SetNumberOfLevels(3); rooster->SetLambdaL0Time(0.1); - rooster->SetL0_iterations(5); + rooster->SetL0_iterations(3); std::cout << "\n\n****** Case 1: Joseph forward projector, voxel-based back projector, positivity, motion mask, " "wavelets spatial denoising, TV temporal denoising, no warping ******" diff --git a/test/rtkfourdsarttest.cxx b/test/rtkfourdsarttest.cxx index 16f3a8f8b..3c2f3c1cc 100644 --- a/test/rtkfourdsarttest.cxx +++ b/test/rtkfourdsarttest.cxx @@ -44,7 +44,7 @@ rtkfourdsarttest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 5; #else - constexpr unsigned int NumberOfProjectionImages = 64; + constexpr unsigned int NumberOfProjectionImages = 32; #endif // Constant image sources @@ -100,7 +100,7 @@ rtkfourdsarttest(int, char *[]) fourDSize[0] = 32; fourDSize[1] = 16; fourDSize[2] = 32; - fourDSize[3] = 8; + fourDSize[3] = 4; fourDSpacing[0] = 4.; fourDSpacing[1] = 4.; fourDSpacing[2] = 4.; @@ -186,7 +186,7 @@ rtkfourdsarttest(int, char *[]) // Ellipse 2 auto e2 = REIType::New(); semiprincipalaxis.Fill(8.); - center[0] = 4 * (std::abs((4 + noProj) % 8 - 4.) - 2.); + center[0] = 4 * (std::abs((fourDSize[3] / 2 + noProj) % fourDSize[3] - fourDSize[3] / 2.) - 2.); center[1] = 0.; center[2] = 0.; e2->SetInput(e1->GetOutput()); @@ -211,7 +211,7 @@ rtkfourdsarttest(int, char *[]) destinationIndex[2]++; // Signal - signalFile << (noProj % 8) / 8. << std::endl; + signalFile << (noProj % fourDSize[3]) / static_cast(fourDSize[3]) << std::endl; } signalFile.close(); @@ -244,7 +244,7 @@ rtkfourdsarttest(int, char *[]) axis2.Fill(8.); de2->SetAxis(axis2); DEType::PointType center2; - center2[0] = 4 * (std::abs((4 + n) % 8 - 4.) - 2.); + center2[0] = 4 * (std::abs((fourDSize[3] / 2 + n) % fourDSize[3] - fourDSize[3] / 2.) - 2.); center2[1] = 0.; center2[2] = 0.; de2->SetCenter(center2); diff --git a/test/rtkgeometryclonetest.cxx b/test/rtkgeometryclonetest.cxx index 3d9364b77..968ea1119 100644 --- a/test/rtkgeometryclonetest.cxx +++ b/test/rtkgeometryclonetest.cxx @@ -8,13 +8,13 @@ rtkgeometryclonetest(int, char *[]) rtk::ThreeDCircularProjectionGeometry::Pointer geometry1, geometry2; geometry1 = rtk::ThreeDCircularProjectionGeometry::New(); - for (double oa = -60.; oa < 65.; oa += 60.) - for (double ia = -60.; ia < 65.; ia += 60.) - for (double ga = 0.; ga < 360.; ga += 45.) - for (double sx = -50.; sx < 55.; sx += 25.) - for (double sy = -50.; sy < 55.; sy += 25.) - for (double px = -50.; px < 55.; px += 25.) - for (double py = -50.; py < 55.; py += 25.) + for (double oa = -60.; oa < 65.; oa += 120.) + for (double ia = -60.; ia < 65.; ia += 120.) + for (double ga = 0.; ga < 360.; ga += 90.) + for (double sx = -50.; sx < 55.; sx += 50.) + for (double sy = -50.; sy < 55.; sy += 50.) + for (double px = -50.; px < 55.; px += 50.) + for (double py = -50.; py < 55.; py += 50.) for (double sdd = 200.; sdd < 405.; sdd += 200.) for (double sid = 100.; sid < 200.; sid += 200.) { diff --git a/test/rtkgeometryfrommatrixtest.cxx b/test/rtkgeometryfrommatrixtest.cxx index dd325e6bd..c77e1ef79 100644 --- a/test/rtkgeometryfrommatrixtest.cxx +++ b/test/rtkgeometryfrommatrixtest.cxx @@ -9,13 +9,13 @@ rtkgeometryfrommatrixtest(int, char *[]) geometry1 = rtk::ThreeDCircularProjectionGeometry::New(); geometry2 = rtk::ThreeDCircularProjectionGeometry::New(); - for (double oa = -60.; oa < 65.; oa += 60.) - for (double ia = -60.; ia < 65.; ia += 60.) - for (double ga = 0.; ga < 360.; ga += 45.) - for (double sx = -50.; sx < 55.; sx += 25.) - for (double sy = -50.; sy < 55.; sy += 25.) - for (double px = -50.; px < 55.; px += 25.) - for (double py = -50.; py < 55.; py += 25.) + for (double oa = -60.; oa < 65.; oa += 120.) + for (double ia = -60.; ia < 65.; ia += 120.) + for (double ga = 0.; ga < 360.; ga += 90.) + for (double sx = -50.; sx < 55.; sx += 50.) + for (double sy = -50.; sy < 55.; sy += 50.) + for (double px = -50.; px < 55.; px += 50.) + for (double py = -50.; py < 55.; py += 50.) for (double sdd = 200.; sdd < 405.; sdd += 200.) for (double sid = 100.; sid < 200.; sid += 200.) { diff --git a/test/rtkosemtest.cxx b/test/rtkosemtest.cxx index 4e9dc29c0..291d0ea87 100644 --- a/test/rtkosemtest.cxx +++ b/test/rtkosemtest.cxx @@ -54,8 +54,8 @@ rtkosemtest(int, char *[]) auto size = itk::MakeSize(2, 2, 2); auto spacing = itk::MakeVector(252., 252., 252.); #else - auto size = itk::MakeSize(34, 34, 34); - auto spacing = itk::MakeVector(4., 4., 4.); + auto size = itk::MakeSize(24, 24, 24); + auto spacing = itk::MakeVector(6., 6., 6.); #endif tomographySource->SetOrigin(origin); tomographySource->SetSpacing(spacing); @@ -76,8 +76,8 @@ rtkosemtest(int, char *[]) size = itk::MakeSize(2, 2, NumberOfProjectionImages); spacing = itk::MakeVector(504., 504., 504.); #else - size = itk::MakeSize(34, 34, NumberOfProjectionImages); - spacing = itk::MakeVector(8., 8., 8.); + size = itk::MakeSize(24, 24, NumberOfProjectionImages); + spacing = itk::MakeVector(12., 12., 8.); #endif projectionsSource->SetOrigin(origin); projectionsSource->SetSpacing(spacing); diff --git a/test/rtkrampfiltertest.cxx b/test/rtkrampfiltertest.cxx index 313183794..e67938514 100644 --- a/test/rtkrampfiltertest.cxx +++ b/test/rtkrampfiltertest.cxx @@ -41,7 +41,7 @@ rtkrampfiltertest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 180; + constexpr unsigned int NumberOfProjectionImages = 160; #endif // Constant image sources diff --git a/test/rtkregularizedconjugategradienttest.cxx b/test/rtkregularizedconjugategradienttest.cxx index a2e38f04e..771645849 100644 --- a/test/rtkregularizedconjugategradienttest.cxx +++ b/test/rtkregularizedconjugategradienttest.cxx @@ -49,8 +49,8 @@ rtkregularizedconjugategradienttest(int, char *[]) auto size = itk::MakeSize(2, 2, 2); auto spacing = itk::MakeVector(252., 252., 252.); #else - auto size = itk::MakeSize(64, 64, 64); - auto spacing = itk::MakeVector(4., 4., 4.); + auto size = itk::MakeSize(48, 48, 48); + auto spacing = itk::MakeVector(5., 5., 5.); #endif tomographySource->SetOrigin(origin); tomographySource->SetSpacing(spacing); @@ -63,8 +63,8 @@ rtkregularizedconjugategradienttest(int, char *[]) size = itk::MakeSize(2, 2, NumberOfProjectionImages); spacing = itk::MakeVector(504., 504., 504.); #else - size = itk::MakeSize(64, 64, NumberOfProjectionImages); - spacing = itk::MakeVector(8., 8., 8.); + size = itk::MakeSize(48, 48, NumberOfProjectionImages); + spacing = itk::MakeVector(10., 10., 8.); #endif projectionsSource->SetOrigin(origin); projectionsSource->SetSpacing(spacing); @@ -113,7 +113,7 @@ rtkregularizedconjugategradienttest(int, char *[]) regularizedConjugateGradient->SetCudaConjugateGradient(false); regularizedConjugateGradient->SetGammaTV(1); - regularizedConjugateGradient->SetTV_iterations(3); + regularizedConjugateGradient->SetTV_iterations(2); regularizedConjugateGradient->SetSoftThresholdWavelets(0.1); regularizedConjugateGradient->SetOrder(3); diff --git a/test/rtksarttest.cxx b/test/rtksarttest.cxx index 45d245020..5856a345f 100644 --- a/test/rtksarttest.cxx +++ b/test/rtksarttest.cxx @@ -38,7 +38,7 @@ rtksarttest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 180; + constexpr unsigned int NumberOfProjectionImages = 120; #endif // Constant image sources diff --git a/test/rtkspectralonesteptest.cxx b/test/rtkspectralonesteptest.cxx index 4fe513af9..95c442c7c 100644 --- a/test/rtkspectralonesteptest.cxx +++ b/test/rtkspectralonesteptest.cxx @@ -82,7 +82,7 @@ rtkspectralonesteptest(int argc, char * argv[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 4; #else - constexpr unsigned int NumberOfProjectionImages = 64; + constexpr unsigned int NumberOfProjectionImages = 32; #endif // Define the material concentrations @@ -240,7 +240,7 @@ rtkspectralonesteptest(int argc, char * argv[]) mechlemOneStep->SetBinnedDetectorResponse(drm); mechlemOneStep->SetMaterialAttenuations(materialAttenuationsMatrix); mechlemOneStep->SetGeometry(geometry); - mechlemOneStep->SetNumberOfIterations(20); + mechlemOneStep->SetNumberOfIterations(5); std::cout << "\n\n****** Case 1: Joseph Backprojector ******" << std::endl; @@ -262,7 +262,7 @@ rtkspectralonesteptest(int argc, char * argv[]) mechlemOneStep->SetBackProjectionFilter(MechlemType::BP_VOXELBASED); mechlemOneStep->SetNumberOfSubsets(4); - mechlemOneStep->SetNumberOfIterations(5); + mechlemOneStep->SetNumberOfIterations(2); MaterialVolumeType::RegionType::SizeType radius; radius.Fill(1); MaterialVolumeType::PixelType weights; diff --git a/test/rtkvariancereconstructiontest.cxx b/test/rtkvariancereconstructiontest.cxx index ead4d424f..1d3c961ac 100644 --- a/test/rtkvariancereconstructiontest.cxx +++ b/test/rtkvariancereconstructiontest.cxx @@ -34,8 +34,8 @@ rtkvariancereconstructiontest(int, char *[]) constexpr unsigned int NumberOfProjectionImages = 3; constexpr unsigned int NumberOfSamples = 2; #else - constexpr unsigned int NumberOfProjectionImages = 50; - constexpr unsigned int NumberOfSamples = 100; + constexpr unsigned int NumberOfProjectionImages = 30; + constexpr unsigned int NumberOfSamples = 30; #endif // Constant image sources diff --git a/test/rtkzengforwardprojectiontest.cxx b/test/rtkzengforwardprojectiontest.cxx index 8b4d46992..da2521f73 100644 --- a/test/rtkzengforwardprojectiontest.cxx +++ b/test/rtkzengforwardprojectiontest.cxx @@ -8,6 +8,10 @@ #include #include +#ifdef USE_CUDA +# include "itkCudaImage.h" +#endif + /** * \file rtkzengforwardprojectiontest.cxx * @@ -35,7 +39,7 @@ rtkzengforwardprojectiontest(int, char *[]) #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else - constexpr unsigned int NumberOfProjectionImages = 40; + constexpr unsigned int NumberOfProjectionImages = 20; #endif // Constant image sources @@ -47,8 +51,8 @@ rtkzengforwardprojectiontest(int, char *[]) auto size = itk::MakeSize(2, 2, 2); auto spacing = itk::MakeVector(252., 252., 252.); #else - auto size = itk::MakeSize(64, 64, 64); - auto spacing = itk::MakeVector(4., 4., 4.); + auto size = itk::MakeSize(33, 33, 33); + auto spacing = itk::MakeVector(8., 8., 8.); #endif auto tomographySource = ConstantImageSourceType::New(); From 4d8b8c5727441e0f2293df74ba4cc3ca32c4bcbd Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 9 Jul 2026 14:24:02 +0200 Subject: [PATCH 5/7] ENH: Add missing cuda variant of existing tests --- test/CMakeLists.txt | 16 ++++++++++++++++ test/rtkdivergencetest.cxx | 4 ++++ test/rtkfbpparalleltest.cxx | 1 + test/rtkfourdadjointoperatorstest.cxx | 6 +++++- test/rtkgradienttest.cxx | 4 ++++ test/rtkhilbertfiltertest.cxx | 5 +++++ test/rtkmotioncompensatedfdktest.cxx | 12 +++++++++++- test/rtkspectralonesteptest.cxx | 6 +++--- test/rtktotalvariationtest.cxx | 4 ++++ 9 files changed, 53 insertions(+), 5 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4fb498491..49cdaa285 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -116,6 +116,7 @@ rtk_add_cuda_test(rtkFDKCudaTest rtkfdktest.cxx) rtk_add_cuda_test(rtkFDKProjWeightCompCudaTest rtkfdkprojweightcompcudatest.cxx) rtk_add_test(rtkFBPParallelTest rtkfbpparalleltest.cxx) +rtk_add_cuda_test(rtkFBPParallelCudaTest rtkfbpparalleltest.cxx) rtk_add_test(rtkImportTest rtkimporttest.cxx) rtk_add_cuda_test(rtkImportCudaTest rtkimporttest.cxx) @@ -123,6 +124,7 @@ rtk_add_cuda_test(rtkImportCudaTest rtkimporttest.cxx) rtk_add_cuda_test(rtkCropFilterCudaTest rtkcroptest.cxx) rtk_add_test(rtkMotionCompensatedFDKTest rtkmotioncompensatedfdktest.cxx) +rtk_add_cuda_test(rtkMotionCompensatedFDKCudaTest rtkmotioncompensatedfdktest.cxx) rtk_add_test(rtkDisplacedDetectorTest rtkdisplaceddetectortest.cxx) rtk_add_cuda_test(rtkDisplacedDetectorCudaTest rtkdisplaceddetectortest.cxx) @@ -141,6 +143,7 @@ rtk_add_test(rtkRampFilterTest2 rtkrampfiltertest2.cxx) rtk_add_cuda_test(rtkRampFilterCudaTest2 rtkrampfiltertest2.cxx) rtk_add_test(rtkHilbertFilterTest rtkhilbertfiltertest.cxx) +rtk_add_cuda_test(rtkHilbertFilterCudaTest rtkhilbertfiltertest.cxx) rtk_add_test(rtkScatterGlareFilterTest rtkscatterglarefiltertest.cxx) rtk_add_cuda_test(rtkScatterGlareFilterCudaTest rtkscatterglarefiltertest.cxx) @@ -154,6 +157,7 @@ rtk_add_test(rtkForwardProjectionTest rtkforwardprojectiontest.cxx) rtk_add_cuda_test(rtkForwardProjectionCudaTest rtkforwardprojectiontest.cxx) rtk_add_test(rtkForwardAttenuatedProjectionTest rtkforwardattenuatedprojectiontest.cxx) rtk_add_test(rtkZengProjectionTest rtkzengforwardprojectiontest.cxx) +rtk_add_cuda_test(rtkZengProjectionCudaTest rtkzengforwardprojectiontest.cxx) rtk_add_test(rtkMaximumIntensityProjectionTest rtkmaximumintensityprojectiontest.cxx) rtk_add_test(rtkGeometryFileTest rtkgeometryfiletest.cxx) @@ -183,6 +187,12 @@ rtk_add_test(rtkSpectralOneStepTest rtkspectralonesteptest.cxx DATA{Input/Spectral/detector_response.mha} DATA{Input/Spectral/material_attenuations.mha} ) +rtk_add_cuda_test(rtkSpectralOneStepCudaTest rtkspectralonesteptest.cxx + DATA{Input/Spectral/OneStep/no_vector_incident_spectrum_64_rows.mha} + DATA{Input/Spectral/OneStep/incident_spectrum_64_rows.mha} + DATA{Input/Spectral/detector_response.mha} + DATA{Input/Spectral/material_attenuations.mha} +) rtk_add_test(rtkVectorImageConvertersTest rtkvectorimageconverterstest.cxx) @@ -298,6 +308,9 @@ rtk_add_cuda_test(rtkAdjointOperatorsCudaTest rtkadjointoperatorstest.cxx) rtk_add_test(rtkFourDAdjointOperatorsTest rtkfourdadjointoperatorstest.cxx DATA{Input/Phases/phases_slow.txt} ) +rtk_add_cuda_test(rtkFourDAdjointOperatorsCudaTest rtkfourdadjointoperatorstest.cxx + DATA{Input/Phases/phases_slow.txt} +) rtk_add_test(rtkInterpolateSplatAdjointTest rtkinterpolatesplatadjointtest.cxx DATA{Input/Phases/phases_slow.txt} @@ -311,10 +324,13 @@ rtk_add_cuda_test(rtkLaplacianCudaTest rtklaplaciantest.cxx ) rtk_add_test(rtkTotalVariationTest rtktotalvariationtest.cxx) +rtk_add_cuda_test(rtkTotalVariationCudaTest rtktotalvariationtest.cxx) rtk_add_test(rtkGradientTest rtkgradienttest.cxx) +rtk_add_cuda_test(rtkGradientCudaTest rtkgradienttest.cxx) rtk_add_test(rtkDivergenceTest rtkdivergencetest.cxx) +rtk_add_cuda_test(rtkDivergenceCudaTest rtkdivergencetest.cxx) rtk_add_test(rtkLagCorrectionTest rtklagcorrectiontest.cxx) rtk_add_cuda_test(rtkLagCorrectionCudaTest rtklagcorrectiontest.cxx) diff --git a/test/rtkdivergencetest.cxx b/test/rtkdivergencetest.cxx index 83e7f2c8b..bb859c50b 100644 --- a/test/rtkdivergencetest.cxx +++ b/test/rtkdivergencetest.cxx @@ -10,6 +10,10 @@ #include "rtkBackwardDifferenceDivergenceImageFilter.h" #include "rtkForwardDifferenceGradientImageFilter.h" +#ifdef USE_CUDA +# include "itkCudaImage.h" +#endif + /** * \file rtkdivergencetest.cxx * diff --git a/test/rtkfbpparalleltest.cxx b/test/rtkfbpparalleltest.cxx index 79fcdcd95..0782a64bc 100644 --- a/test/rtkfbpparalleltest.cxx +++ b/test/rtkfbpparalleltest.cxx @@ -9,6 +9,7 @@ #include "rtkTest.h" #ifdef USE_CUDA +# include "itkCudaImage.h" # include "rtkCudaFDKConeBeamReconstructionFilter.h" #else # include "rtkFDKConeBeamReconstructionFilter.h" diff --git a/test/rtkfourdadjointoperatorstest.cxx b/test/rtkfourdadjointoperatorstest.cxx index 3ddcca519..ff38da02b 100644 --- a/test/rtkfourdadjointoperatorstest.cxx +++ b/test/rtkfourdadjointoperatorstest.cxx @@ -10,6 +10,10 @@ #include +#ifdef USE_CUDA +# include "itkCudaImage.h" +#endif + /** * \file rtkfourdadjointoperatorstest.cxx * @@ -37,7 +41,7 @@ rtkfourdadjointoperatorstest(int argc, char * argv[]) constexpr unsigned int Dimension = 3; using OutputPixelType = float; -#ifdef RTK_USE_CUDA +#ifdef USE_CUDA using ProjectionStackType = itk::CudaImage; using VolumeSeriesType = itk::CudaImage; #else diff --git a/test/rtkgradienttest.cxx b/test/rtkgradienttest.cxx index e736ed624..8844bf21f 100644 --- a/test/rtkgradienttest.cxx +++ b/test/rtkgradienttest.cxx @@ -7,6 +7,10 @@ #include "rtkForwardDifferenceGradientImageFilter.h" #include "rtkTestConfiguration.h" +#ifdef USE_CUDA +# include "itkCudaImage.h" +#endif + template #if FAST_TESTS_NO_CHECKS void diff --git a/test/rtkhilbertfiltertest.cxx b/test/rtkhilbertfiltertest.cxx index 2a1fadd3f..01ceb82bc 100644 --- a/test/rtkhilbertfiltertest.cxx +++ b/test/rtkhilbertfiltertest.cxx @@ -1,6 +1,11 @@ #include "rtkTest.h" #include "rtkFFTHilbertImageFilter.h" + +#ifdef USE_CUDA +# include "itkCudaImage.h" +#endif + /** * \file rtkhilbertfiltertest.cxx * diff --git a/test/rtkmotioncompensatedfdktest.cxx b/test/rtkmotioncompensatedfdktest.cxx index ed8f60b3a..5fdf4b6c5 100644 --- a/test/rtkmotioncompensatedfdktest.cxx +++ b/test/rtkmotioncompensatedfdktest.cxx @@ -4,12 +4,18 @@ #include "rtkConstantImageSource.h" #include "rtkCyclicDeformationImageFilter.h" #include "rtkDrawEllipsoidImageFilter.h" -#include "rtkFDKConeBeamReconstructionFilter.h" #include "rtkFDKWarpBackProjectionImageFilter.h" #include "rtkFieldOfViewImageFilter.h" #include "rtkRayEllipsoidIntersectionImageFilter.h" #include "rtkTest.h" +#ifdef USE_CUDA +# include "itkCudaImage.h" +# include "rtkCudaFDKConeBeamReconstructionFilter.h" +#else +# include "rtkFDKConeBeamReconstructionFilter.h" +#endif + /** * \file rtkmotioncompensatedfdktest.cxx * @@ -28,7 +34,11 @@ int rtkmotioncompensatedfdktest(int, char *[]) { constexpr unsigned int Dimension = 3; +#ifdef USE_CUDA + using OutputImageType = itk::CudaImage; +#else using OutputImageType = itk::Image; +#endif #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; #else diff --git a/test/rtkspectralonesteptest.cxx b/test/rtkspectralonesteptest.cxx index 95c442c7c..8e985a4e7 100644 --- a/test/rtkspectralonesteptest.cxx +++ b/test/rtkspectralonesteptest.cxx @@ -46,7 +46,7 @@ rtkspectralonesteptest(int argc, char * argv[]) using VectorImageType = typename itk::VectorImage; -#ifdef RTK_USE_CUDA +#ifdef USE_CUDA using MaterialVolumeType = itk::CudaImage; using MeasuredProjectionsType = itk::CudaImage; using IncidentSpectrumImageType = itk::CudaImage; @@ -274,7 +274,7 @@ rtkspectralonesteptest(int argc, char * argv[]) CheckVectorImageQuality(mechlemOneStep->GetOutput(), composeVols->GetOutput(), 0.08, 23, 2.0); std::cout << "\n\nTest PASSED! " << std::endl; -#ifdef RTK_USE_CUDA +#ifdef USE_CUDA std::cout << "\n\n****** Case 4: CUDA voxel-based Backprojector, 4 subsets, with regularization ******" << std::endl; mechlemOneStep->SetForwardProjectionFilter(MechlemType::FP_CUDARAYCAST); @@ -285,7 +285,7 @@ rtkspectralonesteptest(int argc, char * argv[]) std::cout << "\n\nTest PASSED! " << std::endl; #endif -#ifdef RTK_USE_CUDA +#ifdef USE_CUDA std::cout << "\n\n****** Case 5: CUDA voxel-based Backprojector, 4 subsets, with regularization, itkVectorImage inputs ******" << std::endl; diff --git a/test/rtktotalvariationtest.cxx b/test/rtktotalvariationtest.cxx index 7948a39f1..9b02c7207 100644 --- a/test/rtktotalvariationtest.cxx +++ b/test/rtktotalvariationtest.cxx @@ -5,6 +5,10 @@ #include "rtkTotalVariationDenoisingBPDQImageFilter.h" #include "rtkTotalVariationImageFilter.h" +#ifdef USE_CUDA +# include "itkCudaImage.h" +#endif + template void CheckTotalVariation(const typename TImage::Pointer & before, const typename TImage::Pointer & after) From 82429b361a909692198592d6170c3e6702a5615a Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 9 Jul 2026 14:26:33 +0200 Subject: [PATCH 6/7] BUG: Fix CUDA Weidinger Hessian initialization size Initialize the CUDA Weidinger forward model outputs with their actual sizes: the gradient output has one value per material, while the Hessian output has one value per material pair. This avoids leaving part of the Hessian buffer uninitialized. --- src/rtkCudaWeidingerForwardModelImageFilter.cu | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/rtkCudaWeidingerForwardModelImageFilter.cu b/src/rtkCudaWeidingerForwardModelImageFilter.cu index 3113721ff..6fee318d2 100644 --- a/src/rtkCudaWeidingerForwardModelImageFilter.cu +++ b/src/rtkCudaWeidingerForwardModelImageFilter.cu @@ -186,10 +186,11 @@ CUDA_WeidingerForwardModel(int projectionSize[3], cudaMemcpyToSymbol(c_binnedDetectorResponse, &(binnedDetectorResponse[0]), nBins * nEnergies * sizeof(float)); cudaMemcpyToSymbol(c_materialAttenuations, &(materialAttenuations[0]), nMaterials * nEnergies * sizeof(float)); - // Set both outputs to zeros - size_t memSize = sizeof(float) * projectionSize[0] * projectionSize[1] * projectionSize[2] * nMaterials; - cudaMemset((void *)pOut1, 0, memSize); - cudaMemset((void *)pOut2, 0, memSize); + // Set both outputs to zeros. Output 1 is a material gradient vector, + // output 2 is a material Hessian matrix. + const size_t numberOfPixels = projectionSize[0] * projectionSize[1] * projectionSize[2]; + cudaMemset((void *)pOut1, 0, sizeof(float) * numberOfPixels * nMaterials); + cudaMemset((void *)pOut2, 0, sizeof(float) * numberOfPixels * nMaterials * nMaterials); dim3 dimBlock = dim3(4, 4, 4); dim3 dimGrid = dim3(iDivUp(projectionSize[0], dimBlock.x), From 036369257987745d6429eddec2b75055bc5c059a Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 9 Jul 2026 14:32:24 +0200 Subject: [PATCH 7/7] ENH: Remove cuda compability in rtkforwardattenuatedprojectiontest The cuda part was likely an error made with a copy/paste --- test/rtkforwardattenuatedprojectiontest.cxx | 22 +++------------------ 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/test/rtkforwardattenuatedprojectiontest.cxx b/test/rtkforwardattenuatedprojectiontest.cxx index bfe771420..15ebcfa58 100644 --- a/test/rtkforwardattenuatedprojectiontest.cxx +++ b/test/rtkforwardattenuatedprojectiontest.cxx @@ -12,11 +12,7 @@ #include -#ifdef USE_CUDA -# include "rtkCudaForwardProjectionImageFilter.h" -#else -# include "rtkJosephForwardAttenuatedProjectionImageFilter.h" -#endif +#include "rtkJosephForwardAttenuatedProjectionImageFilter.h" /** * \file rtkforwardattenuatedprojectiontest.cxx @@ -36,11 +32,7 @@ rtkforwardattenuatedprojectiontest(int, char *[]) constexpr unsigned int Dimension = 3; using OutputPixelType = float; -#ifdef USE_CUDA - using OutputImageType = itk::CudaImage; -#else - using OutputImageType = itk::Image; -#endif +using OutputImageType = itk::Image; #if FAST_TESTS_NO_CHECKS constexpr unsigned int NumberOfProjectionImages = 3; @@ -113,21 +105,13 @@ rtkforwardattenuatedprojectiontest(int, char *[]) projTotal->SetConstant(0.); projTotal->Update(); - // Joseph Forward Projection filter -#ifdef USE_CUDA - using JFPType = rtk::CudaForwardProjectionImageFilter; -#else + // Joseph attenuated forward projection filter using JFPType = rtk::JosephForwardAttenuatedProjectionImageFilter; -#endif auto jfp = JFPType::New(); jfp->InPlaceOff(); jfp->SetInput(projTotal->GetOutput()); jfp->SetInput(1, volumeSource); jfp->SetInput(2, attenuationMap); - -#ifdef USE_CUDA - jfp->SetStepSize(10); -#endif // Geometry using GeometryType = rtk::ThreeDCircularProjectionGeometry; auto geometry_projection = GeometryType::New();