Poisson noise generation script for PETRIC2

CHANGES.md

@@ -2,10 +2,11 @@
 ## vx.x.x
 
 * SIRF/STIR
-  - `ScatterEstimation` has extra methods that allow setting masks for the tail-fitting
+  - `ScatterEstimation` has extra methods that allow setting masks for the tail-fitting.
   - `ImageData` has extra method to zoom image using information from a template image, `zoom_image_as_template`.
   - Error raised in `AcquisitionSensitivityModel.[un]normalise` methods applied to a read-only object.
   - `DataContainer.supports_array_view` to test for zero-copy compatibility.
+  - SIRF interfaces (C++ and Python) for STIR Poisson noise generation utilities provided.
 
 ## v3.8.1
 

examples/Python/PET/generate_noisy_data.py

@@ -0,0 +1,110 @@
+'''Poisson noise generation demo.
+
+Usage:
+  generate_noisy_data [--help | options]
+
+Options:
+  -f <file>, --file=<file>     prompts file [default: my_forward_projection.hs]
+  -p <path>, --path=<path>     path to data files, defaults to data/examples/PET
+                               subfolder of SIRF root folder
+  -o <file>, --output=<file>   output file for Poisson noisy data
+  -r <seed>, --seed=<seed>     random generator seed [default: 1]
+  -F <sf>, --sf=<sf>           scaling factor [default: 1.0] (use a higher value for lower relative noise)
+  -m, --pm                     preserve mean
+  -e <engn>, --engine=<engn>   reconstruction engine [default: STIR]
+  -s <stsc>, --storage=<stsc>  acquisition data storage scheme [default: file]
+  --non-interactive            do not show plots
+'''
+
+## SyneRBI Synergistic Image Reconstruction Framework (SIRF)
+## Copyright 2015 - 2025 Rutherford Appleton Laboratory STFC
+## Copyright 2015 - 2025 University College London.
+##
+## This is software developed for the Collaborative Computational
+## Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
+## (http://www.ccpsynerbi.ac.uk/).
+##
+## Licensed under the Apache License, Version 2.0 (the "License");
+##   you may not use this file except in compliance with the License.
+##   You may obtain a copy of the License at
+##       http://www.apache.org/licenses/LICENSE-2.0
+##   Unless required by applicable law or agreed to in writing, software
+##   distributed under the License is distributed on an "AS IS" BASIS,
+##   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+##   See the License for the specific language governing permissions and
+##   limitations under the License.
+
+__version__ = '0.1.0'
+from docopt import docopt
+args = docopt(__doc__, version=__version__)
+
+from sirf.Utilities import error, examples_data_path, existing_filepath
+
+# import engine module
+import importlib
+engine = args['--engine']
+pet = importlib.import_module('sirf.' + engine)
+
+# process command-line options
+data_file = args['--file']
+data_path = args['--path']
+if data_path is None:
+    data_path = examples_data_path('PET')
+output_file = args['--output']
+seed = int(args['--seed'])
+sf = float(args['--sf'])
+pm = bool(args['--pm'])
+storage = args['--storage']
+show_plot = not args['--non-interactive']
+
+# select acquisition data storage scheme
+# storage = 'file' (default):
+#     all acquisition data generated by the script is kept in
+#     scratch files deleted after the script terminates
+# storage = 'memory':
+#     all acquisition data generated by the script is kept in RAM
+#     (avoid if data is very large)
+scheme = pet.AcquisitionData.get_storage_scheme()
+if scheme != storage:
+    print('default storage scheme is %s' % repr(scheme))
+    print('setting storage scheme to %s' % repr(storage))
+    pet.AcquisitionData.set_storage_scheme(storage)
+else:
+    print('using default storage scheme %s' % repr(scheme))
+
+
+def main():
+    engine_version = pet.get_engine_version_string()
+    print('using %s version %s as the reconstruction engine' % (engine, engine_version))
+
+    # direct all engine's messages to files
+    _ = pet.MessageRedirector('info.txt', 'warn.txt', 'errr.txt')
+
+    # PET acquisition data to be read from this file
+    prompts_file = existing_filepath(data_path, data_file)
+    print('reading prompts from %s' % prompts_file)
+    acq_data = pet.AcquisitionData(prompts_file)
+    dim = acq_data.dimensions()
+    print(f'acquisition data norm: {acq_data.norm()}')
+
+    png = pet.PoissonNoiseGenerator(sf, pm)
+    png.set_seed(seed)
+    png.process(acq_data)
+    noisy_data = png.get_output()
+    print(f'noisy data norm: {noisy_data.norm()}')
+    if show_plot:
+        noisy_data.show(range(dim[1]//4), title='selected noisy sinograms')
+
+    if output_file is not None:
+        noisy_data.write(output_file)
+
+
+try:
+    main()
+    print('\n=== done with %s' % __file__)
+
+except error as err:
+    print('%s' % err.value)
+
+if scheme != storage:
+    pet.AcquisitionData.set_storage_scheme(scheme)

src/xSTIR/cSTIR/cstir.cpp

@@ -204,6 +204,8 @@ void* cSTIR_newObject(const char* name)
                   return NEW_OBJECT_HANDLE(PETScatterEstimator);
 		if (sirf::iequals(name, "SeparableGaussianImageFilter"))
 			return NEW_OBJECT_HANDLE(xSTIR_SeparableGaussianImageFilter);
+                if (sirf::iequals(name, "PoissonNoiseGenerator"))
+                  return NEW_OBJECT_HANDLE(PoissonNoiseGenerator);
 		return unknownObject("object", name, __FILE__, __LINE__);
 	}
 	CATCH;
@@ -290,6 +292,8 @@ void* cSTIR_setParameter
                         return cSTIR_setScatterSimulatorParameter(hs, name, hv);
                 else if(sirf::iequals(obj, "PETScatterEstimator"))
                         return cSTIR_setScatterEstimatorParameter(hs, name, hv);
+                else if(sirf::iequals(obj, "PoissonNoiseGenerator"))
+                        return cSTIR_setPoissonNoiseGeneratorParameter(hs, name, hv);
 		else
 			return unknownObject("object", obj, __FILE__, __LINE__);
 	}
@@ -360,6 +364,8 @@ void* cSTIR_parameter(const void* ptr, const char* obj, const char* name)
 			return cSTIR_FBP2DParameter(handle, name);
                 else if(sirf::iequals(obj, "PETScatterEstimator"))
                         return cSTIR_ScatterEstimatorParameter(handle, name);
+//                else if(sirf::iequals(obj, "PoissonNoiseGenerator"))
+//                        return cSTIR_PoissonNoiseGeneratorParameter(handle, name);
 		return unknownObject("object", obj, __FILE__, __LINE__);
 	}
 	CATCH;
@@ -613,6 +619,30 @@ void* cSTIR_applyImageDataProcessor(const void* ptr_p, void* ptr_i)
 	CATCH;
 }
 
+extern "C"
+void* cSTIR_createPoissonNoiseGenerator(const float scaling_factor, const bool preserve_mean)
+{
+	try {
+		shared_ptr<PoissonNoiseGenerator> 
+			sptr(new PoissonNoiseGenerator(scaling_factor, preserve_mean));
+		return newObjectHandle(sptr);
+	}
+	CATCH;
+}
+
+extern "C"
+void* cSTIR_generatePoissonNoise(const void* ptr_gen, const void* ptr_input)
+{
+	try {
+		auto& generator = objectFromHandle<PoissonNoiseGenerator>(ptr_gen);
+		auto& input = objectFromHandle<STIRAcquisitionData>(ptr_input);
+		auto sptr_output = input.new_acquisition_data();
+		generator.generate_random(*sptr_output, input);
+		return newObjectHandle(sptr_output);
+	}
+	CATCH;
+}
+
 extern "C"
 void* cSTIR_createPETAcquisitionSensitivityModel
 	(const void* ptr_src, const char* src)

src/xSTIR/cSTIR/cstir_p.cpp

@@ -919,6 +919,16 @@ sirf::cSTIR_ScatterEstimatorParameter(DataHandle* hp, const char* name)
 	return parameterNotFound(name, __FILE__, __LINE__);
 }
 
+void*
+sirf::cSTIR_setPoissonNoiseGeneratorParameter
+(const DataHandle *hp, const char* name, const DataHandle* hv)
+{
+    auto& obj = objectFromHandle<PoissonNoiseGenerator>(hp);
+    if (sirf::iequals(name, "seed"))
+        obj.seed(dataFromHandle<int>(hv));
+    return new DataHandle;
+}
+
 void*
 sirf::cSTIR_setGeneralisedObjectiveFunctionParameter
 (DataHandle* hp, const char* name, const DataHandle* hv)

src/xSTIR/cSTIR/include/sirf/STIR/cstir.h

@@ -76,6 +76,13 @@ extern "C" {
 		(const void* ptr_src, const char* src);
 	void* cSTIR_createPETAttenuationModel
 		(const void* ptr_img, const void* ptr_am);
+
+	// Poisson noise generator
+	void* cSTIR_createPoissonNoiseGenerator
+		(const float scaling_factor, const bool preserve_mean);
+	void* cSTIR_generatePoissonNoise
+		(const void* ptr_gen, const void* ptr_input);
+
 	void* cSTIR_computeACF
 		(const void* ptr_sino, const void* ptr_att, void* ptr_acf, void* ptr_iacf);
 	void* cSTIR_chainPETAcquisitionSensitivityModels

src/xSTIR/cSTIR/include/sirf/STIR/cstir_p.h

@@ -168,6 +168,10 @@ namespace sirf {
                 cSTIR_ScatterEstimatorParameter
                 (DataHandle* hp, const char* name);
 
+        void*
+                cSTIR_setPoissonNoiseGeneratorParameter
+                (const DataHandle *hp, const char* name, const DataHandle* hv);
+
 	void*
 		cSTIR_setGeneralisedObjectiveFunctionParameter
 		(DataHandle* hp, const char* name, const DataHandle* hv);

src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h

@@ -41,6 +41,7 @@ limitations under the License.
 #include "sirf/STIR/stir_data_containers.h"
 #include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h"
 #include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h"
+#include "stir/GeneralisedPoissonNoiseGenerator.h"
 
 #define MIN_BIN_EFFICIENCY 1.0e-20f
 //#define MIN_BIN_EFFICIENCY 1.0e-6f
@@ -107,6 +108,43 @@ namespace sirf {
 		//shared_ptr<stir::ChainedBinNormalisation> norm_;
 	};
 
+	/*!
+	\ingroup PET
+	\brief Class for Poisson noise generation.
+
+	Generates noise realisations according to Poisson statistics but allowing for scaling.
+
+	A scaling_factor is used to multiply the input data before generating
+	the Poisson random number. This means that a scaling_factor larger than 1
+	will result in data with lower relative noise.
+
+	If preserve_mean=false, the mean of the output data will
+	be equal to scaling_factor*mean_of_input, otherwise it
+	will be equal to mean_of_input, but then the output is no longer Poisson
+	distributed.
+	*/
+
+	class PoissonNoiseGenerator {
+	public:
+		//! Constructor intialises the random number generator with a fixed seed
+		PoissonNoiseGenerator(const float scaling_factor = 1.0F, const bool preserve_mean = false)
+		{
+			gpng_ = stir::shared_ptr<stir::GeneralisedPoissonNoiseGenerator>
+				(new GeneralisedPoissonNoiseGenerator(scaling_factor, preserve_mean));
+		}
+		//! Sets the seed value for the random number generator
+		void seed(unsigned int s)
+		{
+			gpng_->seed(s);
+		}
+		void generate_random(STIRAcquisitionData& output, const STIRAcquisitionData& input)
+		{
+			gpng_->generate_random(*output.data(), *input.data());
+		}
+
+	protected:
+		stir::shared_ptr<stir::GeneralisedPoissonNoiseGenerator> gpng_;
+	};
 
         /*!
 	\ingroup PET

src/xSTIR/pSTIR/STIR.py

@@ -1625,6 +1625,47 @@ def get_time_at_which_num_prompts_exceeds_threshold(self, threshold):
         return v
 
 
+class PoissonNoiseGenerator(object):
+    """
+    Generates noise realisations according to Poisson statistics but allowing for scaling.
+
+    A scaling_factor is used to multiply the input data before generating
+    the Poisson random number. This means that a scaling_factor larger than 1
+    will result in data with lower relative noise.
+
+    If preserve_mean=false, the mean of the output data will
+    be equal to scaling_factor*mean_of_input, otherwise it
+    will be equal to mean_of_input, but then the output is no longer Poisson
+    distributed.
+    """
+
+    def __init__(self, scaling_factor=1.0, preserve_mean=False):
+        self.name = "PoissonNoiseGenerator"
+        self.handle = pystir.cSTIR_createPoissonNoiseGenerator(scaling_factor, preserve_mean)
+        check_status(self.handle)
+        self.output_handle = None
+
+    def set_seed(self, s):
+        parms.set_int_par(self.handle, self.name, 'seed', s)
+
+    def process(self, acq_data):
+        self.output_handle = pystir.cSTIR_generatePoissonNoise(self.handle, acq_data.handle)
+        check_status(self.output_handle)
+
+    def get_output(self):
+        if self.output_handle is None:
+            raise error('Noise generating not done')
+        output = AcquisitionData()
+        output.handle = self.output_handle
+        return output
+
+    def generate_noisy_data(self, acq_data):
+        noisy_data = AcquisitionData()
+        noisy_data.handle = pystir.cSTIR_generatePoissonNoise(self.handle, acq_data.handle)
+        check_status(noisy_data.handle)
+        return noisy_data
+
+
 class AcquisitionSensitivityModel(object):
     """
     Class that handles PET scanner detector efficiencies and attenuation.