diff --git a/applications/rtkforwardprojections/CMakeLists.txt b/applications/rtkforwardprojections/CMakeLists.txt index 77467f483..a4613b409 100644 --- a/applications/rtkforwardprojections/CMakeLists.txt +++ b/applications/rtkforwardprojections/CMakeLists.txt @@ -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 diff --git a/applications/rtkforwardprojections/rtkforwardprojections.cxx b/applications/rtkforwardprojections/rtkforwardprojections.cxx index 486e8d1a7..effc21f16 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.cxx +++ b/applications/rtkforwardprojections/rtkforwardprojections.cxx @@ -191,12 +191,16 @@ main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(forwardProjection->Update()) } + // Add noise + OutputImageType::Pointer output = + rtk::SetNoiseFromGgo(forwardProjection->GetOutput(), args_info); + // Write if (args_info.verbose_flag) std::cout << "Writing... " << std::endl; auto writer = itk::ImageFileWriter::New(); writer->SetFileName(args_info.output_arg); - writer->SetInput(forwardProjection->GetOutput()); + writer->SetInput(output); if (args_info.lowmem_flag) { writer->SetNumberOfStreamDivisions(sizeOutput[2]); diff --git a/applications/rtkforwardprojections/rtkforwardprojections.py b/applications/rtkforwardprojections/rtkforwardprojections.py index 266510712..fa26fc497 100644 --- a/applications/rtkforwardprojections/rtkforwardprojections.py +++ b/applications/rtkforwardprojections/rtkforwardprojections.py @@ -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 @@ -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: diff --git a/applications/rtknoise_group.py b/applications/rtknoise_group.py new file mode 100644 index 000000000..90193a7a9 --- /dev/null +++ b/applications/rtknoise_group.py @@ -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: Noise level and Noise standard deviation", + type=float, + nargs="+", + ) + rtknoise_group.add_argument( + "--poisson", + help=( + "Poisson noise parameters: Number of impinging photons per pixel " + "and 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 diff --git a/applications/rtknoise_section.ggo b/applications/rtknoise_section.ggo new file mode 100644 index 000000000..2f112f00c --- /dev/null +++ b/applications/rtknoise_section.ggo @@ -0,0 +1,5 @@ +section "Gaussian noise" +option "gaussian" - "Gaussian noise parameters: Noise level and Noise standard deviation" double multiple no + +section "Poisson noise" +option "poisson" - "Poisson noise parameters: Number of impinging photons per pixel and reference linear attenuation coefficient" double multiple no diff --git a/applications/rtkprojectgeometricphantom/CMakeLists.txt b/applications/rtkprojectgeometricphantom/CMakeLists.txt index aabfe8157..f104e3351 100644 --- a/applications/rtkprojectgeometricphantom/CMakeLists.txt +++ b/applications/rtkprojectgeometricphantom/CMakeLists.txt @@ -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 diff --git a/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.cxx b/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.cxx index 868129f33..97c342a3b 100644 --- a/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.cxx +++ b/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.cxx @@ -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(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; } diff --git a/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.py b/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.py index f9a043d72..508fab501 100644 --- a/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.py +++ b/applications/rtkprojectgeometricphantom/rtkprojectgeometricphantom.py @@ -45,6 +45,7 @@ def build_parser(): ) rtk.add_rtk3Doutputimage_group(parser) + rtk.add_rtknoise_group(parser) # Parse the command line arguments return parser @@ -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): diff --git a/applications/rtkprojections/CMakeLists.txt b/applications/rtkprojections/CMakeLists.txt index 6d06e0af4..15649b5a3 100644 --- a/applications/rtkprojections/CMakeLists.txt +++ b/applications/rtkprojections/CMakeLists.txt @@ -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 diff --git a/applications/rtkprojections/rtkprojections.cxx b/applications/rtkprojections/rtkprojections.cxx index 5175f59a9..acf819bd4 100644 --- a/applications/rtkprojections/rtkprojections.cxx +++ b/applications/rtkprojections/rtkprojections.cxx @@ -35,10 +35,13 @@ main(int argc, char * argv[]) auto reader = ReaderType::New(); rtk::SetProjectionsReaderFromGgo(reader, args_info); + OutputImageType::Pointer output = + rtk::SetNoiseFromGgo(reader->GetOutput(), args_info); + // Write auto writer = itk::ImageFileWriter::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)); diff --git a/applications/rtkprojections/rtkprojections.py b/applications/rtkprojections/rtkprojections.py index 07c7ae3b3..b3b01f5a1 100644 --- a/applications/rtkprojections/rtkprojections.py +++ b/applications/rtkprojections/rtkprojections.py @@ -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 @@ -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() diff --git a/applications/rtkprojectshepploganphantom/CMakeLists.txt b/applications/rtkprojectshepploganphantom/CMakeLists.txt index fd961792c..54e8c502a 100644 --- a/applications/rtkprojectshepploganphantom/CMakeLists.txt +++ b/applications/rtkprojectshepploganphantom/CMakeLists.txt @@ -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 diff --git a/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.cxx b/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.cxx index 90e922d19..d35df71a5 100644 --- a/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.cxx +++ b/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.cxx @@ -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::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(slp->GetOutput(), args_info); // Write if (args_info.verbose_flag) diff --git a/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.ggo b/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.ggo index 06721156d..0581c258a 100644 --- a/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.ggo +++ b/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.ggo @@ -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 diff --git a/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.py b/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.py index 1fc336507..d85b86128 100644 --- a/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.py +++ b/applications/rtkprojectshepploganphantom/rtkprojectshepploganphantom.py @@ -32,6 +32,7 @@ def build_parser(): ) rtk.add_rtk3Doutputimage_group(parser) + rtk.add_rtknoise_group(parser) # Parse the command line arguments return parser @@ -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: diff --git a/documentation/docs/rtk_3_migration_guide.md b/documentation/docs/rtk_3_migration_guide.md index 3c9d27b5c..2cde2e30d 100644 --- a/documentation/docs/rtk_3_migration_guide.md +++ b/documentation/docs/rtk_3_migration_guide.md @@ -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 ,`, `--gaussian ,`). The legacy `--noise` option in `rtkprojectshepploganphantom` has been removed. diff --git a/include/rtkGgoFunctions.h b/include/rtkGgoFunctions.h index 6febe9f4b..cde0f66c2 100644 --- a/include/rtkGgoFunctions.h +++ b/include/rtkGgoFunctions.h @@ -28,6 +28,12 @@ #include "rtkProjectionsReader.h" #include #include +#include +#include +#include +#include +#include +#include "rtkAdditiveGaussianNoiseImageFilter.h" namespace rtk { @@ -282,6 +288,65 @@ SetProjectionsReaderFromGgo(TProjectionsReaderType * reader, const TArgsInfo & a TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->UpdateOutputInformation()); } + +template +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::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::New(); + multiply->SetInput(projections); + multiply->SetConstant(-args_info.poisson_arg[1]); // -muref + + auto expf = itk::ExpImageFilter::New(); + expf->SetInput(multiply->GetOutput()); + + auto multiply2 = itk::MultiplyImageFilter::New(); + multiply2->SetInput(expf->GetOutput()); + multiply2->SetConstant(args_info.poisson_arg[0]); // I0 + + auto poisson = itk::ShotNoiseImageFilter::New(); + poisson->SetInput(multiply2->GetOutput()); + + auto threshold = itk::ThresholdImageFilter::New(); + threshold->SetInput(poisson->GetOutput()); + threshold->SetLower(1.); + threshold->SetOutsideValue(1.); + + auto multiply3 = itk::MultiplyImageFilter::New(); + multiply3->SetInput(threshold->GetOutput()); + multiply3->SetConstant(1. / args_info.poisson_arg[0]); + + auto logf = itk::LogImageFilter::New(); + logf->SetInput(multiply3->GetOutput()); + + auto multiply4 = itk::MultiplyImageFilter::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 diff --git a/wrapping/__init_rtk__.py b/wrapping/__init_rtk__.py index f130f0a87..df9813251 100644 --- a/wrapping/__init_rtk__.py +++ b/wrapping/__init_rtk__.py @@ -13,6 +13,7 @@ "itk.rtk3Doutputimage_group", "itk.rtkprojectors_group", "itk.rtkiterations_group", + "itk.rtknoise_group", "itk.rtkExtras", ] for mod_name in rtk_submodules: