# Image Deconvolution using Sparse Regularization

This numerical tour explores the use of sparsity regularization to perform image deconvolution.

## Contents

- Installing toolboxes and setting up the path.
- Sparse Regularization
- Image Blurring
- Soft Thresholding in a Basis
- Deconvolution using Orthogonal Wavelet Sparsity
- Deconvolution using Translation Invariant Wavelet Sparsity

## Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox and general toolbox.

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/');

## Sparse Regularization

This tour consider measurements y=Φ f 0+w
where Φ
is a convolution Φ f=h⋆h
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)m
in a frame Ψ=(ψ m)m
that solves

a⋆∈argmin a 1 2∥y−Φ Ψ a∥2+λ J(a)

where λ
should be adapted to the noise level ∥w∥

Here we used the notation

Ψ a=∑m a m ψ m
to indicate the reconstruction operator, and J(a)
is the ℓ 1
sparsity prior
J(a)=∑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.

    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.

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

Display it.

    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.

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

    hF = real(fft2(h));

Display the kernel h
and its transform h ˆ
. We use `fftshift` to center the filter for display.

    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.

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

Apply the filter.

    y0 = Phi(f0);

Display the filtered observation.

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

Add some noise of variance σ 2
, to obtain y=Φ f 0+w=f 0⋆h+w
.

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

Display.

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

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

Display a curve of the 1D soft thresholding.

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

We set the parameters of the wavelet transform.

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

Shortcut for Ψ
and Ψ∗
in the orthogonal case.

    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 Ψ
is defined as

S Ψ T(f)=∑m s T(⟨f,ψ m⟩)ψ m

It thus corresponds to applying the transform Ψ∗
, thresholding the coefficients using S T
and then undoing the transform using Ψ
.

S Ψ T(f)=Ψ∘S T∘Ψ∗

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

This soft thresholding corresponds to a denoising operator.

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

## Deconvolution using Orthogonal Wavelet Sparsity

If Ψ
is an orthogonal basis, a change of variable shows that the synthesis prior is also an analysis prior, that reads

f⋆∈argmin f E(f)=1 2∥y−Φ f∥2+λ∑m∥⟨f,ψ m⟩∥.

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(ℓ)
defined as

f(ℓ+1)=S Ψ τ λ(f(ℓ)−τ Φ∗(Φ f(ℓ)−y))

Set up the value of the threshold.

    lambda = .02;

In our setting, since h
is symmetric, one has Φ∗f=Φ f=f⋆h
.

For f(ℓ)
to converge to a solution of the problem, the gradient step size should be chosen as

τ<2∥Φ∗Φ∥

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

    tau = 1.5;

Number of iterations.

    niter = 100;

Initialize the solution.

    fSpars = y;

First step: perform one step of gradient descent of the energy ∥y−f⋆h∥2
.

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

Second step: denoise the solution by thresholding.

    fSpars = SoftThreshPsi( fSpars, lambda*tau );

Exercice 1: (the solution is exo1.m) Perform the iterative soft thresholding. Monitor the decay of the energy E
you are minimizing.

    exo1;

Display the result.

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

Exercice 2: (the solution is exo2.m) Try to find the best threshold λ
. To this end, perform a lot of iterations, and progressively decay the threshold λ
during the iterations. Record the best result in `fBestOrtho`.

    exo2;

Display the result.

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

## Deconvolution using Translation Invariant Wavelet Sparsity

Orthogonal sparsity performs a poor regularization because of the lack of translation invariance. This regularization is enhanced by considering Ψ
as a redundant tight frame of translation invariant wavelets.

One thus looks for optimal coefficients a⋆
that solves

a⋆∈argmin a 1 2∥y−Φ Ψ a∥2+λ J(a)

One should be careful that because of the redundancy of the wavelet tight frame, one should use a weighted ℓ 1
norm, where each coefficient is divided by the number of redundancy at each scale.

J(a)=∑j,k u j∥a j,k∥

where the wavelet coefficients a m=a j,k
are indexed by the scale (and orientation) j
and the location k
.

Compute the scaling factor (inverse of the redundancy).

    J = Jmax-Jmin+1;
    u = [4^(-J) 4.^(-floor(J+2/3:-1/3:1)) ];

Value of the regularization parameter.

    lambda = .01;

Shortcut for the wavelet transform. Important: Note that `PsiS` is the shortcut for Ψ∗
, but `Psi` is the shortcut for Ξ=(Ψ∗)+
that satisfy Ξ Φ∗f=f
.

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

The forward-backward algorithm now compute a series of wavelet coefficients a(ℓ)
computed as

a(ℓ+1)=S τ λ(a(ℓ)+Ψ∗Φ(y−Φ Ξ a(ℓ))).

The soft thresholding is defined as:

∀m,S T(a)m=max(0,1−T/∥a m∥)a m.

The step size should satisfy:

τ<2∥Φ Ψ∥2.

    tau = 1.5;

Initialize the wavelet coefficients with those of the observations.

    a = PsiS(y);

Gradient descent.

    a = a + tau * PsiS( Phi( y-Phi(Psi(a)) ) );

Soft threshold.

    a = SoftThresh( a, lambda*tau );

Important: keep in mind that the prior J(a)
is a weighted ℓ 1
norm, it should thus be computed this way:

    U = repmat( reshape(u,[1 1 length(u)]), [n n 1] );
    Ja = sum(sum(sum( abs(a.*U) )));

Exercice 3: (the solution is exo3.m) Perform the iterative soft thresholding. Monitor the decay of the energy.

    exo3;

Perform the reconstruction.

    fTI = Psi(a);

Display the result.

    clf;
    imageplot(fTI);

Exercice 4: (the solution is exo4.m) Compute the optimal value of λ
, and record the optimal reconstruction `fBestTI`.

    exo4;

Display the result.

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

Exercice 5: (the solution is exo5.m) Compare with the result of TV regularization, record the optimal TV result in `fBestTV`.

    exo5;

Display the result.

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

