From 53212f6ddd24ab2020cc7bad1755bc74ad27ca9b Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 23 Jan 2026 11:04:25 +0100 Subject: [PATCH 1/3] ENH: Add rtkspectraldenoiseprojections python application --- .../rtkspectraldenoiseprojections.py | 68 +++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_rtk__.py | 1 + 3 files changed, 70 insertions(+) create mode 100644 applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.py diff --git a/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.py b/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.py new file mode 100644 index 000000000..666f7b402 --- /dev/null +++ b/applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python +import argparse +import itk +from itk import RTK as rtk + + +def build_parser(): + parser = rtk.RTKArgumentParser( + description="Replaces aberrant pixels by the median in a small neighborhood around them. Pixels are aberrant if the difference between their value and the median is larger that threshold multiplier * the standard deviation in the neighborhood" + ) + parser.add_argument( + "--input", "-i", help="Input file name", type=str, required=True + ) + parser.add_argument( + "--output", "-o", help="Output file name", type=str, required=True + ) + parser.add_argument( + "--multiplier", + "-m", + help="Threshold multiplier (actual threshold is obtained by multiplying by standard dev. of neighborhood)", + type=float, + default=1.0, + ) + parser.add_argument( + "--radius", + "-r", + help="Radius of neighborhood in each direction (actual radius is 2r+1)", + type=int, + nargs="+", + ) + return parser + + +def process(args_info: argparse.Namespace): + OutputImageType = itk.VectorImage[itk.F, 3] + + # Reader + if args_info.verbose: + print(f"Reading input from {args_info.input}...") + input_image = itk.imread(args_info.input) + + # Remove aberrant pixels + median = rtk.ConditionalMedianImageFilter[OutputImageType].New() + median.SetThresholdMultiplier(args_info.multiplier) + radius = itk.Size[3]() + if args_info.radius is None: + radius.Fill(1) + else: + radius.Fill(args_info.radius[0]) + for i in range(len(args_info.radius)): + radius[i] = args_info.radius[i] + median.SetRadius(radius) + median.SetInput(input_image) + + # Write + if args_info.verbose: + print(f"Writing output to {args_info.output}...") + itk.imwrite(median.GetOutput(), args_info.output) + + +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..893877297 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -68,6 +68,7 @@ rtkprojectshepploganphantom= "itk.rtkprojectshepploganphantom:main" rtkshowgeometry = "itk.rtkshowgeometry:main" rtksart = "itk.rtksart:main" rtksimulatedgeometry = "itk.rtksimulatedgeometry:main" +rtkspectraldenoiseprojections = "itk.rtkspectraldenoiseprojections:main" rtksubselect = "itk.rtksubselect:main" rtktotalvariationdenoising = "itk.rtktotalvariationdenoising:main" rtkvarianobigeometry = "itk.rtkvarianobigeometry:main" diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index fa3befc93..3c8509de8 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -50,6 +50,7 @@ "rtkshowgeometry", "rtksart", "rtksimulatedgeometry", + "rtkspectraldenoiseprojections", "rtksubselect", "rtktotalvariationdenoising", "rtkvarianobigeometry", From e5859cae7af33bc3e043a7c49249dbcfb4854296 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Tue, 27 Jan 2026 13:30:30 +0100 Subject: [PATCH 2/3] ENH: Add rtkspectralrooster python application --- .../rtkspectralrooster/rtkspectralrooster.py | 318 ++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_rtk__.py | 1 + wrapping/itkCSVFileReaderBaseRTK.wrap | 3 + wrapping/itkImageToImageFilterRTK.wrap | 1 + wrapping/rtkImageToVectorImageFilter.wrap | 1 + wrapping/rtkSignalToInterpolationWeights.wrap | 3 + 7 files changed, 328 insertions(+) create mode 100644 applications/rtkspectralrooster/rtkspectralrooster.py create mode 100644 wrapping/itkCSVFileReaderBaseRTK.wrap create mode 100644 wrapping/rtkSignalToInterpolationWeights.wrap diff --git a/applications/rtkspectralrooster/rtkspectralrooster.py b/applications/rtkspectralrooster/rtkspectralrooster.py new file mode 100644 index 000000000..62c6c0803 --- /dev/null +++ b/applications/rtkspectralrooster/rtkspectralrooster.py @@ -0,0 +1,318 @@ +#!/usr/bin/env python +import argparse +import itk +from itk import RTK as rtk + + +def build_parser(): + parser = rtk.RTKArgumentParser( + description=( + "Reconstructs a 3D + material vector volume from a vector projection stack," + " alternating between conjugate gradient optimization and regularization." + ) + ) + + parser.add_argument( + "--geometry", "-g", help="XML geometry file name", required=True + ) + parser.add_argument("--output", "-o", help="Output file name", 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("--input", "-i", help="Input volume (materials) file") + parser.add_argument( + "--projection", "-p", help="Vector projections file", required=True + ) + parser.add_argument( + "--nodisplaced", + help="Disable the displaced detector filter", + action="store_true", + ) + + # Regularization + parser.add_argument( + "--nopositivity", help="Do not enforce positivity", action="store_true" + ) + parser.add_argument("--tviter", help="TV iterations", type=int, default=10) + parser.add_argument( + "--gamma_space", help="Spatial TV regularization parameter", type=float + ) + parser.add_argument("--threshold", help="Wavelets soft threshold", type=float) + parser.add_argument("--order", help="Wavelets order", type=int, default=5) + parser.add_argument("--levels", help="Wavelets levels", type=int, default=3) + parser.add_argument( + "--gamma_time", help="Temporal TV regularization parameter", type=float + ) + parser.add_argument( + "--lambda_time", help="Temporal L0 regularization parameter", type=float + ) + parser.add_argument("--l0iter", help="L0 iterations", type=int, default=5) + parser.add_argument("--gamma_tnv", help="TNV regularization parameter", type=float) + + # Projector choices + rtk.add_rtkprojectors_group(parser) + rtk.add_rtk3Doutputimage_group(parser) + return parser + + +def process(args_info: argparse.Namespace): + PixelValueType = itk.F + Dimension = 3 + + DecomposedProjectionType = itk.VectorImage[PixelValueType, Dimension] + MaterialsVolumeType = itk.VectorImage[PixelValueType, Dimension] + VolumeSeriesType = itk.Image[PixelValueType, Dimension + 1] + ProjectionStackType = itk.Image[PixelValueType, Dimension] + + # Projections reader + if args_info.verbose: + print(f"Reading decomposed projections from {args_info.projection}...") + proj_reader = itk.ImageFileReader[DecomposedProjectionType].New() + proj_reader.SetFileName(args_info.projection) + proj_reader.Update() + decomposedProjection = proj_reader.GetOutput() + + NumberOfMaterials = decomposedProjection.GetVectorLength() + + # Geometry + if args_info.verbose: + print(f"Reading geometry from {args_info.geometry}...") + geometry = rtk.read_geometry(args_info.geometry) + + # Create 4D input. Fill it either with an existing materials volume + # read from a file or a blank image + vecVol2VolSeries = rtk.VectorImageToImageFilter[ + MaterialsVolumeType, VolumeSeriesType + ].New() + + if args_info.input is not None: + if args_info.like is not None: + print("WARNING: --like ignored since --input was given") + if args_info.verbose: + print(f"Reading input volume {args_info.input}...") + reference = itk.imread(args_info.input) + vecVol2VolSeries.SetInput(reference) + vecVol2VolSeries.Update() + input = vecVol2VolSeries.GetOutput() + elif args_info.like is not None: + if args_info.verbose: + print(f"Reading reference volume {args_info.like} to infer geometry...") + reference = itk.imread(args_info.like) + vecVol2VolSeries.SetInput(reference) + vecVol2VolSeries.UpdateOutputInformation() + constantImageSource = rtk.ConstantImageSource[VolumeSeriesType].New() + constantImageSource.SetInformationFromImage(vecVol2VolSeries.GetOutput()) + constantImageSource.Update() + input = constantImageSource.GetOutput() + else: + # Create new empty volume + constantImageSource = rtk.ConstantImageSource[VolumeSeriesType].New() + + imageSize = itk.Size[4]() + imageSize.Fill(args_info.size[0]) + for i in range(min(len(args_info.size), Dimension)): + imageSize[i] = args_info.size[i] + imageSize[Dimension] = NumberOfMaterials + + imageSpacing = itk.Vector[itk.D, 4]() + imageSpacing.Fill(args_info.spacing[0]) + for i in range(min(len(args_info.spacing), Dimension)): + imageSpacing[i] = args_info.spacing[i] + imageSpacing[Dimension] = 1.0 + + imageOrigin = itk.Point[itk.D, 4]() + for i in range(Dimension): + imageOrigin[i] = imageSpacing[i] * (imageSize[i] - 1) * -0.5 + if args_info.origin is not None: + for i in range(min(len(args_info.origin), Dimension)): + imageOrigin[i] = args_info.origin[i] + imageOrigin[Dimension] = 0.0 + + imageDirection = itk.Matrix[itk.D, 4, 4]() + imageDirection.SetIdentity() + if args_info.direction is not None: + for i in range(Dimension): + for j in range(Dimension): + imageDirection[i][j] = args_info.direction[i * Dimension + j] + + constantImageSource.SetOrigin(imageOrigin) + constantImageSource.SetSpacing(imageSpacing) + constantImageSource.SetDirection(imageDirection) + constantImageSource.SetSize(imageSize) + constantImageSource.SetConstant(0.0) + constantImageSource.Update() + input = constantImageSource.GetOutput() + + # Duplicate geometry and transform the N M-vector projections into N*M scalar projections + # Each material will occupy one frame of the 4D reconstruction, therefore all projections + # of one material need to have the same phase. + # Note : the 4D CG filter is optimized when projections with identical phases are packed together + + # Geometry + initialNumberOfProjections = ( + decomposedProjection.GetLargestPossibleRegion().GetSize()[Dimension - 1] + ) + for material in range(1, NumberOfMaterials): + for proj in range(initialNumberOfProjections): + geometry.AddProjectionInRadians( + geometry.GetSourceToIsocenterDistances()[proj], + geometry.GetSourceToDetectorDistances()[proj], + geometry.GetGantryAngles()[proj], + geometry.GetProjectionOffsetsX()[proj], + geometry.GetProjectionOffsetsY()[proj], + geometry.GetOutOfPlaneAngles()[proj], + geometry.GetInPlaneAngles()[proj], + geometry.GetSourceOffsetsX()[proj], + geometry.GetSourceOffsetsY()[proj], + ) + geometry.SetCollimationOfLastProjection( + geometry.GetCollimationUInf()[proj], + geometry.GetCollimationUSup()[proj], + geometry.GetCollimationVInf()[proj], + geometry.GetCollimationVSup()[proj], + ) + + # Signal + fakeSignal = [] + for material in range(NumberOfMaterials): + v = round(material / NumberOfMaterials * 1000.0) / 1000.0 + for proj in range(initialNumberOfProjections): + fakeSignal.append(v) + + # Projections + vproj2proj = rtk.VectorImageToImageFilter[ + DecomposedProjectionType, ProjectionStackType + ].New() + vproj2proj.SetInput(decomposedProjection) + vproj2proj.Update() + + # Release the memory holding the stack of original projections + decomposedProjection.ReleaseData() + + # Compute the interpolation weights + signalToInterpolationWeights = rtk.SignalToInterpolationWeights.New() + signalToInterpolationWeights.SetSignal(fakeSignal) + signalToInterpolationWeights.SetNumberOfReconstructedFrames(NumberOfMaterials) + signalToInterpolationWeights.Update() + + # Set the forward and back projection filters to be used + # Instantiate ROOSTER with CUDA image types if available, otherwise CPU types + if hasattr(itk, "CudaImage"): + cudaVolumeSeriesType = itk.CudaImage[PixelValueType, Dimension + 1] + cudaProjectionStackType = itk.CudaImage[PixelValueType, Dimension] + rooster = rtk.FourDROOSTERConeBeamReconstructionFilter[ + cudaVolumeSeriesType, cudaProjectionStackType + ].New() + rooster.SetInputVolumeSeries(itk.cuda_image_from_image(input)) + rooster.SetInputProjectionStack( + itk.cuda_image_from_image(vproj2proj.GetOutput()) + ) + else: + rooster = rtk.FourDROOSTERConeBeamReconstructionFilter[ + VolumeSeriesType, ProjectionStackType + ].New() + rooster.SetInputVolumeSeries(input) + rooster.SetInputProjectionStack(vproj2proj.GetOutput()) + + # Configure projectors from args + rtk.SetForwardProjectionFromArgParse(args_info, rooster) + rtk.SetBackProjectionFromArgParse(args_info, rooster) + + rooster.SetCG_iterations(args_info.cgiter) + rooster.SetMainLoop_iterations(args_info.niter) + rooster.SetCudaConjugateGradient(args_info.cudacg) + rooster.SetDisableDisplacedDetectorFilter(args_info.nodisplaced) + + rooster.SetGeometry(geometry) + rooster.SetWeights(signalToInterpolationWeights.GetOutput()) + rooster.SetSignal(fakeSignal) + + # 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) + + # No motion mask is used, since there is no motion + 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) + + if args_info.verbose: + print("Running ROOSTER reconstruction...") + rooster.Update() + + # Convert to result to a vector image + volSeries2VecVol = rtk.ImageToVectorImageFilter[ + VolumeSeriesType, MaterialsVolumeType + ].New() + volSeries2VecVol.SetInput(rooster.GetOutput()) + volSeries2VecVol.Update() + + # Write + if args_info.verbose: + print(f"Writing output to {args_info.output}...") + writer = itk.ImageFileWriter[MaterialsVolumeType].New() + writer.SetFileName(args_info.output) + writer.SetInput(volSeries2VecVol.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 893877297..6a4582414 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,6 +69,7 @@ rtkshowgeometry = "itk.rtkshowgeometry:main" rtksart = "itk.rtksart:main" rtksimulatedgeometry = "itk.rtksimulatedgeometry:main" rtkspectraldenoiseprojections = "itk.rtkspectraldenoiseprojections:main" +rtkspectralrooster = "itk.rtkspectralrooster:main" rtksubselect = "itk.rtksubselect:main" rtktotalvariationdenoising = "itk.rtktotalvariationdenoising:main" rtkvarianobigeometry = "itk.rtkvarianobigeometry:main" diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index 3c8509de8..99a83745b 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -51,6 +51,7 @@ "rtksart", "rtksimulatedgeometry", "rtkspectraldenoiseprojections", + "rtkspectralrooster", "rtksubselect", "rtktotalvariationdenoising", "rtkvarianobigeometry", diff --git a/wrapping/itkCSVFileReaderBaseRTK.wrap b/wrapping/itkCSVFileReaderBaseRTK.wrap new file mode 100644 index 000000000..db604e98e --- /dev/null +++ b/wrapping/itkCSVFileReaderBaseRTK.wrap @@ -0,0 +1,3 @@ +itk_wrap_include(itkCSVFileReaderBase.h) +itk_wrap_simple_class("itk::CSVFileReaderBase" POINTER) +itk_end_wrap_class() diff --git a/wrapping/itkImageToImageFilterRTK.wrap b/wrapping/itkImageToImageFilterRTK.wrap index fd27f2d71..c87bc0689 100644 --- a/wrapping/itkImageToImageFilterRTK.wrap +++ b/wrapping/itkImageToImageFilterRTK.wrap @@ -139,6 +139,7 @@ itk_wrap_class("itk::ImageToImageFilter" POINTER) itk_wrap_template("I${ITKM_${t}}2I${ITKM_${t}}1" "itk::Image<${ITKT_${t}}, 2>, itk::Image<${ITKT_${t}}, 1>") endif() itk_wrap_template("I${ITKM_${t}}3VI${ITKM_${t}}2" "itk::Image<${ITKT_${t}}, 3>, itk::VectorImage<${ITKT_${t}}, 2>") + itk_wrap_template("I${ITKM_${t}}4VI${ITKM_${t}}3" "itk::Image<${ITKT_${t}}, 4>, itk::VectorImage<${ITKT_${t}}, 3>") itk_wrap_template("VI${ITKM_${t}}2I${ITKM_${t}}3" "itk::VectorImage<${ITKT_${t}}, 2>, itk::Image<${ITKT_${t}}, 3>") itk_wrap_template("VI${ITKM_${t}}3I${ITKM_${t}}4" "itk::VectorImage<${ITKT_${t}}, 3>, itk::Image<${ITKT_${t}}, 4>") endforeach() diff --git a/wrapping/rtkImageToVectorImageFilter.wrap b/wrapping/rtkImageToVectorImageFilter.wrap index 5b510c544..94a8011de 100644 --- a/wrapping/rtkImageToVectorImageFilter.wrap +++ b/wrapping/rtkImageToVectorImageFilter.wrap @@ -2,5 +2,6 @@ itk_wrap_class("rtk::ImageToVectorImageFilter" POINTER) foreach(t ${WRAP_ITK_REAL}) itk_wrap_template("I${ITKM_${t}}2VI${ITKM_${t}}2" "itk::Image<${ITKT_${t}}, 2>, itk::VectorImage<${ITKT_${t}}, 2>") itk_wrap_template("I${ITKM_${t}}3VI${ITKM_${t}}2" "itk::Image<${ITKT_${t}}, 3>, itk::VectorImage<${ITKT_${t}}, 2>") + itk_wrap_template("I${ITKM_${t}}4VI${ITKM_${t}}3" "itk::Image<${ITKT_${t}}, 4>, itk::VectorImage<${ITKT_${t}}, 3>") endforeach() itk_end_wrap_class() diff --git a/wrapping/rtkSignalToInterpolationWeights.wrap b/wrapping/rtkSignalToInterpolationWeights.wrap new file mode 100644 index 000000000..980a7a9ea --- /dev/null +++ b/wrapping/rtkSignalToInterpolationWeights.wrap @@ -0,0 +1,3 @@ +itk_wrap_include(rtkSignalToInterpolationWeights.h) +itk_wrap_simple_class("rtk::SignalToInterpolationWeights" POINTER) +itk_end_wrap_class() From f33156547e8804810bb824a897b590cd22e398e6 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Fri, 30 Jan 2026 10:33:35 +0100 Subject: [PATCH 3/3] ENH: Add rtkspectralforwardmodel python application --- .../rtkspectralforwardmodel.py | 160 ++++++++++++++++++ pyproject.toml | 1 + wrapping/__init_rtk__.py | 1 + 3 files changed, 162 insertions(+) create mode 100644 applications/rtkspectralforwardmodel/rtkspectralforwardmodel.py diff --git a/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.py b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.py new file mode 100644 index 000000000..0e0f625e8 --- /dev/null +++ b/applications/rtkspectralforwardmodel/rtkspectralforwardmodel.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python +import argparse +import itk +from itk import RTK as rtk + + +def build_parser(): + parser = rtk.RTKArgumentParser( + description=( + "Computes expected photon counts from incident spectrum, material " + "attenuations, detector response and material-decomposed projections" + ) + ) + parser.add_argument( + "--output", + "-o", + help="Output file name (photon counts)", + type=str, + required=True, + ) + parser.add_argument( + "--input", "-i", help="Material-decomposed projections", type=str, required=True + ) + parser.add_argument( + "--detector", "-d", help="Detector response file", type=str, required=True + ) + parser.add_argument( + "--incident", help="Incident spectrum file", type=str, required=True + ) + parser.add_argument( + "--attenuations", + "-a", + help="Material attenuations file", + type=str, + required=True, + ) + parser.add_argument( + "--thresholds", + "-t", + help="Lower threshold of bins, expressed in pulse height", + type=int, + nargs="+", + required=True, + ) + parser.add_argument( + "--cramer_rao", help="Output Cramer-Rao lower bound file", type=str + ) + parser.add_argument( + "--variances", help="Output variances of photon counts", type=str + ) + return parser + + +def process(args_info: argparse.Namespace): + PixelValueType = itk.F + Dimension = 3 + + DecomposedProjectionType = itk.VectorImage[PixelValueType, Dimension] + MeasuredProjectionsType = itk.VectorImage[PixelValueType, Dimension] + IncidentSpectrumImageType = itk.Image[PixelValueType, Dimension] + DetectorResponseImageType = itk.Image[PixelValueType, Dimension - 1] + MaterialAttenuationsImageType = itk.Image[PixelValueType, Dimension - 1] + + # Read all inputs + if args_info.verbose: + print(f"Reading decomposed projections from {args_info.input}...") + decomposedProjection = itk.imread(args_info.input) + + if args_info.verbose: + print(f"Reading incident spectrum from {args_info.incident}...") + incidentSpectrum = itk.imread(args_info.incident) + + if args_info.verbose: + print(f"Reading detector response from {args_info.detector}...") + detectorResponse = itk.imread(args_info.detector) + + if args_info.verbose: + print(f"Reading material attenuations from {args_info.attenuations}...") + materialAttenuations = itk.imread(args_info.attenuations) + + # Get parameters from the images + NumberOfMaterials = materialAttenuations.GetLargestPossibleRegion().GetSize()[0] + NumberOfSpectralBins = len(args_info.thresholds) + MaximumEnergy = incidentSpectrum.GetLargestPossibleRegion().GetSize()[0] + + # Generate a set of zero-filled photon count projections + measuredProjections = MeasuredProjectionsType.New() + measuredProjections.CopyInformation(decomposedProjection) + measuredProjections.SetVectorLength(NumberOfSpectralBins) + measuredProjections.Allocate() + + # Read the thresholds on command line + thresholds = [int(t) for t in args_info.thresholds] + MaximumPulseHeight = detectorResponse.GetLargestPossibleRegion().GetSize()[1] + # Add the maximum pulse height at the end + thresholds.append(MaximumPulseHeight) + + # Check that the inputs have the expected size + idx = itk.Index[3]() + idx.Fill(0) + if decomposedProjection.GetPixel(idx).Size() != NumberOfMaterials: + raise RuntimeError( + f"Decomposed projections vector size {decomposedProjection.GetPixel(idx).Size()} != {NumberOfMaterials}" + ) + + if measuredProjections.GetPixel(idx).Size() != NumberOfSpectralBins: + raise RuntimeError( + f"Spectral projections vector size {measuredProjections.GetPixel(idx).Size()} != {NumberOfSpectralBins}" + ) + + if detectorResponse.GetLargestPossibleRegion().GetSize()[0] != MaximumEnergy: + raise RuntimeError( + f"Detector response energies {detectorResponse.GetLargestPossibleRegion().GetSize()[0]} != {MaximumEnergy}" + ) + + # Create and set the filter + forward = rtk.SpectralForwardModelImageFilter[ + DecomposedProjectionType, MeasuredProjectionsType, IncidentSpectrumImageType + ].New() + forward.SetInputDecomposedProjections(decomposedProjection) + forward.SetInputMeasuredProjections(measuredProjections) + forward.SetInputIncidentSpectrum(incidentSpectrum) + forward.SetDetectorResponse(detectorResponse) + forward.SetMaterialAttenuations(materialAttenuations) + forward.SetThresholds(thresholds) + if args_info.cramer_rao: + forward.SetComputeCramerRaoLowerBound(True) + if args_info.variances: + forward.SetComputeVariances(True) + + if args_info.verbose: + print("Running spectral forward model...") + forward.Update() + + # Write output + if args_info.verbose: + print(f"Writing output photon counts to {args_info.output}...") + itk.imwrite(forward.GetOutput(), args_info.output) + + # If requested, write the Cramer-Rao lower bound + if args_info.cramer_rao: + if args_info.verbose: + print(f"Writing Cramer-Rao lower bound to {args_info.cramer_rao}...") + itk.imwrite(forward.GetOutputCramerRaoLowerBound(), args_info.cramer_rao) + + # If requested, write the variance + if args_info.variances: + if args_info.verbose: + print(f"Writing variances to {args_info.variances}...") + itk.imwrite(forward.GetOutputVariances(), args_info.variances) + + +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 6a4582414..8ac3ef8dd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -68,6 +68,7 @@ rtkprojectshepploganphantom= "itk.rtkprojectshepploganphantom:main" rtkshowgeometry = "itk.rtkshowgeometry:main" rtksart = "itk.rtksart:main" rtksimulatedgeometry = "itk.rtksimulatedgeometry:main" +rtkspectralforwardmodel = "itk.rtkspectralforwardmodel:main" rtkspectraldenoiseprojections = "itk.rtkspectraldenoiseprojections:main" rtkspectralrooster = "itk.rtkspectralrooster:main" rtksubselect = "itk.rtksubselect:main" diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index 99a83745b..3c54be40d 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -50,6 +50,7 @@ "rtkshowgeometry", "rtksart", "rtksimulatedgeometry", + "rtkspectralforwardmodel", "rtkspectraldenoiseprojections", "rtkspectralrooster", "rtksubselect",