%NOISESTD   Estimate noise standard deviation
%
% SYNOPSIS:
%  function [snoise,tmp,mask] = noisestd(in,mask,frame)
%
% Assumption: White Gaussian noise (since use neighbour differences)
% Limitation: may fail if noise correlated (e.g. result of filtering)
% fail under complex texture.
%
% PARAMETERS:
%  mask = binary regions over which noise is computed ([] means all)
%  frame = the frame over which noise is computed (for 3D input only) 
% 
% EXAMPLE: 
%  a = readim
%  b = noise(a,'gaussian',20)
%  c = hybridf(b)             % better than Lee80
%  d = dip_image(wiener2(double(b),[3 3],noisestd(b)^2))    % Lee80
%  [noisestd(a) noisestd(b) noisestd(c) noisestd(d)] 
%
% LITERATURE:
%  John Immerkær. "Fast Noise Variance Estimation", 
%  Computer Vision and Image Understanding, Vol. 64, No. 2, p. 300-302, 1996. 

% (C) Copyright 1999-2003               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Tuan Pham, December 2003.

function [snoise,tmp,mask] = noisestd(in,mask,frame)

if nargin == 1
   if ischar(in) & strcmp(in,'DIP_GetParamList')
      snoise = struct('menu','none');
      return
   end
end


if length(size(in))==3
   if nargin<3 computed_at_frame = floor(size(in,3)/2), end
   % could do separable convolution in 3D, but costly & 
   % inaccurate because sampling rate in z direction is not the same
   in = squeeze(in(:,:,computed_at_frame));
   if nargin>1&~isempty(mask) mask = squeeze(mask(:,:,frame)); end
end
   
% avoid estimating noise @ pixels near edges
if nargin<2|isempty(mask) mask = ~threshold(gaussf(gradmag(in),3)); end

if length(size(in))==1
   tmp = dip_image(conv([1 -2 1], double(in)));
else
   tmp = dip_image(conv2([1 -2 1], [1 -2 1], double(in)));
end

% three methods listed at
% http://www.mip.sdu.dk/~jimm/ipcv/noiseestimation.html
snoise = 1/6 * sqrt(mean(tmp(mask)^2));
% perform poorly in low-noise cases so discard
% s2     = 1/6 * sqrt(pi/2) * mean(abs(tmp(mask)));
% perform so poorly in high-noise cases so discard
% s3     = 1/6 * sqrt(pi/2) * median(abs(tmp(mask)));
