Rtkfdk python

applications/rtk3Doutputimage_group.py

@@ -1,4 +1,5 @@
 import itk
+from itk import RTK as rtk
 
 __all__ = [
     'add_rtk3Doutputimage_group',
@@ -8,10 +9,10 @@
 # Mimicks rtk3Doutputimage_section.ggo
 def add_rtk3Doutputimage_group(parser):
   rtk3Doutputimage_group = parser.add_argument_group('Output 3D image properties')
-  rtk3Doutputimage_group.add_argument('--origin', help='Origin (default=centered)', type=float, nargs='+')
-  rtk3Doutputimage_group.add_argument('--dimension', help='Dimension', type=int, nargs='+', default=[256])
-  rtk3Doutputimage_group.add_argument('--spacing', help='Spacing', type=float, nargs='+', default=[1])
-  rtk3Doutputimage_group.add_argument('--direction', help='Direction', type=float, nargs='+')
+  rtk3Doutputimage_group.add_argument('--origin', help='Origin (default=centered)', type=rtk.comma_separated_args(float))
+  rtk3Doutputimage_group.add_argument('--dimension', help='Dimension', type=rtk.comma_separated_args(int), default=[256])
+  rtk3Doutputimage_group.add_argument('--spacing', help='Spacing', type=rtk.comma_separated_args(float), default=[1])
+  rtk3Doutputimage_group.add_argument('--direction', help='Direction', type=rtk.comma_separated_args(float))
   rtk3Doutputimage_group.add_argument('--like', help='Copy information from this image (origin, dimension, spacing, direction)')
 
 # Mimicks SetConstantImageSourceFromGgo

applications/rtkconjugategradient/rtkconjugategradient.py

@@ -4,19 +4,6 @@
 import itk
 from itk import RTK as rtk
 
-def GetCudaImageFromImage(img):
-  if hasattr(itk, 'CudaImage'):
-    # TODO: avoid this Update by implementing an itk::ImageToCudaImageFilter
-    img.Update()
-    cuda_img = itk.CudaImage[itk.itkCType.GetCTypeForDType(img.dtype),img.ndim].New()
-    cuda_img.SetPixelContainer(img.GetPixelContainer())
-    cuda_img.CopyInformation(img)
-    cuda_img.SetBufferedRegion(img.GetBufferedRegion())
-    cuda_img.SetRequestedRegion(img.GetRequestedRegion())
-    return cuda_img
-  else:
-    return img
-
 def main():
   # argument parsing
   parser = argparse.ArgumentParser(description=
@@ -106,12 +93,12 @@ def main():
   rtk.SetForwardProjectionFromArgParse(args_info, conjugategradient)
   rtk.SetBackProjectionFromArgParse(args_info, conjugategradient)
   rtk.SetIterationsReportFromArgParse(args_info, conjugategradient)
-  conjugategradient.SetInput(GetCudaImageFromImage(inputFilter.GetOutput()))
-  conjugategradient.SetInput(1, GetCudaImageFromImage(reader.GetOutput()))
-  conjugategradient.SetInput(2, GetCudaImageFromImage(weightsSource.GetOutput()))
+  conjugategradient.SetInput(rtk.CudaImageFromImage(inputFilter.GetOutput()))
+  conjugategradient.SetInput(1, rtk.CudaImageFromImage(reader.GetOutput()))
+  conjugategradient.SetInput(2, rtk.CudaImageFromImage(weightsSource.GetOutput()))
   conjugategradient.SetCudaConjugateGradient(not args_info.nocudacg)
   if args_info.mask is not None:
-    conjugategradient.SetSupportMask(GetCudaImageFromImage(supportmask))
+    conjugategradient.SetSupportMask(rtk.CudaImageFromImage(supportmask))
 
   if args_info.gamma is not None:
     conjugategradient.SetGamma(args_info.gamma)

applications/rtkfdk/rtkfdk.py

@@ -0,0 +1,165 @@
+#!/usr/bin/env python
+import argparse
+import sys
+import itk
+import math
+from itk import RTK as rtk
+
+def main():
+    # Argument parsing
+    parser = argparse.ArgumentParser(
+        description="Reconstructs a 3D volume from a sequence of projections [Feldkamp, David, Kress, 1984].",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
+
+    # General options
+    parser.add_argument('--verbose', '-v', help="Verbose execution", action='store_true')
+    parser.add_argument('--config', help="Config file", type=str)
+    parser.add_argument('--geometry', '-g', help="XML geometry file name", type=str, required=True)
+    parser.add_argument('--output', '-o', help="Output file name", type=str, required=True)
+    parser.add_argument('--hardware', help="Hardware used for computation", choices=['cpu', 'cuda'], default="cpu")
+    parser.add_argument('--lowmem', '-l', help="Load only one projection per thread in memory", action='store_true')
+    parser.add_argument('--divisions', '-d', help="Streaming option: number of stream divisions of the CT", type=int, default=1)
+    parser.add_argument('--subsetsize', help="Streaming option: number of projections processed at a time", type=int, default=16)
+    parser.add_argument('--nodisplaced', help="Disable the displaced detector filter", action='store_true')
+    parser.add_argument('--short', help="Minimum angular gap to detect a short scan (converted to radians)", type=float, default=20.0)
+
+    # Ramp filter options
+    ramp_group = parser.add_argument_group("Ramp filter")
+    ramp_group.add_argument('--pad', help="Data padding parameter to correct for truncation", type=float, default=0.0)
+    ramp_group.add_argument('--hann', help="Cut frequency for Hann window in ]0,1] (0.0 disables it)", type=float, default=0.0)
+    ramp_group.add_argument('--hannY', help="Cut frequency for Hann window in ]0,1] (0.0 disables it)", type=float, default=0.0)
+
+    # Motion compensation options
+    motion_group = parser.add_argument_group(
+        "Motion compensation described in [Rit et al, TMI, 2009] and [Rit et al, Med Phys, 2009]"
+    )
+    motion_group.add_argument('--signal', help="Signal file name", type=str)
+    motion_group.add_argument('--dvf', help="Input 4D DVF", type=str)
+
+    # RTK specific groups
+    rtk.add_rtkinputprojections_group(parser)
+    rtk.add_rtk3Doutputimage_group(parser)
+
+    # Parse the command line arguments
+    args_info = parser.parse_args()
+
+    # Convert `short` from degrees to radians
+    args_info.short = args_info.short * math.pi / 180.0
+
+    # Define output pixel type and dimension
+    OutputPixelType = itk.F
+    Dimension = 3
+    OutputImageType = itk.Image[OutputPixelType, Dimension]
+
+    # Projections reader
+    reader = rtk.ProjectionsReader[OutputImageType].New()
+    rtk.SetProjectionsReaderFromArgParse(reader, args_info)
+
+    # Geometry handling
+    if args_info.verbose:
+        print(f'Reading geometry information from {args_info.geometry}')
+
+    # Read geometry information from the specified XML file
+    geometryReader = rtk.ThreeDCircularProjectionGeometryXMLFileReader.New()
+    geometryReader.SetFilename(args_info.geometry)
+    geometryReader.GenerateOutputInformation()
+    geometry = geometryReader.GetOutputObject()
+
+    # Check on hardware parameter to determine appropriate filter usage
+    # CUDA classes are non-templated, so they do not require a type specification.
+    if args_info.hardware == "cuda":
+        ddf = rtk.CudaDisplacedDetectorImageFilter.New()
+        pssf = rtk.CudaParkerShortScanImageFilter.New()
+    else:
+        ddf = rtk.DisplacedDetectorForOffsetFieldOfViewImageFilter[OutputImageType].New()
+        pssf = rtk.ParkerShortScanImageFilter[OutputImageType].New()
+
+    # Set up the displaced detector filter with geometry and input
+    ddf.SetInput(rtk.GetCudaImageFromImage(reader.GetOutput()))
+    ddf.SetGeometry(geometry)
+    ddf.SetDisable(args_info.nodisplaced)
+
+    # Short scan image filter
+    pssf.SetInput(ddf.GetOutput())
+    pssf.SetGeometry(geometry)
+    pssf.InPlaceOff()
+    pssf.SetAngularGapThreshold(args_info.short)
+
+    # Create reconstructed image
+    constantImageSource = rtk.ConstantImageSource[OutputImageType].New()
+    rtk.SetConstantImageSourceFromArgParse(constantImageSource, args_info)
+
+    # Motion-compensated objects for the compensation of a cyclic deformation.
+    # Although these will only be used if the command line options for motion
+    # compensation are set, we still create the object beforehand to avoid auto
+    # destruction.
+    DVFPixelType = itk.Vector[itk.F, 3]
+    DVFImageType = itk.Image[DVFPixelType, 3]
+    DVFImageSequenceType = itk.Image[DVFPixelType, 4]
+
+    # Create the deformation filter and reader for DVF
+    deformation = rtk.CyclicDeformationImageFilter[DVFImageSequenceType, DVFImageType].New()
+    dvfReader = itk.ImageFileReader[DVFImageSequenceType].New()
+    deformation.SetInput(dvfReader.GetOutput())
+
+    # Set up the back projection filter for motion compensation
+    bp = rtk.FDKWarpBackProjectionImageFilter[OutputImageType, OutputImageType, type(deformation)].New()
+    bp.SetDeformation(deformation)
+    bp.SetGeometry(geometry)
+
+    # FDK reconstruction filtering
+    if (args_info.hardware == "cpu"):
+        feldkamp = rtk.FDKConeBeamReconstructionFilter[OutputImageType].New()
+    else :
+        feldkamp = rtk.CudaFDKConeBeamReconstructionFilter.New()
+
+    # Progress reporting
+    progressCommand = rtk.PercentageProgressCommand(feldkamp)
+    if args_info.verbose:
+        print("Registering progress observer...")
+        feldkamp.AddObserver(itk.ProgressEvent(), progressCommand.callback)  # Register the callback for progress
+        feldkamp.AddObserver(itk.EndEvent(), progressCommand.End)  # Register end notification
+
+    # Set inputs and options for the FDK filter
+    feldkamp.SetInput(0, rtk.GetCudaImageFromImage(constantImageSource.GetOutput()))
+    feldkamp.SetInput(1, pssf.GetOutput())
+    feldkamp.SetGeometry(geometry)
+    feldkamp.GetRampFilter().SetTruncationCorrection(args_info.pad)
+    feldkamp.GetRampFilter().SetHannCutFrequency(args_info.hann)
+    feldkamp.GetRampFilter().SetHannCutFrequencyY(args_info.hannY)
+    feldkamp.SetProjectionSubsetSize(args_info.subsetsize)
+
+
+    # Motion compensation settings
+    if args_info.signal and args_info.dvf:
+        if args_info.hardware == "cuda":
+            print("Motion compensation is not supported in CUDA. Aborting")
+            sys.exit(1)  # Exit if CUDA is selected with motion compensation
+        dvfReader.SetFileName(args_info.dvf)
+        deformation.SetSignalFilename(args_info.signal)
+        feldkamp.SetBackProjectionFilter(bp)
+
+    # Streaming depending on streaming capability of writer
+    streamerBP = itk.StreamingImageFilter[OutputImageType, OutputImageType].New()
+    streamerBP.SetInput(feldkamp.GetOutput())
+    streamerBP.SetNumberOfStreamDivisions(args_info.divisions)
+
+    # Create a splitter to control how the image region is divided during streaming
+    splitter = itk.ImageRegionSplitterDirection.New()
+    splitter.SetDirection(2)
+    streamerBP.SetRegionSplitter(splitter)
+
+    # Write the reconstructed output image
+    if args_info.verbose:
+        print("Reconstructing and writing...")
+
+    # Attempt to write the image and handle exceptions
+    try:
+        itk.imwrite(streamerBP.GetOutput(), args_info.output)
+    except Exception as e:
+        print(f"Error while writing image: {e}")
+        sys.exit(1)  # Exit the program with an error status
+
+if __name__ == '__main__':
+    main()

applications/rtkinputprojections_group.py

@@ -1,4 +1,5 @@
 import itk
+from itk import RTK as rtk
 
 __all__ = [
     'add_rtkinputprojections_group',
@@ -13,20 +14,20 @@ def add_rtkinputprojections_group(parser):
   rtkinputprojections_group.add_argument('--nsort', help='Numeric sort for regular expression matches', action='store_true')
   rtkinputprojections_group.add_argument('--submatch', help='Index of the submatch that will be used to sort matches', type=int, default=0)
   rtkinputprojections_group.add_argument('--nolineint', help='Disable raw to line integral conversion, just casts to float', type=bool, default=False)
-  rtkinputprojections_group.add_argument('--newdirection', help='New value of input projections (before pre-processing)', type=float, nargs='+')
-  rtkinputprojections_group.add_argument('--neworigin', help='New origin of input projections (before pre-processing)', type=float, nargs='+')
-  rtkinputprojections_group.add_argument('--newspacing', help='New spacing of input projections (before pre-processing)', type=float, nargs='+')
-  rtkinputprojections_group.add_argument('--lowercrop', help='Lower boundary crop size', type=int, nargs='+', default=[0])
-  rtkinputprojections_group.add_argument('--uppercrop', help='Upper boundary crop size', type=int, nargs='+', default=[0])
-  rtkinputprojections_group.add_argument('--binning', help='Shrink / Binning factos in each direction', type=int, nargs='+', default=[1])
-  rtkinputprojections_group.add_argument('--wpc', help='Water precorrection coefficients (default is no correction)', type=float, nargs='+')
+  rtkinputprojections_group.add_argument('--newdirection', help='New value of input projections (before pre-processing)', type=rtk.comma_separated_args(float))
+  rtkinputprojections_group.add_argument('--neworigin', help='New origin of input projections (before pre-processing)', type=rtk.comma_separated_args(float))
+  rtkinputprojections_group.add_argument('--newspacing', help='New spacing of input projections (before pre-processing)', type=rtk.comma_separated_args(float))
+  rtkinputprojections_group.add_argument('--lowercrop', help='Lower boundary crop size', type=rtk.comma_separated_args(int))
+  rtkinputprojections_group.add_argument('--uppercrop', help='Upper boundary crop size', type=rtk.comma_separated_args(int))
+  rtkinputprojections_group.add_argument('--binning', help='Shrink / Binning factos in each direction', type=rtk.comma_separated_args(int), default=[1])
+  rtkinputprojections_group.add_argument('--wpc', help='Water precorrection coefficients (default is no correction)', type=rtk.comma_separated_args(float))
   rtkinputprojections_group.add_argument('--spr', help='Boellaard scatter correction: scatter-to-primary ratio', type=float, default=0)
   rtkinputprojections_group.add_argument('--nonneg', help='Boellaard scatter correction: non-negativity threshold', type=float)
   rtkinputprojections_group.add_argument('--airthres', help='Boellaard scatter correction: air threshold', type=float)
   rtkinputprojections_group.add_argument('--i0', help='I0 value (when assumed constant per projection), 0 means auto', type=float)
   rtkinputprojections_group.add_argument('--idark', help='IDark value, i.e., value when beam is off', type=float, default=0)
   rtkinputprojections_group.add_argument('--component', help='Vector component to extract, for multi-material projections', type=int, default=0)
-  rtkinputprojections_group.add_argument('--radius', help='Radius of neighborhood for conditional median filtering', type=int, nargs='+', default=[0])
+  rtkinputprojections_group.add_argument('--radius', help='Radius of neighborhood for conditional median filtering', type=rtk.comma_separated_args(int), default=[0])
   rtkinputprojections_group.add_argument('--multiplier', help='Threshold multiplier for conditional median filtering', type=float, default=0)
 
 # Mimicks GetProjectionsFileNamesFromGgo

pyproject.toml

@@ -46,6 +46,7 @@ rtkelektasynergygeometry = "itk.rtkelektasynergygeometry:main"
 rtkorageometry = "itk.rtkorageometry:main"
 rtksimulatedgeometry = "itk.rtksimulatedgeometry:main"
 rtkvarianobigeometry = "itk.rtkvarianobigeometry:main"
+rtkfdk = "itk.rtkfdk:main"
 
 [project.urls]
 Download = "https://github.com/RTKConsortium/RTK"

wrapping/CMakeLists.txt

@@ -29,13 +29,16 @@ itk_end_wrap_module()
 
 wrap_itk_python_bindings_install(/itk "RTK"
     __init_rtk__.py
+
     ${RTK_SOURCE_DIR}/applications/rtk3Doutputimage_group.py
     ${RTK_SOURCE_DIR}/applications/rtkinputprojections_group.py
     ${RTK_SOURCE_DIR}/applications/rtkiterations_group.py
     ${RTK_SOURCE_DIR}/applications/rtkprojectors_group.py
     ${RTK_SOURCE_DIR}/applications/rtkconjugategradient/rtkconjugategradient.py
     ${RTK_SOURCE_DIR}/applications/rtkelektasynergygeometry/rtkelektasynergygeometry.py
+    ${RTK_SOURCE_DIR}/applications/rtkfdk/rtkfdk.py
     ${RTK_SOURCE_DIR}/applications/rtkorageometry/rtkorageometry.py
     ${RTK_SOURCE_DIR}/applications/rtksimulatedgeometry/rtksimulatedgeometry.py
     ${RTK_SOURCE_DIR}/applications/rtkvarianobigeometry/rtkvarianobigeometry.py
+    ${RTK_SOURCE_DIR}/wrapping/rtkExtra.py
 )

wrapping/rtkExtra.py

@@ -0,0 +1,48 @@
+import itk
+from itk import RTK as rtk
+
+# Write a 3D circular projection geometry to a file.
+def WriteGeometry(geometry, filename):
+    # Create and configure the writer
+    writer = rtk.ThreeDCircularProjectionGeometryXMLFileWriter.New()
+    writer.SetObject(geometry)
+    writer.SetFilename(filename)
+
+    # Write the geometry to file
+    writer.WriteFile()
+
+# Read a 3D circular projection geometry from a file.
+def ReadGeometry(filename):
+    # Create and configure the reader
+    reader = rtk.ThreeDCircularProjectionGeometryXMLFileReader.New()
+    reader.SetFilename(filename)
+    reader.GenerateOutputInformation()
+
+    # Return the geometry object
+    return reader.GetOutputObject()
+
+# Convert an ITK image to a CUDA image if CUDA is available.
+def CudaImageFromImage(img):
+    if hasattr(itk, 'CudaImage'):
+        img.Update()  # Ensure the image is up to date
+        cuda_img = itk.CudaImage[itk.itkCType.GetCTypeForDType(img.dtype), img.ndim].New()
+        cuda_img.SetPixelContainer(img.GetPixelContainer())
+        cuda_img.CopyInformation(img)
+        cuda_img.SetBufferedRegion(img.GetBufferedRegion())
+        cuda_img.SetRequestedRegion(img.GetRequestedRegion())
+        return cuda_img
+    return img  # Return the original image if CUDA is not available
+
+class PercentageProgressCommand:
+    def __init__(self,caller):
+        self.percentage = -1
+        self.caller=caller
+
+    def callback(self):
+        new_percentage = int(self.caller.GetProgress() * 100)
+        if new_percentage > self.percentage:
+            print(f"\r{self.caller.GetNameOfClass()} {new_percentage}% completed.", end='', flush=True)
+            self.percentage = new_percentage
+
+    def End(self):
+        print()  # Print newline when execution ends