4D applications

applications/rtk3Doutputimage_group.py

@@ -46,7 +46,7 @@ def SetConstantImageSourceFromArgParse(source, args_info):
     ImageType = type(source.GetOutput())
 
     # Handle deprecated --dimension argument
-    if args_info.dimension is not None and args_info.Size is None:
+    if args_info.dimension is not None and args_info.size is None:
         print(
             "Warning: '--dimension' is deprecated and will be removed in a future release. "
             "Please use '--size' instead."

applications/rtk4Doutputimage_group.py

@@ -0,0 +1,26 @@
+def add_rtk4Doutputimage_group(parser):
+    group = parser.add_argument_group("Output 4D image properties")
+    group.add_argument(
+        "--origin", help="Origin (default=centered)", type=float, nargs="+"
+    )
+    group.add_argument(
+        "--dimension",
+        help="Dimension (Deprecated) Use --size instead.",
+        type=int,
+        nargs="+",
+        default=[256],
+    )
+    group.add_argument("--size", help="Size", type=int, nargs="+", default=[256])
+    group.add_argument("--spacing", help="Spacing", type=float, nargs="+", default=[1])
+    group.add_argument("--direction", help="Direction", type=float, nargs="+")
+    group.add_argument(
+        "--frames",
+        help="Number of time frames in the 4D reconstruction",
+        type=int,
+        default=10,
+    )
+    group.add_argument(
+        "--like",
+        help="Copy information from this image (origin, size, spacing, direction, frames)",
+        type=str,
+    )

applications/rtkelektasynergygeometry/rtkelektasynergygeometry.py

@@ -4,26 +4,55 @@
 
 
 def build_parser():
-    # Argument parsing
     parser = rtk.RTKArgumentParser(
         description="Creates an RTK geometry file from an Elekta Synergy acquisition."
     )
 
-    parser.add_argument(
-        "--xml", "-x", help="XML file name (starting with XVI5)", required=True
-    )
+    parser.add_argument("--image-db", "-i", help="Image table filename (prior to XVI5)")
+    parser.add_argument("--frame-db", "-f", help="Frame table filename (prior to XVI5)")
+    parser.add_argument("--xml", "-x", help="XML file name (starting with XVI5)")
+    parser.add_argument("--dicom-uid", "-u", help="Dicom uid of the acquisition")
     parser.add_argument("--output", "-o", help="Output file name", required=True)
 
     return parser
 
 
 def process(args_info: argparse.Namespace):
 
-    reader = rtk.ElektaXVI5GeometryXMLFileReader.New()
-    reader.SetFilename(args_info.xml)
-    reader.GenerateOutputInformation()
+    if (
+        args_info.image_db is not None
+        and args_info.frame_db is not None
+        and args_info.dicom_uid is not None
+        and args_info.xml is None
+    ):
+        if args_info.verbose:
+            print("Building geometry from DBF tables...")
+        reader = rtk.ElektaSynergyGeometryReader.New()
+        reader.SetDicomUID(args_info.dicom_uid)
+        reader.SetImageDbfFileName(args_info.image_db)
+        reader.SetFrameDbfFileName(args_info.frame_db)
+        reader.UpdateOutputData()
+        geometry = reader.GetGeometry()
+    elif (
+        args_info.image_db
+        and args_info.frame_db
+        and args_info.dicom_uid
+        and args_info.xml is not None
+    ):
+        if args_info.verbose:
+            print(f"Reading geometry information from {args_info.xml}...")
+        geometryReader = rtk.ElektaXVI5GeometryXMLFileReader.New()
+        geometryReader.SetFilename(args_info.xml)
+        geometryReader.GenerateOutputInformation()
+        geometry = geometryReader.GetGeometry()
+    else:
+        raise ValueError(
+            "You must either provide image_db, frame_db and dicom_uid for versions up to v4 or xml starting with v5."
+        )
 
-    rtk.write_geometry(reader.GetGeometry(), args_info.output)
+    if args_info.verbose:
+        print(f"Writing geometry to {args_info.output}...")
+    rtk.write_geometry(geometry, args_info.output)
 
 
 def main(argv=None):

applications/rtkfourdconjugategradient/rtkfourdconjugategradient.py

@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+import argparse
+import itk
+from itk import RTK as rtk
+import numpy as np
+
+
+def build_parser():
+    parser = rtk.RTKArgumentParser(
+        description="Reconstructs a 3D + time sequence of volumes from a projection stack and a respiratory/cardiac signal, with a conjugate gradient technique."
+    )
+    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(
+        "--niterations", "-n", help="Number of iterations", type=int, default=5
+    )
+    parser.add_argument(
+        "--cudacg",
+        help="Perform conjugate gradient calculations on GPU",
+        action="store_true",
+    )
+    parser.add_argument("--input", "-i", help="Input volume", type=str)
+    parser.add_argument(
+        "--nodisplaced",
+        help="Disable the displaced detector filter",
+        action="store_true",
+    )
+    parser.add_argument(
+        "--signal",
+        help="File containing the phase of each projection",
+        type=str,
+        required=True,
+    )
+    rtk.add_rtkinputprojections_group(parser)
+    rtk.add_rtkprojectors_group(parser)
+    rtk.add_rtkiterations_group(parser)
+    rtk.add_rtk4Doutputimage_group(parser)
+    return parser
+
+
+def process(args_info: argparse.Namespace):
+    OutputPixelType = itk.F
+    VolumeSeriesType = itk.Image[OutputPixelType, 4]
+    ProjectionStackType = itk.Image[OutputPixelType, 3]
+
+    # Projections reader
+    reader = rtk.ProjectionsReader[ProjectionStackType].New()
+    rtk.SetProjectionsReaderFromArgParse(reader, args_info)
+    reader.UpdateLargestPossibleRegion()
+
+    # Geometry
+    if args_info.verbose:
+        print(f"Reading geometry information from {args_info.geometry}...")
+    geometry = rtk.read_geometry(args_info.geometry)
+
+    # Create input: either an existing volume read from a file or a blank image
+    if args_info.input:
+        # Read an existing image to initialize the volume
+        inputReader = itk.ImageFileReader[VolumeSeriesType].New()
+        inputReader.SetFileName(args_info.input)
+        inputFilter = inputReader
+    else:
+        constantImageSource = rtk.ConstantImageSource[VolumeSeriesType].New()
+        rtk.SetConstantImageSourceFromArgParse(constantImageSource, args_info)
+        # Set 4th dimension (frames)
+        size = constantImageSource.GetSize()
+        size[3] = args_info.frames
+        constantImageSource.SetSize(size)
+        inputFilter = constantImageSource
+
+    inputFilter.Update()
+    inputFilter.ReleaseDataFlagOn()
+
+    # Reorder projections
+    signal = np.loadtxt(args_info.signal).tolist()
+    reorder = rtk.ReorderProjectionsImageFilter[
+        ProjectionStackType, ProjectionStackType
+    ].New()
+    reorder.SetInput(reader.GetOutput())
+    reorder.SetInputGeometry(geometry)
+    reorder.SetInputSignal(signal)
+    reorder.Update()
+    # Release the memory holding the stack of original projections
+    reader.GetOutput().ReleaseData()
+
+    # Interpolation weights
+    signalToInterpolationWeights = rtk.SignalToInterpolationWeights.New()
+    signalToInterpolationWeights.SetSignal(reorder.GetOutputSignal())
+    signalToInterpolationWeights.SetNumberOfReconstructedFrames(
+        inputFilter.GetOutput().GetLargestPossibleRegion().GetSize()[3]
+    )
+    signalToInterpolationWeights.Update()
+
+    if args_info.cudacg and not hasattr(itk, "CudaImage"):
+        raise RuntimeError(
+            "CUDA conjugate gradient requested but RTK was not built with CUDA support."
+        )
+
+    if args_info.cudacg:
+        CudaVolumeSeriesType = itk.CudaImage[OutputPixelType, 4]
+        CudaProjectionStackType = itk.CudaImage[OutputPixelType, 3]
+        # Conjugate gradient filter
+        conjugategradient = rtk.FourDConjugateGradientConeBeamReconstructionFilter[
+            CudaVolumeSeriesType, CudaProjectionStackType
+        ].New()
+        conjugategradient.SetInputVolumeSeries(
+            itk.cuda_image_from_image(inputFilter.GetOutput())
+        )
+        conjugategradient.SetWeights(signalToInterpolationWeights.GetOutput())
+        conjugategradient.SetInputProjectionStack(
+            itk.cuda_image_from_image(reorder.GetOutput())
+        )
+    else:
+        # Conjugate gradient filter
+        conjugategradient = rtk.FourDConjugateGradientConeBeamReconstructionFilter[
+            VolumeSeriesType, ProjectionStackType
+        ].New()
+        conjugategradient.SetInputVolumeSeries(inputFilter.GetOutput())
+        conjugategradient.SetWeights(signalToInterpolationWeights.GetOutput())
+        conjugategradient.SetInputProjectionStack(reorder.GetOutput())
+
+    rtk.SetForwardProjectionFromArgParse(args_info, conjugategradient)
+    rtk.SetBackProjectionFromArgParse(args_info, conjugategradient)
+    conjugategradient.SetNumberOfIterations(args_info.niterations)
+    conjugategradient.SetCudaConjugateGradient(args_info.cudacg)
+    conjugategradient.SetDisableDisplacedDetectorFilter(args_info.nodisplaced)
+    conjugategradient.SetGeometry(reorder.GetOutputGeometry())
+    conjugategradient.SetSignal(reorder.GetOutputSignal())
+
+    rtk.SetIterationsReportFromArgParse(args_info, conjugategradient)
+
+    if args_info.verbose:
+        print("Running 4D conjugate gradient reconstruction...")
+    conjugategradient.Update()
+
+    if args_info.verbose:
+        print(f"Writing output to {args_info.output}...")
+    writer = itk.ImageFileWriter[VolumeSeriesType].New()
+    writer.SetFileName(args_info.output)
+    writer.SetInput(conjugategradient.GetOutput())
+    writer.Update()
+
+
+def main(argv=None):
+    parser = build_parser()
+    args_info = parser.parse_args(argv)
+    process(args_info)
+
+
+if __name__ == "__main__":
+    main()