Poisson noise generation script for PETRIC2

[evgueni-ovtchinnikov]

Changes in this pull request

Poisson noise generation script for PETRIC2 to be added.

Testing performed

A simple example of Poisson noise generation added (temporarily?) to examples/Python/PET/acquisition_data.py

Related issues

Checklist before requesting a review

• I have performed a self-review of my code
• I have added docstrings/doxygen in line with the guidance in the developer guide
• I have implemented unit tests that cover any new or modified functionality
• The code builds and runs on my machine
• CHANGES.md has been updated with any functionality change

Contribution Notes

Please read and adhere to the contribution guidelines.

Please tick the following:

• The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in SIRF (the Work) under the terms and conditions of the Apache-2.0 License.

[evgueni-ovtchinnikov]

@KrisThielemans I am afraid I find your instructions re Poisson noise generation script a bit cryptic:

• Read in prompts (EO: ok)
• Scale (EO: scale prompts? What scale factor to use? Or is this merely an alternative to using scaling_factor?)
• Generate Poisson realisation (EO: ok, but with what scaling_factor and preserve_mean? Use defaults if prompts are already scaled?)
• Also scale additive term (EO: what scale? Same as with prompts? And how the additive term participates in noise generation?)

I still hope that what I have done so far based on poisson_noise.cxx and common sense does make sense (please have a look at PR #1321 - your comments and more developer-friendly instructions would be very much appreciated. :-))

[KrisThielemans]

ah, this is far more than I thought. https://github.com/SyneRBI/SIRF-Exercises/blob/84eee654ad2ee2dee0387b783e5477174e1c1dc7/notebooks/PET/ML_reconstruction.ipynb has an example where numpy.random.poisson is used. That could be used directly in PETRIC repo without any SIRF modifications.

However, exposing the STIR Poisson generator is of course a good strategy. It all looks fine to me (aside from missing documentation).

• Read in prompts (EO: ok)
• Scale (EO: scale prompts? What scale factor to use? Or is this merely an alternative to using scaling_factor?)

• for PETRIC, we'll have to scale prompts and additive term, I haven't thought through yet what we'd do with the OSEM_image. In principle, people would have to recompute the OSEM image (otherwise it's cheating), but in practice for initialisation, they could experiment with just scaling the OSEM image in the same way.
• it is the same as scaling_factor (numpy.random.poisson doesn't have it)

• Generate Poisson realisation (EO: ok, but with what scaling_factor and preserve_mean? Use defaults if prompts are already scaled?)

we'll have to find the scaling_factor for each data-set as discussed during the meeting. I would set preserve_mean=false. (It's a trick to be able to reconstruct the image in the same scale, but messes up Poisson stats, penalties etc)

• Also scale additive term (EO: what scale? Same as with prompts? And how the additive term participates in noise generation?)

yes. same scale. The additive term is part of the acq_model is NOT noisy. So, no noise generation for it.

The idea is to provide a function in PETRIC repo that does all this (with the scale factor as input). I wouldn't put it in SIRF.

[KrisThielemans]

some comments.

Shouldn't we follow our standard scheme with process() and get_output, or whatever we do in other places?

[KrisThielemans]

let's be explicit (STIR isn't) and say ""will result in data with lower relative noise"

[KrisThielemans]

I would rather have had process() and get_output() in the C++ version as well for 2 reasons:

• It's better for the developer if C++ and Python SIRF are similar
• I still hope to be able to get rid of the C layer by using SWIG in the future

but this has taken too much of our time already, so fine to merge.