From 2a4491eb4c0316c9e7c411cfa1cccc525c8292e9 Mon Sep 17 00:00:00 2001 From: leoca Date: Tue, 15 Sep 2026 05:31:27 +0200 Subject: [PATCH] CTnoise.add: use MATLAB's default photon count and clamp non-positive readings With Poisson left to its default, the Python port set the photon count to ceil(log2(max(|proj|))), which is 9 for projections peaking at 300. That turns the requested realistic noise into noise as large as the signal, and the Gaussian term then pushes many readings to zero or below, so -log() returns NaN and inf. MATLAB's addCTnoise uses I0 = 60000 (max(proj)/5 if the projections exceed it) and sets readings <= 0 to 1e-6 before the log; do the same here. Fixes CERN/TIGRE#780 --- Python/tigre/utilities/CTnoise.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/Python/tigre/utilities/CTnoise.py b/Python/tigre/utilities/CTnoise.py index a009132b..8fd4a553 100644 --- a/Python/tigre/utilities/CTnoise.py +++ b/Python/tigre/utilities/CTnoise.py @@ -9,7 +9,10 @@ def add(projections, Gaussian=None, Poisson=None): "Poisson value should be an scalar, is " + str(type(Poisson)) + " instead." ) else: - Poisson = np.ceil(np.log2(np.max(np.abs(projections)))) # nextpow2 + # same default photon count as MATLAB's addCTnoise + Poisson = 60000 + if np.max(projections) > Poisson: + Poisson = np.max(projections) / 5 if Gaussian is not None: if not isinstance(Gaussian, np.ndarray): raise ValueError( @@ -23,6 +26,8 @@ def add(projections, Gaussian=None, Poisson=None): projections = Poisson * np.exp(-projections / max_proj) projections = RNG.add_noise(projections, Gaussian[0], Gaussian[1]) + # the Gaussian term can push a reading to zero or below, where the log is not finite + projections[projections <= 0] = 1e-6 projections = -np.log(projections / Poisson) * max_proj projections = np.float32(projections)