Implemented new way of noise simulation

simulatedmicroscopy/image.py

@@ -22,8 +22,11 @@ class Image:
     is_downsampled = False
     """Whether the image has undergone downsampling"""
 
-    has_noise = False
-    """Whether the image has undergone noise addition"""
+    has_shot_noise = False
+    """Whether the image has undergone shot noise addition"""
+
+    has_read_noise = False
+    """Whether the image has undergone read noise addition"""
 
     def __init__(
         self,
@@ -382,7 +385,7 @@ def convolve(self, other: type[Image]) -> type[Image]:
         return self
 
     def noisify(self, lam: float = 1.0) -> type[Image]:
-        """Add Poisson noise to the image
+        """Add Poisson noise to the image [DEPRECATED]
 
         Parameters
         ----------
@@ -394,14 +397,73 @@ def noisify(self, lam: float = 1.0) -> type[Image]:
         type[Image]
             The image with noise added
         """
-        if self.has_noise:
-            warnings.warn("Image has already undergone noisification once")
+        raise DeprecationWarning(
+            "This method is deprecated, please use either add_shot_noise or add_read_noise"
+        )
+
+    def add_shot_noise(self, SNR: float = 30.0) -> type[Image]:
+        """Add shot noise to the image, realised by sampling the intensities of the image from a Poisson distribution.
+        At each pixel, the current value (after rescaling) dictates the mean value of the Poisson distribution, yielding
+        signal-dependent noise.
+
+        Parameters
+        ----------
+        SNR : float, optional
+            The signal-to-noise ratio of the resulting image, by default 30.0.
+
+        Returns
+        -------
+        type[Image]
+            The image with shot noise
+        """
+        if self.has_shot_noise:
+            warnings.warn("This image already has shot noise, proceeding anyway")
+
+        scaling_factor = (
+            SNR**2 / self.image.mean()
+        )  # scaling factor to get a proper Poisson sampling
+
+        self.image = (
+            np.random.poisson(scaling_factor * self.image.mean()) / scaling_factor
+        )
+        # division by `scaling_factor` to revert pixel values back to their (roughly) original range
+
+        self.has_shot_noise = True  # set shot noise flag
+
+        return self
+
+    def add_read_noise(self, SNR: float = 50.0, background: float = 0.0) -> type[Image]:
+        """Add read noise to the image, realised by adding noise sampled from a normal distribution. The standard deviation
+        of the distribution is constant across the entire image and determined by the SNR parameter, it is therefore
+        signal-independent.
+
+        Parameters
+        ----------
+        SNR : float, optional
+            The signal-to-noise ratio of the resulting image, by default 50.0
+        background : float, optional
+            Optional offset of the background, this will be the mean of the normal distribution, by default 0.0
+
+        Returns
+        -------
+        type[Image]
+            The image with read noise
+        """
+        if self.has_read_noise:
+            warnings.warn("This image already has read noise, proceeding anyway")
+
+        peak = np.percentile(
+            self.image, 99.9
+        )  # get peak intensity at the 99.9 percentile
+
+        # As SNR = peak / sigma; this yields that sigma = peak/SNR
+        sigma = peak / SNR
 
-        rng = np.random.default_rng()
-        self.image = self.image * rng.poisson(lam, size=self.image.shape)
+        self.image += np.random.normal(
+            loc=background, scale=sigma
+        )  # centred around value of `background`
 
-        self.has_noise = True
-        self.metadata["has_noise"] = True
+        self.has_read_noise = True  # set read noise flag
 
         return self
 

simulatedmicroscopy/psf.py

@@ -54,9 +54,7 @@ def __init__(
         mu = [0.0] * 3
 
         # generate Gaussian with given mu and sigmas
-        gauss_psf = scipy.stats.multivariate_normal.pdf(
-            zyx, mu, np.diag(sigmas_nm**2)
-        )
+        gauss_psf = scipy.stats.multivariate_normal.pdf(zyx, mu, np.diag(sigmas_nm**2))
 
         # reshape to 3D image
         gauss_psf = gauss_psf.reshape(z.shape)

tests/test_image.py

@@ -154,18 +154,34 @@ def test_convolution_wrongpixelsize():
         im.convolve(psf)
 
 
-def test_noise():
+def test_noise_deprecated():
+    with pytest.raises(DeprecationWarning):
+        im = create_demo_image()
+        im.noisify()
+
+def test_shot_noise():
     im = create_demo_image()
 
     # check if returns correct image
-    assert im.noisify() != create_demo_image()
+    assert im.add_shot_noise(SNR = 10.0) != create_demo_image()
 
     # check if original image was also changed
     assert im != create_demo_image()
 
     # should be set to True
-    assert im.has_noise
+    assert im.has_shot_noise
+
+def test_read_noise():
+    im = create_demo_image()
+
+    # check if returns correct image
+    assert im.add_read_noise(SNR = 10.0, background = 1e-3) != create_demo_image()
 
+    # check if original image was also changed
+    assert im != create_demo_image()
+
+    # should be set to True
+    assert im.has_read_noise
 
 def test_point_image(tmp_path):
     from simulatedmicroscopy import Coordinates