Poisson noise generation script for PETRIC2

CHANGES.md

@@ -6,7 +6,7 @@
   - `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 Python interface for STIR Poisson noise generation provided.
+  - SIRF interfaces (C++ and Python) for STIR Poisson noise generation utilities provided.
 
 ## v3.8.1
 

examples/Python/PET/generate_noisy_data.py

@@ -9,7 +9,7 @@
                                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]
+  -F <sf>, --sf=<sf>           scaling factor [default: 1.0] (use a higher value for higher relative noise)
   -m, --pm                     preserve mean
   -e <engn>, --engine=<engn>   reconstruction engine [default: STIR]
   -s <stsc>, --storage=<stsc>  acquisition data storage scheme [default: file]
@@ -89,7 +89,8 @@ def main():
 
     png = pet.PoissonNoiseGenerator(sf, pm)
     png.set_seed(seed)
-    noisy_data = png.generate_noisy_data(acq_data)
+    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')

src/xSTIR/cSTIR/cstir.cpp

@@ -634,10 +634,9 @@ extern "C"
 void* cSTIR_generatePoissonNoise(const void* ptr_gen, const void* ptr_input)
 {
 	try {
-		PoissonNoiseGenerator& generator =
-			objectFromHandle<PoissonNoiseGenerator>(ptr_gen);
-		STIRAcquisitionData& input = objectFromHandle<STIRAcquisitionData>(ptr_input);
-		std::shared_ptr<STIRAcquisitionData> sptr_output = input.new_acquisition_data();
+		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);
 	}

src/xSTIR/cSTIR/cstir_p.cpp

@@ -923,8 +923,7 @@ void*
 sirf::cSTIR_setPoissonNoiseGeneratorParameter
 (const DataHandle *hp, const char* name, const DataHandle* hv)
 {
-    PoissonNoiseGenerator& obj =
-            objectFromHandle<PoissonNoiseGenerator>(hp);
+    auto& obj = objectFromHandle<PoissonNoiseGenerator>(hp);
     if (sirf::iequals(name, "seed"))
         obj.seed(dataFromHandle<int>(hv));
     return new DataHandle;

src/xSTIR/pSTIR/STIR.py

@@ -1627,22 +1627,40 @@ def get_time_at_which_num_prompts_exceeds_threshold(self, threshold):
 
 class PoissonNoiseGenerator(object):
     """
-    Class that generates Poisson noise.
+    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 less noisy data.
+
+    If preserve_mean=false, the mean of the output data will
+    be equal to <tt>scaling_factor*mean_of_input</tt>, 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.scale = scaling_factor
-
-    @property
-    def scaling_factor(self):
-        return self.scale
+        self.scaling_factor = scaling_factor
+        self.preserve_mean = preserve_mean
+        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)