Skip to content

WIP : Spectral python applications #889

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
axel-grc wants to merge 3 commits into
base: main
Choose a base branch
from
Open

WIP : Spectral python applications #889

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
53212f6
ENH: Add rtkspectraldenoiseprojections python application
axel-grc Jan 23, 2026
e5859ca
ENH: Add rtkspectralrooster python application
axel-grc Jan 27, 2026
f331565
ENH: Add rtkspectralforwardmodel python application
axel-grc Jan 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 68 additions & 0 deletions applications/rtkspectraldenoiseprojections/rtkspectraldenoiseprojections.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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()
160 changes: 160 additions & 0 deletions applications/rtkspectralforwardmodel/rtkspectralforwardmodel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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()
Loading
Loading