Skip to content

Add gaussian and poisson noise options to projection applications #860

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 1 commit into
base: main
Choose a base branch
from
Open

Add gaussian and poisson noise options to projection applications #860

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 1 addition & 1 deletion applications/rtkforwardprojections/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
wrap_ggo(rtkforwardprojections_GGO_C rtkforwardprojections.ggo ../rtk3Doutputimage_section.ggo)
wrap_ggo(rtkforwardprojections_GGO_C rtkforwardprojections.ggo ../rtk3Doutputimage_section.ggo ../rtknoise_section.ggo)
add_executable(
rtkforwardprojections
rtkforwardprojections.cxx
Expand Down
6 changes: 5 additions & 1 deletion applications/rtkforwardprojections/rtkforwardprojections.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,12 +191,16 @@ main(int argc, char * argv[])
TRY_AND_EXIT_ON_ITK_EXCEPTION(forwardProjection->Update())
}

// Add noise
OutputImageType::Pointer output =
rtk::SetNoiseFromGgo<OutputImageType, args_info_rtkforwardprojections>(forwardProjection->GetOutput(), args_info);

// Write
if (args_info.verbose_flag)
std::cout << "Writing... " << std::endl;
auto writer = itk::ImageFileWriter<OutputImageType>::New();
writer->SetFileName(args_info.output_arg);
writer->SetInput(forwardProjection->GetOutput());
writer->SetInput(output);
if (args_info.lowmem_flag)
{
writer->SetNumberOfStreamDivisions(sizeOutput[2]);
Expand Down
6 changes: 5 additions & 1 deletion applications/rtkforwardprojections/rtkforwardprojections.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,7 @@ def build_parser():
warp_forwardprojection_group.add_argument("--dvf", help="Input 4D DVF", type=str)

rtk.add_rtk3Doutputimage_group(parser)
rtk.add_rtknoise_group(parser)

return parser

Expand Down Expand Up @@ -227,13 +228,16 @@ def process(args_info):
if not args_info.lowmem:
forwardProjection.Update()

# Add noise
output = rtk.SetNoiseFromArgParse(forwardProjection.GetOutput(), args_info)

# Write
if args_info.verbose:
print("Writing...")

writer = itk.ImageFileWriter[OutputImageType].New()
writer.SetFileName(args_info.output)
writer.SetInput(forwardProjection.GetOutput())
writer.SetInput(output)
writer.Update()

if args_info.verbose:
Expand Down
88 changes: 88 additions & 0 deletions applications/rtknoise_group.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@
import itk

__all__ = [
"add_rtknoise_group",
"SetNoiseFromArgParse",
]


# Mimics rtknoise_section.ggo
def add_rtknoise_group(parser):
rtknoise_group = parser.add_argument_group("Noise")
rtknoise_group.add_argument(
"--gaussian",
help="Gaussian noise parameters: <mean> Noise level and <std> Noise standard deviation",
type=float,
nargs="+",
)
rtknoise_group.add_argument(
"--poisson",
help=(
"Poisson noise parameters: <I0> Number of impinging photons per pixel "
"and <muref> reference linear attenuation coefficient"
),
type=float,
nargs="+",
)


def SetNoiseFromArgParse(projections, args_info):
OutputImageType = type(projections)
output = projections

if args_info.gaussian is not None:
noisy = itk.AdditiveGaussianNoiseImageFilter[
OutputImageType, OutputImageType
].New()
mean, std = args_info.gaussian
noisy.SetInput(output)
noisy.SetMean(mean)
noisy.SetStandardDeviation(std)
noisy.UpdateOutputInformation()
output = noisy.GetOutput()

if args_info.poisson is not None:
I0, muref = args_info.poisson

multiply = itk.MultiplyImageFilter[
OutputImageType, OutputImageType, OutputImageType
].New()
multiply.SetInput(output)
multiply.SetConstant(-muref)

expf = itk.ExpImageFilter[OutputImageType, OutputImageType].New()
expf.SetInput(multiply.GetOutput())

multiply2 = itk.MultiplyImageFilter[
OutputImageType, OutputImageType, OutputImageType
].New()
multiply2.SetInput(expf.GetOutput())
multiply2.SetConstant(I0)

poisson = itk.ShotNoiseImageFilter[OutputImageType].New()
poisson.SetInput(multiply2.GetOutput())

threshold = itk.ThresholdImageFilter[OutputImageType].New()
threshold.SetInput(poisson.GetOutput())
threshold.SetLower(1.0)
threshold.SetOutsideValue(1.0)

multiply3 = itk.MultiplyImageFilter[
OutputImageType, OutputImageType, OutputImageType
].New()
multiply3.SetInput(threshold.GetOutput())
multiply3.SetConstant(1.0 / I0)

logf = itk.LogImageFilter[OutputImageType, OutputImageType].New()
logf.SetInput(multiply3.GetOutput())

multiply4 = itk.MultiplyImageFilter[
OutputImageType, OutputImageType, OutputImageType
].New()
multiply4.SetInput(logf.GetOutput())
multiply4.SetConstant(-1.0 / muref)

multiply4.UpdateOutputInformation()
output = multiply4.GetOutput()

return output
5 changes: 5 additions & 0 deletions applications/rtknoise_section.ggo
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
section "Gaussian noise"
option "gaussian" - "Gaussian noise parameters: <mean> Noise level and <std> Noise standard deviation" double multiple no

section "Poisson noise"
option "poisson" - "Poisson noise parameters: <I0> Number of impinging photons per pixel and <muref> reference linear attenuation coefficient" double multiple no
Comment on lines +1 to +5
Copy link
Collaborator

@lesaintjerome lesaintjerome Nov 28, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems to make the mu_ref parameter mandatory. Why not. Though we could mention something like "e.g. set 0.01879, which is the attenuation coefficient of water at 75 keV".

Copy link
Collaborator

@lesaintjerome lesaintjerome Nov 28, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same could be said about the mean parameter for gaussian noise. Default=0. would make sense, to me.

Copy link
Collaborator Author

@axel-grc axel-grc Dec 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The problem is that we can't use different default values for lists with gengetopt, so the same default value will be used for mean and std for example.

We can maybe add theses default values to the implementation part ?

Copy link
Collaborator

@lesaintjerome lesaintjerome Dec 15, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or have slightly more complex Gaussian and Poisson ggo sections, with 2 parameters each?

2 changes: 1 addition & 1 deletion applications/rtkprojectgeometricphantom/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
wrap_ggo(rtkprojectgeometricphantom_GGO_C rtkprojectgeometricphantom.ggo ../rtk3Doutputimage_section.ggo)
wrap_ggo(rtkprojectgeometricphantom_GGO_C rtkprojectgeometricphantom.ggo ../rtk3Doutputimage_section.ggo ../rtknoise_section.ggo)
add_executable(
rtkprojectgeometricphantom
rtkprojectgeometricphantom.cxx
Expand Down
6 changes: 5 additions & 1 deletion applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,10 +102,14 @@ main(int argc, char * argv[])

TRY_AND_EXIT_ON_ITK_EXCEPTION(ppc->Update())

// Add noise
OutputImageType::Pointer output = ppc->GetOutput();
output = rtk::SetNoiseFromGgo<OutputImageType, args_info_rtkprojectgeometricphantom>(output, args_info);

// Write
if (args_info.verbose_flag)
std::cout << "Projecting and writing... " << std::flush;
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(ppc->GetOutput(), args_info.output_arg))
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(output, args_info.output_arg))

return EXIT_SUCCESS;
}
7 changes: 5 additions & 2 deletions applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ def build_parser():
)

rtk.add_rtk3Doutputimage_group(parser)
rtk.add_rtknoise_group(parser)

# Parse the command line arguments
return parser
Expand Down Expand Up @@ -104,11 +105,13 @@ def process(args_info: argparse.Namespace):
ppc.SetConfigFile(args_info.phantomfile)
ppc.Update()

# Write
# Add noise
output = rtk.SetNoiseFromArgParse(ppc.GetOutput(), args_info)

if args_info.verbose:
print("Projecting and writing... ")

itk.imwrite(ppc.GetOutput(), args_info.output)
itk.imwrite(output, args_info.output)


def main(argv=None):
Expand Down
2 changes: 1 addition & 1 deletion applications/rtkprojections/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
wrap_ggo(rtkprojections_GGO_C rtkprojections.ggo ../rtkinputprojections_section.ggo)
wrap_ggo(rtkprojections_GGO_C rtkprojections.ggo ../rtkinputprojections_section.ggo ../rtknoise_section.ggo)
add_executable(
rtkprojections
rtkprojections.cxx
Expand Down
5 changes: 4 additions & 1 deletion applications/rtkprojections/rtkprojections.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,10 +35,13 @@ main(int argc, char * argv[])
auto reader = ReaderType::New();
rtk::SetProjectionsReaderFromGgo<ReaderType, args_info_rtkprojections>(reader, args_info);

OutputImageType::Pointer output =
rtk::SetNoiseFromGgo<OutputImageType, args_info_rtkprojections>(reader->GetOutput(), args_info);

// Write
auto writer = itk::ImageFileWriter<OutputImageType>::New();
writer->SetFileName(args_info.output_arg);
writer->SetInput(reader->GetOutput());
writer->SetInput(output);
TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->UpdateOutputInformation())
writer->SetNumberOfStreamDivisions(1 + reader->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels() /
(1024 * 1024 * 4));
Expand Down
6 changes: 5 additions & 1 deletion applications/rtkprojections/rtkprojections.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ def build_parser():
"--output", "-o", help="Output file name", type=str, required=True
)
rtk.add_rtkinputprojections_group(parser)
rtk.add_rtknoise_group(parser)

# Parse the command line arguments
return parser
Expand All @@ -29,10 +30,13 @@ def process(args_info: argparse.Namespace):
if args_info.verbose:
print(f"Reading projections...")

# Add noise
output = rtk.SetNoiseFromArgParse(reader.GetOutput(), args_info)

# Write
writer = itk.ImageFileWriter[OutputImageType].New()
writer.SetFileName(args_info.output)
writer.SetInput(reader.GetOutput())
writer.SetInput(output)
if args_info.verbose:
print(f"Writing output to: {args_info.output}")
writer.Update()
Expand Down
2 changes: 1 addition & 1 deletion applications/rtkprojectshepploganphantom/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
wrap_ggo(rtkprojectshepploganphantom_GGO_C rtkprojectshepploganphantom.ggo ../rtk3Doutputimage_section.ggo)
wrap_ggo(rtkprojectshepploganphantom_GGO_C rtkprojectshepploganphantom.ggo ../rtk3Doutputimage_section.ggo ../rtknoise_section.ggo)
add_executable(
rtkprojectshepploganphantom
rtkprojectshepploganphantom.cxx
Expand Down
12 changes: 2 additions & 10 deletions applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,16 +75,8 @@ main(int argc, char * argv[])
TRY_AND_EXIT_ON_ITK_EXCEPTION(slp->Update())

// Add noise
OutputImageType::Pointer output = slp->GetOutput();
if (args_info.noise_given)
{
auto noisy = rtk::AdditiveGaussianNoiseImageFilter<OutputImageType>::New();
noisy->SetInput(slp->GetOutput());
noisy->SetMean(0.0);
noisy->SetStandardDeviation(args_info.noise_arg);
TRY_AND_EXIT_ON_ITK_EXCEPTION(noisy->Update())
output = noisy->GetOutput();
}
OutputImageType::Pointer output =
rtk::SetNoiseFromGgo<OutputImageType, args_info_rtkprojectshepploganphantom>(slp->GetOutput(), args_info);

// Write
if (args_info.verbose_flag)
Expand Down
1 change: 0 additions & 1 deletion applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.ggo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,5 +3,4 @@ purpose "Computes projections through a 3D Shepp & Logan phantom, according to a
option "geometry" g "XML geometry file name" string yes
option "output" o "Output projections file name" string yes
option "phantomscale" - "Scaling factor for the phantom dimensions" double multiple no default="128"
option "noise" - "Gaussian noise parameter (SD)" double no
option "offset" - "3D spatial offset of the phantom center" double multiple no
13 changes: 2 additions & 11 deletions applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ def build_parser():
)

rtk.add_rtk3Doutputimage_group(parser)
rtk.add_rtknoise_group(parser)

# Parse the command line arguments
return parser
Expand Down Expand Up @@ -74,18 +75,8 @@ def process(args_info: argparse.Namespace):
slp.SetOriginOffset(offset)
slp.Update()

output = slp.GetOutput()

# Add noise
if args_info.noise:
noisy = rtk.AdditiveGaussianNoiseImageFilter[
OutputImageType, OutputImageType
].New()
noisy.SetInput(slp.GetOutput())
noisy.SetMean(0.0)
noisy.SetStandardDeviation(args_info.noise)
noisy.Update()
output = noisy.GetOutput()
output = rtk.SetNoiseFromArgParse(slp.GetOutput(), args_info)

# Write
if args_info.verbose:
Expand Down
4 changes: 4 additions & 0 deletions documentation/docs/rtk_3_migration_guide.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,3 +43,7 @@ Most users intuitively expect this for the backprojector matched to `rtk::CudaFo
## CUDA forward / ray-cast step size default

The CUDA forward projector and the CUDA ray-cast backprojector now use the minimum voxel spacing of the input volume as the ray step when StepSize is not set. The previous default value was 1.

## Centralized noise handling

Noise options are now provided centrally via `rtknoise_section.ggo` and the Python helper `rtknoise_group.py` (for example: `--poisson <I0>,<muref>`, `--gaussian <mean>,<std>`). The legacy `--noise` option in `rtkprojectshepploganphantom` has been removed.
65 changes: 65 additions & 0 deletions include/rtkGgoFunctions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,12 @@
#include "rtkProjectionsReader.h"
#include <itkRegularExpressionSeriesFileNames.h>
#include <itksys/RegularExpression.hxx>
#include <itkMultiplyImageFilter.h>
#include <itkExpImageFilter.h>
#include <itkShotNoiseImageFilter.h>
#include <itkThresholdImageFilter.h>
#include <itkLogImageFilter.h>
#include "rtkAdditiveGaussianNoiseImageFilter.h"

namespace rtk
{
Expand Down Expand Up @@ -282,6 +288,65 @@ SetProjectionsReaderFromGgo(TProjectionsReaderType * reader, const TArgsInfo & a
TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->UpdateOutputInformation());
}


template <class TImageType, class TArgsInfo>
typename TImageType::Pointer
SetNoiseFromGgo(typename TImageType::Pointer projections, const TArgsInfo & args_info)
{
using ImageType = TImageType;

// Gaussian noise
if (args_info.gaussian_given)
{
auto noisy = rtk::AdditiveGaussianNoiseImageFilter<ImageType>::New();
noisy->SetInput(projections);
noisy->SetMean(args_info.gaussian_arg[0]);
noisy->SetStandardDeviation(args_info.gaussian_arg[1]);

TRY_AND_EXIT_ON_ITK_EXCEPTION(noisy->Update());
projections = noisy->GetOutput();
}

// Poisson noise
if (args_info.poisson_given)
{
auto multiply = itk::MultiplyImageFilter<ImageType>::New();
multiply->SetInput(projections);
multiply->SetConstant(-args_info.poisson_arg[1]); // -muref

auto expf = itk::ExpImageFilter<ImageType, ImageType>::New();
expf->SetInput(multiply->GetOutput());

auto multiply2 = itk::MultiplyImageFilter<ImageType>::New();
multiply2->SetInput(expf->GetOutput());
multiply2->SetConstant(args_info.poisson_arg[0]); // I0

auto poisson = itk::ShotNoiseImageFilter<ImageType>::New();
poisson->SetInput(multiply2->GetOutput());

auto threshold = itk::ThresholdImageFilter<ImageType>::New();
threshold->SetInput(poisson->GetOutput());
threshold->SetLower(1.);
threshold->SetOutsideValue(1.);

auto multiply3 = itk::MultiplyImageFilter<ImageType>::New();
multiply3->SetInput(threshold->GetOutput());
multiply3->SetConstant(1. / args_info.poisson_arg[0]);

auto logf = itk::LogImageFilter<ImageType, ImageType>::New();
logf->SetInput(multiply3->GetOutput());

auto multiply4 = itk::MultiplyImageFilter<ImageType>::New();
multiply4->SetInput(logf->GetOutput());
multiply4->SetConstant(-1. / args_info.poisson_arg[1]);

TRY_AND_EXIT_ON_ITK_EXCEPTION(multiply4->Update());
projections = multiply4->GetOutput();
}

return projections;
}

/** \brief Set the correct RTK backprojection type from gengetopt
* Given a gengetopt bp_arg option from the rtkprojectors_section.ggo, will
* set the corresponding enum value from rtk::IterativeConeBeamReconstructionFilter
Expand Down
1 change: 1 addition & 0 deletions wrapping/__init_rtk__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
"itk.rtk3Doutputimage_group",
"itk.rtkprojectors_group",
"itk.rtkiterations_group",
"itk.rtknoise_group",
"itk.rtkExtras",
]
for mod_name in rtk_submodules:
Expand Down
Loading