📄 index.html
字号:
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"><html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Total Variation Minimization</title> <meta name="generator" content="MATLAB 7.4"> <meta name="date" content="2008-10-17"> <meta name="m-file" content="index"> <LINK REL="stylesheet" HREF="style.css" TYPE="text/css"> </head> <body> <div class="content"> <h1>Total Variation Minimization</h1> <introduction> <p>This numerical tour explores total variation minimization and regularization.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">Installing toolboxes and setting up the path.</a></li> <li><a href="#8">Total Variant in Images</a></li> <li><a href="#10">Total Variation Minizing Flow</a></li> <li><a href="#14">Total Variation Flow for Denoising</a></li> <li><a href="#18">Total Variation Regularization</a></li> </ul> </div> <h2>Installing toolboxes and setting up the path.<a name="1"></a></h2> <p>You need to download the <a href="../toolbox_general.zip">general purpose toolbox</a> and the <a href="../toolbox_signal.zip">signal toolbox</a>. </p> <p>You need to unzip these toolboxes in your working directory, so that you have <tt>toolbox_general/</tt> and <tt>toolbox_signal/</tt> in your directory. </p> <p><b>For Scilab user:</b> you must replace the Matlab comment '%' by its Scilab counterpart '//'. </p> <p><b>Recommandation:</b> You should create a text file named for instance <tt>numericaltour.sce</tt> (in Scilabe) or <tt>numericaltour.m</tt> to write all the Scilab/Matlab command you want to execute. Then, simply run <tt>exec('numericaltour.sce');</tt> (in Scilab) or <tt>numericaltour;</tt> (in Matlab) to run the commands. </p> <p>Execute this line only if you are using Matlab.</p><pre class="codeinput">getd = @(p)path(path,p); <span class="comment">% scilab users must *not* execute this</span></pre><p>Then you can add these toolboxes to the path.</p><pre class="codeinput"><span class="comment">% Add some directories to the path</span>getd(<span class="string">'toolbox_signal/'</span>);getd(<span class="string">'toolbox_general/'</span>);</pre><h2>Total Variant in Images<a name="8"></a></h2> <p>The total variation of an image is the total length of its level sets. It is computed as the L1 norme of the gradient, viewed as a complex operator, which is the sum of the length of all the gradient vectors. </p> <p><i>Exercice 1:</i> (the solution is <a href="../private/tv_minimization/exo1.m">exo1.m</a>) Compute the total variation of several images. </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_01.png"> <h2>Total Variation Minizing Flow<a name="10"></a></h2> <p>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. </p> <p>Load an image</p><pre class="codeinput">n = 256;M = load_image(<span class="string">'lena'</span>, n);M = rescale(M);</pre><p>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. </p><pre class="codeinput"><span class="comment">% gradient</span>G = grad(M);<span class="comment">% norm of gradient</span>d = sqrt(sum(G.^2,3));<span class="comment">% avoid instabilities</span>epsilon = 1e-5;d = max(d,epsilon);<span class="comment">% normalize</span>G = G ./ repmat(d, [1 1 2]);<span class="comment">% the gradient of the TV norm</span>G = div(G);<span class="comment">% display</span>clf;imageplot(d, <span class="string">'Gradient magnitude'</span>, 1,2,1)imageplot(G, <span class="string">'TV Gradient'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_02.png"> <p><i>Exercice 2:</i> (the solution is <a href="../private/tv_minimization/exo2.m">exo2.m</a>) Compute a TV flow by performing a gradient descent of the TV norm. You can use a step size <tt>tau = .003</tt> for instance (it needs to be small enough) and the number of iterations is <tt>niter=T/tau</tt>, where <tt>T</tt> is the final time of the flow (should be like 1/2). </p><pre class="codeinput">exo2;</pre><img vspace="5" hspace="5" src="index_03.png"> <h2>Total Variation Flow for Denoising<a name="14"></a></h2> <p>The TV flow can also be used to perform enoising. The difficulty is to choose the correct time <tt>T</tt>. It should be large enoug to remove lots of noise, but not too large to avoid destroying the edges in the image. </p> <p>Load an image and add some noise.</p><pre class="codeinput">n = 256;M0 = load_image(<span class="string">'lena'</span>, n);M0 = rescale(M0,.05,.95);sigma = .1;M = M0 + randn(n,n)*sigma;</pre><p><i>Exercice 3:</i> (the solution is <a href="../private/tv_minimization/exo3.m">exo3.m</a>) Compute the total variation minimization flow, and record the time T optimal that gives the smallest denoising error, and record the optimal denoising <tt>Mopt</tt>. </p><pre class="codeinput">exo3;</pre><img vspace="5" hspace="5" src="index_04.png"> <p>Display the result.</p><pre class="codeinput">clf;imageplot(clamp(M), strcat([<span class="string">'Noisy, SNR='</span> num2str(snr(M0,M),3) <span class="string">'dB'</span>]), 1,2,1 );imageplot(clamp(Mopt), strcat([<span class="string">'TV Flow Denoised, SNR='</span> num2str(snr(M0,Mopt),3) <span class="string">'dB'</span>]), 1,2,2 );</pre><img vspace="5" hspace="5" src="index_05.png"> <h2>Total Variation Regularization<a name="18"></a></h2> <p>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 <tt>M0</tt> from obervations <tt>M=M0+W</tt></p><pre> |min_Mtv 1/2*norm(M-Mtv,'fro')^2 + lambda*normTV(Mtv)| (*)</pre><p>Where <tt>normTV</tt> is the TV norm of the image. </p> <p>The parameter <tt>lambda</tt> is related to the time <tt>T</tt> of the TV flow. It should be adjusted to match the noise level <tt>W</tt> (a larger <tt>lambda</tt> implies a stronger denoising). </p> <p>The minimization is performed using the algorithm proposed by Antonin Chambolle in</p> <div> <ul> <li><i>An Algorithm for Total Variation Minimization and Applications</i>, Journal of Mathematical Imaging and Vision, 20(1-2), 2004 </li> </ul> </div> <p>The idea is to replace the optimization of the image <tt>Mtv</tt> by the optimization of a vector field <tt>G</tt> that is related to <tt>Mtv</tt> by <tt>Mtv=M-lambda*div(G)</tt>. The vector field is the one that minimize </p><pre>|min_G norm(M-lambda*div(G),'fro')| (**)</pre><p>Subject to the constraints that all the entries <tt>G(i,j,k)</tt> of G are smaller than 1 in magnitude. </p> <p>Chambolle proposed to perform this minimization (*) with a fixed point iteration. Another approach is simply to perform a projected gradient descent on <tt>G</tt>. </p> <p>We initialize the iteration with a zero vector field.</p><pre class="codeinput">lambda = .2;G = zeros(n,n,2);</pre><p>Do a step of gradient descent of <tt>norm(M-lambda*div(G))^2</tt> 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. </p><pre class="codeinput"><span class="comment">% step size</span>tau = 1/4;<span class="comment">% gradient of the energy</span>dG = grad( div(G) - M/lambda );<span class="comment">% gradient descent</span>G = G + tau*dG;<span class="comment">% display</span>clf;imageplot( G, <span class="string">'Vector field before projection'</span>, 1, 2, 1 );<span class="comment">% projection on Linfty constraints</span>d = repmat( sqrt(sum(G.^2,3)), [1 1 2] ); <span class="comment">% norm of the vectors</span>G = G ./ max(d,1);imageplot( G, <span class="string">'Vector field after projection'</span>, 1, 2, 2 );</pre><img vspace="5" hspace="5" src="index_06.png"> <p><i>Exercice 4:</i> (the solution is <a href="../private/tv_minimization/exo4.m">exo4.m</a>) Perform Chambolle algorithm to find the solution of (<b>*) and thus (</b>) by iterating this gradient descent. Reconstruct the denoise image according to the <tt>Mtv=M-lambda*div(G)</tt></p><pre class="codeinput">exo4;</pre><img vspace="5" hspace="5" src="index_07.png"> <p>Display the denoised image</p><pre class="codeinput">clf;imageplot(M, <span class="string">'Noisy image'</span>, 1,2,1);imageplot(Mtv, strcat([<span class="string">'TV regularized using \lambda='</span> num2str(lambda)]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_08.png"> <p>The issue with the minimization (*) is that it is difficult to set the value of <tt>lambda</tt> for denoising. It is often simpler to solve the following constraint minimization. </p><pre> |min_Mtv normTV(Mtv)| subject to |norm(M-Mtv,'fro') < epsilon| (***)</pre><p>Where <tt>epsilon</tt> is the noise level, set to <tt>epsilon = n*sigma where |sigma</tt> is the standard deviation of the image. </p> <p>The constraint minimization (<b>*</b>) can be solved by modifying the value of <tt>lambda</tt> at each iteration of Chambolle's algorithm. If the solution, at a given step, is <tt>Mtv = M-lambda*div(G)</tt>, then the <tt>lambda</tt> is modified according to </p><pre> |lambda <- lambda * n*sigma / norm( M-Mtv, 'fro' )|</pre><p><i>Exercice 5:</i> (the solution is <a href="../private/tv_minimization/exo5.m">exo5.m</a>) Implement Chambolle's algorithm with this update rule for <tt>lambda</tt>. Display the evolution of <tt>lambda</tt> and <tt>norm( M-Mtv,'fro')-n*sigma</tt> during the iterations. </p><pre class="codeinput">exo5;</pre><img vspace="5" hspace="5" src="index_09.png"> <p>Display the result</p><pre class="codeinput">clf;imageplot(clamp(M), strcat([<span class="string">'Noisy, SNR='</span> num2str(snr(M0,M),3) <span class="string">'dB'</span>]), 1,2,1 );imageplot(clamp(Mtv), strcat([<span class="string">'TV Regularization Denoised, SNR='</span> num2str(snr(M0,Mtv),3) <span class="string">'dB'</span>]), 1,2,2 );</pre><img vspace="5" hspace="5" src="index_10.png"> <p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Total Variation Minimization% This numerical tour explores total variation minimization and% regularization.%% Installing toolboxes and setting up the path.%%% You need to download the % <../toolbox_general.zip general purpose toolbox>% and the <../toolbox_signal.zip 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 Scilabe) or |numericaltour.m| 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(path,p); % scilab users must *not* execute this%%% Then you can add these toolboxes to the path.% Add some directories to the pathgetd('toolbox_signal/');getd('toolbox_general/');%% Total Variant in Images% The total variation of an image is the total length of its level sets. It% is computed as the L1 norme of the gradient, viewed as a complex% operator, which is the sum of the length of all the gradient vectors.%%% _Exercice 1:_ (the solution is <../private/tv_minimization/exo1.m exo1.m>)% Compute the total variation of several images.exo1;⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -