📄 ricedemo.html
字号:
<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>Rice/Rician Distribution Demo</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2007-03-12"> <meta name="m-file" content="ricedemo"><style>body { background-color: white; margin:10px;}h1 { color: #990000; font-size: x-large;}h2 { color: #990000; font-size: medium;}/* Make the text shrink to fit narrow windows, but not stretch too far in wide windows. On Gecko-based browsers, the shrink-to-fit doesn't work. */ p,h1,h2,div.content div { /* for MATLAB's browser */ width: 600px; /* for Mozilla, but the "width" tag overrides it anyway */ max-width: 600px; /* for IE */ width:expression(document.body.clientWidth > 620 ? "600px": "auto" );}pre.codeinput { background: #EEEEEE; padding: 10px;}span.keyword {color: #0000FF}span.comment {color: #228B22}span.string {color: #A020F0}span.untermstring {color: #B20000}span.syscmd {color: #B28C00}pre.codeoutput { color: #666666; padding: 10px;}pre.error { color: red;}p.footer { text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;} </style></head> <body> <div class="content"> <h1>Rice/Rician Distribution Demo</h1> <introduction> <p>Demonstrate ricernd, ricepdf, and ricestat, in the context of simulating Rician distributed noise for Magnetic Resonance Imaging magnitude data. </p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">Example: Rician noise vs (additive) Gaussian noise</a></li> <li><a href="#2">The Rician PDF</a></li> <li><a href="#3">Approximations to the Rician PDF</a></li> <li><a href="#6">A pitfal: Adding approximate Rician noise</a></li> </ul> </div> <h2>Example: Rician noise vs (additive) Gaussian noise<a name="1"></a></h2> <p>"Rician noise" depends on the data itself, it is not additive, so to "add" Rician noise to data, what we really mean is make the data Rician distributed. </p><pre class="codeinput">clear, close(<span class="string">'all'</span>), clcdata = load(<span class="string">'mri'</span>); <span class="comment">% MRI dataset (magnitude image) distributed with MATLAB</span>im = double(data.D(:,:,1,16)); <span class="comment">% (particular slice through the brain)</span>s=5; <span class="comment">% noise level (NB actual Rician stdev depends on signal, see ricestat)</span>im_g = im + 5 * randn(size(im)); <span class="comment">% *Add* Gaussian noise</span>im_r = ricernd(im, s); <span class="comment">% "Add" Rician noise (make Rician distributed)</span><span class="comment">% Note that the Gaussian noise has made some data go negative.</span>min_o = round(min(im(:))); max_o = round(max(im(:)));min_g = round(min(im_g(:))); max_g = round(max(im_g(:)));min_r = round(min(im_r(:))); max_r = round(max(im_r(:)));<span class="comment">% Show each image with the same color scaling limits</span>clim = [min_g max(max_g, max_r)];figure(<span class="string">'Position'</span>, [30 500 800 300]); colormap(data.map);subplot(1,3,1); imagesc(im, clim); axis <span class="string">image</span>;title(<span class="string">'Original'</span>);xlabel([<span class="string">'Range: ('</span> num2str(min_o) <span class="string">', '</span> num2str(max_o) <span class="string">')'</span>])subplot(1,3,2); imagesc(im_g, clim); axis <span class="string">image</span>;title(<span class="string">'Gaussian noise'</span>);xlabel([<span class="string">'Range: ('</span> num2str(min_g) <span class="string">', '</span> num2str(max_g) <span class="string">')'</span>])subplot(1,3,3); imagesc(im_r, clim); axis <span class="string">image</span>;title(<span class="string">'Rician noise'</span>)xlabel([<span class="string">'Range: ('</span> num2str(min_r) <span class="string">', '</span> num2str(max_r) <span class="string">')'</span>])</pre><img vspace="5" hspace="5" src="ricedemo_01.png"> <h2>The Rician PDF<a name="2"></a></h2><pre class="codeinput">x = linspace(0, 8, 100);figure;subplot(3, 1, 1)plot(x, ricepdf(x, 0, 1), x, ricepdf(x, 1, 1),<span class="keyword">...</span> x, ricepdf(x, 2, 1), x, ricepdf(x, 4, 1))title(<span class="string">'Rice PDF with s = 1'</span>)legend(<span class="string">'v=0'</span>, <span class="string">'v=1'</span>, <span class="string">'v=2'</span>, <span class="string">'v=4'</span>)subplot(3,1,2)plot(x, ricepdf(x, 1, 0.25), x, ricepdf(x, 1, 0.50),<span class="keyword">...</span> x, ricepdf(x, 1, 1.00), x, ricepdf(x, 1, 2.00))title(<span class="string">'Rice PDF with v = 1'</span>)legend(<span class="string">'s=0.25'</span>, <span class="string">'s=0.50'</span>, <span class="string">'s=1.00'</span>, <span class="string">'s=2.00'</span>)subplot(3,1,3)plot(x, ricepdf(x, 0.5, 0.25), x, ricepdf(x, 1, 0.50),<span class="keyword">...</span> x, ricepdf(x, 2.0, 1.00), x, ricepdf(x, 4, 2.00))title(<span class="string">'Rice PDF with v/s = 2'</span>)legend(<span class="string">'s=0.25'</span>, <span class="string">'s=0.50'</span>, <span class="string">'s=1.00'</span>, <span class="string">'s=2.00'</span>)</pre><img vspace="5" hspace="5" src="ricedemo_02.png"> <h2>Approximations to the Rician PDF<a name="3"></a></h2> <p>At high SNR (e.g. v/s > 3), Rician data is approximately Gaussian, though note that v is not the mean of the Rician distribution, nor of the Gaussian approximation (though for very high SNR it tends towards it). </p><pre class="codeinput">v = 1.75; s = 0.5; x = linspace(0, 4, 100);mu1 = v;gpdf1 = (2*pi*s^2)^(-0.5) * exp(-0.5*(x-mu1).^2/s^2);mu2 = ricestat(v, s); <span class="comment">% (mu2 is approximately sqrt(v^2 + s^2))</span>gpdf2 = (2*pi*s^2)^(-0.5) * exp(-0.5*(x-mu2).^2/s^2);figure; plot(x, ricepdf(x, v, s), x, gpdf1, x, gpdf2);title(<span class="string">'High SNR'</span>);legend(<span class="string">'Rician'</span>, <span class="string">'Gaussian, naive mean'</span>, <span class="string">'Gaussian, corrected mean'</span>)</pre><img vspace="5" hspace="5" src="ricedemo_03.png"> <p>At very low SNR, Rician data is approximately Rayleigh distributed (with no signal, the distributions are exactly equivalent)</p><pre class="codeinput">v = 0.25; s = 0.5; x = linspace(0, 3, 100);raypdf = (x / s^2) .* exp(-x.^2 / (2*s^2));figure; plot(x, ricepdf(x, v, s), x, raypdf);title(<span class="string">'Low SNR'</span>); legend(<span class="string">'Rician'</span>, <span class="string">'Rayleigh'</span>)</pre><img vspace="5" hspace="5" src="ricedemo_04.png"> <p>At low-medium SNR, neither Gaussian nor Rayleigh is a great approximation</p><pre class="codeinput">v = 0.75; s = 0.5; x = linspace(-1, 3, 100);mu = ricestat(v, s);gpdf = (2*pi*s^2)^(-0.5) * exp(-0.5*(x-mu).^2/s^2);raypdf = (x / s^2) .* exp(-x.^2 / (2*s^2)); raypdf(x <= 0) = 0;figure; plot(x, ricepdf(x, v, s), x, gpdf, x, raypdf);title(<span class="string">'Medium SNR'</span>); legend(<span class="string">'Rician'</span>, <span class="string">'Gaussian'</span>, <span class="string">'Rayleigh'</span>)</pre><img vspace="5" hspace="5" src="ricedemo_05.png"> <h2>A pitfal: Adding approximate Rician noise<a name="6"></a></h2> <p>Since the Rician distribution with zero signal is equivalent to the Rayleigh, and with high SNR is approximated by a Gaussian, it is tempting to add Rayleigh or Gaussian noise (depending on SNR) to existing data, to simulate Rician distributed data (aka "adding Rician noise"). </p><pre class="codeinput">v = 75; s = 50; <span class="comment">% medium SNR => Add Rayleigh noise</span><span class="comment">% Generate Rayleigh samples</span>N = 10000;gsamp = s * randn(2, N);rsamp = sqrt(sum(gsamp .^ 2));</pre><p>Compare data with added Rayleigh samples to expected Rician distribution</p><pre class="codeinput">r = v + rsamp; <span class="comment">% Add Rayleigh noise to (constant) signal</span>c = linspace(0, 250, 20);w = c(2); <span class="comment">% histogram bin-width</span>h = histc(r, c);figure; hold <span class="string">on</span>; bar(c, h, <span class="string">'histc'</span>);xl = xlim; x = linspace(xl(1), xl(2), 100);yl = ylim; <span class="comment">% (for second plot)</span>plot(x, N*w*ricepdf(x, v, s), <span class="string">'r'</span>);title(<span class="string">'Histogram of data, vs desired Rician PDF'</span>)legend(<span class="string">'Data with added Rayleigh noise'</span>, <span class="string">'Rician distribution'</span>)</pre><img vspace="5" hspace="5" src="ricedemo_06.png"> <p>Compare Rician samples to expected Rician distribution</p><pre class="codeinput">r = ricernd(v*ones(1, N), s);c = linspace(0, 250, 20);w = c(2); <span class="comment">% histogram bin-width</span>h = histc(r, c);figure; hold <span class="string">on</span>; bar(c, h, <span class="string">'histc'</span>);xl = xlim; x = linspace(xl(1), xl(2), 100);ylim(yl);plot(x, N*w*ricepdf(x, v, s), <span class="string">'r'</span>);title(<span class="string">'Histogram of data, vs desired Rician PDF'</span>)legend(<span class="string">'Rician sampled data'</span>, <span class="string">'Rician distribution'</span>)</pre><img vspace="5" hspace="5" src="ricedemo_07.png"> <p>We see that the additive Rayleigh result is very poor (with medium SNR).</p> <p>With high noise levels the appearance of results can differ dramatically:</p><pre class="codeinput">s = 0.15 * max_o;<span class="comment">% Approximate additive Rayleigh/Gaussian noise:</span>im_r = im;⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -