index.html

信号处理系列导航

%% Total Variation Minizing Flow
% A edge respecting flow is obtained by performing a gradient descent of
% the total variation of the image. The TV norm is not differentiable at an
% image with a point of vanighing gradient. It means one needs to slightly
% regularize the TV norm to avoid numerical instability.

%%
% Load an image

n = 256;
M = load_image('lena', n);
M = rescale(M);

%%
% The sub-gradient of the total variation is minus the divergence of the
% normalized gradient of the image. You must be carefull to avoid division by zero during the normalization,
% which corresponds to smoothing a little the TV normmm. 

% gradient
G = grad(M);
% norm of gradient
d = sqrt(sum(G.^2,3));
% avoid instabilities
epsilon = 1e-5;
d = max(d,epsilon);
% normalize
G = G ./ repmat(d, [1 1 2]);
% the gradient of the TV norm
G = div(G);
% display
clf;
imageplot(d, 'Gradient magnitude', 1,2,1)
imageplot(G, 'TV Gradient', 1,2,2);

%%
% _Exercice 2:_ (the solution is <../private/tv_minimization/exo2.m exo2.m>)
% Compute a TV flow by performing a gradient descent of the TV norm.
% You can use a step size |tau = .003| for instance (it needs to be small
% enough) and the number of iterations is |niter=T/tau|, where |T| is the
% final time of the flow (should be like 1/2).

exo2;

%% Total Variation Flow for Denoising
% The TV flow can also be used to perform enoising. The difficulty is to
% choose the correct time |T|. It should be large enoug to remove lots of
% noise, but not too large to avoid destroying the edges in the image.

%% 
% Load an image and add some noise.

n = 256;
M0 = load_image('lena', n);
M0 = rescale(M0,.05,.95);
sigma = .1;
M = M0 + randn(n,n)*sigma;

%%
% _Exercice 3:_ (the solution is <../private/tv_minimization/exo3.m exo3.m>)
% Compute the total variation minimization flow, and record the time T
% optimal that gives the smallest denoising error, and record the optimal denoising |Mopt|.

exo3;

%%
% Display the result.

clf;
imageplot(clamp(M), strcat(['Noisy, SNR=' num2str(snr(M0,M),3) 'dB']), 1,2,1 );
imageplot(clamp(Mopt), strcat(['TV Flow Denoised, SNR=' num2str(snr(M0,Mopt),3) 'dB']), 1,2,2 );


%% Total Variation Regularization
% To avoid the difficulty of regularizing the TV norm (which is not differentiable), one can instead
% compute the solution of a variational regularization to denoise an image
% |M0| from obervations |M=M0+W|

%%
%   |min_Mtv 1/2*norm(M-Mtv,'fro')^2 + lambda*normTV(Mtv)|    (*)

%%
% Where |normTV| is the TV norm of the image.

%%
% The parameter |lambda| is related to the time |T| of the TV flow. It
% should be adjusted to match the noise level |W| (a larger |lambda|
% implies a stronger denoising).


%%
% The minimization is performed using the algorithm proposed by Antonin
% Chambolle in

%%
% * _An Algorithm for Total Variation Minimization and Applications_, Journal of Mathematical Imaging and Vision, 20(1-2), 2004
% 

%%
% The idea is to replace the optimization of the image |Mtv| by the
% optimization of a vector field |G| that is related to |Mtv| by
% |Mtv=M-lambda*div(G)|. The vector field is the one that minimize 

%%
%  |min_G norm(M-lambda*div(G),'fro')| (**)

%%
% Subject to the constraints that all the entries |G(i,j,k)| of G are
% smaller than 1 in magnitude. 


%%
% Chambolle proposed to perform this minimization (*) with a fixed point
% iteration. Another approach is simply to perform a projected gradient
% descent on |G|.

%%
% We initialize the iteration with a zero vector field.

lambda = .2;
G = zeros(n,n,2);

%%
% Do a step of gradient descent of |norm(M-lambda*div(G))^2| followed by a
% projection to ensure that the entries are smaller than 1. The gradient
% descent step should be smaller than 1/4 for the iterations to converge.

% step size
tau = 1/4;
% gradient of the energy
dG = grad( div(G) - M/lambda );
% gradient descent
G = G + tau*dG;
% display
clf;
imageplot( G, 'Vector field before projection', 1, 2, 1 );
% projection on Linfty constraints
d = repmat( sqrt(sum(G.^2,3)), [1 1 2] ); % norm of the vectors
G = G ./ max(d,1);
imageplot( G, 'Vector field after projection', 1, 2, 2 );

%%
% _Exercice 4:_ (the solution is <../private/tv_minimization/exo4.m exo4.m>)
% Perform Chambolle algorithm to find the solution of (**) and thus (*) by iterating
% this gradient descent. Reconstruct the denoise image according to the
% |Mtv=M-lambda*div(G)|

exo4;

%%
% Display the denoised image
clf;
imageplot(M, 'Noisy image', 1,2,1);
imageplot(Mtv, strcat(['TV regularized using \lambda=' num2str(lambda)]), 1,2,2);


%%
% The issue with the minimization (*) is that it is difficult to set the
% value of |lambda| for denoising. It is often simpler to solve the
% following constraint minimization.

%%
%   |min_Mtv normTV(Mtv)| subject to |norm(M-Mtv,'fro') < epsilon|    (***)

%%
% Where |epsilon| is the noise level, set to |epsilon = n*sigma where |sigma| is 
% the standard deviation of the image. 

%%
% The constraint minimization (***) can be solved by modifying the value of
% |lambda| at each iteration of Chambolle's algorithm. If the solution, at
% a given step, is |Mtv = M-lambda*div(G)|, then the |lambda| is modified according to 

%%
%   |lambda <- lambda * n*sigma / norm( M-Mtv, 'fro' )|

%%
% _Exercice 5:_ (the solution is <../private/tv_minimization/exo5.m exo5.m>)
% Implement Chambolle's algorithm with this update rule for |lambda|.
% Display the evolution of |lambda| and |norm( M-Mtv,'fro')-n*sigma| during the
% iterations.

exo5;



%%
% Display the result

clf;
imageplot(clamp(M), strcat(['Noisy, SNR=' num2str(snr(M0,M),3) 'dB']), 1,2,1 );
imageplot(clamp(Mtv), strcat(['TV Regularization Denoised, SNR=' num2str(snr(M0,Mtv),3) 'dB']), 1,2,2 );
##### SOURCE END #####
-->
   </body>
</html>