From 2598cc91cc76fd0664560f5f750a79d78d60d13d Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 12 Nov 2025 10:28:46 +0100 Subject: [PATCH 1/7] BUG: Corrected variable name from 'Size' to 'size' --- applications/rtk3Doutputimage_group.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/rtk3Doutputimage_group.py b/applications/rtk3Doutputimage_group.py index fa506023e..1b4853801 100644 --- a/applications/rtk3Doutputimage_group.py +++ b/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." From 4538ff2f4feca2f8e25b75b1d4adb2ba7659b1fc Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 12 Nov 2025 10:33:09 +0100 Subject: [PATCH 2/7] ENH: Correct wrong version of rtkelektasynergygeometry The version of rtkelektasynergygeometry added in commit d0b215f93cccc90d245aefd708da8a0850a5b208 was incomplete --- .../rtkelektasynergygeometry.py | 45 +++++++++++++++---- 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/applications/rtkelektasynergygeometry/rtkelektasynergygeometry.py b/applications/rtkelektasynergygeometry/rtkelektasynergygeometry.py index 6de5e1aa9..83422d9f7 100755 --- a/applications/rtkelektasynergygeometry/rtkelektasynergygeometry.py +++ b/applications/rtkelektasynergygeometry/rtkelektasynergygeometry.py @@ -4,14 +4,14 @@ 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 @@ -19,11 +19,40 @@ def build_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): From efbd550c1c323fd6e0dd03f9413e99be32c9396e Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 12 Nov 2025 10:40:35 +0100 Subject: [PATCH 3/7] ENH: Add rtk4Doutputimage_group --- applications/rtk4Doutputimage_group.py | 26 ++++++++++++++++++++++++++ wrapping/__init_rtk__.py | 2 +- 2 files changed, 27 insertions(+), 1 deletion(-) create mode 100644 applications/rtk4Doutputimage_group.py diff --git a/applications/rtk4Doutputimage_group.py b/applications/rtk4Doutputimage_group.py new file mode 100644 index 000000000..6daf21e37 --- /dev/null +++ b/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, + ) diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index fa3befc93..5f5919a01 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -1,6 +1,5 @@ import sys import importlib -import os itk_module = sys.modules["itk"] rtk_module = getattr(itk_module, "RTK") @@ -11,6 +10,7 @@ "itk.rtkargumentparser", "itk.rtkinputprojections_group", "itk.rtk3Doutputimage_group", + "itk.rtk4Doutputimage_group", "itk.rtkprojectors_group", "itk.rtkiterations_group", "itk.rtkExtras", From 05e87797edbbf45726b01be5e9b05904006e20c9 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 12 Nov 2025 10:41:50 +0100 Subject: [PATCH 4/7] ENH: Add rtkfourdconjugategradient python application --- .../rtkfourdconjugategradient.py | 155 ++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_rtk__.py | 1 + wrapping/itkCSVFileReaderBaseRTK.wrap | 1 + wrapping/rtkSignalToInterpolationWeights.wrap | 1 + 5 files changed, 159 insertions(+) create mode 100644 applications/rtkfourdconjugategradient/rtkfourdconjugategradient.py create mode 100644 wrapping/itkCSVFileReaderBaseRTK.wrap create mode 100644 wrapping/rtkSignalToInterpolationWeights.wrap diff --git a/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.py b/applications/rtkfourdconjugategradient/rtkfourdconjugategradient.py new file mode 100644 index 000000000..f25c85a07 --- /dev/null +++ b/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() diff --git a/pyproject.toml b/pyproject.toml index 7db108f2f..d0112dc71 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,6 +56,7 @@ rtkextractphasesignal = "itk.rtkextractphasesignal:main" rtkfdk = "itk.rtkfdk:main" rtkfieldofview = "itk.rtkfieldofview:main" rtkforwardprojections = "itk.rtkforwardprojections:main" +rtkfourdconjugategradient = "itk.rtkfourdconjugategradient:main" rtki0estimation = "itk.rtki0estimation:main" rtkimagexgeometry = "itk.rtkimagxgeometry:main" rtkmaskcollimation = "itk.rtkmaskcollimation:main" diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index 5f5919a01..22957a734 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -38,6 +38,7 @@ "rtkfdk", "rtkfieldofview", "rtkforwardprojections", + "rtkfourdconjugategradient", "rtki0estimation", "rtkimagxgeometry", "rtkmaskcollimation", diff --git a/wrapping/itkCSVFileReaderBaseRTK.wrap b/wrapping/itkCSVFileReaderBaseRTK.wrap new file mode 100644 index 000000000..d6d190505 --- /dev/null +++ b/wrapping/itkCSVFileReaderBaseRTK.wrap @@ -0,0 +1 @@ +itk_wrap_simple_class("itk::CSVFileReaderBase" POINTER) diff --git a/wrapping/rtkSignalToInterpolationWeights.wrap b/wrapping/rtkSignalToInterpolationWeights.wrap new file mode 100644 index 000000000..67ffde3e3 --- /dev/null +++ b/wrapping/rtkSignalToInterpolationWeights.wrap @@ -0,0 +1 @@ +itk_wrap_simple_class("rtk::SignalToInterpolationWeights" POINTER) From eafea844d3694a3428b8ac400882e893bb14de2e Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 14 Nov 2025 14:37:52 +0100 Subject: [PATCH 5/7] ENH: Add rtkfourdfdk python application --- applications/rtkfourdfdk/rtkfourdfdk.py | 192 ++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_rtk__.py | 1 + ...electOneProjectionPerCycleImageFilter.wrap | 8 + wrapping/rtkSubSelectImageFilter.wrap | 6 + 5 files changed, 208 insertions(+) create mode 100644 applications/rtkfourdfdk/rtkfourdfdk.py diff --git a/applications/rtkfourdfdk/rtkfourdfdk.py b/applications/rtkfourdfdk/rtkfourdfdk.py new file mode 100644 index 000000000..2b5d59f35 --- /dev/null +++ b/applications/rtkfourdfdk/rtkfourdfdk.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python +import argparse +import sys +import itk +from itk import RTK as rtk + + +def build_parser(): + parser = rtk.RTKArgumentParser( + description="Reconstructs a 3D + time (4D) sequence using FDK with one projection per respiratory cycle in each frame." + ) + 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( + "--signal", + help="File containing the phase of each projection", + type=str, + required=True, + ) + parser.add_argument( + "--hardware", + help="Hardware used for computation", + choices=["cpu", "cuda"], + default="cpu", + ) + parser.add_argument( + "--divisions", + "-d", + help="Streaming option: number of stream divisions of the 3D reconstruction per frame", + 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", + ) + # Ramp filter + ramp = parser.add_argument_group("Ramp filter") + ramp.add_argument( + "--pad", + help="Data padding parameter to correct for truncation", + type=float, + default=0.0, + ) + ramp.add_argument( + "--hann", + help="Cut frequency for Hann window in ]0,1] (0.0 disables it)", + type=float, + default=0.0, + ) + ramp.add_argument( + "--hannY", + help="Cut frequency for Hann window in ]0,1] (0.0 disables it)", + type=float, + default=0.0, + ) + + # RTK common groups + rtk.add_rtkinputprojections_group(parser) + rtk.add_rtk4Doutputimage_group(parser) + return parser + + +def process(args_info: argparse.Namespace): + OutputPixelType = itk.F + OutputImageType = itk.Image[OutputPixelType, 3] + + # Check on hardware parameter + if args_info.hardware == "cuda" and not hasattr(itk, "CudaImage"): + print("The program has not been compiled with cuda option") + sys.exit(1) + + # Projections reader + reader = rtk.ProjectionsReader[OutputImageType].New() + rtk.SetProjectionsReaderFromArgParse(reader, args_info) + + # Geometry + if args_info.verbose: + print(f"Reading geometry information from {args_info.geometry}...") + geometry = rtk.read_geometry(args_info.geometry) + + # Part specific to 4D + selector = rtk.SelectOneProjectionPerCycleImageFilter[OutputImageType].New() + selector.SetInput(reader.GetOutput()) + selector.SetInputGeometry(geometry) + selector.SetSignalFilename(args_info.signal) + + # Displaced detector weighting + if args_info.hardware == "cuda": + ddf = rtk.CudaDisplacedDetectorImageFilter.New() + ddf.SetInput(itk.cuda_image_from_image(selector.GetOutput())) + else: + ddf = rtk.DisplacedDetectorForOffsetFieldOfViewImageFilter[ + OutputImageType + ].New() + ddf.SetInput(selector.GetOutput()) + ddf.SetGeometry(selector.GetOutputGeometry()) + ddf.SetDisable(args_info.nodisplaced) + + # Short scan image filter + if args_info.hardware == "cuda": + pssf = rtk.CudaParkerShortScanImageFilter.New() + else: + pssf = rtk.ParkerShortScanImageFilter[OutputImageType].New() + pssf.SetInput(ddf.GetOutput()) + pssf.SetGeometry(selector.GetOutputGeometry()) + pssf.InPlaceOff() + + # Create one frame of the reconstructed image + constantImageSource = rtk.ConstantImageSource[OutputImageType].New() + rtk.SetConstantImageSourceFromArgParse(constantImageSource, args_info) + constantImageSource.Update() + + # FDK reconstruction filtering + if args_info.hardware == "cuda": + feldkamp = rtk.CudaFDKConeBeamReconstructionFilter.New() + feldkamp.SetInput(0, itk.cuda_image_from_image(constantImageSource.GetOutput())) + feldkamp.SetInput(1, pssf.GetOutput()) + else: + feldkamp = rtk.FDKConeBeamReconstructionFilter[OutputImageType].New() + feldkamp.SetInput(0, constantImageSource.GetOutput()) + feldkamp.SetInput(1, pssf.GetOutput()) + feldkamp.SetGeometry(selector.GetOutputGeometry()) + feldkamp.GetRampFilter().SetTruncationCorrection(args_info.pad) + feldkamp.GetRampFilter().SetHannCutFrequency(args_info.hann) + feldkamp.GetRampFilter().SetHannCutFrequencyY(args_info.hannY) + feldkamp.SetProjectionSubsetSize(args_info.subsetsize) + + # Streaming depending on streaming capability of writer + streamerBP = itk.StreamingImageFilter[OutputImageType, OutputImageType].New() + streamerBP.SetInput(feldkamp.GetOutput()) + streamerBP.SetNumberOfStreamDivisions(args_info.divisions) + # Prefer streaming splits along Z for cache friendliness (like rtkfdk) + splitter = itk.ImageRegionSplitterDirection.New() + splitter.SetDirection(2) + streamerBP.SetRegionSplitter(splitter) + # Initialize meta-information once + streamerBP.UpdateOutputInformation() + + # Create empty 4D image + fourDConstantImageSource = rtk.ConstantImageSource[ + itk.Image[OutputPixelType, 4] + ].New() + rtk.SetConstantImageSourceFromArgParse(fourDConstantImageSource, args_info) + fourDInputSize = fourDConstantImageSource.GetSize() + + fourDInputSize[3] = args_info.frames + fourDConstantImageSource.SetSize(fourDInputSize) + fourDConstantImageSource.Update() + + # Go over each frame, reconstruct 3D frame and paste with NumPy views into the 4D image + volOut = fourDConstantImageSource.GetOutput() + arr4d = itk.GetArrayViewFromImage(volOut) + num_frames = args_info.frames + for f in range(num_frames): + if args_info.verbose: + print(f"Reconstructing frame #{f}...") + selector.SetPhase(f / float(num_frames)) + streamerBP.UpdateLargestPossibleRegion() + + arr3d = itk.GetArrayViewFromImage(streamerBP.GetOutput()) + arr4d[f] = arr3d + + # Write + if args_info.verbose: + print(f"Writing output to {args_info.output}...") + writer = itk.ImageFileWriter[itk.Image[OutputPixelType, 4]].New() + writer.SetFileName(args_info.output) + writer.SetInput(volOut) + writer.SetNumberOfStreamDivisions(args_info.divisions) + writer.Update() + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index d0112dc71..7fc9eaa10 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,6 +57,7 @@ rtkfdk = "itk.rtkfdk:main" rtkfieldofview = "itk.rtkfieldofview:main" rtkforwardprojections = "itk.rtkforwardprojections:main" rtkfourdconjugategradient = "itk.rtkfourdconjugategradient:main" +rtkfourdfdk = "itk.rtkfourdfdk:main" rtki0estimation = "itk.rtki0estimation:main" rtkimagexgeometry = "itk.rtkimagxgeometry:main" rtkmaskcollimation = "itk.rtkmaskcollimation:main" diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index 22957a734..09ef9400d 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -39,6 +39,7 @@ "rtkfieldofview", "rtkforwardprojections", "rtkfourdconjugategradient", + "rtkfourdfdk", "rtki0estimation", "rtkimagxgeometry", "rtkmaskcollimation", diff --git a/wrapping/rtkSelectOneProjectionPerCycleImageFilter.wrap b/wrapping/rtkSelectOneProjectionPerCycleImageFilter.wrap index ab48905e2..70b1d943d 100644 --- a/wrapping/rtkSelectOneProjectionPerCycleImageFilter.wrap +++ b/wrapping/rtkSelectOneProjectionPerCycleImageFilter.wrap @@ -1,3 +1,11 @@ itk_wrap_class("rtk::SelectOneProjectionPerCycleImageFilter" POINTER) + # Ensure SWIG sees base classes for inheritance to avoid 'Nothing known about base class' warnings + itk_wrap_include(rtkSubSelectImageFilter.h) + itk_wrap_include(itkCSVFileReaderBase.h) itk_wrap_image_filter("${WRAP_ITK_REAL}" 1 3) + + if(RTK_USE_CUDA) + itk_wrap_include(itkCudaImage.h) + itk_wrap_template("CIF3" "itk::CudaImage") + endif() itk_end_wrap_class() diff --git a/wrapping/rtkSubSelectImageFilter.wrap b/wrapping/rtkSubSelectImageFilter.wrap index da98a0129..b8ae93d18 100644 --- a/wrapping/rtkSubSelectImageFilter.wrap +++ b/wrapping/rtkSubSelectImageFilter.wrap @@ -1,3 +1,9 @@ itk_wrap_class("rtk::SubSelectImageFilter" POINTER) + # Base class for projection selection filters; ensure both CPU and CUDA instantiations + itk_wrap_include(rtkSubSelectImageFilter.h) itk_wrap_image_filter("${WRAP_ITK_REAL}" 1 3) + if(RTK_USE_CUDA) + itk_wrap_include(itkCudaImage.h) + itk_wrap_template("CIF3" "itk::CudaImage") + endif() itk_end_wrap_class() From 063995082e55b34161c405cc9d7c66d0544b9108 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Mon, 17 Nov 2025 11:30:25 +0100 Subject: [PATCH 6/7] ENH: Add rtkfourdrooster python application --- .../rtkfourdrooster/rtkfourdrooster.py | 358 ++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_rtk__.py | 1 + wrapping/itkImageRTK.wrap | 8 +- wrapping/rtkExtras.py | 20 + 5 files changed, 384 insertions(+), 4 deletions(-) create mode 100644 applications/rtkfourdrooster/rtkfourdrooster.py diff --git a/applications/rtkfourdrooster/rtkfourdrooster.py b/applications/rtkfourdrooster/rtkfourdrooster.py new file mode 100644 index 000000000..963bec023 --- /dev/null +++ b/applications/rtkfourdrooster/rtkfourdrooster.py @@ -0,0 +1,358 @@ +#!/usr/bin/env python +import argparse +import sys +import itk +from itk import RTK as rtk + + +def build_parser(): + parser = rtk.RTKArgumentParser( + description=( + "Reconstructs a 3D + time sequence of volumes from a projection stack " + "and a respiratory/cardiac signal, applying TV regularization in space " + "and time, and restricting motion to a region of interest." + ) + ) + + parser.add_argument( + "--geometry", "-g", help="XML geometry file name", type=str, required=True + ) + parser.add_argument("--input", "-i", help="Input volume", type=str) + parser.add_argument( + "--output", "-o", help="Output file name", type=str, required=True + ) + parser.add_argument( + "--niter", + "-n", + help="Number of main loop iterations", + type=int, + default=5, + ) + parser.add_argument( + "--cgiter", + help="Number of conjugate gradient nested iterations", + type=int, + default=4, + ) + parser.add_argument( + "--cudacg", + help="Perform conjugate gradient calculations on GPU", + action="store_true", + ) + parser.add_argument( + "--cudadvfinterpolation", + help="Perform DVF interpolation calculations on GPU", + action="store_true", + ) + parser.add_argument( + "--nodisplaced", + help="Disable the displaced detector filter", + action="store_true", + ) + + # Phase gating + parser.add_argument( + "--signal", + help="File containing the phase of each projection", + type=str, + required=True, + ) + parser.add_argument( + "--shift", + help="Phase shift applied on the DVFs to simulate phase estimation errors", + type=float, + default=0.0, + ) + + # Regularization + parser.add_argument( + "--nopositivity", + help="Do not enforce positivity", + action="store_true", + ) + parser.add_argument( + "--motionmask", + help="Motion mask file: binary image with ones where movement can occur and zeros elsewhere", + type=str, + ) + parser.add_argument( + "--tviter", + help="Total variation (spatial, temporal and nuclear) regularization: number of iterations", + type=int, + default=10, + ) + parser.add_argument( + "--gamma_space", + help="Total variation spatial regularization parameter. The larger, the smoother", + type=float, + ) + parser.add_argument( + "--threshold", + help="Daubechies wavelets spatial regularization: soft threshold", + type=float, + ) + parser.add_argument( + "--order", + help="Daubechies wavelets spatial regularization: order of the wavelets", + type=int, + default=5, + ) + parser.add_argument( + "--levels", + help="Daubechies wavelets spatial regularization: number of decomposition levels", + type=int, + default=3, + ) + parser.add_argument( + "--gamma_time", + help="Total variation temporal regularization parameter. The larger, the smoother", + type=float, + ) + parser.add_argument( + "--lambda_time", + help="Temporal gradient's L0 norm regularization parameter. The larger, the stronger", + type=float, + ) + parser.add_argument( + "--l0iter", + help="Temporal gradient's L0 norm regularization: number of iterations", + type=int, + default=5, + ) + parser.add_argument( + "--gamma_tnv", + help="Total nuclear variation regularization parameter. The larger, the smoother", + type=float, + ) + + # Motion-compensation + parser.add_argument( + "--dvf", + help="Input 4D DVF", + type=str, + ) + parser.add_argument( + "--idvf", + help=( + "Input 4D inverse DVF. Inverse transform computed by conjugate gradient if not provided" + ), + type=str, + ) + parser.add_argument( + "--nn", + help="Nearest neighbor interpolation (default is trilinear)", + action="store_true", + ) + + # RTK common groups (projections input) + rtk.add_rtkinputprojections_group(parser) + rtk.add_rtkprojectors_group(parser) + rtk.add_rtk4Doutputimage_group(parser) + rtk.add_rtkiterations_group(parser) + + return parser + + +def process(args_info: argparse.Namespace): + OutputPixelType = itk.F + DVFVectorType = itk.CovariantVector[OutputPixelType, 3] + VolumeSeriesType = itk.Image[OutputPixelType, 4] + ProjectionStackType = itk.Image[OutputPixelType, 3] + DVFSequenceImageType = itk.Image[ + DVFVectorType, VolumeSeriesType.GetImageDimension() + ] + + VolumeType = ProjectionStackType + + # Check CUDA-related options vs. build capabilities + if (args_info.cudacg or args_info.cudadvfinterpolation) and not hasattr( + itk, "CudaImage" + ): + print( + "Error: CUDA options (--cudacg, --cudadvfinterpolation) were requested but ITK/RTK was " + "not built with CUDA support.", + ) + sys.exit(1) + + # Projections reader + reader = rtk.ProjectionsReader[ProjectionStackType].New() + rtk.SetProjectionsReaderFromArgParse(reader, args_info) + + # 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 existing 4D volume + inputFilter = itk.ImageFileReader[VolumeSeriesType].New() + inputFilter.SetFileName(args_info.input) + else: + # Constant 4D volume + 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() + + # Re-order geometry and projections + # In the new order, projections with identical phases are packed together + signal = rtk.read_signal_file(args_info.signal) + 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() + + # Compute the interpolation weights + signalToInterpolationWeights = rtk.SignalToInterpolationWeights.New() + signalToInterpolationWeights.SetSignal(reorder.GetOutputSignal()) + nb_frames = inputFilter.GetOutput().GetLargestPossibleRegion().GetSize()[3] + signalToInterpolationWeights.SetNumberOfReconstructedFrames(nb_frames) + signalToInterpolationWeights.Update() + + # Create the 4DROOSTER filter, connect the basic inputs, and set the basic parameters + # Also set the forward and back projection filters to be used + if hasattr(itk, "CudaImage"): + CudaVolumeSeriesType = itk.CudaImage[OutputPixelType, 4] + CudaProjectionStackType = itk.CudaImage[OutputPixelType, 3] + rooster = rtk.FourDROOSTERConeBeamReconstructionFilter[ + CudaVolumeSeriesType, CudaProjectionStackType + ].New() + rooster.SetInputVolumeSeries(itk.cuda_image_from_image(inputFilter.GetOutput())) + rooster.SetInputProjectionStack(itk.cuda_image_from_image(reorder.GetOutput())) + else: + rooster = rtk.FourDROOSTERConeBeamReconstructionFilter[ + VolumeSeriesType, ProjectionStackType + ].New() + rooster.SetInputVolumeSeries(inputFilter.GetOutput()) + rooster.SetInputProjectionStack(reorder.GetOutput()) + + rooster.SetCG_iterations(args_info.cgiter) + rooster.SetMainLoop_iterations(args_info.niter) + rooster.SetPhaseShift(args_info.shift) + rooster.SetCudaConjugateGradient(args_info.cudacg) + rooster.SetUseCudaCyclicDeformation(args_info.cudadvfinterpolation) + rooster.SetDisableDisplacedDetectorFilter(args_info.nodisplaced) + rooster.SetGeometry(reorder.GetOutputGeometry()) + rooster.SetWeights(signalToInterpolationWeights.GetOutput()) + rooster.SetSignal(reorder.GetOutputSignal()) + + rtk.SetIterationsReportFromArgParse(args_info, rooster) + + # For each optional regularization step, set whether or not + # it should be performed, and provide the necessary inputs + + # Positivity + rooster.SetPerformPositivity(not args_info.nopositivity) + + rtk.SetForwardProjectionFromArgParse(args_info, rooster) + rtk.SetBackProjectionFromArgParse(args_info, rooster) + # Motion mask + if args_info.motionmask: + reader = itk.ImageFileReader[VolumeType].New() + reader.SetFileName(args_info.motionmask) + if hasattr(itk, "CudaImage"): + motionMask = itk.cuda_image_from_image(reader.GetOutput()) + else: + motionMask = reader.GetOutput() + rooster.SetMotionMask(motionMask) + rooster.SetPerformMotionMask(True) + else: + rooster.SetPerformMotionMask(False) + + # Spatial TV + if args_info.gamma_space is not None: + rooster.SetGammaTVSpace(args_info.gamma_space) + rooster.SetTV_iterations(args_info.tviter) + rooster.SetPerformTVSpatialDenoising(True) + else: + rooster.SetPerformTVSpatialDenoising(False) + + # Spatial wavelets + if args_info.threshold is not None: + rooster.SetSoftThresholdWavelets(args_info.threshold) + rooster.SetOrder(args_info.order) + rooster.SetNumberOfLevels(args_info.levels) + rooster.SetPerformWaveletsSpatialDenoising(True) + else: + rooster.SetPerformWaveletsSpatialDenoising(False) + + # Temporal TV + if args_info.gamma_time is not None: + rooster.SetGammaTVTime(args_info.gamma_time) + rooster.SetTV_iterations(args_info.tviter) + rooster.SetPerformTVTemporalDenoising(True) + else: + rooster.SetPerformTVTemporalDenoising(False) + + # Temporal L0 + if args_info.lambda_time is not None: + rooster.SetLambdaL0Time(args_info.lambda_time) + rooster.SetL0_iterations(args_info.l0iter) + rooster.SetPerformL0TemporalDenoising(True) + else: + rooster.SetPerformL0TemporalDenoising(False) + + # Total nuclear variation + if args_info.gamma_tnv is not None: + rooster.SetGammaTNV(args_info.gamma_tnv) + rooster.SetTV_iterations(args_info.tviter) + rooster.SetPerformTNVDenoising(True) + else: + rooster.SetPerformTNVDenoising(False) + + # Warping + if args_info.dvf: + rooster.SetPerformWarping(True) + + if args_info.nn: + rooster.SetUseNearestNeighborInterpolationInWarping(True) + + reader = itk.ImageFileReader[DVFSequenceImageType].New() + reader.SetFileName(args_info.dvf) + dvf = reader.GetOutput() + if hasattr(itk, "CudaImage"): + dvf = itk.cuda_image_from_image(dvf) + rooster.SetDisplacementField(dvf) + + if args_info.idvf: + rooster.SetComputeInverseWarpingByConjugateGradient(False) + reader.SetFileName(args_info.idvf) + idvf = reader.GetOutput() + if hasattr(itk, "CudaImage"): + idvf = itk.cuda_image_from_image(idvf) + rooster.SetInverseDisplacementField(idvf) + else: + rooster.SetPerformWarping(False) + + rooster.Update() + + # Write + if args_info.verbose: + print(f"Writing output to {args_info.output}...") + writer = itk.ImageFileWriter[itk.Image[OutputPixelType, 4]].New() + writer.SetFileName(args_info.output) + writer.SetInput(rooster.GetOutput()) + writer.Update() + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 7fc9eaa10..f1400f8aa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,6 +58,7 @@ rtkfieldofview = "itk.rtkfieldofview:main" rtkforwardprojections = "itk.rtkforwardprojections:main" rtkfourdconjugategradient = "itk.rtkfourdconjugategradient:main" rtkfourdfdk = "itk.rtkfourdfdk:main" +rtkfourdrooster = "itk.rtkfourdrooster:main" rtki0estimation = "itk.rtki0estimation:main" rtkimagexgeometry = "itk.rtkimagxgeometry:main" rtkmaskcollimation = "itk.rtkmaskcollimation:main" diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index 09ef9400d..a47b2569f 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -40,6 +40,7 @@ "rtkforwardprojections", "rtkfourdconjugategradient", "rtkfourdfdk", + "rtkfourdrooster", "rtki0estimation", "rtkimagxgeometry", "rtkmaskcollimation", diff --git a/wrapping/itkImageRTK.wrap b/wrapping/itkImageRTK.wrap index 5a37e9ccf..533a1fc46 100644 --- a/wrapping/itkImageRTK.wrap +++ b/wrapping/itkImageRTK.wrap @@ -90,8 +90,8 @@ itk_wrap_simple_type("itk::Image< float, 3 >::ConstPointer" "itkImageF3_ConstPoi itk_wrap_simple_type_swig_interface("itk::Image< float, 4 >::ConstPointer" "itkImageF4_ConstPointer") itk_wrap_simple_type("itk::Image< float, 4 >::ConstPointer" "itkImageF4_ConstPointer") -#itk_wrap_simple_type_swig_interface("itk::Image< itk::CovariantVector< float, 3 >, 4 >::ConstPointer" "itkImageCVF34_ConstPointer") -#itk_wrap_simple_type("itk::Image< itk::CovariantVector< float, 3 >, 4 >::ConstPointer" "itkImageCVF34_ConstPointer") +itk_wrap_simple_type_swig_interface("itk::Image< itk::CovariantVector< float, 3 >, 4 >::ConstPointer" "itkImageCVF34_ConstPointer") +itk_wrap_simple_type("itk::Image< itk::CovariantVector< float, 3 >, 4 >::ConstPointer" "itkImageCVF34_ConstPointer") if (ITK_WRAP_double) itk_wrap_simple_type_swig_interface("itk::Image< double, 3 >::ConstPointer" "itkImageD3_ConstPointer") @@ -100,8 +100,8 @@ if (ITK_WRAP_double) itk_wrap_simple_type_swig_interface("itk::Image< double, 4 >::ConstPointer" "itkImageD4_ConstPointer") itk_wrap_simple_type("itk::Image< double, 4 >::ConstPointer" "itkImageD4_ConstPointer") - #itk_wrap_simple_type_swig_interface("itk::Image< itk::CovariantVector< double, 3 >, 4 >::ConstPointer" "itkImageCVD34_ConstPointer") - #itk_wrap_simple_type("itk::Image< itk::CovariantVector< double, 3 >, 4 >::ConstPointer" "itkImageCVD34_ConstPointer") + itk_wrap_simple_type_swig_interface("itk::Image< itk::CovariantVector< double, 3 >, 4 >::ConstPointer" "itkImageCVD34_ConstPointer") + itk_wrap_simple_type("itk::Image< itk::CovariantVector< double, 3 >, 4 >::ConstPointer" "itkImageCVD34_ConstPointer") endif() set(vectorComponents 2 3 4 5) diff --git a/wrapping/rtkExtras.py b/wrapping/rtkExtras.py index e16c07237..0ecccb23d 100644 --- a/wrapping/rtkExtras.py +++ b/wrapping/rtkExtras.py @@ -2,6 +2,7 @@ from itk import RTK as rtk import importlib import shlex +import math # Write a 3D circular projection geometry to a file. @@ -26,6 +27,25 @@ def read_geometry(filename): return reader.GetOutputObject() +# Read a signal file +def read_signal_file(filename): + signal_vector = [] + try: + with open(filename) as f: + for line in f: + line = line.strip() + if line: + value = float(line) + rounded = round(value * 100) / 100 + if rounded == 1.0: + signal_vector.append(0.0) + else: + signal_vector.append(rounded) + except OSError: + raise RuntimeError(f"Could not open signal file {filename}") + return signal_vector + + # Returns the progress percentage class PercentageProgressCommand: def __init__(self, caller): From f059aec31287a6ab6eeb13c471d0549982654925 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 23 Jan 2026 16:02:33 +0100 Subject: [PATCH 7/7] ENH: Add rtkfourdsart python application --- applications/rtkfourdsart/rtkfourdsart.py | 168 ++++++++++++++++++ wrapping/rtkPhasesToInterpolationWeights.wrap | 1 + 2 files changed, 169 insertions(+) create mode 100644 applications/rtkfourdsart/rtkfourdsart.py create mode 100644 wrapping/rtkPhasesToInterpolationWeights.wrap diff --git a/applications/rtkfourdsart/rtkfourdsart.py b/applications/rtkfourdsart/rtkfourdsart.py new file mode 100644 index 000000000..dfdf85572 --- /dev/null +++ b/applications/rtkfourdsart/rtkfourdsart.py @@ -0,0 +1,168 @@ +#!/usr/bin/env python +import argparse +import sys +import itk +from itk import RTK as rtk + + +def build_parser(): + parser = rtk.RTKArgumentParser( + description=( + "Reconstructs a 4D sequence of volumes from a sequence of projections " + "with a 4D version of the Simultaneous Algebraic Reconstruction Technique [Andersen, 1984]." + ) + ) + + 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( + "--lambda", + "-l", + help="Convergence factor", + type=float, + default=0.3, + dest="lambda_val", + ) + parser.add_argument( + "--positivity", + help="Enforces positivity during the reconstruction", + action="store_true", + ) + parser.add_argument("--input", "-i", help="Input volume", type=str) + parser.add_argument( + "--nprojpersubset", + help="Number of projections processed between each update of the reconstructed volume (1 for SART, several for OSSART, all for SIRT)", + type=int, + default=1, + ) + parser.add_argument( + "--nodisplaced", + help="Disable the displaced detector filter", + action="store_true", + ) + + # Phase gating + parser.add_argument( + "--signal", + help="File containing the phase of each projection", + type=str, + required=True, + ) + + # RTK common groups + rtk.add_rtkinputprojections_group(parser) + rtk.add_rtk4Doutputimage_group(parser) + rtk.add_rtkprojectors_group(parser) + rtk.add_rtkiterations_group(parser) + + return parser + + +def process(args_info: argparse.Namespace): + OutputPixelType = itk.F + VolumeSeriesType = itk.Image[OutputPixelType, 4] + ProjectionStackType = itk.Image[OutputPixelType, 3] + + # Check CUDA-related options vs. build capabilities + if hasattr(itk, "CudaImage"): + CudaVolumeSeriesType = itk.CudaImage[OutputPixelType, 4] + CudaProjectionStackType = itk.CudaImage[OutputPixelType, 3] + + # Projections reader + reader = rtk.ProjectionsReader[ProjectionStackType].New() + rtk.SetProjectionsReaderFromArgParse(reader, args_info) + + # 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 existing 4D volume + inputFilter = itk.ImageFileReader[VolumeSeriesType].New() + inputFilter.SetFileName(args_info.input) + else: + # Constant 4D volume + 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() + + # Read the phases file + phaseReader = rtk.PhasesToInterpolationWeights.New() + phaseReader.SetFileName(args_info.signal) + phaseReader.SetNumberOfReconstructedFrames( + inputFilter.GetOutput().GetLargestPossibleRegion().GetSize()[3] + ) + phaseReader.Update() + + # 4D SART reconstruction filter + if hasattr(itk, "CudaImage"): + fourdsart = rtk.CudaFourDSARTConeBeamReconstructionFilter[ + CudaVolumeSeriesType, CudaProjectionStackType + ].New() + fourdsart.SetInputVolumeSeries( + itk.cuda_image_from_image(inputFilter.GetOutput()) + ) + fourdsart.SetInputProjectionStack(itk.cuda_image_from_image(reader.GetOutput())) + else: + fourdsart = rtk.FourDSARTConeBeamReconstructionFilter[ + VolumeSeriesType, ProjectionStackType + ].New() + fourdsart.SetInputVolumeSeries(inputFilter.GetOutput()) + fourdsart.SetInputProjectionStack(reader.GetOutput()) + + # Set the forward and back projection filters + rtk.SetForwardProjectionFromArgParse(args_info, fourdsart) + rtk.SetBackProjectionFromArgParse(args_info, fourdsart) + + fourdsart.SetGeometry(geometry) + fourdsart.SetNumberOfIterations(args_info.niterations) + fourdsart.SetNumberOfProjectionsPerSubset(args_info.nprojpersubset) + fourdsart.SetWeights(phaseReader.GetOutput()) + fourdsart.SetSignal(rtk.read_signal_file(args_info.signal)) + fourdsart.SetLambda(args_info.lambda_val) + fourdsart.SetDisableDisplacedDetectorFilter(args_info.nodisplaced) + + if args_info.positivity: + fourdsart.SetEnforcePositivity(True) + + rtk.SetIterationsReportFromArgParse(args_info, fourdsart) + + fourdsart.Update() + + # Write + if args_info.verbose: + print(f"Writing output to {args_info.output}...") + writer = itk.ImageFileWriter[itk.Image[OutputPixelType, 4]].New() + writer.SetFileName(args_info.output) + writer.SetInput(fourdsart.GetOutput()) + writer.Update() + + +def main(argv=None): + parser = build_parser() + args_info = parser.parse_args(argv) + process(args_info) + + +if __name__ == "__main__": + main() diff --git a/wrapping/rtkPhasesToInterpolationWeights.wrap b/wrapping/rtkPhasesToInterpolationWeights.wrap new file mode 100644 index 000000000..e92a3f4eb --- /dev/null +++ b/wrapping/rtkPhasesToInterpolationWeights.wrap @@ -0,0 +1 @@ +itk_wrap_simple_class("rtk::PhasesToInterpolationWeights" POINTER)