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>

documentation/docs/ExternalData/FourDConjugateGradient.gif.sha512

@@ -0,0 +1 @@
+430dd919f1ad889d78a04cc70c32f1851b790f0423d6ee63ae3a9c626786a3fe9617a38b5655cecfb1890f56e6ffe298a01fe2d39210a7f7a16e539074b3f259

documentation/docs/ExternalData/FourDFDK.gif.sha512

@@ -0,0 +1 @@
+6db1332af81b0ac92d9fbbf3ad42ce8d3f7c839923e979cd4aa7bb1285e7b1f482221d77c7983547020317ab652b1b05d998916f0952f96c741b58d37a0b2c80

documentation/docs/ExternalData/FourDRooster.gif.sha512

@@ -0,0 +1 @@
+6db1332af81b0ac92d9fbbf3ad42ce8d3f7c839923e979cd4aa7bb1285e7b1f482221d77c7983547020317ab652b1b05d998916f0952f96c741b58d37a0b2c80

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/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()

examples/FourDConjugateGradient/README.md

@@ -0,0 +1,23 @@
+#   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).
+
+![img_4D](../../documentation/docs/ExternalData/FourDConjugateGradient.gif){w=200px alt="conjugate gradient 4D reconstruction"}
+
+`````{tab-set}
+
+````{tab-item} C++
+
+```{literalinclude} ./FourDConjugateGradient.cxx
+:language: c++
+```
+````
+````{tab-item} Python
+
+```{literalinclude} ./FourDConjugateGradient.py
+:language: python
+```
+````
+
+`````

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