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.

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) \in \mathbb{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 \[ \frac{\partial f_t}{\partial t} = -\text{Grad}J(f_t) \]

where \(\text{Grad}J(f)\) is the gradient of the functional \(f \mapsto 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^{(\ell)}\) indexed by \(\ell \in \mathbb{N}\) as \[ f^{(\ell+1)} = f^{(\ell)} + \tau \text{Grad}J(f^{(\ell)}). \]

One can sees that \(f^{(\ell)} \approx f_t\) for \(t = \ell \tau\), and the agreement is better if \(\tau\) is small, and hence the number of steps \(\ell\) is large.

Note that for \(f^{(\ell)}\) to converge with \(\ell \rightarrow +\infty\) toward a minimizer of \(J\), \(\tau\) needs to be small enough, more precisely, if the functional \(J\) is twice differentiable, \[ \tau < \frac{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) = \int \|\nabla f(x)\|^2 d x \]

Where the gradient vector at point \(x\) is defined as \[ \nabla f(x) = \left( \frac{\partial f(x)}{\partial x_1},\frac{\partial f(x)}{\partial x_2} \right) \]

For a discrete pixelized image \(f \in \mathbb{R}^N\) where \(N=n \times n\) is the number of pixel, \(\nabla f(x) \in \mathbb{R}^2\) is computed using finite difference.

Gr = grad(f0);

One can compute the norm of gradient, \(d(x) = \|\nabla 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 \(\nabla f \in \mathbb{R}^{N \times 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. \[\text{Grad}J(f) = -\Delta f = -\text{div}(\nabla f)\]

The Heat diffusion PDE thus reads: \[ \frac{\partial f_t}{\partial t} = \Delta 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: \[ \tau < \frac{2}{\|\Delta\|} = \frac{1}{4} \]

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)=\int \|\nabla f(x)\| d x\]

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)=\int_{-\infty}^{+\infty} L( S_t(f) ) dt\]

Where \(S_t(f)\) is the level set at \(t\) of the image \(f\) \[S_t(f)=\{ x \backslash f(x)=t \}\]

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

exo2;

The Gradient of the TV norm is \[ \text{Grad}J(f) = \text{div}\left( \frac{\nabla f}{\|\nabla f\|} \right) . \]

The gradient of the TV norm is not defined if at a pixel \(x\) one has \(\nabla 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_\epsilon(f) = \int \sqrt{ \varepsilon^2 + \| \nabla f(x) \|^2 } d x\]

This corresponds to replacing \(\|u\|\) by \( \sqrt{\varepsilon^2 + \|u\|^2} \) 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 \[ \text{Grad}J(f) = \text{div}\left( \frac{\nabla f}{\sqrt{\varepsilon^2 + \|\nabla f\|^2}} \right) . \]

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

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: \[ \tau \leq \frac{2 \varepsilon}{\|\Delta\|} = \frac{\varepsilon}{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');