diff --git a/applications/rtk3Doutputimage_group.py b/applications/rtk3Doutputimage_group.py index 333266857..e7ade3004 100644 --- a/applications/rtk3Doutputimage_group.py +++ b/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 diff --git a/applications/rtkconjugategradient/rtkconjugategradient.py b/applications/rtkconjugategradient/rtkconjugategradient.py index 01bcd003e..354380379 100755 --- a/applications/rtkconjugategradient/rtkconjugategradient.py +++ b/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) diff --git a/applications/rtkfdk/rtkfdk.py b/applications/rtkfdk/rtkfdk.py new file mode 100644 index 000000000..80e68067d --- /dev/null +++ b/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() \ No newline at end of file diff --git a/applications/rtkinputprojections_group.py b/applications/rtkinputprojections_group.py index 29be147dd..3f8d820cf 100644 --- a/applications/rtkinputprojections_group.py +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index db42c8645..fd18d4bba 100644 --- a/pyproject.toml +++ b/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" diff --git a/wrapping/CMakeLists.txt b/wrapping/CMakeLists.txt index 83afdc330..68596dc28 100644 --- a/wrapping/CMakeLists.txt +++ b/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 ) diff --git a/wrapping/rtkExtra.py b/wrapping/rtkExtra.py new file mode 100644 index 000000000..5acd83aca --- /dev/null +++ b/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 \ No newline at end of file