ENH: Add FourDConjugateGradient example and corresponding doc page

Author: cyrilmory

applications/rtkfourdconjugategradient/rtkfourdconjugategradient.cxx

@@ -26,7 +26,6 @@
 
 #ifdef RTK_USE_CUDA
 #  include "itkCudaImage.h"
-#  include "rtkCudaConstantVolumeSeriesSource.h"
 #endif
 
 #include <itkImageFileWriter.h>

examples/FourDConjugateGradient/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(FourDConjugateGradient)
+
+# Find ITK with RTK
+find_package(ITK REQUIRED COMPONENTS RTK)
+include(${ITK_USE_FILE})
+
+# Executable(s)
+add_executable(FourDConjugateGradient FourDConjugateGradient.cxx)
+target_link_libraries(FourDConjugateGradient ${ITK_LIBRARIES})

examples/FourDConjugateGradient/FourDConjugateGradient.cxx

@@ -0,0 +1,79 @@
+// #include "rtkThreeDCircularProjectionGeometryXMLFile.h"
+#include "rtkFourDConjugateGradientConeBeamReconstructionFilter.h"
+#include "rtkIterationCommands.h"
+#include "rtkSignalToInterpolationWeights.h"
+#include "rtkReorderProjectionsImageFilter.h"
+#include "../../test/rtkFourDTestHelper.h"
+
+#ifdef RTK_USE_CUDA
+#  include <itkCudaImage.h>
+#endif
+#include <itkImageFileWriter.h>
+
+int
+main(int argc, char * argv[])
+{
+  using OutputPixelType = float;
+
+#ifdef RTK_USE_CUDA
+  using VolumeSeriesType = itk::CudaImage<OutputPixelType, 4>;
+  using ProjectionStackType = itk::CudaImage<OutputPixelType, 3>;
+  using VolumeType = itk::CudaImage<OutputPixelType, 3>;
+#else
+  using VolumeSeriesType = itk::Image<OutputPixelType, 4>;
+  using ProjectionStackType = itk::Image<OutputPixelType, 3>;
+  using VolumeType = itk::Image<OutputPixelType, 3>;
+#endif
+
+  auto data = rtk::GenerateFourDTestData<OutputPixelType>(FAST_TESTS_NO_CHECKS);
+
+  // Re-order geometry and projections
+  // In the new order, projections with identical phases are packed together
+  auto reorder = rtk::ReorderProjectionsImageFilter<ProjectionStackType>::New();
+  reorder->SetInput(data.Projections);
+  reorder->SetInputGeometry(data.Geometry);
+  reorder->SetInputSignal(data.Signal);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update())
+
+  // Release the memory holding the stack of original projections
+  data.Projections->ReleaseData();
+
+  // Compute the interpolation weights
+  auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New();
+  signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal());
+  signalToInterpolationWeights->SetNumberOfReconstructedFrames(
+    data.InitialVolumeSeries->GetLargestPossibleRegion().GetSize(3));
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update())
+
+  // Set the forward and back projection filters to be used
+  using ConjugateGradientFilterType =
+    rtk::FourDConjugateGradientConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>;
+  auto fourdconjugategradient = ConjugateGradientFilterType::New();
+
+  fourdconjugategradient->SetInputVolumeSeries(data.InitialVolumeSeries);
+  fourdconjugategradient->SetNumberOfIterations(2);
+#ifdef RTK_USE_CUDA
+  fourdconjugategradient->SetCudaConjugateGradient(true);
+  fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_CUDARAYCAST);
+  fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_CUDAVOXELBASED);
+#else
+  fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_JOSEPH);
+  fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_VOXELBASED);
+#endif
+
+  // Set the newly ordered arguments
+  fourdconjugategradient->SetInputProjectionStack(reorder->GetOutput());
+  fourdconjugategradient->SetGeometry(reorder->GetOutputGeometry());
+  fourdconjugategradient->SetWeights(signalToInterpolationWeights->GetOutput());
+  fourdconjugategradient->SetSignal(reorder->GetOutputSignal());
+
+  auto verboseIterationCommand = rtk::VerboseIterationCommand<ConjugateGradientFilterType>::New();
+  fourdconjugategradient->AddObserver(itk::AnyEvent(), verboseIterationCommand);
+
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(fourdconjugategradient->Update())
+
+  // Write
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(fourdconjugategradient->GetOutput(), "fourdconjugategradient.mha"));
+
+  return EXIT_SUCCESS;
+}

examples/FourDConjugateGradient/README.md

@@ -0,0 +1,17 @@
+#   4D Conjugate Gradient
+
+FourDConjugateGradient shows how to perform iterative cone-beam CT reconstruction using either CPU or GPU resources.
+You can refer to the [projectors documentation](../../documentation/docs/Projectors.md) to see all options available for the back and forwardprojections.
+
+This example generates its own input data: a phantom made of two ellipsoids, one of which is moving. Projections are computed through this moving phantom.
+
+`````{tab-set}
+
+````{tab-item} C++
+
+```{literalinclude} ./FourDConjugateGradient.cxx
+:language: c++
+```
+````
+
+`````

test/CMakeLists.txt

@@ -376,6 +376,9 @@ rtk_add_test(rtkConjugateGradientExampleTest
   ../examples/ConjugateGradient/ConjugateGradient.cxx
   DATA{Input/Forbild/Thorax}
 )
+rtk_add_test(rtkFourDConjugateGradientExampleTest
+  ../examples/FourDConjugateGradient/FourDConjugateGradient.cxx
+)
 
 if(ITK_WRAP_PYTHON)
   itk_python_add_test(NAME rtkMaximumIntensityPythonTest COMMAND rtkMaximumIntensity.py)