Image Deconvolution using Variational Method¶

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 Sobolev and TV regularization to perform image deconvolution.

In [2]:

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

This tour is concerned with the deconvolution problem. The measurement are assumed to be blurry and noisy: $$y=\Phi f_0 + w = h \star f_0 + w$$

Where here |h| is the filter (low pass) and |w| some noise (here assumed to be white Gaussian).

We consider variational deconvolution methods, that finds a regularizer through a convex optimization: $$f^\star \in \text{argmin}_f \frac{1}{2}\|y-\Phi f\|^2 + \lambda J(f)$$

Where $J(f)$ is a prior energy. In this tour we consider a simple L2 prior (the image is assumed to have a bounded energy), a Sobolev prior (the image is uniformly smooth) and an approximate total variation (the image has edges of bounded perimeter).

Note that the parameter $\lambda$ should be carefully chosen to fit the noise level.

Image Blurring¶

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

First we load the image to be inpainted.

In [3]:

n = 256;
name = 'lena';
name = 'mri';
name = 'boat';
f0 = load_image(name);
f0 = rescale(crop(f0,n));

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.

Width of the kernel, in pixel.

In [4]:

s = 3;

Kernel.

In [5]:

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 [6]:

hF = real(fft2(h));

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

In [7]:

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

We use this short hand for the filtering. Note that this is a symmetric operator.

In [8]:

if using_matlab()
    Phi = @(x,h)real(ifft2(fft2(x).*fft2(h)));
end

Important Scilab user should define a function |Phi| in a separate file |Phi.sci| to perform this.

Apply the filter.

In [9]:

y0 = Phi(f0,h);

Display the filtered observation.

In [10]:

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

Variance $\sigma^2$ of the noise $w$.

In [11]:

sigma = .02;

Add some noise to obtain the measurements $y = \Phi f_0 + w$.

In [12]:

y = y0 + randn(n)*sigma;

Display.

In [13]:

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

Deconvolution with L2 Regularization¶

Deconvolution is obtained by dividing the Fourier transform of $y$ by $\hat h$. $$f^\star(\omega) = \frac{\hat y(\omega)}{\hat h(\omega)} = \hat f_0(\omega) + \hat w(\omega)/{\hat h(\omega)}$$

Unfortunately this creates an explosion of the Noise.

To avoid this explosion, we consider a simple regularization. $$f^{\star} = \text{argmin}_f \: \|y-\Phi f\|^2 + \lambda \|f\|^2$$

Since the filtering is diagonalized over Fourier, the solution is simply computed over the Fourier domain as: $$\hat f^\star(\omega) = \frac{\hat y(\omega) \hat h(\omega)}{ \|\hat h(\omega)\|^2 + \lambda }$$

Useful for later: Fourier transform of the observations.

In [14]:

yF = fft2(y);

Select a value for the regularization parameter.

In [15]:

lambda = 0.02;

Perform the inversion.

In [16]:

fL2 = real( ifft2( yF .* hF ./ ( abs(hF).^2 + lambda) ) );

Display.

In [17]:

clf;
imageplot(y, strcat(['Observation, SNR=' num2str(snr(f0,y),3) 'dB']), 1,2,1);
imageplot(clamp(fL2), strcat(['L2 deconvolution, SNR=' num2str(snr(f0,fL2),3) 'dB']), 1,2,2);

Exercise 1

Find the optimal solution |fL2| by testing several value of |lambda|.

In [18]:

exo1()

In [19]:

%% Insert your code here.

Display optimal result.

In [20]:

clf;
imageplot(y, strcat(['Observation, SNR=' num2str(snr(f0,y),3) 'dB']), 1,2,1);
imageplot(clamp(fL2), strcat(['L2 deconvolution, SNR=' num2str(snr(f0,fL2),3) 'dB']), 1,2,2);