📄 index.html

📁 信号处理系列导航
💻 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>Outliers and Median Denoiser</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>Outliers and Median Denoiser</h1>         <introduction>            <p>This numerical tour explores non-linear denoisers based on local median computations.</p>         </introduction>         <h2>Contents</h2>         <div>            <ul>               <li><a href="#1">Installing toolboxes and setting up the path.</a></li>               <li><a href="#8">Robust Statistics and Median Estimator</a></li>               <li><a href="#13">Impulse Noise Removal with Local Median</a></li>               <li><a href="#20">Iterated median.</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>Robust Statistics and Median Estimator<a name="8"></a></h2>         <p>Robust statistic are estimators that can compute accurately means and standard deviations of process contaminated with outliers            (few exceptional events like spikes). This is related to replacing the L2 norm for which traditional estimator (empirical            mean and empirical deviation) are optimal by L1 norm which is robust to outliers.         </p>         <p>First we generate a white noise of mean and deviation 1 with a few outliers that have a much larger deviation and that are            positive.         </p><pre class="codeinput"><span class="comment">% dimension</span>n = 10000;<span class="comment">% deviation for base signal and outliers</span>sigma0 = 1;sigma1 = 20;<span class="comment">% ratio of outliers</span>rho = .02;<span class="comment">% number of outliers</span>p = round(rho*n);<span class="comment">% noisy base signal</span>x = randn(n,1)*sigma0;<span class="comment">% positions of the outliers</span>sel = randperm(n); sel = sel(1:p);<span class="comment">% add outliers</span>x(sel) = x(sel) + rand(p,1)*sigma1;<span class="comment">% display</span>clf;plot(x); axis(<span class="string">'tight'</span>);</pre><img vspace="5" hspace="5" src="index_01.png"> <p>We can compute the mean and the median. Note how the mean is shifted because of the positive outliers.</p><pre class="codeinput">disp(strcat([<span class="string">'Mean  ='</span> num2str(mean(x))] ));disp(strcat([<span class="string">'Median='</span> num2str(median(x))] ));</pre><pre class="codeoutput">Mean  =0.1774Median=-0.0025397</pre><p>We can compute the standard deviation and the L1 norm of absolute deviation (<tt>mad</tt>). The mad estimator is normalized by         </p><pre class="codeinput">disp(strcat([<span class="string">'Mean ='</span> num2str(std(x))] ));disp(strcat([<span class="string">'MAD  ='</span> num2str(mad(x,1)/0.6745)] ));</pre><pre class="codeoutput">Mean =1.8435MAD  =1.0329</pre><p>Another options to stabilize the estimation is to threshold the signal to detect outliers. In wavelet denoising, such a technic            is used because the outlier are in fact the data one is interested in (the signal content). The Gaussian variable is with            high probabilitiy below <tt>T=sqrt(2*log(n))</tt>, so that this universal threshold can be used. In practice, a smaller, less conservative threshold can be used.         </p><pre class="codeinput">T = sqrt(2*log(n));xT = x .* (abs(x)&lt;=T);plot(xT); axis(<span class="string">'tight'</span>);title(<span class="string">'Signal below the threshold'</span>);<span class="comment">% compute the thresholded mean and deviation</span>disp(strcat([<span class="string">'Mean ='</span> num2str(mean(xT))] ));disp(strcat([<span class="string">'Std  ='</span> num2str(std(xT))] ));</pre><pre class="codeoutput">Mean =0.0028432Std  =1.0116</pre><img vspace="5" hspace="5" src="index_02.png"> <h2>Impulse Noise Removal with Local Median<a name="13"></a></h2>         <p>A median is able to corretly estimate the mean of a signal. The mean of a smooth image is a good approximation of the image            itself. For a piecewise smooth image, this mean needs to be computed locally. Applying local mean (average) leads to blurring            and also does not works to remove impulse noise. Both issues are solved by computing a local median.         </p>         <p>Load a clean image</p><pre class="codeinput">n = 128;M0 = load_image(<span class="string">'lena'</span>, n*2);M0 = rescale(crop(M0,n));</pre><p>Parameters for the noise</p><pre class="codeinput">rho = .4;sigma = 2;</pre><p><i>Exercice 1:</i> (the solution is <a href="../private/tv_median/exo1.m">exo1.m</a>) Compute an image <tt>M</tt> corrupted by adding impulse gaussian noise to <tt>M0</tt> for a fraction <tt>rho</tt> of its pixels. For the display you need to clamp the images.         </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_03.png"> <p>For Scilab, one might need to extend a little the size of the stack.</p><pre class="codeinput"><span class="comment">% execute |extend_stack_size(4);| if needed.</span></pre><p>A local median applied over a window of <tt>w x w</tt> pixels removes the outlier of the image. The parameter <tt>w</tt> should be large enough for the median to be well estimated (the more noise, the larger <tt>w</tt>) and small enough not too blur too much the edges.         </p><pre class="codeinput">k = 2;      <span class="comment">% half window</span>w = 2*k+1;  <span class="comment">% window size</span>Mmed = perform_median_filtering(M,k);clf;imageplot(clamp(M), <span class="string">'Noisy'</span>, 1,2,1);imageplot(clamp(Mmed), <span class="string">'Median filtered'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_04.png"> <p><i>Exercice 2:</i> (the solution is <a href="../private/tv_median/exo2.m">exo2.m</a>) Try for various patch size <tt>w</tt> to find the best result.         </p><pre class="codeinput">exo2;</pre><img vspace="5" hspace="5" src="index_05.png"> <h2>Iterated median.<a name="20"></a></h2>         <p>Instead of applying the median only once, one can iteratively filter an image to simulate a time evolution. One can prove            that asymptotically, this flow converges (up to renormalizaiton) to the Total Variation (TV) flow of Rudin/Osher/Fatemi.         </p>         <p><i>Exercice 3:</i> (the solution is <a href="../private/tv_median/exo3.m">exo3.m</a>) Implemented an iterated median filtering, starting from the image <tt>M0</tt>. What do you notice. Compare this median with a total variation flow.         </p><pre class="codeinput">exo3;</pre><img vspace="5" hspace="5" src="index_06.png"> <p class="footer"><br>            Copyright  &reg; 2008 Gabriel Peyre<br></p>      </div>      <!--##### SOURCE BEGIN #####%% Outliers and Median Denoiser% This numerical tour explores non-linear denoisers based on local median% computations.%% 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/');%% Robust Statistics and Median Estimator% Robust statistic are estimators that can compute accurately means and% standard deviations of process contaminated with outliers (few exceptional% events like spikes). This is related to replacing the L2 norm for which% traditional estimator (empirical mean and empirical deviation) are% optimal by L1 norm which is robust to outliers.%%% First we generate a white noise of mean and deviation 1 with a few outliers that have a% much larger deviation and that are positive.% dimensionn = 10000;% deviation for base signal and outlierssigma0 = 1;sigma1 = 20;% ratio of outliersrho = .02; % number of outliersp = round(rho*n);% noisy base signalx = randn(n,1)*sigma0;% positions of the outlierssel = randperm(n); sel = sel(1:p);% add outliersx(sel) = x(sel) + rand(p,1)*sigma1;% displayclf;plot(x); axis('tight');%%% We can compute the mean and the median. Note how the mean is shifted% because of the positive outliers.disp(strcat(['Mean  =' num2str(mean(x))] ));disp(strcat(['Median=' num2str(median(x))] ));%%% We can compute the standard deviation and the L1 norm of absolute% deviation (|mad|). The mad estimator is normalized by disp(strcat(['Mean =' num2str(std(x))] ));disp(strcat(['MAD  =' num2str(mad(x,1)/0.6745)] ));%%% Another options to stabilize the estimation is to threshold the signal to% detect outliers. In wavelet denoising, such a technic is used% because the outlier are in fact the data one is interested in (the signal content).% The Gaussian variable is with high probabilitiy below |T=sqrt(2*log(n))|,% so that this universal threshold can be used. In practice, a smaller,% less conservative threshold can be used.T = sqrt(2*log(n));xT = x .* (abs(x)<=T);plot(xT); axis('tight'); title('Signal below the threshold');% compute the thresholded mean and deviationdisp(strcat(['Mean =' num2str(mean(xT))] ));disp(strcat(['Std  =' num2str(std(xT))] ));%% Impulse Noise Removal with Local Median% A median is able to corretly estimate the mean of a signal.% The mean of a smooth image is a good approximation of the image itself.% For a piecewise smooth image, this mean needs to be computed locally.% Applying local mean (average) leads to blurring and also does not works% to remove impulse noise. Both issues are solved by computing a local% median.%%% Load a clean imagen = 128;M0 = load_image('lena', n*2);M0 = rescale(crop(M0,n));%%% Parameters for the noiserho = .4;sigma = 2;%%% _Exercice 1:_ (the solution is <../private/tv_median/exo1.m exo1.m>)% Compute an image |M| corrupted by adding impulse gaussian noise to |M0|% for a fraction |rho| of its pixels. For the display you need to clamp% the images.exo1;%%% For Scilab, one might need to extend a little the size of the stack.% execute |extend_stack_size(4);| if needed.%%% A local median applied over a window of |w x w| pixels removes the% outlier of the image. The parameter |w| should be large enough for the% median to be well estimated (the more noise, the larger |w|) and small% enough not too blur too much the edges.k = 2;      % half windoww = 2*k+1;  % window sizeMmed = perform_median_filtering(M,k);clf;imageplot(clamp(M), 'Noisy', 1,2,1);imageplot(clamp(Mmed), 'Median filtered', 1,2,2);%%% _Exercice 2:_ (the solution is <../private/tv_median/exo2.m exo2.m>)% Try for various patch size |w| to find the best result.exo2;%% Iterated median.% Instead of applying the median only once, one can iteratively filter an% image to simulate a time evolution. One can prove that asymptotically,% this flow converges (up to renormalizaiton) to the Total Variation (TV) flow of Rudin/Osher/Fatemi.%%% _Exercice 3:_ (the solution is <../private/tv_median/exo3.m exo3.m>)% Implemented an iterated median filtering, starting from the image |M0|. What do you notice. Compare% this median with a total variation flow.exo3;##### SOURCE END #####-->   </body></html>
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -