02. Noise operators

This tutorial shows how to use noise operators using the spyrit.core.noise submodule.

Load a batch of images

We load a batch of images from the /images/ folder. Using the transform_gray_norm() function with the normalize=False argument returns images with values in (0,1).

import os

import torch
import torchvision
import matplotlib.pyplot as plt

from spyrit.misc.disp import imagesc
from spyrit.misc.statistics import transform_gray_norm

spyritPath = os.getcwd()
imgs_path = os.path.join(spyritPath, "images/")

# Grayscale images of size 64 x 64, no normalization to keep values in (0,1)
transform = transform_gray_norm(img_size=64, normalize=False)

# Create dataset and loader (expects class folder 'images/test/')
dataset = torchvision.datasets.ImageFolder(root=imgs_path, transform=transform)
dataloader = torch.utils.data.DataLoader(dataset, batch_size=7)

x, _ = next(iter(dataloader))
print(f"Shape of input images: {x.shape}")

Shape of input images: torch.Size([7, 1, 64, 64])

We select the first image in the batch and plot it.

i_plot = 1
imagesc(x[i_plot, 0, :, :], r"$x$ in (0, 1)")

/home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.1.1/lib/python3.11/site-packages/matplotlib/cbook.py:684: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments. To learn more, see the migration guide https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword
  x = np.array(x, subok=True, copy=copy)

Gaussian noise

We consider additive Gaussian noise,

\[y \sim z + \mathcal{N}(0,\sigma^2),\]

where \(\mathcal{N}(\mu, \sigma^2)\) is a Gaussian distribution with mean \(\mu\) and variance \(\sigma^2\), and \(z\) is the noiseless image. The larger \(\sigma\), the lower the signal-to-noise ratio.

To add 10% Gaussian noise, we instantiate a spyrit.core.noise operator with sigma=0.1.

from spyrit.core.noise import Gaussian

noise_op = Gaussian(sigma=0.1)
x_noisy = noise_op(x)

imagesc(x_noisy[1, 0, :, :], r"10% Gaussian noise")
# sphinx_gallery_thumbnail_number = 2

/home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.1.1/lib/python3.11/site-packages/matplotlib/cbook.py:684: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments. To learn more, see the migration guide https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword
  x = np.array(x, subok=True, copy=copy)

To add 2% Gaussian noise, we update the class attribute sigma.

noise_op.sigma = 0.02
x_noisy = noise_op(x)

imagesc(x_noisy[1, 0, :, :], r"2% Gaussian noise")

/home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.1.1/lib/python3.11/site-packages/matplotlib/cbook.py:684: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments. To learn more, see the migration guide https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword
  x = np.array(x, subok=True, copy=copy)

Poisson noise

We now consider Poisson noise,

\[y \sim \mathcal{P}(\alpha z), \quad z \ge 0,\]

where \(\alpha \ge 0\) is a scalar value that represents the maximum image intensity (in photons). The larger \(\alpha\), the higher the signal-to-noise ratio.

We consider the spyrit.core.noise.Poisson class and set \(\alpha\) to 100 photons (which corresponds to the Poisson parameter).

from spyrit.core.noise import Poisson
from spyrit.misc.disp import add_colorbar, noaxis

alpha = 100  # number of photons
noise_op = Poisson(alpha)

We simulate two noisy versions of the same images

y1 = noise_op(x)  # first sample
y2 = noise_op(x)  # another sample

We now consider the case \(\alpha = 1000\) photons.

noise_op.alpha = 1000
y3 = noise_op(x)  # noisy measurement vector

We finally plot the noisy images

# plot
f, axs = plt.subplots(1, 3, figsize=(10, 5))
axs[0].set_title("100 photons")
im = axs[0].imshow(y1[1, 0].reshape(64, 64), cmap="gray")
add_colorbar(im, "bottom")

axs[1].set_title("100 photons")
im = axs[1].imshow(y2[1, 0].reshape(64, 64), cmap="gray")
add_colorbar(im, "bottom")

axs[2].set_title("1000 photons")
im = axs[2].imshow(
    y3[
        1,
        0,
    ].reshape(64, 64),
    cmap="gray",
)
add_colorbar(im, "bottom")

noaxis(axs)

/home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.1.1/lib/python3.11/site-packages/matplotlib/cbook.py:684: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments. To learn more, see the migration guide https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword
  x = np.array(x, subok=True, copy=copy)
/home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.1.1/lib/python3.11/site-packages/matplotlib/cbook.py:684: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments. To learn more, see the migration guide https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword
  x = np.array(x, subok=True, copy=copy)
/home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.1.1/lib/python3.11/site-packages/matplotlib/cbook.py:684: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments. To learn more, see the migration guide https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword
  x = np.array(x, subok=True, copy=copy)

As expected the signal-to-noise ratio of the measurement vector is higher for 1,000 photons than for 100 photons

Note

Not only the signal-to-noise, but also the scale of the measurements depends on \(\alpha\), which motivates the introduction of the preprocessing operator.