matlab-imnoise分析及泊松噪声

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. 噪音生成（泊松噪音生成方法）