.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tuto_03_pseudoinverse_linear.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_03_pseudoinverse_linear.py: 03. Pseudoinverse solution from linear measurements =================================================== .. _tuto_pseudoinverse_linear: This tutorial shows how to simulate measurements and perform image reconstruction using the :class:`spyrit.core.inverse.PseudoInverse` class of the :mod:`spyrit.core.inverse` submodule. .. image:: ../fig/tuto3_pinv.png :width: 600 :align: center :alt: Reconstruction architecture sketch | .. GENERATED FROM PYTHON SOURCE LINES 17-19 Loads images ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 21-24 We load a batch of images from the :attr:`/images/` folder. Using the :func:`spyrit.misc.statistics.transform_gray_norm` function with the :attr:`normalize=False` argument returns images with values in (0,1). .. GENERATED FROM PYTHON SOURCE LINES 24-43 .. code-block:: Python import os import torchvision import torch.nn from spyrit.misc.statistics import transform_gray_norm spyritPath = os.getcwd() imgs_path = os.path.join(spyritPath, "images/") # Grayscale images of size 32 x 32, 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"Ground-truth images: {x.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Ground-truth images: torch.Size([7, 1, 64, 64]) .. GENERATED FROM PYTHON SOURCE LINES 44-46 Linear measurements without noise ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 48-49 We consider a Hadamard matrix in "2D". The matrix has a shape of (64*64, 64*64)and values in {-1, 1}. .. GENERATED FROM PYTHON SOURCE LINES 49-56 .. code-block:: Python from spyrit.core.torch import walsh_matrix_2d H = walsh_matrix_2d(64) print(f"Acquisition matrix: {H.shape}", end=" ") print(rf"with values in {{{H.min()}, {H.max()}}}") .. rst-class:: sphx-glr-script-out .. code-block:: none Acquisition matrix: torch.Size([4096, 4096]) with values in {-1.0, 1.0} .. GENERATED FROM PYTHON SOURCE LINES 57-58 We instantiate a :class:`spyrit.core.meas.Linear` operator. To indicate that the operator works in 2D, on images with shape (64, 64), we use the :attr:`meas_shape` argument. .. GENERATED FROM PYTHON SOURCE LINES 58-62 .. code-block:: Python from spyrit.core.meas import Linear meas_op = Linear(H, (64, 64)) .. GENERATED FROM PYTHON SOURCE LINES 63-64 We simulate the measurement vectors, which have a shape of (7, 1, 4096). .. GENERATED FROM PYTHON SOURCE LINES 64-68 .. code-block:: Python y = meas_op(x) print(f"Measurement vectors: {y.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Measurement vectors: torch.Size([7, 1, 4096]) .. GENERATED FROM PYTHON SOURCE LINES 69-70 We now compute the pseudo inverse solutions, which have a shape of (7, 1, 64, 64). .. GENERATED FROM PYTHON SOURCE LINES 70-77 .. code-block:: Python from spyrit.core.inverse import PseudoInverse pinv = PseudoInverse(meas_op) x_rec = pinv(y) print(f"Reconstructed images: {x_rec.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Reconstructed images: torch.Size([7, 1, 64, 64]) .. GENERATED FROM PYTHON SOURCE LINES 78-79 We plot the reconstruction of the second image in the batch .. GENERATED FROM PYTHON SOURCE LINES 79-84 .. code-block:: Python from spyrit.misc.disp import imagesc, add_colorbar imagesc(x_rec[1, 0]) # sphinx_gallery_thumbnail_number = 1 .. image-sg:: /gallery/images/sphx_glr_tuto_03_pseudoinverse_linear_001.png :alt: tuto 03 pseudoinverse linear :srcset: /gallery/images/sphx_glr_tuto_03_pseudoinverse_linear_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.0.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) .. GENERATED FROM PYTHON SOURCE LINES 85-87 .. note:: The measurement operator is chosen as a Hadamard matrix with positive but this matrix can be replaced by any other matrix. .. GENERATED FROM PYTHON SOURCE LINES 89-91 LinearSplit measurements with Gaussian noise ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 93-94 We consider a linear operator where the positive and negative components are split, i.e. acquired separately. To do so, we instantiate a :class:`spyrit.core.meas.LinearSplit` operator. .. GENERATED FROM PYTHON SOURCE LINES 94-98 .. code-block:: Python from spyrit.core.meas import LinearSplit meas_op = LinearSplit(H, (64, 64)) .. GENERATED FROM PYTHON SOURCE LINES 99-100 We consider additive Gaussian noise with standard deviation 2. .. GENERATED FROM PYTHON SOURCE LINES 100-104 .. code-block:: Python from spyrit.core.noise import Gaussian meas_op.noise_model = Gaussian(2) .. GENERATED FROM PYTHON SOURCE LINES 105-106 We simulate the measurement vectors, which have shape (7, 1, 8192). .. GENERATED FROM PYTHON SOURCE LINES 106-110 .. code-block:: Python y = meas_op(x) print(f"Measurement vectors: {y.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Measurement vectors: torch.Size([7, 1, 8192]) .. GENERATED FROM PYTHON SOURCE LINES 111-112 We preprocess measurement vectors by computing the difference of the positive and negative components of the measurement vectors. To do so, we use the :class:`spyrit.core.prep.Unsplit` class. The preprocess measurements have a shape of (7, 1, 4096). .. GENERATED FROM PYTHON SOURCE LINES 112-120 .. code-block:: Python from spyrit.core.prep import Unsplit prep = Unsplit() m = prep(y) print(f"Preprocessed measurement vectors: {m.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Preprocessed measurement vectors: torch.Size([7, 1, 4096]) .. GENERATED FROM PYTHON SOURCE LINES 121-122 We now compute the pseudo inverse solutions, which have a shape of (7, 1, 64, 64). .. GENERATED FROM PYTHON SOURCE LINES 122-129 .. code-block:: Python from spyrit.core.inverse import PseudoInverse pinv = PseudoInverse(meas_op) x_rec = pinv(m) print(f"Reconstructed images: {x_rec.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Reconstructed images: torch.Size([7, 1, 64, 64]) .. GENERATED FROM PYTHON SOURCE LINES 130-131 We plot the reconstruction .. GENERATED FROM PYTHON SOURCE LINES 131-135 .. code-block:: Python from spyrit.misc.disp import imagesc, add_colorbar imagesc(x_rec[1, 0]) .. image-sg:: /gallery/images/sphx_glr_tuto_03_pseudoinverse_linear_002.png :alt: tuto 03 pseudoinverse linear :srcset: /gallery/images/sphx_glr_tuto_03_pseudoinverse_linear_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/spyrit/envs/3.0.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) .. GENERATED FROM PYTHON SOURCE LINES 136-138 HadamSplit2d with x4 subsampling with Poisson noise ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 140-141 We consider the acquisition of the 2D Hadamard transform of an image, where the positive and negative components of acquisition matrix are acquired separately. To do so, we use the dedicated :class:`spyrit.core.meas.HadamSplit2d` operator. It also allows for subsampling the rows the Hadamard matrix, using a sampling map. .. GENERATED FROM PYTHON SOURCE LINES 141-152 .. code-block:: Python from spyrit.core.meas import HadamSplit2d # Sampling map with ones in the top left corner and zeros elsewhere (low-frequency subsampling) sampling_map = torch.ones((64, 64)) sampling_map[:, 64 // 2 :] = 0 sampling_map[64 // 2 :, :] = 0 # Linear operator with HadamSplit2d meas_op = HadamSplit2d(64, 64**2 // 4, order=sampling_map, reshape_output=True) .. GENERATED FROM PYTHON SOURCE LINES 153-154 We consider additive Poisson noise with an intensity of 100 photons. .. GENERATED FROM PYTHON SOURCE LINES 154-159 .. code-block:: Python from spyrit.core.noise import Poisson meas_op.noise_model = Poisson(100) .. GENERATED FROM PYTHON SOURCE LINES 160-161 We simulate the measurement vectors, which have a shape of (7, 1, 2048) .. GENERATED FROM PYTHON SOURCE LINES 163-165 .. note:: The :class:`spyrit.core.noise.Poisson` class noise assumes that the images are in the range [0, 1] .. GENERATED FROM PYTHON SOURCE LINES 165-170 .. code-block:: Python y = meas_op(x) print(rf"Reference images with values in {{{x.min()}, {x.max()}}}") print(f"Measurement vectors: {y.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Reference images with values in {0.003921568859368563, 0.9960784316062927} Measurement vectors: torch.Size([7, 1, 2048]) .. GENERATED FROM PYTHON SOURCE LINES 171-172 We preprocess measurement vectors by i) computing the difference of the positive and negative components, and ii) normalizing the intensity. To do so, we use the :class:`spyrit.core.prep.UnsplitRescale` class. The preprocessed measurements have a shape of (7, 1, 1024). .. GENERATED FROM PYTHON SOURCE LINES 172-180 .. code-block:: Python from spyrit.core.prep import UnsplitRescale prep = UnsplitRescale(100) m = prep(y) # (y+ - y-)/alpha print(f"Preprocessed measurement vectors: {m.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Preprocessed measurement vectors: torch.Size([7, 1, 1024]) .. GENERATED FROM PYTHON SOURCE LINES 181-182 We compute the pseudo inverse solution, which has a shape of (7, 1, 64, 64). .. GENERATED FROM PYTHON SOURCE LINES 182-187 .. code-block:: Python x_rec = meas_op.fast_pinv(m) print(f"Reconstructed images: {x_rec.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none Reconstructed images: torch.Size([7, 1, 64, 64]) .. GENERATED FROM PYTHON SOURCE LINES 188-190 .. note:: There is no need to use the :class:`spyrit.core.inverse.PseudoInverse` class here, as the :class:`spyrit.core.meas.HadamSplit2d` class includes a method that returns the pseudo inverse solution. .. GENERATED FROM PYTHON SOURCE LINES 192-193 We plot the reconstruction .. GENERATED FROM PYTHON SOURCE LINES 193-196 .. code-block:: Python from spyrit.misc.disp import imagesc, add_colorbar imagesc(x_rec[1, 0]) .. image-sg:: /gallery/images/sphx_glr_tuto_03_pseudoinverse_linear_003.png :alt: tuto 03 pseudoinverse linear :srcset: /gallery/images/sphx_glr_tuto_03_pseudoinverse_linear_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/3.0.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) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 18.806 seconds) .. _sphx_glr_download_gallery_tuto_03_pseudoinverse_linear.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tuto_03_pseudoinverse_linear.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tuto_03_pseudoinverse_linear.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: tuto_03_pseudoinverse_linear.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_