Examples 3D+time

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,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})

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 <itkCudaImage.h>
+#endif
+#include <itkImageFileWriter.h>
+
+int
+main(int argc, char * argv[])
+{
+  using OutputPixelType = float;
+
+#ifdef RTK_USE_CUDA
+  using VolumeSeriesType = itk::CudaImage<OutputPixelType, 4>;
+  using ProjectionStackType = itk::CudaImage<OutputPixelType, 3>;
+  using VolumeType = itk::CudaImage<OutputPixelType, 3>;
+#else
+  using VolumeSeriesType = itk::Image<OutputPixelType, 4>;
+  using ProjectionStackType = itk::Image<OutputPixelType, 3>;
+  using VolumeType = itk::Image<OutputPixelType, 3>;
+#endif
+
+  // Generate the input volume series, used as initial estimate by 4D conjugate gradient
+  auto fourDOrigin = itk::MakePoint(-63., -31., -63., 0.);
+  auto fourDSpacing = itk::MakeVector(4., 4., 4., 1.);
+  auto fourDSize = itk::MakeSize(32, 16, 32, 8);
+
+  using ConstantFourDSourceType = rtk::ConstantImageSource<VolumeSeriesType>;
+  auto fourDSource = ConstantFourDSourceType::New();
+
+  fourDSource->SetOrigin(fourDOrigin);
+  fourDSource->SetSpacing(fourDSpacing);
+  fourDSource->SetSize(fourDSize);
+  fourDSource->Update();
+
+  // Read geometry, projections and signal
+  rtk::ThreeDCircularProjectionGeometry::Pointer geometry;
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml"));
+
+  using ReaderType = rtk::ProjectionsReader<ProjectionStackType>;
+  auto                     projectionsReader = ReaderType::New();
+  std::vector<std::string> fileNames = std::vector<std::string>();
+  fileNames.push_back("four_d_projections.mha");
+  projectionsReader->SetFileNames(fileNames);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update());
+
+  // Re-order geometry and projections
+  // In the new order, projections with identical phases are packed together
+  std::vector<double> signal = rtk::ReadSignalFile("four_d_signal.txt");
+  auto                reorder = rtk::ReorderProjectionsImageFilter<ProjectionStackType>::New();
+  reorder->SetInput(projectionsReader->GetOutput());
+  reorder->SetInputGeometry(geometry);
+  reorder->SetInputSignal(signal);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(reorder->Update())
+
+  // Release the memory holding the stack of original projections
+  projectionsReader->GetOutput()->ReleaseData();
+
+  // Compute the interpolation weights
+  auto signalToInterpolationWeights = rtk::SignalToInterpolationWeights::New();
+  signalToInterpolationWeights->SetSignal(reorder->GetOutputSignal());
+  signalToInterpolationWeights->SetNumberOfReconstructedFrames(fourDSize[3]);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(signalToInterpolationWeights->Update())
+
+  // Set the forward and back projection filters to be used
+  using ConjugateGradientFilterType =
+    rtk::FourDConjugateGradientConeBeamReconstructionFilter<VolumeSeriesType, ProjectionStackType>;
+  auto fourdconjugategradient = ConjugateGradientFilterType::New();
+
+  fourdconjugategradient->SetInputVolumeSeries(fourDSource->GetOutput());
+  fourdconjugategradient->SetNumberOfIterations(2);
+#ifdef RTK_USE_CUDA
+  fourdconjugategradient->SetCudaConjugateGradient(true);
+  fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_CUDARAYCAST);
+  fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_CUDAVOXELBASED);
+#else
+  fourdconjugategradient->SetForwardProjectionFilter(ConjugateGradientFilterType::FP_JOSEPH);
+  fourdconjugategradient->SetBackProjectionFilter(ConjugateGradientFilterType::BP_VOXELBASED);
+#endif
+
+  // Set the newly ordered arguments
+  fourdconjugategradient->SetInputProjectionStack(reorder->GetOutput());
+  fourdconjugategradient->SetGeometry(reorder->GetOutputGeometry());
+  fourdconjugategradient->SetWeights(signalToInterpolationWeights->GetOutput());
+  fourdconjugategradient->SetSignal(reorder->GetOutputSignal());
+
+  auto verboseIterationCommand = rtk::VerboseIterationCommand<ConjugateGradientFilterType>::New();
+  fourdconjugategradient->AddObserver(itk::AnyEvent(), verboseIterationCommand);
+
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(fourdconjugategradient->Update())
+
+  // Write
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(fourdconjugategradient->GetOutput(), "fourdconjugategradient.mha"));
+
+  return EXIT_SUCCESS;
+}

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++
+```
+````
+
+`````

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})

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 <itkCudaImage.h>
+#endif
+#include <itkImageFileWriter.h>
+
+int
+main(int argc, char * argv[])
+{
+  using OutputPixelType = float;
+  constexpr unsigned int Dimension = 3;
+
+  using CPUOutputImageType = itk::Image<OutputPixelType, Dimension>;
+#ifdef RTK_USE_CUDA
+  using OutputImageType = itk::CudaImage<OutputPixelType, Dimension>;
+#else
+  using OutputImageType = CPUOutputImageType;
+#endif
+
+  // Read geometry, projections and signal
+  rtk::ThreeDCircularProjectionGeometry::Pointer geometry;
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("four_d_geometry.xml"));
+
+  using ReaderType = rtk::ProjectionsReader<OutputImageType>;
+  auto                     projectionsReader = ReaderType::New();
+  std::vector<std::string> fileNames = std::vector<std::string>();
+  fileNames.push_back("four_d_projections.mha");
+  projectionsReader->SetFileNames(fileNames);
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(projectionsReader->Update());
+
+  // Part specific to 4D
+  auto selector = rtk::SelectOneProjectionPerCycleImageFilter<OutputImageType>::New();
+  selector->SetInput(projectionsReader->GetOutput());
+  selector->SetInputGeometry(geometry);
+  selector->SetSignalFilename("four_d_signal.txt");
+
+  // Create one frame of the reconstructed image
+  using ConstantImageSourceType = rtk::ConstantImageSource<OutputImageType>;
+  auto constantImageSource = ConstantImageSourceType::New();
+  constantImageSource->SetOrigin(itk::MakePoint(-63., -31., -63.));
+  constantImageSource->SetSpacing(itk::MakeVector(4., 4., 4.));
+  constantImageSource->SetSize(itk::MakeSize(32, 16, 32));
+  constantImageSource->SetConstant(0.);
+  constantImageSource->Update();
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(constantImageSource->Update())
+
+  // FDK reconstruction filtering
+#ifdef RTK_USE_CUDA
+  using FDKType = rtk::CudaFDKConeBeamReconstructionFilter;
+#else
+  using FDKType = rtk::FDKConeBeamReconstructionFilter<OutputImageType>;
+#endif
+  auto feldkamp = FDKType::New();
+  feldkamp->SetInput(0, constantImageSource->GetOutput());
+  feldkamp->SetInput(1, selector->GetOutput());
+  feldkamp->SetGeometry(selector->GetOutputGeometry());
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->Update());
+
+  // Create empty 4D image
+  using FourDOutputImageType = itk::Image<OutputPixelType, Dimension + 1>;
+  using ConstantFourDSourceType = rtk::ConstantImageSource<FourDOutputImageType>;
+  auto fourDSource = ConstantFourDSourceType::New();
+  fourDSource->SetOrigin(itk::MakePoint(-63., -31., -63., 0.));
+  fourDSource->SetSpacing(itk::MakeVector(4., 4., 4., 1.));
+  fourDSource->SetSize(itk::MakeSize(32, 16, 32, 4));
+  fourDSource->SetConstant(0.);
+  fourDSource->Update();
+
+  // Go over each frame, reconstruct 3D frame and paste with iterators in 4D image
+  for (int f = 0; f < 4; f++)
+  {
+    selector->SetPhase(f / (double)4);
+    TRY_AND_EXIT_ON_ITK_EXCEPTION(feldkamp->UpdateLargestPossibleRegion())
+
+    ConstantFourDSourceType::OutputImageRegionType region;
+    region = fourDSource->GetOutput()->GetLargestPossibleRegion();
+    region.SetIndex(3, f);
+    region.SetSize(3, 1);
+
+    itk::ImageRegionIterator<FourDOutputImageType> it4D(fourDSource->GetOutput(), region);
+    itk::ImageRegionIterator<CPUOutputImageType>   it3D(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;
+}

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++
+```
+````
+
+`````

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})