Compressed Sensing of Images¶

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This tour explores compressed sensing of natural images, using different sparsity priors over a wavelet basis.

In [2]:

addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/sparsity_2_cs_images')

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.

In [3]:

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.

In [4]:

k0 = 2;
J = log2(n)-k0;
Wav  = @(f)perform_wavelet_transf(f,J,+1);
WavI = @(x)perform_wavelet_transf(x,J,-1);

Compute the wavelet transform.

In [5]:

fw = Wav(f);

Display the coefficients.

In [6]:

clf;
plot_wavelet(fw, J);

Exercise 1

Compute an approximation |fLow| using the $ P=2^{2J}=(n/k_0)^2 $ low pass coefficients.

In [7]:

exo1()

In [8]:

%% Insert your code here.

Randomized Orthogonal Measurements¶

We consider a compressed sensing operator that corresponds to randomized orthogonal projections.

Extract the high pass wavelet coefficients, $x_0 = \{ \dotp{f}{\psi_m} \}_{m \in I_0}$.

In [9]:

A = ones(n,n); A(1:2^J,1:2^J) = 0;
I0 = find(A==1);
x0 = fw(I0);

Number of coefficients.

In [10]:

N = length(x0);

Number $ P_0 = 2^{2J}=(n/k_0)^2 $ of low pass measurements.

In [11]:

P0 = (n/2^k0)^2;

Number of CS measurements.

In [12]:

P = 4 * P0;

Generate random permutation operators $S_1,S_2 : \RR^N \rightarrow \RR^N$ so that $S_k(x)_i = x_{\sigma_k(i)}$ where $ \sigma_k \in \Sigma_N $ is a random permutation of $\{1,\ldots,N\}$.

In [13]:

sigma1 = randperm(N)'; 
sigma2 = randperm(N)'; 
S1 = @(x)x(sigma1);
S2 = @(x)x(sigma2);

The adjoint (and also inverse) operators $S_1^*,S_2^*$ (denoted |S1S,S2S|) corresponds to the inverse permutation $\sigma_k^*$ such that $\sigma_k^* \circ \sigma_k(i)=i$.

In [14]:

sigma1S = 1:N; sigma1S(sigma1) = 1:N; 
sigma2S = 1:N; sigma2S(sigma2) = 1:N; 
S1S = @(x)x(sigma1S);
S2S = @(x)x(sigma2S);

We consider a CS operator $ \Phi : \RR^N \rightarrow \RR^P $ that corresponds to a projection on randomized atoms $$ (\Phi x)_i = \dotp{x}{ \phi_{\sigma_2(i)}} $$ where $ \phi_i $ is a scrambled orthogonal basis $$ \phi_i(x) = c_i( \sigma_1(x) ) $$ where $\{ c_i \}_i$ is the orthogonal DCT basis.

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.

In [15]:

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.

In [16]:

uparrow = @(x)[x; zeros(N-P,1)];
PhiS = @(x)S1S(idct(S2S(uparrow(x))));

Perform the CS (noiseless) measurements.

In [17]:

y = Phi(x0);

Exercise 2

Reconstruct an image using the pseudo inverse coefficients $\Phi^+ y = \Phi^* y$.

In [18]:

exo2()

In [19]:

%% Insert your code here.

Compressed Sensing Recovery using Douglas Rachford Scheme¶

We consider the minimum $\ell^1$ recovery from the measurements $y = \Phi x_0 \in \RR^P$ $$ \umin{\Phi x = y} \normu{x}. $$ This can be written as $$ \umin{ x } F(x) + G(x) \qwhereq \choice{ F(x) = i_{\Cc}(x), \\ G(x) = \normu{x}. }$ $$ where $\Cc = \enscond{x}{\Phi x =y}$.

One can solve this problem using the Douglas-Rachford iterations $$ \tilde x_{k+1} = \pa{1-\frac{\mu}{2}} \tilde x_k + \frac{\mu}{2} \text{rPox}_{\gamma G}( \text{rProx}_{\gamma F}(\tilde x_k) ) \qandq x_{k+1} = \text{Prox}_{\gamma F}(\tilde x_{k+1},) $$

We have use the following definition for the proximal and reversed-proximal mappings: $$ \text{rProx}_{\gamma F}(x) = 2\text{Prox}_{\gamma F}(x)-x $$ $$ \text{Prox}_{\gamma F}(x) = \uargmin{y} \frac{1}{2}\norm{x-y}^2 + \ga F(y). $$

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$.

Exercise 3

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.

In [20]:

exo3()

In [21]:

%% Insert your code here.

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

In [22]:

mu = 1;
gamma = 1;

Exercise 4

Implement the DR iterative algorithm. Keep track of the evolution of the $\ell^1$ norm $G(x_k)$.

In [23]:

exo4()

In [24]:

%% Insert your code here.

Exercise 5

Display the image reconstructed using the $P_0$ linear and $P$ CS measurements. The total number of used measurements is thus $P+P_0$.

In [25]:

exo5()