Note

Go to the endto download the full example code.

This tutorial is an extension of the Tutorials 01 and 02 where:

we introduce split measurements to handle a Hadamard measurements,

we discuss subsampling for accelerated acquisitions.

## Load a batch of images

Images \(x\) for training neural networks expect values in [-1,1]. The images are normalizedusing the `transform_gray_norm()`

function.

import osimport torchimport torchvisionimport numpy as npimport matplotlib.pyplot as pltfrom spyrit.misc.disp import imagescfrom spyrit.misc.statistics import transform_gray_normh = 64 # image size hxhi = 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 tensorstransform = 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 imagex = x[i : i + 1, :, :, :]x = x.detach().clone()b, c, h, w = x.shape# plotx_plot = x.view(-1, h, h).cpu().numpy()imagesc(x_plot[0, :, :], r"$x$ in [-1, 1]")

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

## The measurement and noise operators

Noise operators are defined in the noise module. A noiseoperator computes the following three steps sequentially:

Normalization of the image \(x\) with values in [-1,1] to get animage \(\tilde{x}=\frac{x+1}{2}\) in [0,1], as it is required for measurement simulation

Application of the measurement model, i.e., computation of \(P\tilde{x}\)

Application of the noise model

\[y \sim \texttt{Noise}(P\tilde{x}) = \texttt{Noise}\left(\frac{P(x+1)}{2}\right).\]

The normalization is useful when considering distributions suchas the Poisson distribution that are defined on positive values.

## Split measurement operator and no noise

\[\begin{split}y = P\tilde{x}= \begin{bmatrix} H_{+} \\ H_{-} \end{bmatrix} \tilde{x}.\end{split}\]

Hadamard split measurement operator is defined in the spyrit.core.meas.HadamSplit class.It computes linear measurements from incoming images, where \(P\) is alinear operator (matrix) with positive entries and \(\tilde{x}\) is a vectorized image.The class relies on a matrix \(H\) withshape \((M,N)\) where \(N\) represents the number of pixels in theimage and \(M \le N\) the number of measurements. The matrix \(P\)is obtained by splitting the matrix \(H\) as \(H = H_{+}-H_{-}\) where\(H_{+} = \max(0,H)\) and \(H_{-} = \max(0,-H)\).

## Subsampling

Subsampling is done by retaining the first \(M\) rows of a permuted Hadamard matrix \(H=GF\), where \(G\) is a subsampled permutation matrix with shape \((M,N)\) and \(F\) is a “full” Hadamard matrix with shape \((N,N)\).

We consider two subsampling strategies:

“Naive subsampling” by retaining only the first \(M\) rows of the measurement matrix.

“Variance subsampling” by retaining only the first \(M\) rows of a permuted measurement matrix where the first rows corresponds to the coefficients with largest variance and the last ones to the coefficients that are close to constant. The motivation is that almost constant coefficients are less informative than the others. This can be supported by principal component analysis, which states that preserving the components with largest variance leads to the best linear predictor.

First, we download the covariance matrix from our warehouse and load it. The covariance matrixhas been computed from ImageNet 2012 dataset.

from spyrit.misc.load_data import download_girder# api Rest url of the warehouseurl = "https://pilot-warehouse.creatis.insa-lyon.fr/api/v1"dataId = "63935b624d15dd536f0484a5" # for reconstruction (imageNet, 64)data_folder = "./stat/"cov_name = "Cov_64x64.npy"# download the covariance matrix and get the file pathfile_abs_path = download_girder(url, dataId, data_folder, cov_name)try: # Load covariance matrix for "variance subsampling" Cov = np.load(file_abs_path) print(f"Cov matrix {cov_name} loaded")except: # Set to the identity if not found for "naive subsampling" Cov = np.eye(h * h) print(f"Cov matrix {cov_name} not found! Set to the identity")

Local folder not found, creating it... done.Downloading file... done.Cov matrix Cov_64x64.npy loaded

The permutation matrix is defined from a sampling matrix with shape \((\sqrt{N},\sqrt{N})\) (see the sampling submodule).

We compute the sampling matrix for the “naive” subsampling

from spyrit.misc.statistics import Cov2Varfrom spyrit.misc.disp import add_colorbar, noaxisM = 64 * 64 // 4 # number of measurements (here, 1/4 of the pixels)Cov_eye = np.eye(h * h)Ord_nai = Cov2Var(Cov_eye)

And for the “variance” subsampling

Ord_var = Cov2Var(Cov)

Further insight on the two strategies can be gained by plotting the masks corresponding to the sampling matrices.

# sphinx_gallery_thumbnail_number = 2from spyrit.misc.sampling import sort_by_significancemask_basis = np.zeros((h, h))mask_basis.flat[:M] = 1# Mask for "naive subsampling"mask_nai = sort_by_significance(mask_basis, Ord_nai, axis="flatten")# Mask for "variance subsampling"mask_var = sort_by_significance(mask_basis, Ord_var, axis="flatten")# Plot the masksf, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))im1 = ax1.imshow(mask_nai, vmin=0, vmax=1)ax1.set_title("Mask \n'naive subsampling'", fontsize=20)noaxis(ax1)add_colorbar(im1, "bottom", size="20%")im2 = ax2.imshow(mask_var, vmin=0, vmax=1)ax2.set_title("Mask \n'variance subsampling'", fontsize=20)noaxis(ax2)add_colorbar(im2, "bottom", size="20%")

<matplotlib.colorbar.Colorbar object at 0x7f9120c56150>

Note

Note that in this tutorial the covariance matrix is used only for choosing the subsampling strategy. Although the covariance matrix can be also exploited to improve the reconstruction, this will be considered in a future tutorial.

## Measurement and noise operators

We compute the measurement and noise operators and thensimulate the measurement vector \(y\).

We consider Poisson noise, i.e., a noisy measurement vector given by

\[y \sim \mathcal{P}(\alpha P \tilde{x}),\]

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

We use the spyrit.core.noise.Poisson class, set \(\alpha\)to 100 photons, and simulate a noisy measurement vector for the two sampling strategies. Subsampling is handled internally by the HadamSplit class.

from spyrit.core.noise import Poissonfrom spyrit.core.meas import HadamSplitfrom spyrit.core.noise import Poissonalpha = 100.0 # number of photons# "Naive subsampling"# Measurement and noise operatorsmeas_nai_op = HadamSplit(M, h, torch.from_numpy(Ord_nai))noise_nai_op = Poisson(meas_nai_op, alpha)# Measurement operatorx = x.view(b * c, h * w) # vectorized imagey_nai = noise_nai_op(x) # a noisy measurement vector# "Variance subsampling"meas_var_op = HadamSplit(M, h, torch.from_numpy(Ord_var))noise_var_op = Poisson(meas_var_op, alpha)y_var = noise_var_op(x) # a noisy measurement vectorx = x.view(b * c, h * w) # vectorized imageprint(f"Shape of vectorized image: {x.shape}")print(f"Shape of simulated measurements y: {y_var.shape}")

Shape of vectorized image: torch.Size([1, 4096])Shape of simulated measurements y: torch.Size([1, 2048])

## The preprocessing operator measurements for split measurements

We compute the preprocessing operators for the three cases considered above,using the spyrit.core.prep module. As previously introduced,a preprocessing operator applies to the noisy measurements in order tocompensate for the scaling factors that appear in the measurement or noise operators:

\[m = \texttt{Prep}(y),\]

We consider the spyrit.core.prep.SplitPoisson class that intendsto “undo” the spyrit.core.noise.Poisson class, for split measurements, 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

\[m = \frac{2(y_+-y_-)}{\alpha} - P\mathbb{1},\]

where \(y_+=H_+\tilde{x}\) and \(y_-=H_-\tilde{x}\).This is handled internally by the spyrit.core.prep.SplitPoisson class.

We compute the preprocessing operator and the measurements vectors forthe two sampling strategies.

from spyrit.core.prep import SplitPoisson# "Naive subsampling"## Preprocessing operatorprep_nai_op = SplitPoisson(alpha, meas_nai_op)# Preprocessed measurementsm_nai = prep_nai_op(y_nai)# "Variance subsampling"prep_var_op = SplitPoisson(alpha, meas_var_op)m_var = prep_var_op(y_var)

## Noiseless measurements

We consider now noiseless measurements for the “naive subsampling” strategy.We compute the required operators and the noiseless measurement vector.For this we use the spyrit.core.noise.NoNoise class, which normalizesthe input vector to get an image in [0,1], as explained inacquisition operators tutorial.For the preprocessing operator, we assign the number of photons equal to one.

from spyrit.core.noise import NoNoisenonoise_nai_op = NoNoise(meas_nai_op)y_nai_nonoise = nonoise_nai_op(x) # a noisy measurement vectorprep_nonoise_op = SplitPoisson(1.0, meas_nai_op)m_nai_nonoise = prep_nonoise_op(y_nai_nonoise)

We can now plot the three measurement vectors

from spyrit.misc.sampling import meas2img# Plot the three measurement vectorsm_plot = m_nai_nonoise.numpy()m_plot = meas2img(m_plot, Ord_nai)m_plot_max = np.max(m_plot[0, :, :])m_plot_min = np.min(m_plot[0, :, :])m_plot2 = m_nai.numpy()m_plot2 = meas2img(m_plot2, Ord_nai)m_plot3 = m_var.numpy()m_plot3 = meas2img(m_plot3, Ord_var)f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 7))im1 = ax1.imshow(m_plot[0, :, :], cmap="gray")ax1.set_title("Noiseless measurements $m$ \n 'Naive' subsampling", fontsize=20)noaxis(ax1)add_colorbar(im1, "bottom", size="20%")im2 = ax2.imshow(m_plot2[0, :, :], cmap="gray", vmin=m_plot_min, vmax=m_plot_max)ax2.set_title("Measurements $m$ \n 'Naive' subsampling", fontsize=20)noaxis(ax2)add_colorbar(im2, "bottom", size="20%")im3 = ax3.imshow(m_plot3[0, :, :], cmap="gray", vmin=m_plot_min, vmax=m_plot_max)ax3.set_title("Measurements $m$ \n 'Variance' subsampling", fontsize=20)noaxis(ax3)add_colorbar(im3, "bottom", size="20%")

<matplotlib.colorbar.Colorbar object at 0x7f9120b468d0>

## PinvNet network

We use the spyrit.core.recon.PinvNet class wherethe pseudo inverse reconstruction is performed by a neural network

from spyrit.core.recon import PinvNet# PinvNet(meas_op, prep_op, denoi=torch.nn.Identity())pinvnet_nai_nonoise = PinvNet(nonoise_nai_op, prep_nonoise_op)pinvnet_nai = PinvNet(noise_nai_op, prep_nai_op)pinvnet_var = PinvNet(noise_var_op, prep_var_op)# Reconstructionz_nai_nonoise = pinvnet_nai_nonoise.reconstruct(y_nai_nonoise)z_nai = pinvnet_nai.reconstruct(y_nai)z_var = pinvnet_var.reconstruct(y_var)

We can now plot the three reconstructed images

from spyrit.misc.disp import add_colorbar, noaxis# Plotx_plot = x.view(-1, h, h).numpy()z_plot_nai_nonoise = z_nai_nonoise.view(-1, h, h).numpy()z_plot_nai = z_nai.view(-1, h, h).numpy()z_plot_var = z_var.view(-1, h, h).numpy()f, axs = plt.subplots(2, 2, figsize=(10, 10))im1 = axs[0, 0].imshow(x_plot[0, :, :], cmap="gray")axs[0, 0].set_title("Ground-truth image")noaxis(axs[0, 0])add_colorbar(im1, "bottom")im2 = axs[0, 1].imshow(z_plot_nai_nonoise[0, :, :], cmap="gray")axs[0, 1].set_title("Reconstruction noiseless")noaxis(axs[0, 1])add_colorbar(im2, "bottom")im3 = axs[1, 0].imshow(z_plot_nai[0, :, :], cmap="gray")axs[1, 0].set_title("Reconstruction \n 'Naive' subsampling")noaxis(axs[1, 0])add_colorbar(im3, "bottom")im4 = axs[1, 1].imshow(z_plot_var[0, :, :], cmap="gray")axs[1, 1].set_title("Reconstruction \n 'Variance' subsampling")noaxis(axs[1, 1])add_colorbar(im4, "bottom")plt.show()

Note

Note that reconstructed images are pixelized when using the “naive subsampling”,while they are smoother and more similar to the ground-truth image when using the“variance subsampling”.

Another way to further improve results is to include a nonlinear post-processing step,which we will consider in a future tutorial.

**Total running time of the script:** (2 minutes 55.224 seconds)

Download Jupyter notebook: tuto_05_acquisition_split_measurements.ipynb

Download Python source code: tuto_05_acquisition_split_measurements.py