?? 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>Image Denoising with Wavelets</title> <meta name="generator" content="MATLAB 7.4"> <meta name="date" content="2008-10-15"> <meta name="m-file" content="index"> <LINK REL="stylesheet" HREF="style.css" TYPE="text/css"> </head> <body> <div class="content"> <h1>Image Denoising with Wavelets</h1> <introduction> <p>This numerical tour uses wavelets to perform both linear and non-linear image denoising.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">Installing toolboxes and setting up the path.</a></li> <li><a href="#8">Thresholding Estimator and Sparsity</a></li> <li><a href="#14">Image loading and adding Gaussian noise.</a></li> <li><a href="#17">Hard Thresholding vs. Soft Thresholding</a></li> <li><a href="#19">Orthogonal Wavelet Denoising</a></li> <li><a href="#26">Estimating the noise level</a></li> <li><a href="#29">Translation Invariant Wavelet Transform</a></li> <li><a href="#33">Translation Invariant Wavelet Denoising</a></li> <li><a href="#39">Wavelet Block Thresholding</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>Thresholding Estimator and Sparsity<a name="8"></a></h2> <p>The idea of non-linear denoising is to use an orthogonal basis in which the coefficients <tt>x</tt> of the signal or image <tt>M0</tt> is sparse (a few large coefficients). In this case, the noisy coefficients <tt>x</tt> of the noisy data <tt>M</tt> (perturbated with Gaussian noise) are <tt>x0+noise</tt> where <tt>noise</tt> is Gaussian. A thresholding set to 0 the noise coefficients that are below <tt>T</tt>. The threshold level <tt>T</tt> should be chosen judiciously to be just above the noise level. </p> <p>First we generate a spiky signal.</p><pre class="codeinput"><span class="comment">% dimension</span>n = 4096;<span class="comment">% probability of spiking</span>rho = .05;<span class="comment">% location of the spike</span>x0 = rand(n,1)<rho;<span class="comment">% random amplitude in [-1 1]</span>x0 = 2 * x0 .* ( rand(n,1)-.5 );</pre><p>We add some gaussian noise</p><pre class="codeinput">sigma = .1;x = x0 + randn(size(x0))*sigma;<span class="comment">% display</span>clf;subplot(2,1,1);plot(x0); axis([1 n -1 1]);set_graphic_sizes([], 20);title(<span class="string">'Original signal'</span>);subplot(2,1,2);plot(x); axis([1 n -1 1]);set_graphic_sizes([], 20);title(<span class="string">'Noisy signal'</span>);</pre><img vspace="5" hspace="5" src="index_01.png"> <p><i>Exercice 1:</i> (the solution is <a href="../private/denoising_wavelet/exo1.m">exo1.m</a>) What is the optimal threshold <tt>T</tt> to remove as much as possible of noise ? Try several values of <tt>T</tt>. </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_02.png"> <p>In order to be optimal without knowing in advance the amplitude of the coefficients of <tt>x0</tt>, one needs to set <tt>T</tt> just above the noise level. This means that <tt>T</tt> should be roughly equal to the maximum value of a Gaussian white noise of size <tt>n</tt>. </p> <p><i>Exercice 2:</i> (the solution is <a href="../private/denoising_wavelet/exo2.m">exo2.m</a>) The theory predicts that the maximum of <tt>n</tt> Gaussian variable of variance <tt>sigma^2</tt> is smaller than <tt>sqrt(2*log(n))</tt> with large probability (that tends to 1 when <tt>n</tt> increases). This is also a sharp result. Check this numerically by computing with Monte Carlo sampling the maximum with <tt>n</tt> increasing (in power of 2). Check also the deviation of the maximum when you perform several trial with <tt>n</tt> fixed. </p><pre class="codeinput">exo2;</pre><img vspace="5" hspace="5" src="index_03.png"> <h2>Image loading and adding Gaussian noise.<a name="14"></a></h2> <p>A simple noise model is additive Gaussian noise.</p> <p>First we load an image.</p><pre class="codeinput">name = <span class="string">'boat'</span>;n = 256;M0 = load_image(name,n);M0 = rescale( M0, .05, .95 );</pre><p>Then we add some gaussian noise to it.</p><pre class="codeinput">sigma = .08; <span class="comment">% noise level</span>M = M0 + sigma*randn(size(M0));clf;imageplot(M0, <span class="string">'Original'</span>, 1,2,1);imageplot(clamp(M), <span class="string">'Noisy'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_04.png"> <h2>Hard Thresholding vs. Soft Thresholding<a name="17"></a></h2> <p>A thresholding is a 1D non-linear function applied to each wavelet coefficients. The most important thresholding are the hard thresholding (related to L0 minimization) and the soft thresholding (related to L1 minimization). </p><pre class="codeinput"><span class="comment">% threshold value</span>T = 1;v = -linspace(-3,3,2000);<span class="comment">% hard thresholding of the t values</span>v_hard = v.*(abs(v)>T);<span class="comment">% soft thresholding of the t values</span>v_soft = max(1-T./abs(v), 0).*v;<span class="comment">% display</span>clf;hold(<span class="string">'on'</span>);plot(v, v_hard);plot(v, v_soft, <span class="string">'r--'</span>);axis(<span class="string">'equal'</span>); axis(<span class="string">'tight'</span>);legend(<span class="string">'Hard thresholding'</span>, <span class="string">'Soft thresholding'</span>);hold(<span class="string">'off'</span>);</pre><img vspace="5" hspace="5" src="index_05.png"> <h2>Orthogonal Wavelet Denoising<a name="19"></a></h2> <p>It is possible to perform non linear denoising by thresholding the wavelet coefficients. This allows to better respect the sharp features of the image. </p> <p>First we compute the wavelet coefficients of the noisy image.</p><pre class="codeinput">options.ti = 0;Jmin = 4;MW = perform_wavelet_transf(M,Jmin,+1,options);</pre><p>Then we hard threshold the coefficients below the noise level. In practice a threshold of <tt>3*sigma</tt> is close to optimal for natural images. </p><pre class="codeinput">T = 3*sigma;MWT = perform_thresholding(MW,T,<span class="string">'hard'</span>);clf;subplot(1,2,1);plot_wavelet(MW,Jmin);title(<span class="string">'Noisy coefficients'</span>);set_axis(0);subplot(1,2,2);plot_wavelet(MWT,Jmin);title(<span class="string">'Thresholded coefficients'</span>);set_axis(0);</pre><img vspace="5" hspace="5" src="index_06.png"> <p>One can then reconstruct from these noisy coefficients.</p><pre class="codeinput">Mhard = perform_wavelet_transf(MWT,Jmin,-1,options);<span class="comment">% display</span>clf;imageplot(clamp(M), <span class="string">'Noisy'</span>, 1,2,1);imageplot(clamp(Mhard), strcat([<span class="string">'Hard denoising, SNR='</span> num2str(snr(M0,Mhard))]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_07.png"> <p>The image suffers from many artifacts (wavelets poping arround). It is possible to improve the result by using soft thresholding. Two important remark should be made </p> <div> <ul> <li>First one must use a lower threshold because soft thresholding also lower the value of non-thresholded coeffcients. In practice, a threshold of <tt>3/2*sigma</tt> works well. </li> <li>The low frequency part of the coefficients must not be thresholded.</li> </ul> </div><pre class="codeinput">T = 3/2*sigma;MWT = perform_thresholding(MW,T,<span class="string">'soft'</span>);<span class="comment">% re-inject the low frequencies</span>MWT(1:2^Jmin,1:2^Jmin) = MW(1:2^Jmin,1:2^Jmin);<span class="comment">% re-construct</span>Msoft = perform_wavelet_transf(MWT,Jmin,-1,options);<span class="comment">% display</span>clf;imageplot(clamp(Mhard), strcat([<span class="string">'Hard denoising, SNR='</span> num2str(snr(M0,Mhard))]), 1,2,1);imageplot(clamp(Msoft), strcat([<span class="string">'Soft denoising, SNR='</span> num2str(snr(M0,Msoft))]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_08.png"> <p><i>Exercice 3:</i> (the solution is <a href="../private/denoising_wavelet/exo3.m">exo3.m</a>) Determine the best threshold <tt>T</tt> for both hard and soft thresholding. To that end, check for <tt>T=alpha*sigma</tt> (for hard) and <tt>T=alpha*sigma/2</tt> (for hard) and compute the denoising error. What can you conclude from these results ? Test with another image. </p><pre class="codeinput">exo3;</pre><img vspace="5" hspace="5" src="index_09.png"> <h2>Estimating the noise level<a name="26"></a></h2> <p>In practice, the noise level <tt>sigma</tt> is unknown. A good estimator is given by the median of the wavelet coefficients at the finer scale. An even simple estimator is given by the normalized derivate along X or Y direction </p> <p>First we extract the high frequency residual.</p><pre class="codeinput">H = (M(1:n-1,:) - M(2:n,:))/sqrt(2);<span class="comment">% histograms</span>[h,t] = hist(H(:), 100);h = h/sum(h);<span class="comment">% display</span>clf;imageplot(H, <span class="string">'High freq. coefficients'</span>, 2,1,1);subplot(2,1,2);bar(t, h);axis([-.5 .5 0 max(h)]);</pre><img vspace="5" hspace="5" src="index_10.png"> <p>The mad estimator (median of median) must be rescaled so that it gives the correct variance for gaussian noise.</p><pre class="codeinput">sigma_est = mad(H(:),1)/0.6745;disp( strcat([<span class="string">'Estimated noise level='</span> num2str(sigma_est), <span class="string">', true='</span> num2str(sigma)]) );</pre><pre class="codeoutput">Estimated noise level=0.090125, true=0.08</pre><h2>Translation Invariant Wavelet Transform<a name="29"></a></h2> <p>Orthogonal wavelet transforms are not translation invariant. It means that the processing of an image and of a translated version of the image give different results. A translation invariant wavelet transform is implemented by ommitting the sub-sampling at each stage of the transform. This correspond to the decomposition of the image in a redundant familly of N*(J+1) atoms where N is the number of pixel and J is the number of scales of the transforms. </p> <p>For Scilab, we need to extend a little the available memory.</p><pre class="codeinput">extend_stack_size(4);</pre><p>The invariant transform is obtained using the same function, by activating the switch <tt>options.ti=1</tt>. </p><pre class="codeinput">options.ti = 1;MW = perform_wavelet_transf(M0,Jmin,+1,options);</pre><p><tt>MW(:,:,1)</tt> corresponds to the low scale residual. Each <tt>MW(:,:,3*j+k+1)</tt> for k=1:3 (orientation) corresponds to a scale of wavelet coefficient, and has the same size as the original image. </p><pre class="codeinput">clf;i = 0;<span class="keyword">for</span> j=1:2 <span class="keyword">for</span> k=1:3 i = i+1; imageplot(MW(:,:,i+1), strcat([<span class="string">'Scale='</span> num2str(j) <span class="string">' Orientation='</span> num2str(k)]), 2,3,i ); <span class="keyword">end</span><span class="keyword">end</span></pre><img vspace="5" hspace="5" src="index_11.png"> <h2>Translation Invariant Wavelet Denoising<a name="33"></a></h2> <p>Orthogonal wavelet denoising does not performs very well because of its lack of translation invariance. A much better result is obtained by not sub-sampling the wavelet transform, which leads to a redundant tight-frame. </p> <p>First we compute the translation invariant wavelet transform</p><pre class="codeinput">options.ti = 1;MW = perform_wavelet_transf(M,Jmin,+1,options);</pre><p>Then we threshold the set of coefficients.</p><pre class="codeinput">T = 3.5*sigma;MWT = perform_thresholding(MW,T,<span class="string">'hard'</span>);</pre><p>We can display some wavelets coefficients</p><pre class="codeinput">J = size(MW,3)-5;clf;imageplot(MW(:,:,J), <span class="string">'Noisy coefficients'</span>, 1,2,1);imageplot(MWT(:,:,J), <span class="string">'Thresholded coefficients'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_12.png"> <p>We can now reconstruct</p><pre class="codeinput">Mti = perform_wavelet_transf(MWT,Jmin,-1,options);<span class="comment">% display</span>clf;imageplot(clamp(Msoft), strcat([<span class="string">'Soft orthogonal, SNR='</span> num2str(snr(M0,Msoft))]), 1,2,1);imageplot(clamp(Mti), strcat([<span class="string">'Hard invariant, SNR='</span> num2str(snr(M0,Mti))]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_13.png"> <p><i>Exercice 4:</i> (the solution is <a href="../private/denoising_wavelet/exo4.m">exo4.m</a>) Determine the best threshold <tt>T</tt> for both hard and soft thresholding, but now in the translation invariant case. What can you conclude ? </p><pre class="codeinput">exo4;</pre><img vspace="5" hspace="5" src="index_14.png"> <h2>Wavelet Block Thresholding<a name="39"></a></h2> <p>Wavelets coefficients of natural images are not independant one from each other. One can thus improve the denoising results by thresholding block of coefficients togethers. Block thresholding is only efficient when used as a soft thresholder. </p> <p>You can perform the block thresholding for an arbitrary block size.</p><pre class="codeinput">options.ti = 0;MW = perform_wavelet_transf(M,Jmin,+1,options);<span class="comment">% soft block thresholding</span>T = 2.5*sigma/2;options.block_size = 4;MWT = perform_thresholding(MW,T,<span class="string">'block'</span>,options);<span class="comment">% display</span>plot_wavelet(MWT,Jmin);</pre><img vspace="5" hspace="5" src="index_15.png"> <p>You can reconstruct the image. Test with several values for <tt>T</tt> in order to determine the best threshold. </p><pre class="codeinput">Mblock = perform_wavelet_transf(MWT,Jmin,-1,options);<span class="comment">% display</span>clf;imageplot(clamp(Msoft), strcat([<span class="string">'Soft orthogonal, SNR='</span> num2str(snr(M0,Msoft))]), 1,2,1);imageplot(clamp(Mblock), strcat([<span class="string">'Block thresholding, SNR='</span> num2str(snr(M0,Mblock))]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_16.png"> <p><i>Exercice 5:</i> (the solution is <a href="../private/denoising_wavelet/exo5.m">exo5.m</a>) Try block thresholding for a variety of block size and determine the best SNR. </p><pre class="codeinput">exo5;</pre><img vspace="5" hspace="5" src="index_17.png"> <p>Block thresholding can also be applied to a translation invariant wavelet transform. It gives state of the art denoising results.</p><pre class="codeinput"><span class="comment">% transform</span>options.ti = 1;MW = perform_wavelet_transf(M,Jmin,+1,options);<span class="comment">% threshold</span>T = 2.5*sigma/2;options.block_size = 5;MWT = perform_thresholding(MW,T,<span class="string">'block'</span>,options);<span class="comment">% transform back</span>Mblockti = perform_wavelet_transf(MWT,Jmin,-1,options);<span class="comment">% display</span>clf;imageplot(clamp(Mti), strcat([<span class="string">'Hard TI, SNR='</span> num2str(snr(M0,Mti))]), 1,2,1);imageplot(clamp(Mblockti), strcat([<span class="string">'Block TI, SNR='</span> num2str(snr(M0,Mblockti))]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_18.png"> <p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Image Denoising with Wavelets% This numerical tour uses wavelets to perform both linear and non-linear% image denoising.%% 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/');%% Thresholding Estimator and Sparsity% The idea of non-linear denoising is to use an orthogonal basis in which% the coefficients |x| of the signal or image |M0| is sparse (a few large% coefficients). In this case, the noisy coefficients |x| of the noisy% data |M| (perturbated with Gaussian noise) are |x0+noise| where |noise|% is Gaussian. A thresholding set to 0 the noise coefficients that are% below |T|. The threshold level |T| should be chosen judiciously to be% just above the noise level.%%% First we generate a spiky signal.% dimensionn = 4096;% probability of spikingrho = .05;% location of the spikex0 = rand(n,1)<rho;% random amplitude in [-1 1]x0 = 2 * x0 .* ( rand(n,1)-.5 );%%% We add some gaussian noisesigma = .1;x = x0 + randn(size(x0))*sigma;% display?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -