From 95c6923764a217145c3f6615c4eae7dab263a41d Mon Sep 17 00:00:00 2001 From: Cyril Mory Date: Thu, 8 Jan 2026 13:44:40 +0100 Subject: [PATCH] ENH: Add 4D C++ reconstruction examples and corresponding doc pages --- .../rtkfourdconjugategradient.cxx | 1 - .../FourDConjugateGradient/CMakeLists.txt | 30 +++ .../FourDConjugateGradient.cxx | 101 +++++++++ examples/FourDConjugateGradient/README.md | 15 ++ examples/FourDFDK/CMakeLists.txt | 30 +++ examples/FourDFDK/FourDFDK.cxx | 103 +++++++++ examples/FourDFDK/README.md | 15 ++ examples/FourDROOSTER/CMakeLists.txt | 30 +++ examples/FourDROOSTER/FourDROOSTER.cxx | 143 ++++++++++++ examples/FourDROOSTER/README.md | 15 ++ examples/GenerateFourDData/CMakeLists.txt | 12 + .../GenerateFourDData/GenerateFourDData.cxx | 214 ++++++++++++++++++ examples/GenerateFourDData/README.md | 21 ++ examples/README.md | 4 + include/rtkCudaExternTemplates.h | 4 + src/rtkCudaExternTemplates.cxx | 2 + test/CMakeLists.txt | 26 +++ 17 files changed, 765 insertions(+), 1 deletion(-) create mode 100644 examples/FourDConjugateGradient/CMakeLists.txt create mode 100644 examples/FourDConjugateGradient/FourDConjugateGradient.cxx create mode 100644 examples/FourDConjugateGradient/README.md create mode 100644 examples/FourDFDK/CMakeLists.txt create mode 100644 examples/FourDFDK/FourDFDK.cxx create mode 100644 examples/FourDFDK/README.md create mode 100644 examples/FourDROOSTER/CMakeLists.txt create mode 100644 examples/FourDROOSTER/FourDROOSTER.cxx create mode 100644 examples/FourDROOSTER/README.md create mode 100644 examples/GenerateFourDData/CMakeLists.txt create mode 100644 examples/GenerateFourDData/GenerateFourDData.cxx create mode 100644 examples/GenerateFourDData/README.md diff --git a/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx b/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx index fc57e78a0..f66820fbb 100644 --- a/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx +++ b/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx @@ -26,7 +26,6 @@ #ifdef RTK_USE_CUDA # include "itkCudaImage.h" -# include "rtkCudaConstantVolumeSeriesSource.h" #endif #include diff --git a/examples/FourDConjugateGradient/CMakeLists.txt b/examples/FourDConjugateGradient/CMakeLists.txt new file mode 100644 index 000000000..b3f694b08 --- /dev/null +++ b/examples/FourDConjugateGradient/CMakeLists.txt @@ -0,0 +1,30 @@ +cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR) + +# This project is designed to be built outside the RTK source tree. +project(FourDConjugateGradient) + +# Find ITK with RTK +find_package( + ITK + REQUIRED + COMPONENTS + RTK + ITKImageIO +) +if("${ITK_VERSION_MAJOR}" VERSION_LESS "6") + include(${ITK_USE_FILE}) + set(RTK_EXAMPLE_TARGETS ${ITK_LIBRARIES}) +else() + itk_generate_factory_registration() + set( + RTK_EXAMPLE_TARGETS + ${ITK_INTERFACE_LIBRARIES} + ITK::ITKImageIO + ITK::ITKFFTImageFilterInit + ) + message(ITK_INTERFACE_LIBRARIES=${ITK_INTERFACE_LIBRARIES}) +endif() + +# Executable(s) +add_executable(FourDConjugateGradient FourDConjugateGradient.cxx) +target_link_libraries(FourDConjugateGradient ${RTK_EXAMPLE_TARGETS}) diff --git a/examples/FourDConjugateGradient/FourDConjugateGradient.cxx b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx new file mode 100644 index 000000000..95c2c77af --- /dev/null +++ b/examples/FourDConjugateGradient/FourDConjugateGradient.cxx @@ -0,0 +1,101 @@ +#include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h" +#include "rtkIterationCommands.h" +#include "rtkSignalToInterpolationWeights.h" +#include "rtkReorderProjectionsImageFilter.h" +#include "rtkProjectionsReader.h" +#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h" + +#ifdef RTK_USE_CUDA +# include +#endif +#include + +int +main(int argc, char * argv[]) +{ + using OutputPixelType = float; + +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; +#endif + + // Generate the input volume series, used as initial estimate by 4D conjugate gradient + auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.); + auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.); + auto fourDSize = itk::MakeSize(32, 16, 32, 8); + + using ConstantFourDSourceType = rtk::ConstantImageSource; + auto fourDSource = ConstantFourDSourceType::New(); + + fourDSource->SetOrigin(fourDOrigin); + fourDSource->SetSpacing(fourDSpacing); + fourDSource->SetSize(fourDSize); + fourDSource->Update(); + + // Read geometry, projections and signal + rtk::ThreeDCircularProjectionGeometry::Pointer geometry; + TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml")); + + using ReaderType = rtk::ProjectionsReader; + auto projectionsReader = ReaderType::New(); + std::vector fileNames = std::vector(); + fileNames.push_back("four_d_projections.mha"); + projectionsReader->SetFileNames(fileNames); + TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update()); + + // Re-order geometry and projections + // In the new order, projections with identical phases are packed together + std::vector signal = rtk::ReadSignalFile("four_d_signal.txt"); + auto reorder = rtk::ReorderProjectionsImageFilter::New(); + reorder->SetInput(projectionsReader->GetOutput()); + reorder->SetInputGeometry(geometry); + reorder->SetInputSignal(signal); + TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update()) + + // Release the memory holding the stack of original projections + projectionsReader->GetOutput()->ReleaseData(); + + // Compute the interpolation weights + auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New(); + signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal()); + signalToInterpolationWeights->SetNumberOfReconstructedFrames(fourDSize[3]); + TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update()) + + // Set the forward and back projection filters to be used + using ConjugateGradientFilterType = + rtk::FourDConjugateGradientConeBeamReconstructionFilter; + auto fourdconjugategradient = ConjugateGradientFilterType::New(); + + fourdconjugategradient->SetInputVolumeSeries(fourDSource->GetOutput()); + fourdconjugategradient->SetNumberOfIterations(2); +#ifdef RTK_USE_CUDA + fourdconjugategradient->SetCudaConjugateGradient(true); + fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_CUDARAYCAST); + fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_CUDAVOXELBASED); +#else + fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_JOSEPH); + fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_VOXELBASED); +#endif + + // Set the newly ordered arguments + fourdconjugategradient->SetInputProjectionStack(reorder->GetOutput()); + fourdconjugategradient->SetGeometry(reorder->GetOutputGeometry()); + fourdconjugategradient->SetWeights(signalToInterpolationWeights->GetOutput()); + fourdconjugategradient->SetSignal(reorder->GetOutputSignal()); + + auto verboseIterationCommand = rtk::VerboseIterationCommand::New(); + fourdconjugategradient->AddObserver(itk::AnyEvent(), verboseIterationCommand); + + TRY_AND_EXIT_ON_ITK_EXCEPTION(fourdconjugategradient->Update()) + + // Write + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(fourdconjugategradient->GetOutput(), "fourdconjugategradient.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/FourDConjugateGradient/README.md b/examples/FourDConjugateGradient/README.md new file mode 100644 index 000000000..77c1a644d --- /dev/null +++ b/examples/FourDConjugateGradient/README.md @@ -0,0 +1,15 @@ +# 4D Conjugate Gradient + +This example shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources. +It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./FourDConjugateGradient.cxx +:language: c++ +``` +```` + +````` diff --git a/examples/FourDFDK/CMakeLists.txt b/examples/FourDFDK/CMakeLists.txt new file mode 100644 index 000000000..0b4cd4703 --- /dev/null +++ b/examples/FourDFDK/CMakeLists.txt @@ -0,0 +1,30 @@ +cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR) + +# This project is designed to be built outside the RTK source tree. +project(FourDFDK) + +# Find ITK with RTK +find_package( + ITK + REQUIRED + COMPONENTS + RTK + ITKImageIO +) +if("${ITK_VERSION_MAJOR}" VERSION_LESS "6") + include(${ITK_USE_FILE}) + set(RTK_EXAMPLE_TARGETS ${ITK_LIBRARIES}) +else() + itk_generate_factory_registration() + set( + RTK_EXAMPLE_TARGETS + ${ITK_INTERFACE_LIBRARIES} + ITK::ITKImageIO + ITK::ITKFFTImageFilterInit + ) + message(ITK_INTERFACE_LIBRARIES=${ITK_INTERFACE_LIBRARIES}) +endif() + +# Executable(s) +add_executable(FourDFDK FourDFDK.cxx) +target_link_libraries(FourDFDK ${RTK_EXAMPLE_TARGETS}) diff --git a/examples/FourDFDK/FourDFDK.cxx b/examples/FourDFDK/FourDFDK.cxx new file mode 100644 index 000000000..7211622a4 --- /dev/null +++ b/examples/FourDFDK/FourDFDK.cxx @@ -0,0 +1,103 @@ +#include "rtkConfiguration.h" +#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h" +#include "rtkFDKConeBeamReconstructionFilter.h" +#ifdef RTK_USE_CUDA +# include "rtkCudaFDKConeBeamReconstructionFilter.h" +#endif +#include "rtkSelectOneProjectionPerCycleImageFilter.h" +#include "rtkProjectionsReader.h" + +#ifdef RTK_USE_CUDA +# include +#endif +#include + +int +main(int argc, char * argv[]) +{ + using OutputPixelType = float; + constexpr unsigned int Dimension = 3; + + using CPUOutputImageType = itk::Image; +#ifdef RTK_USE_CUDA + using OutputImageType = itk::CudaImage; +#else + using OutputImageType = CPUOutputImageType; +#endif + + // Read geometry, projections and signal + rtk::ThreeDCircularProjectionGeometry::Pointer geometry; + TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml")); + + using ReaderType = rtk::ProjectionsReader; + auto projectionsReader = ReaderType::New(); + std::vector fileNames = std::vector(); + fileNames.push_back("four_d_projections.mha"); + projectionsReader->SetFileNames(fileNames); + TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update()); + + // Part specific to 4D + auto selector = rtk::SelectOneProjectionPerCycleImageFilter::New(); + selector->SetInput(projectionsReader->GetOutput()); + selector->SetInputGeometry(geometry); + selector->SetSignalFilename("four_d_signal.txt"); + + // Create one frame of the reconstructed image + using ConstantImageSourceType = rtk::ConstantImageSource; + auto constantImageSource = ConstantImageSourceType::New(); + constantImageSource->SetOrigin(itk::MakePoint(-63., -31., -63.)); + constantImageSource->SetSpacing(itk::MakeVector(4., 4., 4.)); + constantImageSource->SetSize(itk::MakeSize(32, 16, 32)); + constantImageSource->SetConstant(0.); + constantImageSource->Update(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(constantImageSource->Update()) + + // FDK reconstruction filtering +#ifdef RTK_USE_CUDA + using FDKType = rtk::CudaFDKConeBeamReconstructionFilter; +#else + using FDKType = rtk::FDKConeBeamReconstructionFilter; +#endif + auto feldkamp = FDKType::New(); + feldkamp->SetInput(0, constantImageSource->GetOutput()); + feldkamp->SetInput(1, selector->GetOutput()); + feldkamp->SetGeometry(selector->GetOutputGeometry()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->Update()); + + // Create empty 4D image + using FourDOutputImageType = itk::Image; + using ConstantFourDSourceType = rtk::ConstantImageSource; + auto fourDSource = ConstantFourDSourceType::New(); + fourDSource->SetOrigin(itk::MakePoint(-63., -31., -63., 0.)); + fourDSource->SetSpacing(itk::MakeVector(4., 4., 4., 1.)); + fourDSource->SetSize(itk::MakeSize(32, 16, 32, 4)); + fourDSource->SetConstant(0.); + fourDSource->Update(); + + // Go over each frame, reconstruct 3D frame and paste with iterators in 4D image + for (int f = 0; f < 4; f++) + { + selector->SetPhase(f / (double)4); + TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->UpdateLargestPossibleRegion()) + + ConstantFourDSourceType::OutputImageRegionType region; + region = fourDSource->GetOutput()->GetLargestPossibleRegion(); + region.SetIndex(3, f); + region.SetSize(3, 1); + + itk::ImageRegionIterator it4D(fourDSource->GetOutput(), region); + itk::ImageRegionIterator it3D(feldkamp->GetOutput(), + feldkamp->GetOutput()->GetLargestPossibleRegion()); + while (!it3D.IsAtEnd()) + { + it4D.Set(it3D.Get()); + ++it4D; + ++it3D; + } + } + + // Write + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(fourDSource->GetOutput(), "fourdfdk.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/FourDFDK/README.md b/examples/FourDFDK/README.md new file mode 100644 index 000000000..0b06ab133 --- /dev/null +++ b/examples/FourDFDK/README.md @@ -0,0 +1,15 @@ +# 4D FDK + +This example shows how to perform filtered backprojection cone-beam CT reconstruction using either CPU or GPU resources depending on how RTK was compiled. +It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./FourDFDK.cxx +:language: c++ +``` +```` + +````` diff --git a/examples/FourDROOSTER/CMakeLists.txt b/examples/FourDROOSTER/CMakeLists.txt new file mode 100644 index 000000000..eaa95e17c --- /dev/null +++ b/examples/FourDROOSTER/CMakeLists.txt @@ -0,0 +1,30 @@ +cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR) + +# This project is designed to be built outside the RTK source tree. +project(FourDROOSTER) + +# Find ITK with RTK +find_package( + ITK + REQUIRED + COMPONENTS + RTK + ITKImageIO +) +if("${ITK_VERSION_MAJOR}" VERSION_LESS "6") + include(${ITK_USE_FILE}) + set(RTK_EXAMPLE_TARGETS ${ITK_LIBRARIES}) +else() + itk_generate_factory_registration() + set( + RTK_EXAMPLE_TARGETS + ${ITK_INTERFACE_LIBRARIES} + ITK::ITKImageIO + ITK::ITKFFTImageFilterInit + ) + message(ITK_INTERFACE_LIBRARIES=${ITK_INTERFACE_LIBRARIES}) +endif() + +# Executable(s) +add_executable(FourDROOSTER FourDROOSTER.cxx) +target_link_libraries(FourDROOSTER ${RTK_EXAMPLE_TARGETS}) diff --git a/examples/FourDROOSTER/FourDROOSTER.cxx b/examples/FourDROOSTER/FourDROOSTER.cxx new file mode 100644 index 000000000..825641211 --- /dev/null +++ b/examples/FourDROOSTER/FourDROOSTER.cxx @@ -0,0 +1,143 @@ +#include "rtkFourDROOSTERConeBeamReconstructionFilter.h" +#include "rtkIterationCommands.h" +#include "rtkSignalToInterpolationWeights.h" +#include "rtkReorderProjectionsImageFilter.h" +#include "rtkProjectionsReader.h" +#include "rtkThreeDCircularProjectionGeometryXMLFileReader.h" + +#ifdef RTK_USE_CUDA +# include +#endif +#include + +int +main(int argc, char * argv[]) +{ + using OutputPixelType = float; + using DVFVectorType = itk::CovariantVector; +#ifdef RTK_USE_CUDA + using VolumeSeriesType = itk::CudaImage; + using ProjectionStackType = itk::CudaImage; + using VolumeType = itk::CudaImage; + using DVFSequenceImageType = itk::CudaImage; +#else + using VolumeSeriesType = itk::Image; + using ProjectionStackType = itk::Image; + using VolumeType = itk::Image; + using DVFSequenceImageType = itk::Image; +#endif + + // Generate the input volume series, used as initial estimate by 4D conjugate gradient + auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.); + auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.); + auto fourDSize = itk::MakeSize(32, 16, 32, 8); + + using ConstantFourDSourceType = rtk::ConstantImageSource; + auto fourDSource = ConstantFourDSourceType::New(); + + fourDSource->SetOrigin(fourDOrigin); + fourDSource->SetSpacing(fourDSpacing); + fourDSource->SetSize(fourDSize); + fourDSource->Update(); + + // Read geometry, projections and signal + rtk::ThreeDCircularProjectionGeometry::Pointer geometry; + TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml")); + + using ReaderType = rtk::ProjectionsReader; + auto projectionsReader = ReaderType::New(); + std::vector fileNames = std::vector(); + fileNames.push_back("four_d_projections.mha"); + projectionsReader->SetFileNames(fileNames); + TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update()); + + // Re-order geometry and projections + // In the new order, projections with identical phases are packed together + std::vector signal = rtk::ReadSignalFile("four_d_signal.txt"); + auto reorder = rtk::ReorderProjectionsImageFilter::New(); + reorder->SetInput(projectionsReader->GetOutput()); + reorder->SetInputGeometry(geometry); + reorder->SetInputSignal(signal); + TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update()) + + // Release the memory holding the stack of original projections + projectionsReader->GetOutput()->ReleaseData(); + + // Compute the interpolation weights + auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New(); + signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal()); + signalToInterpolationWeights->SetNumberOfReconstructedFrames(fourDSize[3]); + TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update()) + + // Set the forward and back projection filters to be used + using ROOSTERFilterType = rtk::FourDROOSTERConeBeamReconstructionFilter; + auto rooster = ROOSTERFilterType::New(); + + rooster->SetInputVolumeSeries(fourDSource->GetOutput()); + rooster->SetCG_iterations(2); + rooster->SetMainLoop_iterations(2); +#ifdef RTK_USE_CUDA + rooster->SetCudaConjugateGradient(true); + rooster->SetUseCudaCyclicDeformation(true); + rooster->SetForwardProjectionFilter(ROOSTERFilterType::FP_CUDARAYCAST); + rooster->SetBackProjectionFilter(ROOSTERFilterType::BP_CUDAVOXELBASED); +#else + rooster->SetForwardProjectionFilter(ROOSTERFilterType::FP_JOSEPH); + rooster->SetBackProjectionFilter(ROOSTERFilterType::BP_VOXELBASED); +#endif + + // Set the newly ordered arguments + rooster->SetInputProjectionStack(reorder->GetOutput()); + rooster->SetGeometry(reorder->GetOutputGeometry()); + rooster->SetWeights(signalToInterpolationWeights->GetOutput()); + rooster->SetSignal(reorder->GetOutputSignal()); + + // For each optional regularization step, set whether or not + // it should be performed, and provide the necessary inputs + + // Positivity + rooster->SetPerformPositivity(true); + + // Motion mask + rooster->SetPerformMotionMask(false); + + // Spatial TV + rooster->SetGammaTVSpace(0.1); + rooster->SetTV_iterations(4); + rooster->SetPerformTVSpatialDenoising(true); + + // Spatial wavelets + rooster->SetPerformWaveletsSpatialDenoising(false); + + // Temporal TV + rooster->SetGammaTVTime(0.1); + rooster->SetTV_iterations(4); + rooster->SetPerformTVTemporalDenoising(true); + + // Temporal L0 + rooster->SetPerformL0TemporalDenoising(false); + + // Total nuclear variation + rooster->SetPerformTNVDenoising(false); + + // Warping + rooster->SetPerformWarping(true); + rooster->SetUseNearestNeighborInterpolationInWarping(false); + DVFSequenceImageType::Pointer dvf; + TRY_AND_EXIT_ON_ITK_EXCEPTION(dvf = itk::ReadImage("four_d_dvf.mha")) + rooster->SetDisplacementField(dvf); + rooster->SetComputeInverseWarpingByConjugateGradient(false); + DVFSequenceImageType::Pointer idvf; + TRY_AND_EXIT_ON_ITK_EXCEPTION(idvf = itk::ReadImage("four_d_idvf.mha")) + rooster->SetInverseDisplacementField(idvf); + + auto verboseIterationCommand = rtk::VerboseIterationCommand::New(); + rooster->AddObserver(itk::AnyEvent(), verboseIterationCommand); + + TRY_AND_EXIT_ON_ITK_EXCEPTION(rooster->Update()) + + // Write + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(rooster->GetOutput(), "fourdrooster.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/FourDROOSTER/README.md b/examples/FourDROOSTER/README.md new file mode 100644 index 000000000..477d17f1e --- /dev/null +++ b/examples/FourDROOSTER/README.md @@ -0,0 +1,15 @@ +# 4D ROOSTER + +This example shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources. +It reads its data from disk. The data can be generated by the [*Generate 4D data* example](../GenerateFourDData/README.md) or [downloaded from Girder](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./FourDROOSTER.cxx +:language: c++ +``` +```` + +````` diff --git a/examples/GenerateFourDData/CMakeLists.txt b/examples/GenerateFourDData/CMakeLists.txt new file mode 100644 index 000000000..04b65011a --- /dev/null +++ b/examples/GenerateFourDData/CMakeLists.txt @@ -0,0 +1,12 @@ +cmake_minimum_required(VERSION 3.9.5 FATAL_ERROR) + +# This project is designed to be built outside the RTK source tree. +project(GenerateFourDData) + +# Find ITK with RTK +find_package(ITK REQUIRED COMPONENTS RTK) +include(${ITK_USE_FILE}) + +# Executable(s) +add_executable(GenerateFourDData GenerateFourDData.cxx) +target_link_libraries(GenerateFourDData ${ITK_LIBRARIES}) diff --git a/examples/GenerateFourDData/GenerateFourDData.cxx b/examples/GenerateFourDData/GenerateFourDData.cxx new file mode 100644 index 000000000..0044accb3 --- /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::MakeVector(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::MakeVector(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::MakeVector(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::MakeVector(4 * (itk::Math::abs((4 + t) % 8 - 4.) - 2.), 0., 0.)); + de2->InPlaceOff(); + de2->Update(); + + using DuplicatorType = itk::ImageDuplicator; + auto duplicator = DuplicatorType::New(); + duplicator->SetInputImage(de2->GetOutput()); + duplicator->Update(); + + join->SetInput(t, duplicator->GetOutput()); + } + + join->Update(); + TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(join->GetOutput(), "four_d_ground_truth.mha")); + + return EXIT_SUCCESS; +} diff --git a/examples/GenerateFourDData/README.md b/examples/GenerateFourDData/README.md new file mode 100644 index 000000000..1cc19be41 --- /dev/null +++ b/examples/GenerateFourDData/README.md @@ -0,0 +1,21 @@ +# Generate 4D data + +This example shows how to generate a full set of input data to run 4D examples. It contains: +- a moving phantom, made of two ellipsoids, one of which moves along the first dimension +- a geometry +- a phase signal +- the set of projections of the moving phantom +- a deformation vector field describing the motion of the moving ellipsoid +- the corresponding inverse deformation vector field + +You can also skip this part and [download the data](https://data.kitware.com/#collection/5a7706878d777f0649e04776/folder/69611c373b9bcc32c3188531). + +`````{tab-set} + +````{tab-item} C++ + +```{literalinclude} ./GenerateFourDData.cxx +:language: c++ +``` +```` +````` diff --git a/examples/README.md b/examples/README.md index d6d775f18..f67e085b3 100644 --- a/examples/README.md +++ b/examples/README.md @@ -22,4 +22,8 @@ providing efficient implementations for high-performance applications. ./InlineReconstruction/README.md ./AddNoise/README.md ./GeometricPhantom/README.md +./GenerateFourDData/README.md +./FourDFDK/README.md +./FourDConjugateGradient/README.md +./FourDROOSTER/README.md ``` diff --git a/include/rtkCudaExternTemplates.h b/include/rtkCudaExternTemplates.h index 1d7704ba1..1638acb50 100644 --- a/include/rtkCudaExternTemplates.h +++ b/include/rtkCudaExternTemplates.h @@ -22,7 +22,9 @@ #include "rtkConfiguration.h" #ifdef RTK_USE_CUDA +# include # include +# include # include # include # include @@ -49,6 +51,8 @@ extern template class RTK_EXPORT_EXPLICIT itk::ImageSource, itk::CudaImage>; extern template class RTK_EXPORT_EXPLICIT itk::InPlaceImageFilter, itk::CudaImage>; extern template class RTK_EXPORT_EXPLICIT itk::ImageSource, 3>>; +extern template class RTK_EXPORT_EXPLICIT itk::ExtractImageFilter, itk::CudaImage>; +extern template class RTK_EXPORT_EXPLICIT itk::CropImageFilter, itk::CudaImage>; ITK_GCC_PRAGMA_DIAG_POP() # endif diff --git a/src/rtkCudaExternTemplates.cxx b/src/rtkCudaExternTemplates.cxx index 594c9b396..9d6d8fd6f 100644 --- a/src/rtkCudaExternTemplates.cxx +++ b/src/rtkCudaExternTemplates.cxx @@ -25,6 +25,8 @@ template class itk::ImageSource>; template class itk::ImageToImageFilter, itk::CudaImage>; template class itk::InPlaceImageFilter, itk::CudaImage>; template class itk::ImageSource, 3>>; +template class itk::ExtractImageFilter, itk::CudaImage>; +template class itk::CropImageFilter, itk::CudaImage>; // RTK base class explicit instantiation definitions # include "rtkBackProjectionImageFilter.h" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fb744596c..c066b14ad 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -409,6 +409,32 @@ rtk_add_test(rtkConjugateGradientExampleTest ../examples/ConjugateGradient/ConjugateGradient.cxx DATA{Input/Forbild/Thorax} ) +rtk_add_test(rtkGenerateFourDDataExampleTest + ../examples/GenerateFourDData/GenerateFourDData.cxx +) +set_tests_properties( + rtkGenerateFourDDataExampleTest + PROPERTIES + FIXTURES_SETUP + FourDData +) +rtk_add_test(rtkFourDConjugateGradientExampleTest + ../examples/FourDConjugateGradient/FourDConjugateGradient.cxx +) +rtk_add_test(rtkFourDFDKExampleTest + ../examples/FourDFDK/FourDFDK.cxx +) +rtk_add_test(rtkFourDROOSTERExampleTest + ../examples/FourDROOSTER/FourDROOSTER.cxx +) +set_tests_properties( + rtkFourDConjugateGradientExampleTest + rtkFourDFDKExampleTest + rtkFourDROOSTERExampleTest + PROPERTIES + FIXTURES_REQUIRED + FourDData +) if(ITK_WRAP_PYTHON) itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py)