To print higher-resolution math symbols, click the Hi-Res Fonts for Printing button on the jsMath control panel.

No jsMath TeX fonts found -- using image fonts instead. These may be slow and might not print well. Use the jsMath control panel to get additional information. jsMath Control Panel Hide this Message

jsMath
Processing Math: Done

# Denoising by Sobolev and Total Variation Regularization

This numerical tour explores the use of variational minimization to perform denoising. It consider the Sobolev and the Total Variation regularization functional (priors).

## Contents

- Installing toolboxes and setting up the path.
- Prior and Regularization
- Sobolev Prior
- Heat Regularization for Denoising
- Total Variation Prior
- Total Variation Regulariation for Denoising

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

## Prior and Regularization

For a given image f, a prior J(f)R assign a score supposed to be small for the image of interest.

A denoising of some noisy image y is obtained by a variational minimization that mixes a fit to the data (usually using an L 2 norm) and the prior.

min f 2 1 y−f 2+J(f)

If J(f) is a convex function of f, then the minimum exists and is unique.

The parameter 0 should be adapted to the noise level. Increasing its value means a more agressive denoising.

If J(f) is a smooth functional of the image f, a minimizer of this problem can be computed by gradient descent. It defines a series of images f() indexed by N as

f(+1)=f()+f()−y+GradJ(f())

Note that for f() to converge with + toward a solution of the problem, needs to be small enough, more precisely, if the functional J is twice differentiable,

2 1+max f D 2 J(f)

## Sobolev Prior

The Sobolev image prior is a quadratic prior, i.e. an Hilbert (pseudo)-norm.

First we load a clean image.

    n = 256;
    name = 'hibiscus';
    f0 = load_image(name,n);
    f0 = rescale( sum(f0,3) );

For a smooth continuous function f it is defined as

J(f)=f(x)2 dx

Where the gradient vector at point x is defined as

f(x)=x 1 f(x)x 2 f(x)

For a discrete pixelized image f R N where N=n n is the number of pixel, f(x)R 2 is computed using finite difference.

    Gr = grad(f0);

One can compute the norm of gradient, d(x)=f(x).

    d = sqrt(sum3(Gr.^2,3));

Display.

    clf;
    imageplot(Gr, strcat(['grad']), 1,2,1);
    imageplot(d, strcat(['|grad|']), 1,2,2);

The Sobolev norm is the (squared) L 2 norm of f R N 2.

    sob = sum(d(:).^2);

## Heat Regularization for Denoising

Heat regularization smoothes the image using a low pass filter. Increasing the value of `\lambda` increases the amount of smoothing.

Add some noise to the original image.

    sigma = .1;
    y = f0 + randn(n,n)*sigma;

The solution f of the Sobolev regularization can be computed exactly using the Fourier transform.

f()=y()1+S()where S()=2

This shows that f is a filtering of y.

Useful for later: Fourier transform of the observations.

    yF = fft2(y);

Compute the Sobolev prior penalty `S` (rescale to [0,1]).

    x = [0:n/2-1, -n/2:-1];
    [Y,X] = meshgrid(x,x);
    S = (X.^2 + Y.^2)*(2/n)^2;

Regularization parameter:

    lambda = 20;

Perform the denoising by filtering.

    fSob = real( ifft2( yF ./ ( 1 + lambda*S) ) );

Display.

    clf;
    imageplot(clamp(fSob));

Exercice 1: (the solution is exo1.m) Compute the solution for several value of and choose the optimal `lambda` and the corresponding optimal denoising `fSob0`. You can increase progressively lambda and reduce considerably the number of iterations.

    exo1;

Display best "oracle" denoising result.

    esob = snr(f0,fSob0);  enoisy = snr(f0,y);
    clf;
    imageplot(clamp(y), strcat(['Noisy ' num2str(enoisy,3) 'dB']), 1,2,1);
    imageplot(clamp(fSob0), strcat(['Sobolev regularization ' num2str(esob,3) 'dB']), 1,2,2);

## Total Variation Prior

The total variation is a Banach norm. On the contrary to the Sobolev norm, it is able to take into account step edges.

The total variation of a smooth image f is defined as

J(f)=f(x)dx

It is extended to non-smooth images having step discontinuities.

The total variation of an image is also equal to the total length of its level sets.

J(f)=−+L(S t(f))dt

Where S t(f) is the level set at t of the image f

S t(f)=xf(x)=t

Exercice 2: (the solution is exo2.m) Compute the total variation of `f0`.

    exo2;

The Gradient of the TV norm is

GradJ(f)=div f f

The gradient of the TV norm is not defined if at a pixel x one has f(x)=0. This means that the TV norm is difficult to minimize, and its gradient flow is not well defined.

To define a gradient flow, we consider instead a smooth TV norm

J(f)=2+f(x)2 dx

This corresponds to replacing u by 2+u2 which is a smooth function.

We display (in 1D) the smoothing of the absolute value.

    u = linspace(-5,5)';
    clf;
    subplot(2,1,1); hold('on');
    plot(u, abs(u), 'b');
    plot(u, sqrt(.5^2+u.^2), 'r');
    title('\epsilon=1/2'); axis('square');
    subplot(2,1,2); hold('on');
    plot(u, abs(u), 'b');
    plot(u, sqrt(1^2+u.^2), 'r');
    title('\epsilon=1'); axis('square');

The Gradient of the smoothed TV norm is

GradJ(f)=div f 2+f 2

When increases, the smoothed TV gradient approaches the Laplacian (normalized by 1).

    epsilon_list = [1e-9 1e-2 1e-1 .5];
    clf;
    for i=1:length(epsilon_list)
        G = div( Gr ./ repmat( sqrt( epsilon_list(i)^2 + d.^2 ) , [1 1 2]) );
        imageplot(G, strcat(['epsilon=' num2str(epsilon_list(i))]), 2,2,i);
    end

## Total Variation Regulariation for Denoising

Total variation regularization was introduced by Rudin, Osher and Fatemi, to better respect the edge of image than linear filtering.

We set a small enough regularization parameter.

    epsilon = 1e-2;

Define the regularization parameter .

    lambda = .1;

The step size for diffusion should satisfy:

2 1+8

In practice a larger step size can be used.

    tau = 2 / ( 1 + lambda * 8 / epsilon);

Initialization of the minimization.

    fTV = y;

Compute the gradient of the smoothed TV norm.

    Gr = grad(fTV);
    d = sqrt(sum3(Gr.^2,3));
    G = -div( Gr ./ repmat( sqrt( epsilon^2 + d.^2 ) , [1 1 2]) );

One step of descent.

    fTV = fTV - tau*( y-fTV + lambda* G);

Exercice 3: (the solution is exo3.m) Compute the gradient descent and monitor the minimized energy.

    exo3;

Exercice 4: (the solution is exo4.m) Compute the solution for several value of and choose the optimal and the corresponding optimal denoising `fSob0`. You can increase progressively and reduce considerably the number of iterations.

    exo4;

Display best "oracle" denoising result.

    etvr = snr(f0,fTV0);
    clf;
    imageplot(clamp(y), strcat(['Noisy ' num2str(enoisy,3) 'dB']), 1,2,1);
    imageplot(clamp(fTV0), strcat(['TV regularization ' num2str(etvr,3) 'dB']), 1,2,2);

Exercice 5: (the solution is exo5.m) Compare the TV denoising with a hard thresholding in a translation invariant tight frame of wavelets.

    exo5;