Denoising by Sobolev and Total Variation Flows
It is possible to denoise an image using either a variation flow of (linear Sobolev or non-linear TV) or a variational regularization (Sobolev or TV).
Contents
Installing toolboxes and setting up the path.
You need to download the general purpose toolbox and the signal toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ 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 commands 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(path,p); % scilab users must *not* execute this
Then you can add these toolboxes to the path.
% Add some directories to the path getd('toolbox_signal/'); getd('toolbox_general/');
Sobolev and TV Norms
The Sobolev image prior is a linear norm whereas the TV norm is a Banach (non-linear) norm.
First we load an image
n = 256;
M = load_image('lena', n);
M = rescale(M);
We can compute the gradient (and display it using R/G colors) and its norm.
% gradient Gr = grad(M); % norm of gradient 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 L2 norm of grad whereas the Sobolev norm is the L1 norm. They do not have the same scaling so it is difficult to compare them directly.
sob = sum(d(:).^2); tv = sum(d(:));
The gradient of the Sobolev norm is (minus) the Laplacian.
L = div(grad(M)); % display clf; imageplot(M, 'Image', 1,2,1); imageplot(L, 'Laplacian', 1,2,2);
The gradient of the TV norm is not defined if the gradient of M vanishes. One can regularize the TV norm by replacing abs(x) by sqrt(epsilon^2 + x^2), and obtained a well-defined approximate gradient.
t = linspace(-5,5)'; clf; subplot(2,1,1); hold('on'); plot(t, abs(t), 'b'); plot(t, sqrt(.5^2+t.^2), 'r'); title('\epsilon=1/2'); axis('square'); subplot(2,1,2); hold('on'); plot(t, abs(t), 'b'); plot(t, sqrt(1^2+t.^2), 'r'); title('\epsilon=1'); axis('square');
When epsilon increases, the TV gradient approaches the Laplacian (normalized by 1/epsilon).
epsilon_list = [1e-9 1e-2 1e-1 .5]; clf; for i=1:length(epsilon_list) epsilon = epsilon_list(i); G = div( Gr ./ repmat( sqrt( epsilon^2 + d.^2 ) , [1 1 2]) ); imageplot(G, strcat(['epsilon=' num2str(epsilon)]), 2,2,i); end
Heat Diffusion for Denoising
The gradient descent of the Sobolev norm corresponds to the heat diffusion. It can be used to perform denoising.
Load an image and add some noise.
n = 256;
M0 = load_image('cameraman', n);
M0 = rescale(M0,.05,.95);
sigma = .14;
M = M0 + randn(n,n)*sigma;
We set up the parameters for the diffusion.
% maximal smoothing time T = 3; % step size (should be small enough) tau = .05; % number of iterations to reach T. niter = round(T/tau); % initial solution for time t=0. Mheat = M;
An iteration of the diffusion simply corresponds to a gradient step.
% gradient of the Sobolev norm is minus the Laplacian G = -div(grad(Mheat)); % descent Mheat = Mheat - tau*G;
Exercice 1: (the solution is exo1.m) Compute niter iteration of the Heat diffusion. Compute the SNR denoising error err(i) at step i and keep track of the optimal denoising result Mheat0.
exo1;
Display the error err.
clf; plot( linspace(0,T,niter), err ); axis tight; set_label('t', 'SNR');
Display best "oracle" denoising result.
eheat = max(err); enoisy = snr(M0,M); clf; imageplot(clamp(M), strcat(['Noisy ' num2str(enoisy) 'dB']), 1,2,1); imageplot(clamp(Mheat0), strcat(['Heat diffusion ' num2str(eheat) 'dB']), 1,2,2);
Heat Regulariation for Denoising
Sobolev regulariation is very close to heat diffusion, except that a data-fidelity term avoids the diffusion to converge to 0. The solution Msob minimizes norm(M-Msob)^2+lambda*Sobolev(Msob) where Sobolev(Msob) is the Sobolev norm.
First we define a regularization parameter lambda that plays the same role as the time t.
lambda = .5;
Then we set up the paramters for the minimization.
% gradient descent step size, should be small enough tau = .05; % initialization of the minimization Msob = M;
The gradient of the regularization functional includes a drifting term Msob-M that avoid the computed solution Msob to converge to 0.
G = Msob-M-lambda*div(grad(Msob));
% descent step
Msob = Msob - tau*G;
Exercice 2: (the solution is exo2.m) Compute the gradient descent using a small enough step size and monitor the energy norm(M-Msob)^2+lambda*Sobolev(Msob).
exo2;
Exercice 3: (the solution is exo3.m) Compute the solution for several value of lambda and choose the optimal lambda and the corresponding optimal denoising Msob0. You can increase progressively lambda and reduce considerably the number of iterations.
exo3;
Display best "oracle" denoising result.
esob = snr(M0,Msob0); enoisy = snr(M0,M); clf; imageplot(clamp(M), strcat(['Noisy ' num2str(enoisy) 'dB']), 1,2,1); imageplot(clamp(Msob0), strcat(['Sobolev regularization ' num2str(esob) 'dB']), 1,2,2);
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 up the parameters for the diffusion. Not that they are different from the parameter of the heat diffusion because the scalings are different.
% maximal smoothing time T = 1/4; % step size tau = .003; % number of iterations to reach T. niter = round(T/tau); % initial solution for time t=0. Mtv = M;
An iteration of the diffusion simply corresponds to a gradient step. We use a regularized TV norm with a small value of epsilon
% regularization parameter epsilon = 1e-3; % gradient and norm Gr = grad(Mtv); d = sqrt(sum3(Gr.^2,3)); % gradient of the TV norm G = -div( Gr ./ repmat( sqrt( epsilon^2 + d.^2 ) , [1 1 2]) ); % descent Mtv = Mtv - tau*G;
Exercice 4: (the solution is exo4.m) Compute niter iteration of the TV diffusion. Compute the SNR denoising error err(i) at step i and keep track of the optimal denoising result Mtv0.
exo4;
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(M0,M); clf; imageplot(clamp(M), strcat(['Noisy ' num2str(enoisy) 'dB']), 1,2,1); imageplot(clamp(Mtv0), strcat(['TV diffusion ' num2str(etv) 'dB']), 1,2,2);
Total Variation Regulariation for Denoising
TV regulariation is very close to TV diffusion, except that a data-fidelity term avoids the diffusion to converge to 0. The solution Mtvr minimizes norm(M-Mtvr)^2+lambda*TV(Mtvr) where TV(Mtvr) is the TV norm.
First we define a regularization parameter lambda that plays the same role as the time t.
lambda = .1;
Then we set up the paramters for the minimization.
% gradient descent step size, should be small enough tau = .05; % initialization of the minimization Mtvr = M;
The gradient of the regularization functional includes a drifting term Msob-M that avoid the computed solution Msob to converge to 0.
% gradient of TV Gr = grad(Mtvr); d = sqrt(sum3(Gr.^2,3)); G0 = -div( Gr ./ repmat( sqrt( epsilon^2 + d.^2 ) , [1 1 2]) ); % gradient of TV regularization G = Mtvr-M+lambda*G0; % descent step Mtvr = Mtvr - tau*G;
Exercice 5: (the solution is exo5.m) Compute the gradient descent using a small enough step size and monitor the energy 1/2*norm(M-Msob)^2+lambda*Sobolev(Msob).
exo5;
Exercice 6: (the solution is exo6.m) Compute the solution for several value of lambda and choose the optimal lambda and the corresponding optimal denoising Msob0. You can increase progressively lambda and reduce considerably the number of iterations.
exo6;
Display best "oracle" denoising result.
etvr = snr(M0,Mtvr0); clf; imageplot(clamp(M), strcat(['Noisy ' num2str(enoisy) 'dB']), 1,2,1); imageplot(clamp(Mtvr0), strcat(['TV regularization ' num2str(etvr) 'dB']), 1,2,2);
Exercice 7: (the solution is exo7.m) Compare the TV denoising with a hard thresholding in a translation invariant tight frame of wavelets.
exo7;
