Image Deconvolution using Sparse Regularization¶

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 numerical tour explores the use of sparsity regularization to perform image deconvolution.

In [2]:

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

Sparse Regularization¶

This tour consider measurements $y=\Phi f_0 + w$ where $\Phi$ is a convolution $ \Phi f = h \star f $ and $w$ is an additive noise.

This tour is focussed on using sparsity to recover an image from the measurements $y$. It consider a synthesis-based regularization, that compute a sparse set of coefficients $ (a_m^{\star})_m $ in a frame $\Psi = (\psi_m)_m$ that solves $$a^{\star} \in \text{argmin}_a \: \frac{1}{2}\|y-\Phi \Psi a\|^2 + \lambda J(a)$$

where $\lambda$ should be adapted to the noise level $\|w\|$

Here we used the notation $$\Psi a = \sum_m a_m \psi_m$$ to indicate the reconstruction operator, and $J(a)$ is the $\ell^1$ sparsity prior $$J(a)=\sum_m \|a_m\|.$$

Image Blurring¶

Deconvolution corresponds to removing a blur from an image. We use here a Gaussian blur.

Parameters for the tour: width of the kernel (in pixel) and noise level.

In [3]:

setting = 1;
switch setting
    case 1
        % difficult
        s = 3;
        sigma = .02;
    case 2
        % easy
        s = 1.2;
        sigma = .02;
end

First we load the image to be processed.

In [4]:

n = 128*2;
name = 'lena';
name = 'boat';
name = 'mri';
f0 = load_image(name);
f0 = rescale(crop(f0,n));

Display it.

In [5]:

clf;
imageplot(f0);

We build a convolution kernel. Since we are going to use Fourier to compute the convolution, we set the center of the kernel in the (1,1) pixel location.

Kernel.

In [6]:

x = [0:n/2-1, -n/2:-1];
[Y,X] = meshgrid(x,x);
h = exp( (-X.^2-Y.^2)/(2*s^2) );
h = h/sum(h(:));

Useful for later : the Fourier transform (should be real because of symmetry).

In [7]:

hF = real(fft2(h));

Display the kernel $h$ and its transform $\hat h$. We use |fftshift| to center the filter for display.

In [8]:

clf;
imageplot(fftshift(h), 'Filter', 1,2,1);
imageplot(fftshift(hF), 'Fourier transform', 1,2,2);

We use this short hand for the filtering. Scilab user should define a function in a separate file to perform this. Note that this is a symmetric operator.

In [9]:

Phi = @(x)real(ifft2(fft2(x).*hF));

Apply the filter.

In [10]:

y0 = Phi(f0);

Display the filtered observation.

In [11]:

clf;
imageplot(f0, 'Image f0', 1,2,1);
imageplot(y0, 'Observation without noise', 1,2,2);

Add some noise of variance $\sigma^2$, to obtain $y=\Phi f_0 + w = f_0 \star h + w$.

In [12]:

y = y0 + randn(n,n)*sigma;

Display.

In [13]:

clf;
imageplot(y0, 'Observation without noise', 1,2,1);
imageplot(clamp(y), 'Observation with noise', 1,2,2);

Soft Thresholding in a Basis¶

The soft thresholding operator is at the heart of $\ell^1$ minimization schemes. It can be applied to coefficients $a$, or to an image $f$ in an ortho-basis.

The soft thresholding is a 1-D functional that shrinks the value of coefficients. $$ s_T(u)=\max(0,1-T/|u|)u $$

Define a shortcut for this soft thresholding 1-D functional.

In [14]:

SoftThresh = @(x,T)x.*max( 0, 1-T./max(abs(x),1e-10) );

Display a curve of the 1D soft thresholding.

In [15]:

clf;
T = linspace(-1,1,1000);
plot( T, SoftThresh(T,.5) );
axis('equal');

Note that the function |SoftThresh| can also be applied to vector (because of Matlab/Scilab vectorialized computation), which defines an operator on coefficients: $$ S_T(a) = ( s_T(a_m) )_m. $$

In the next section, we use an orthogonal wavelet basis $\Psi$.

We set the parameters of the wavelet transform.

In [16]:

Jmax = log2(n)-1;
Jmin = Jmax-3;

Shortcut for $\Psi$ and $\Psi^*$ in the orthogonal case.

In [17]:

options.ti = 0; % use orthogonality.
Psi = @(a)perform_wavelet_transf(a, Jmin, -1,options);
PsiS = @(f)perform_wavelet_transf(f, Jmin, +1,options);

The soft thresholding opterator in the basis $\Psi$ is defined as $$S_T^\Psi(f) = \sum_m s_T( \langle f,\psi_m \rangle ) \psi_m $$

It thus corresponds to applying the transform $\Psi^*$, thresholding the coefficients using $S_T$ and then undoing the transform using $\Psi$. $$ S_T^\Psi(f) = \Psi \circ S_T \circ \Psi^*$$

In [18]:

SoftThreshPsi = @(f,T)Psi(SoftThresh(PsiS(f),T));

This soft thresholding corresponds to a denoising operator.

In [19]:

clf;
imageplot( clamp(SoftThreshPsi(f0,.1)) );

Deconvolution using Orthogonal Wavelet Sparsity¶

If $\Psi$ is an orthogonal basis, a change of variable shows that the synthesis prior is also an analysis prior, that reads $$f^{\star} \in \text{argmin}_f \: E(f) = \frac{1}{2}\|y-\Phi f\|^2 + \lambda \sum_m \|\langle f,\psi_m \rangle\|. $$

To solve this non-smooth optimization problem, one can use forward-backward splitting, also known as iterative soft thresholding.

It computes a series of images $f^{(\ell)}$ defined as $$ f^{(\ell+1)} = S_{\tau\lambda}^{\Psi}( f^{(\ell)} - \tau \Phi^{*}$ (\Phi f^{(\ell)} - y) ) $$

Set up the value of the threshold.

In [20]:

lambda = .02;

In our setting, since $h$ is symmetric, one has $\Phi^* f = \Phi f = f \star h$.

For $f^{(\ell)}$ to converge to a solution of the problem, the gradient step size should be chosen as $$\tau < \frac{2}{\|\Phi^* \Phi\|}$$

Since the filtering is an operator of norm 1, this must be smaller than 2.

In [21]:

tau = 1.5;

Number of iterations.

In [22]:

niter = 100;

Initialize the solution.

In [23]:

fSpars = y;

First step: perform one step of gradient descent of the energy $ \|y-f\star h\|^2 $.

In [24]:

fSpars = fSpars + tau * Phi( y-Phi(fSpars) );

Second step: denoise the solution by thresholding.

In [25]:

fSpars = SoftThreshPsi( fSpars, lambda*tau );

Exercise 1

Perform the iterative soft thresholding. Monitor the decay of the energy $E$ you are minimizing.

In [26]:

exo1()

In [27]:

%% Insert your code here.

Display the result.

In [28]:

clf;
imageplot(clamp(fSpars), ['Sparsity deconvolution, SNR=' num2str(snr(f0,fSpars),3) 'dB']);