.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tuto_01_b_splitting.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_b_splitting.py: 01.b. Acquisition operators (splitting) ==================================================== .. _tuto_acquisition_operators_splitting: This tutorial shows how to simulate linear measurements by splitting an acquisition matrix :math:`H\in \mathbb{R}^{M\times N}` that contains negative values. It based on the :class:`spyrit.core.meas.LinearSplit` class of the :mod:`spyrit.core.meas` submodule. .. image:: ../fig/tuto1.png :width: 600 :align: center :alt: Reconstruction architecture sketch | In practice, only positive values can be implemented using a digital micromirror device (DMD). Therefore, we acquire .. math:: y =Ax, where :math:`A \colon\, \mathbb{R}_+^{2M\times N}` is the acquisition matrix that contains positive DMD patterns, :math:`x \in \mathbb{R}^N` is the signal of interest, :math:`2M` is the number of DMD patterns, and :math:`N` is the dimension of the signal. .. important:: The vector :math:`x \in \mathbb{R}^N` represents a multi-dimensional array (e.g, an image :math:`X \in \mathbb{R}^{N_1 \times N_2}` with :math:`N = N_1 \times N_2`). Both variables are related through vectorization , i.e., :math:`x = \texttt{vec}(X)`. Given a matrix :math:`H` with negative values, we define the positive DMD patterns :math:`A` from the positive and negative components :math:`H`. In practice, the even rows of :math:`A` contain the positive components of :math:`H`, while odd rows of :math:`A` contain the negative components of :math:`H`. .. math:: \begin{cases} A[0::2, :] = H_{+}, \text{ with } H_{+} = \max(0,H),\\ A[1::2, :] = H_{-}, \text{ with } H_{-} = \max(0,-H). \end{cases} .. GENERATED FROM PYTHON SOURCE LINES 37-39 Splitting in 1D ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 41-42 We instantiate a measurement operator from a matrix of shape (10, 15). .. GENERATED FROM PYTHON SOURCE LINES 42-48 .. code-block:: Python import torch from spyrit.core.meas import LinearSplit H = torch.randn(10, 15) meas_op = LinearSplit(H) .. GENERATED FROM PYTHON SOURCE LINES 49-50 We consider 3 signals of length 15. .. GENERATED FROM PYTHON SOURCE LINES 50-52 .. code-block:: Python x = torch.randn(3, 15) .. GENERATED FROM PYTHON SOURCE LINES 53-55 We apply the operator to the batch of images, which produces 3 measurements of length 10*2 = 20. .. GENERATED FROM PYTHON SOURCE LINES 55-58 .. code-block:: Python y = meas_op(x) print(y.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none torch.Size([3, 20]) .. GENERATED FROM PYTHON SOURCE LINES 59-61 .. note:: The number of measurements is twice the number of rows of the matrix H that contains negative values. .. GENERATED FROM PYTHON SOURCE LINES 63-65 Illustration ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 67-68 We plot the positive and negative components of H that are concatenated in the matrix A. .. GENERATED FROM PYTHON SOURCE LINES 68-100 .. code-block:: Python A = meas_op.A H_pos = meas_op.A[::2, :] # Even rows H_neg = meas_op.A[1::2, :] # Odd rows from spyrit.misc.disp import add_colorbar, noaxis import matplotlib.pyplot as plt fig = plt.figure(figsize=(10, 5)) gs = fig.add_gridspec(2, 2) ax1 = fig.add_subplot(gs[:, 0]) ax2 = fig.add_subplot(gs[0, 1]) ax3 = fig.add_subplot(gs[1, 1]) ax1.set_title("Forward matrix A") im = ax1.imshow(A, cmap="gray") add_colorbar(im) ax2.set_title("Forward matrix H_pos") im = ax2.imshow(H_pos, cmap="gray") add_colorbar(im) ax3.set_title("Measurements H_neg") im = ax3.imshow(H_neg, cmap="gray") add_colorbar(im) noaxis(ax1) noaxis(ax2) noaxis(ax3) # sphinx_gallery_thumbnail_number = 1 .. image-sg:: /gallery/images/sphx_glr_tuto_01_b_splitting_001.png :alt: Forward matrix A, Forward matrix H_pos, Measurements H_neg :srcset: /gallery/images/sphx_glr_tuto_01_b_splitting_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.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/3.0.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/3.0.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 101-102 We can verify numerically that H = H_pos - H_neg .. GENERATED FROM PYTHON SOURCE LINES 102-108 .. code-block:: Python H = meas_op.H diff = torch.linalg.norm(H - (H_pos - H_neg)) print(f"|| H - (H_pos - H_neg) || = {diff}") .. rst-class:: sphx-glr-script-out .. code-block:: none || H - (H_pos - H_neg) || = 0.0 .. GENERATED FROM PYTHON SOURCE LINES 109-110 We now plot the matrix-vector products between A and x. .. GENERATED FROM PYTHON SOURCE LINES 110-126 .. code-block:: Python f, axs = plt.subplots(1, 3, figsize=(10, 5)) axs[0].set_title("Forward matrix A") im = axs[0].imshow(A, cmap="gray") add_colorbar(im, "bottom") axs[1].set_title("Signals x") im = axs[1].imshow(x.T, cmap="gray") add_colorbar(im, "bottom") axs[2].set_title("Split measurements y") im = axs[2].imshow(y.T, cmap="gray") add_colorbar(im, "bottom") noaxis(axs) .. image-sg:: /gallery/images/sphx_glr_tuto_01_b_splitting_002.png :alt: Forward matrix A, Signals x, Split measurements y :srcset: /gallery/images/sphx_glr_tuto_01_b_splitting_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.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/3.0.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/3.0.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 127-129 Simulations with noise and using the matrix H -------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 131-137 The operators in the :mod:`spyrit.core.meas` submodule allow for simulating noisy measurements .. math:: y =\mathcal{N}\left(Ax\right), where :math:`\mathcal{N} \colon\, \mathbb{R}^{2M} \to \mathbb{R}^{2M}` represents a noise operator (e.g., Gaussian). By default, no noise is applied to the measurement, i.e., :math:`\mathcal{N}` is the identity. We can consider noise by setting the :attr:`noise_model` attribute of the :class:`spyrit.core.meas.LinearSplit` class. .. GENERATED FROM PYTHON SOURCE LINES 139-140 For instance, we can consider additive Gaussian noise with standard deviation 2. .. GENERATED FROM PYTHON SOURCE LINES 140-145 .. code-block:: Python from spyrit.core.noise import Gaussian meas_op.noise_model = Gaussian(2) .. GENERATED FROM PYTHON SOURCE LINES 146-148 .. note:: To learn more about noise models, please refer to :ref:`tutorial 2 `. .. GENERATED FROM PYTHON SOURCE LINES 150-151 We simulate the noisy measurement vectors .. GENERATED FROM PYTHON SOURCE LINES 151-153 .. code-block:: Python y_noise = meas_op(x) .. GENERATED FROM PYTHON SOURCE LINES 154-155 Noiseless measurements can be simulated using the :meth:`spyrit.core.LinearSplit.measure` method. .. GENERATED FROM PYTHON SOURCE LINES 155-157 .. code-block:: Python y_nonoise = meas_op.measure(x) .. GENERATED FROM PYTHON SOURCE LINES 158-159 The :meth:`spyrit.core.LinearSplit.measure_H` method simulates noiseless measurements using the matrix H, i.e., :math:`m = Hx`. .. GENERATED FROM PYTHON SOURCE LINES 159-161 .. code-block:: Python m_nonoise = meas_op.measure_H(x) .. GENERATED FROM PYTHON SOURCE LINES 162-163 We now plot the noisy and noiseless measurements .. GENERATED FROM PYTHON SOURCE LINES 163-177 .. code-block:: Python f, axs = plt.subplots(1, 3, figsize=(8, 5)) axs[0].set_title("Split measurements y \n with noise") im = axs[0].imshow(y_noise.mT, cmap="gray") add_colorbar(im) axs[1].set_title("Split measurements y \n without noise") im = axs[1].imshow(y_nonoise.mT, cmap="gray") add_colorbar(im) axs[2].set_title("Measurements m \n without noise") im = axs[2].imshow(m_nonoise.mT, cmap="gray") add_colorbar(im) noaxis(axs) .. image-sg:: /gallery/images/sphx_glr_tuto_01_b_splitting_003.png :alt: Split measurements y with noise, Split measurements y without noise, Measurements m without noise :srcset: /gallery/images/sphx_glr_tuto_01_b_splitting_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.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/3.0.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/3.0.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 0.323 seconds) .. _sphx_glr_download_gallery_tuto_01_b_splitting.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_b_splitting.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tuto_01_b_splitting.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: tuto_01_b_splitting.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_