imnoise
C:\Program Files\MATLAB\R2017a\toolbox\images\images\imnoise.m
function b = imnoise(varargin)
[a, code, classIn, classChanged, p3, p4] = ParseInputs(varargin{:});
clear varargin;
b = images.internal.algimnoise(a, code, classIn, classChanged, p3, p4);
images.internal.algimnoise
C:\Program Files\MATLAB\R2017a\toolbox\images\images+images+internal\algimnoise.m
function b = algimnoise(a, code, classIn, classChanged, p3, p4)
sizeA = size(a);
switch code
%高斯白噪声
case 'gaussian' % Gaussian white noise
b = a + sqrt(p4)*randn(sizeA) + p3;
case 'localvar_1' % Gaussian white noise with variance varying locally
% imnoise(a,'localvar',v)
% v is local variance array
b = a + sqrt(p3).*randn(sizeA); % Use a local variance array
case 'localvar_2' % Gaussian white noise with variance varying locally
% Use an empirical intensity-variance relation
intensity = p3(:);
var = p4(:);
minI = min(intensity);
maxI = max(intensity);
b = min(max(a,minI),maxI);
b = reshape(interp1(intensity,var,b(:)),sizeA);
b = a + sqrt(b).*randn(sizeA);
%泊松噪声
case 'poisson' % Poisson noise
switch classIn
case 'uint8'
a = round(a*255);
case 'uint16'
a = round(a*65535);
case 'single'
a = a * 1e6; % Recalibration
case 'double'
a = a * 1e12; % Recalibration
end
a = a(:);
% (Monte-Carlo Rejection Method) Ref. Numerical
% Recipes in C, 2nd Edition, Press, Teukolsky,
% Vetterling, Flannery (Cambridge Press)
b = zeros(size(a),'like', a);
%对于小于50的像素值,按照泊松分布生成
% 下面的泊松随机数算法参照Donald Knuth的算法
idx1 = find(a<50); % Cases where pixel intensities are less than 50 units
if (~isempty(idx1))
g = exp(-a(idx1));
em = -ones(size(g));
t = ones(size(g));
idx2 = (1:length(idx1))';
while ~isempty(idx2)
em(idx2) = em(idx2) + 1;
t(idx2) = t(idx2) .* rand(size(idx2));
idx2 = idx2(t(idx2) > g(idx2));
end
b(idx1) = em;
end
% 对于大于50的像素值,当作高斯分布
% For large pixel intensities the Poisson pdf becomes
% very similar to a Gaussian pdf of mean and of variance
% equal to the local pixel intensities. Ref. Mathematical Methods
% of Physics, 2nd Edition, Mathews, Walker (Addison Wesley)
idx1 = find(a >= 50); % Cases where pixel intensities are at least 50 units
if (~isempty(idx1))
b(idx1) = round(a(idx1) + sqrt(a(idx1)) .* randn(size(idx1)));
end
b = reshape(b,sizeA);
%椒盐噪声
case 'salt & pepper' % Salt & pepper noise
b = a;
x = rand(sizeA);
b(x < p3/2) = 0; % Minimum value
b(x >= p3/2 & x < p3) = 1; % Maximum (saturated) value
%乘性噪声
case 'speckle' % Speckle (multiplicative) noise
b = a + sqrt(12*p3)*a.*(rand(sizeA)-.5);
end
% Truncate the output array data if necessary
if strcmp(code,{'poisson'})
switch classIn
case 'uint8'
b = uint8(b);
case 'uint16'
b = uint16(b);
case 'single'
b = max(0, min(b / 1e6, 1));
case 'double'
b = max(0, min(b / 1e12, 1));
end
else
b = max(0,min(b,1));
% The output class should be the same as the input class
if classChanged,
b = images.internal.changeClass(classIn, b);
end
end
关于泊松随机的算法
见
Generating Poisson random values
数字图像学笔记——7. 噪音生成(泊松噪音生成方法)