.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tuto_01_acquisition_operators.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_tuto_01_acquisition_operators.py: 01. Acquisition operators ========================== .. _tuto_acquisition_operators: This tutorial shows how to simulate measurements using the :mod:`spyrit.core` submodule. The simulation is based on three modules: 1. **Measurement operators** compute linear measurements :math:`y = Hx` from images :math:`x`, where :math:`H` is a linear operator (matrix) and :math:`x` is an image (see :mod:`spyrit.core.meas`) 2. **Noise operator** corrupts measurements :math:`y` with noise (see :mod:`spyrit.core.noise`) 3. **Preprocessing operators** are typically used to process the noisy measurements prior to reconstruction (see :mod:`spyrit.core.prep`) .. image:: ../fig/tuto1.png :width: 600 :align: center :alt: Measurement, noise, and preprocessing sketches These tutorials load image samples from `/images/`. Please note that as of v.2.4.0, the inputs to the measurement operators are no longer vectorized images, but rather image tensors with shape :math:`(*, H, W)`, where :math:`*` represents any number of additional dimensions, e.g. batch size and number of channels. .. GENERATED FROM PYTHON SOURCE LINES 32-34 Load a batch of images ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 36-42 Images :math:`x` for training neural networks expect values in [-1,1]. The images are normalized using the :func:`transform_gray_norm` function. Spyrit can handle images with the shape :math:`(h, w)` or :math:`(*, h, w)`, where :math:`*` represents any number of additional dimensions, e.g. batch size and number of channels. In this case, we load a batch of black and white images of size :math:`64 \times 64`, and select one image for the tutorial. This results in a tensor of shape :math:`(1, 1, 64, 64)`. .. GENERATED FROM PYTHON SOURCE LINES 42-79 .. code-block:: Python 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 # sphinx_gallery_thumbnail_path = 'fig/tuto1.png' h = 64 # image size hxh i = 1 # Image index (modify to change the image) spyritPath = os.getcwd() imgs_path = os.path.join(spyritPath, "images/") # Create a transform for natural images to normalized grayscale image tensors transform = transform_gray_norm(img_size=h) # 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}") # Select image x = x[i : i + 1, :, :, :] x = x.detach().clone() print(f"Shape of selected image: {x.shape}") b, c, h, w = x.shape # plot imagesc(x[0, 0, :, :], r"$x$ in [-1, 1]") .. image-sg:: /gallery/images/sphx_glr_tuto_01_acquisition_operators_001.png :alt: $x$ in [-1, 1] :srcset: /gallery/images/sphx_glr_tuto_01_acquisition_operators_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Shape of input images: torch.Size([7, 1, 64, 64]) Shape of selected image: torch.Size([1, 1, 64, 64]) /home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/2.4.0/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) .. GENERATED FROM PYTHON SOURCE LINES 80-82 The measurement and noise operators ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 84-104 Noise operators are defined in the :mod:`~spyrit.core.noise` module. A noise operator computes the following three steps sequentially: 1. Normalization of the image :math:`x` with values in [-1,1] to get an image :math:`\tilde{x}=\frac{x+1}{2}` in [0,1], as it is required for measurement simulation 2. Application of the measurement model, i.e., computation of :math:`H\tilde{x}` 3. Application of the noise model .. math:: y \sim \texttt{Noise}(H\tilde{x}) = \texttt{Noise}\left(\frac{H(x+1)}{2}\right), The normalization is usefull when considering distributions such as the Poisson distribution that are defined on positive values. .. note:: The noise operator is constructed from a measurement operator (see the :mod:`~spyrit.core.meas` submodule) in order to compute the measurements :math:`H\tilde{x}`, as given by step #2. .. GENERATED FROM PYTHON SOURCE LINES 107-109 A simple example: identity measurement matrix and no noise ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 111-113 .. math:: y = \tilde{x} .. GENERATED FROM PYTHON SOURCE LINES 115-124 We start with a simple example where the measurement matrix :math:`H` is the identity, which can be handled by the more general :class:`spyrit.core.meas.Linear` class. We consider the noiseless case handled by the :class:`spyrit.core.noise.NoNoise` class. Usually, the measurement tensor is in another space than the image tensor (e.g. Fourier space or Hadamard space), but using the identity matrix results in the measurement vector being (identical and) in the same space as the image tensor. As measurements are always vectorized, the measurement vector is a vectorized image. .. GENERATED FROM PYTHON SOURCE LINES 124-131 .. code-block:: Python from spyrit.core.meas import Linear from spyrit.core.noise import NoNoise meas_op = Linear(torch.eye(h * h)) noise_op = NoNoise(meas_op) .. GENERATED FROM PYTHON SOURCE LINES 132-135 We simulate the measurement vector :math:`y` that we want to visualise as an image. Note that the measurement vector :math:`y` lost a dimension compared to the image :math:`x`, because the measurement operator acts on the last 2 dimensions of the image tensor. .. GENERATED FROM PYTHON SOURCE LINES 135-142 .. code-block:: Python y_eye = noise_op(x) # noisy measurement vector print(f"Shape of simulated measurements y: {y_eye.shape}") # plot imagesc(y_eye[0, 0, :].reshape(h, h), r"$\tilde{x}$ in [0, 1]") .. image-sg:: /gallery/images/sphx_glr_tuto_01_acquisition_operators_002.png :alt: $\tilde{x}$ in [0, 1] :srcset: /gallery/images/sphx_glr_tuto_01_acquisition_operators_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Shape of simulated measurements y: torch.Size([1, 1, 4096]) /home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/2.4.0/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) .. GENERATED FROM PYTHON SOURCE LINES 143-146 .. note:: Note that the image is identical to the original one, except it has been normalized in [0,1]. .. GENERATED FROM PYTHON SOURCE LINES 148-150 Same example with Poisson noise ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 152-159 We now consider Poisson noise, i.e., a noisy measurement vector given by .. math:: y \sim \mathcal{P}(\alpha H \tilde{x}), where :math:`\alpha` is a scalar value that represents the maximum image intensity (in photons). The larger :math:`\alpha`, the higher the signal-to-noise ratio. .. GENERATED FROM PYTHON SOURCE LINES 162-164 We consider the :class:`spyrit.core.noise.Poisson` class and set :math:`\alpha` to 100 photons (which corresponds to the Poisson parameter). .. GENERATED FROM PYTHON SOURCE LINES 164-171 .. code-block:: Python from spyrit.core.noise import Poisson from spyrit.misc.disp import add_colorbar, noaxis alpha = 100 # number of photons noise_op = Poisson(meas_op, alpha) .. GENERATED FROM PYTHON SOURCE LINES 172-173 We simulate two noisy measurement vectors .. GENERATED FROM PYTHON SOURCE LINES 173-177 .. code-block:: Python y1 = noise_op(x) # a noisy measurement vector y2 = noise_op(x) # another noisy measurement vector .. GENERATED FROM PYTHON SOURCE LINES 178-179 We now consider the case :math:`\alpha = 1000` photons. .. GENERATED FROM PYTHON SOURCE LINES 179-183 .. code-block:: Python noise_op.alpha = 1000 y3 = noise_op(x) # noisy measurement vector .. GENERATED FROM PYTHON SOURCE LINES 184-185 We finally plot the measurement vectors as images .. GENERATED FROM PYTHON SOURCE LINES 185-203 .. code-block:: Python # plot f, axs = plt.subplots(1, 3, figsize=(10, 5)) axs[0].set_title("100 photons") im = axs[0].imshow(y1[0, 0, :].reshape(h, h), cmap="gray") add_colorbar(im, "bottom") axs[1].set_title("100 photons") im = axs[1].imshow(y2[0, 0, :].reshape(h, h), cmap="gray") add_colorbar(im, "bottom") axs[2].set_title("1000 photons") im = axs[2].imshow(y3[0, 0, :].reshape(h, h), cmap="gray") add_colorbar(im, "bottom") noaxis(axs) plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_01_acquisition_operators_003.png :alt: 100 photons, 100 photons, 1000 photons :srcset: /gallery/images/sphx_glr_tuto_01_acquisition_operators_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/2.4.0/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/2.4.0/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/2.4.0/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) .. GENERATED FROM PYTHON SOURCE LINES 204-211 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 :math:`\alpha`, which motivates the introduction of the preprocessing operator. .. GENERATED FROM PYTHON SOURCE LINES 213-215 The preprocessing operator ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 217-228 Preprocessing operators are defined in the :mod:`spyrit.core.prep` module. A preprocessing operator applies to the noisy measurements .. math:: m = \texttt{Prep}(y), For instance, a preprocessing operator can be used to compensate for the scaling factors that appear in the measurement or noise operators. In this case, a preprocessing operator is closely linked to its measurement and/or noise operator counterpart. While scaling factors are required to simulate realistic measurements, they are not required for reconstruction. .. GENERATED FROM PYTHON SOURCE LINES 230-232 Preprocessing measurements corrupted by Poisson noise ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 234-246 We consider the :class:`spyrit.core.prep.DirectPoisson` class that intends to "undo" the :class:`spyrit.core.noise.Poisson` class by compensating for: * the scaling that appears when computing Poisson-corrupted measurements * the affine transformation to get images in [0,1] from images in [-1,1] For this, it computes .. math:: m = \frac{2}{\alpha} y - H1 .. GENERATED FROM PYTHON SOURCE LINES 248-250 We consider the :class:`spyrit.core.prep.DirectPoisson` class and set :math:`\alpha` to 100 photons. .. GENERATED FROM PYTHON SOURCE LINES 250-256 .. code-block:: Python from spyrit.core.prep import DirectPoisson alpha = 100 # number of photons prep_op = DirectPoisson(alpha, meas_op) .. GENERATED FROM PYTHON SOURCE LINES 257-258 We preprocess the first two noisy measurement vectors .. GENERATED FROM PYTHON SOURCE LINES 258-262 .. code-block:: Python m1 = prep_op(y1) m2 = prep_op(y2) .. GENERATED FROM PYTHON SOURCE LINES 263-265 We now consider the case :math:`\alpha = 1000` photons to preprocess the third measurement vector .. GENERATED FROM PYTHON SOURCE LINES 265-269 .. code-block:: Python prep_op.alpha = 1000 m3 = prep_op(y3) .. GENERATED FROM PYTHON SOURCE LINES 270-271 We finally plot the preprocessed measurement vectors as images .. GENERATED FROM PYTHON SOURCE LINES 271-289 .. code-block:: Python # plot f, axs = plt.subplots(1, 3, figsize=(10, 5)) axs[0].set_title("100 photons") im = axs[0].imshow(m1[0, 0, :].reshape(h, h), cmap="gray") add_colorbar(im, "bottom") axs[1].set_title("100 photons") im = axs[1].imshow(m2[0, 0, :].reshape(h, h), cmap="gray") add_colorbar(im, "bottom") axs[2].set_title("1000 photons") im = axs[2].imshow(m3[0, 0, :].reshape(h, h), cmap="gray") add_colorbar(im, "bottom") noaxis(axs) plt.show() .. image-sg:: /gallery/images/sphx_glr_tuto_01_acquisition_operators_004.png :alt: 100 photons, 100 photons, 1000 photons :srcset: /gallery/images/sphx_glr_tuto_01_acquisition_operators_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/2.4.0/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/2.4.0/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/2.4.0/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) .. GENERATED FROM PYTHON SOURCE LINES 290-294 .. note:: The preprocessed measurements still have different the signal-to-noise ratios depending on :math:`\alpha`; however, they (approximately) all lie within the same range (here, [-1, 1]). .. GENERATED FROM PYTHON SOURCE LINES 297-298 We show again one of the preprocessed measurement vectors (tutorial thumbnail purpose) .. GENERATED FROM PYTHON SOURCE LINES 298-301 .. code-block:: Python # Plot imagesc(m2[0, 0, :].reshape(h, h), "100 photons", title_fontsize=20) .. image-sg:: /gallery/images/sphx_glr_tuto_01_acquisition_operators_005.png :alt: 100 photons :srcset: /gallery/images/sphx_glr_tuto_01_acquisition_operators_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/2.4.0/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) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.182 seconds) .. _sphx_glr_download_gallery_tuto_01_acquisition_operators.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tuto_01_acquisition_operators.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tuto_01_acquisition_operators.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: tuto_01_acquisition_operators.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_