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
jsMath v3.6e (Image symbol fonts) [help]

• Quick Start
• Missing Fonts
• Printing Issues
• Home Page
• Documentation
• User Community
• Known Bugs
• Report a Bug
• License
Click the jsMath button or ALT-click on mathematics to reopen this panel.

Loading jsMath-controls.html

# Denoising by Sobolev and Total Variation Flows

It is possible to denoise an image using either a variation flow of some regularization prior. This leads to a linear or non-linear PDE. We consider here the linear Sobolev and non-linear TV flow, that corresponds to the heat and mean-curvature differential equations.

## Contents

- Installing toolboxes and setting up the path.
- Prior and Gradient Flows
- Sobolev Prior
- Heat Diffusion for Denoising
- Total Variation Prior
- Total Variation Diffusion 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 Gradient Flows

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

A PDE flow performs a continuous gradient descent of the energy. It defines a continuous series of images f t(x) that satisfies

t f t=−GradJ(f t)

where GradJ(f) is the gradient of the functional f J(f), and hence is also an image.

In a discrete setting, this continuous flow is discretized using an explicit Euler scheme in time, which corresponds to a gradient descent algorithm to minimize the energy J(f).

It defines a series of images f() indexed by N as

f(+1)=f()+GradJ(f())

One can sees that f()f t for t=, and the agreement is better if is small, and hence the number of steps is large.

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

2 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 Diffusion for Denoising

The gradient descent of the Sobolev norm corresponds to the heat diffusion. It can be used to perform denoising.

Add some noise to the original image.

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

The gradient of the Sobolev norm is (minus) the Laplacian.

GradJ(f)=−f=−div(f)

The Heat diffusion PDE thus reads:

t f t=f t

Initialize the PDE at time t=0 using the noisy image, `f_0 = y`.

    fHeat = y;

The Laplacian is computed using a finite differences divergence

    L = div(grad(fHeat));

Display.

    clf;
    imageplot(fHeat, 'Image', 1,2,1);
    imageplot(L, 'Laplacian', 1,2,2);

For the discretized PDE in time to minimize J(f), the step size for the gradient descent should satisfy:

2=4 1

To have a good precision, we use an even smaller step size:

    tau = .05;

Stopping time for the PDE resolution.

    T = 3;

Number of iterations to reach T.

    niter = ceil(T/tau);

Gradient of the Sobolev norm is minus the Laplacian

    G = -div(grad(fHeat));

Descent.

    fHeat = fHeat - tau*G;

Exercice 1: (the solution is exo1.m) Compute `niter` iterations of the Heat diffusion. Keep track of `err(i)=snr(f0,fHeat)`. Keep track of the optimal denoising result `fHeat0`.

    exo1;

Plot the evolution of `snr(f0,fHeat)` during the diffusion.

    clf;
    plot( linspace(0,T,niter), err ); axis tight;
    set_label('t', 'SNR');

Display best "oracle" denoising result.

    eheat = max(err); enoisy = snr(f0,y);
    clf;
    imageplot(clamp(y), strcat(['Noisy ' num2str(enoisy) 'dB']), 1,2,1);
    imageplot(clamp(fHeat0), strcat(['Heat diffusion ' num2str(eheat) '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)2 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 Diffusion for Denoising

The gradient descent of the TV norm corresponds to a non-linear, edge-preservig diffusion. It can be used to perform denoising.

We set a small enough regularization parameter.

    epsilon = 1e-2;

Maximal smoothing time. Not that it is different from the parameter of the heat diffusion because the scalings are different.

    T = 1/4;

The step size for diffusion should satisfy:

2=4

    tau = epsilon/4;

Number of iterations to reach T.

    niter = ceil(T/tau);

Initial solution for time t=0.

    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*G;

Exercice 3: (the solution is exo3.m) Compute `niter` iterations of the smoothed TV diffusion. Keep track of `err(i)=snr(f0,fTV)`. Keep track of the optimal denoising result `fTV0`.

    exo3;

Plot the evolution of `snr(f0,fTV)` during the diffusion. Display the error `err`.

    clf;
    plot( linspace(0,T,niter), err ); axis tight;
    set_label('t', 'SNR');

Display best "oracle" denoising result.

    etv = max(err);
    enoisy = snr(f0,fTV0);
    clf;
    imageplot(clamp(y), strcat(['Noisy ' num2str(enoisy) 'dB']), 1,2,1);
    imageplot(clamp(fTV0), strcat(['TV diffusion ' num2str(etv) 'dB']), 1,2,2);