Installing toolboxes and setting up the path.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');

Low Pass Linear Measures

We first make use of \(P\) low pass linear measurements to remove the low frequency content of the image.

Natural images are not only sparse over a wavelet domain. They also exhibit a fast decay of the coefficient through the scale.
The coarse (low pass) wavelets caries much of the image energy. It thus make sense to measure directly the low pass coefficients.

We load an image \(f \in \RR^{n^2}\) of \(n \times n\) pixels.

name = 'boat';
n = 256;
f = load_image(name, n);
f = rescale(f);

Shortcuts for the wavelet transform \( \{\dotp{f}{\psi_m}\}_m \). We only compute up to a scale \(J\) so that only \(k_0\)
sub-bands are transformed.

This can be rewritten in compact operator form as \[ \Phi x = ( S_2 \circ C \circ S_1 (x) ) \downarrow_P \] where \(S_1,S_2\)
are the permutation operators, and \(\downarrow_P\) selects the \(P\) first entries of a vector.

downarrow = @(x)x(1:P);
Phi = @(x)downarrow(S2(dct(S1(x))));

The adjoint operator is \[ \Phi^* x = S_1^* \circ C^* \circ S_2^* (x\uparrow_P) \] where \(\uparrow_P\) append \(N-P\) zeros
at the end of a vector, and \(C^*\) is the inverse DCT transform.

One can show that for any value of \(\gamma>0\), any \( 0 < \mu < 2 \), and any \(\tilde x_0\), \(x_k \rightarrow x^\star\)
which is a solution of the minimization of \(F+G\).

Exercice 3: (check the solution) Implement the proximal and reversed-proximal mappings of \(F\) (the orthogonal projector on \(\Cc\) and \(G\) (soft thresholding).
In Matlab, use inline function with the @ operator.

exo3;

Value for the \(0 < \mu < 2\) and \(\gamma>0\) parameters. You can use other values, this might speed up the convergence.

mu = 1;
gamma = 1;

Exercice 4: (check the solution) Implement the DR iterative algorithm. Keep track of the evolution of the \(\ell^1\) norm \(G(x_k)\).

exo4;

Exercice 5: (check the solution) Display the image reconstructed using the \(P_0\) linear and \(P\) CS measurements. The total number of used measurements
is thus \(P+P_0\).

exo5;

Compressed Sensing Reconstruction using Block Sparsity

In order to enhance the CS reconstruction, it is possible to use more advanced priors than plain \( \ell^1 \).

One can for instance use a block \( \ell^1 \) norm \[ G(x) = \sum_i \norm{x_{B_i}} \] where \( (B_i)_i \) is a disjoint segmentation
of the index set \( \{1,\ldots,N\} \), where \(x_{B} = \{ x_i \}_{i \in B} \in \RR^{|B|}\) extracts the coefficients within
\(B\), and \( \norm{x_B} \) is the \(\ell^2\) norm.