?? 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 Linear Methods</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 Linear Methods</h1> <introduction> <p>This numerical tour introduces some basics about 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">Noise Distributions</a></li> <li><a href="#13">Noise in Signal and Image</a></li> <li><a href="#18">Linear image denoising</a></li> <li><a href="#22">Wiener filtering</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>Noise Distributions<a name="8"></a></h2> <p>Image and signal processing usually assumes that an random noise vector is added to the data. The simplest noises distribution are uniform and Gaussian. </p> <p>A bounded noise is simulated with the function <tt>rand</tt>. A Gaussian noise is simulated with the function <tt>randn</tt>. We show here there two noises in 1D. </p><pre class="codeinput">n = 64*64;<span class="comment">% uniform distribution in [-a,a] so that the variance is 1</span>a = sqrt(3);wu = 2*(rand(n,1)-.5)*a;<span class="comment">% gaussian distribution with variance 1</span>wn = randn(n,1);<span class="comment">% check the empirical variance</span>disp(strcat([<span class="string">'Empirical variances: uniform='</span> num2str(std(wu)) <span class="string">', Gaussian='</span> num2str(std(wn)) ]));</pre><pre class="codeoutput">Empirical variances: uniform=0.99876, Gaussian=1.0026</pre><p>We can display the histograms</p><pre class="codeinput">nbins = 50;[hu,tu] = hist(wu, nbins); hu = hu/sum(hu);[hn,tn] = hist(wn, nbins); hn = hn/sum(hn);clf;subplot(2,1,1);bar(tu, hu); axis([-tau tau 0 max(hu)]);subplot(2,1,2);bar(tn, hn); axis([-tau tau 0 max(hn)]);</pre><img vspace="5" hspace="5" src="index_01.png"> <p>We can display the noises as 1D signals.</p><pre class="codeinput">tau = 3;clf;subplot(2,1,1);plot(wu); axis([1 n -tau tau]);subplot(2,1,2);plot(wn); axis([1 n -tau tau]);</pre><img vspace="5" hspace="5" src="index_02.png"> <p>Or as 2D images.</p><pre class="codeinput">Wu = reshape(wu, sqrt(n)*[1 1]); Wu(1) = -tau; Wu(2) = tau;Wn = reshape(wn, sqrt(n)*[1 1]); Wn(1) = -tau; Wn(2) = tau;clf;imageplot(Wu, <span class="string">'Uniform noise'</span>, 1,2,1);imageplot(Wn, <span class="string">'Gaussian noise'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_03.png"> <h2>Noise in Signal and Image<a name="13"></a></h2> <p>In these numerical tour, we simulate noisy acquisition by adding some white noise (each pixel is corrupted by adding an independant Gaussian variable). This is helpful since we know the original, clean, image, and can test and compare several algorihtms by computing the recovery error. </p> <p>We load a signal.</p><pre class="codeinput">name = <span class="string">'piece-regular'</span>;n = 1024;x0 = load_signal(name,n);x0 = rescale(x0);</pre><p>We add some noise to it.</p><pre class="codeinput">sigma = .08; <span class="comment">% noise level</span>x = x0 + sigma*randn(size(x0));clf;subplot(2,1,1);plot(x0); axis([1 n -.05 1.05]);subplot(2,1,2);plot(x); axis([1 n -.05 1.05]);</pre><img vspace="5" hspace="5" src="index_04.png"> <p>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_05.png"> <h2>Linear image denoising<a name="18"></a></h2> <p>The simplest way to denoise an image is simply to apply a linear operator. In practice, one uses a Gaussian blur, and the only parameter is the width (variance) of the filter. It is non-trivial to select this parameter, since it should account for both the variance of the noise and the power spectrum of the image. </p> <p>We use a simple gaussian blur to denoise it.</p><pre class="codeinput"><span class="comment">% we use cyclic boundary condition since it is quite faster</span>options.bound = <span class="string">'per'</span>;<span class="comment">% number of pixel of the filter</span>mu = 10;Mh = perform_blurring(M,mu,options);clf;imageplot(clamp(M), <span class="string">'Noisy'</span>, 1,2,1);imageplot(clamp(Mh), <span class="string">'Blurred'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_06.png"> <p><i>Exercice 1:</i> (the solution is <a href="../private/denoising_linear/exo1.m">exo1.m</a>) Try for various Gaussian variance to compute the denoising <tt>Mh</tt>. Compute, in an oracle manner, the best variance <tt>muopt</tt> by computing the residual error <tt>snr(M0,Mh)</tt>. </p><pre class="codeinput">exo1;</pre><pre class="codeoutput">The optimal smoothing width is 2.8947 pixels, SNR=22.3067dB.</pre><img vspace="5" hspace="5" src="index_07.png"> <p>Display the results</p><pre class="codeinput"><span class="comment">% optimal filter</span>Mgauss = perform_blurring(M,muopt,options);<span class="comment">% display</span>clf;imageplot(M, strcat([<span class="string">'Noisy, SNR='</span> num2str(snr(M0,M)) <span class="string">'dB'</span>]), 1,2,1);imageplot(Mgauss, strcat([<span class="string">'Gaussian denoise, SNR='</span> num2str(snr(M0,Mgauss)) <span class="string">'dB'</span>]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_08.png"> <h2>Wiener filtering<a name="22"></a></h2> <p>In a probabilistic setting, for translation invariant signal distributions, the Wiener filtering is the optimal filtering.</p> <p>Perform the wiener filtering</p><pre class="codeinput">[Mwien,Hwien] = peform_wiener_filtering(M0,M,sigma);</pre><p>display the filter</p><pre class="codeinput">k = 5;clf;imageplot(Hwien(n/2-k+2:n/2+k,n/2-k+2:n/2+k), <span class="string">'Wiener filter (zoom)'</span>);</pre><img vspace="5" hspace="5" src="index_09.png"> <p>display the result</p><pre class="codeinput"><span class="comment">% display</span>clf;imageplot( clamp(Mgauss), strcat([<span class="string">'Gaussian denoise, SNR='</span> num2str(snr(M0,Mgauss)) <span class="string">'dB'</span>]), 1,2,1);imageplot( clamp(Mwien), strcat([<span class="string">'Wiener denoise, SNR='</span> num2str(snr(M0,Mwien)) <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 #####%% Image Denoising with Linear Methods% This numerical tour introduces some basics about 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/');%% Noise Distributions% Image and signal processing usually assumes that an random noise vector% is added to the data. The simplest noises distribution are uniform and% Gaussian.%%% A bounded noise is simulated with the function |rand|.% A Gaussian noise is simulated with the function |randn|.% We show here there two noises in 1D.n = 64*64;% uniform distribution in [-a,a] so that the variance is 1a = sqrt(3);wu = 2*(rand(n,1)-.5)*a;% gaussian distribution with variance 1wn = randn(n,1);% check the empirical variancedisp(strcat(['Empirical variances: uniform=' num2str(std(wu)) ', Gaussian=' num2str(std(wn)) ]));%%% We can display the histogramsnbins = 50;[hu,tu] = hist(wu, nbins); hu = hu/sum(hu);[hn,tn] = hist(wn, nbins); hn = hn/sum(hn);clf;subplot(2,1,1);bar(tu, hu); axis([-tau tau 0 max(hu)]);subplot(2,1,2);bar(tn, hn); axis([-tau tau 0 max(hn)]);%%% We can display the noises as 1D signals.tau = 3;clf;subplot(2,1,1);plot(wu); axis([1 n -tau tau]);subplot(2,1,2);plot(wn); axis([1 n -tau tau]);%%% Or as 2D images.Wu = reshape(wu, sqrt(n)*[1 1]); Wu(1) = -tau; Wu(2) = tau;Wn = reshape(wn, sqrt(n)*[1 1]); Wn(1) = -tau; Wn(2) = tau;clf;imageplot(Wu, 'Uniform noise', 1,2,1);imageplot(Wn, 'Gaussian noise', 1,2,2);%% Noise in Signal and Image% In these numerical tour, we simulate noisy acquisition by adding some% white noise (each pixel is corrupted by adding an independant Gaussian% variable). This is helpful since we know the original, clean, image, and% can test and compare several algorihtms by computing the recovery error.%% % We load a signal.name = 'piece-regular';n = 1024;x0 = load_signal(name,n);x0 = rescale(x0);%%% We add some noise to it.sigma = .08; % noise levelx = x0 + sigma*randn(size(x0));clf;subplot(2,1,1);plot(x0); axis([1 n -.05 1.05]);subplot(2,1,2);plot(x); axis([1 n -.05 1.05]);%%% We load an image.name = 'boat';n = 256;M0 = load_image(name,n);M0 = rescale( M0, .05, .95 );%%% Then we add some gaussian noise to it.sigma = .08; % noise levelM = M0 + sigma*randn(size(M0));clf;imageplot(M0, 'Original', 1,2,1);imageplot(clamp(M), 'Noisy', 1,2,2);%% Linear image denoising% The simplest way to denoise an image is simply to apply a linear% operator. In practice, one uses a Gaussian blur, and the only parameter% is the width (variance) of the filter. It is non-trivial to select this% parameter, since it should account for both the variance of the noise and% the power spectrum of the image. %%% We use a simple gaussian blur to denoise it.% we use cyclic boundary condition since it is quite fasteroptions.bound = 'per';% number of pixel of the filtermu = 10;Mh = perform_blurring(M,mu,options);clf;imageplot(clamp(M), 'Noisy', 1,2,1);imageplot(clamp(Mh), 'Blurred', 1,2,2);%%% _Exercice 1:_ (the solution is <../private/denoising_linear/exo1.m exo1.m>)% Try for various Gaussian variance to compute the denoising |Mh|.% Compute, in an oracle manner, the best variance |muopt| by computing the% residual error |snr(M0,Mh)|.exo1;%%% Display the results% optimal filterMgauss = perform_blurring(M,muopt,options);% displayclf;imageplot(M, strcat(['Noisy, SNR=' num2str(snr(M0,M)) 'dB']), 1,2,1);imageplot(Mgauss, strcat(['Gaussian denoise, SNR=' num2str(snr(M0,Mgauss)) 'dB']), 1,2,2);%% Wiener filtering% In a probabilistic setting, for translation invariant signal% distributions, the Wiener filtering is the optimal filtering.%%% Perform the wiener filtering[Mwien,Hwien] = peform_wiener_filtering(M0,M,sigma);%% % display the filterk = 5;clf;imageplot(Hwien(n/2-k+2:n/2+k,n/2-k+2:n/2+k), 'Wiener filter (zoom)');%%% display the result% displayclf;imageplot( clamp(Mgauss), strcat(['Gaussian denoise, SNR=' num2str(snr(M0,Mgauss)) 'dB']), 1,2,1);imageplot( clamp(Mwien), strcat(['Wiener denoise, SNR=' num2str(snr(M0,Mwien)) 'dB']), 1,2,2);##### SOURCE END #####--> </body></html>?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -