function [alpha, sigma2, delta]=pawgn_noise_est(x, blkSize)
% 
% Estimate gain, offset and variance in image data with mixed Poisson-AWGN noise
% 
% SYNTAX:
% 
% [alpha, sigma2, delta]=pawgn_noise_est(x)
% 
% [alpha, sigma2, delta]=pawgn_noise_est(x, blkSize)
% 
% MANDATORY INPUT:
% 
% x       : vector or matrix data (double precision)
% 
% OPTIONAL INPUT:
% 
% blkSize : scalar, integer (default 8)
% 
% OUTPUTS:
% 
% alpha   : the gain
% 
% sigma2  : variance
% 
% delta   : offset
% 
% See also batch_purelet_denoise, purelet_denoise, purelet_theta
% 
% References:
% 
% 
% Luisier, F., T. Blu, and M. Unser. 2011. Image denoising in mixed
% Poisson-Gaussian noise. IEEE Trans. Image Process. 20:696?708.
% 
% Boulanger, J., et al. 2010. Patch-based nonlocal functional for denoising
% fluorescence microscopy image sequences. IEEE Trans. Med. Imag. 29:442?454.
% 
% Copyright 2011-2013 by Cezar M. Tigaret <Cezar.Tigaret@bristol.ac.uk>
% 
 
%{
   ///***************************************************************************
   // *   Copyright (C) 2011-2013 by Cezar M. Tigaret <Cezar.Tigaret@bristol.ac.uk>  *
   // *                                                                         *
   // *                                                                         *
   // *   This program is free software; you can redistribute it and/or modify  *
   // *   it under the terms of the GNU General Public License as published by  *
   // *   the Free Software Foundation; either version 3 of the License, or     *
   // *   (at your option) any later version.                                   *
   // *                                                                         *
   // *   This program is distributed in the hope that it will be useful,       *
   // *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
   // *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
   // *   GNU General Public License for more details.                          *
   // *                                                                         *
   // *   You should have received a copy of the GNU General Public License     *
   // *   along with this program; if not, write to the                         *
   // *   Free Software Foundation, Inc.,                                       *
   // *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
   // ***************************************************************************/
%}

  
  
  if(~isnumeric(x))
    error(generatemsgid('BadArgument'),'Numeric data expected');
  end;
 
  if(ndims(x)>2)
    error(generatemsgid('BadArgument'),'Can only operate on vectors and matrices');
  end;
  
  if(nargin<2 || isempty(blkSize))
    blkSize=8;
  end;
  
  if(isvector(x))
    x=x(:); % prefer column vectors
  end;
  
  sz=size(x);
  
  % number of blkSize by blkSize blocks
  nblocks=fix(sz./blkSize);
  
  if(any(rem(sz,blkSize)>0))
    newSize = nblocks.*blkSize;
    if(isvector(x))
      x = x(1:newSize(1));
    else
      x = x(1:newSize(1),1:newSize(2));
    end;
  end;
    
  % note: 
  % we calculate var_x as the standard deviation, then we square it below
  if(isvector(x))
    xblocks=reshape(x,[blkSize,nblocks(1)]);
    mu_x  = nanmean(xblocks); mu_x  = mu_x(:);
    var_x = nanstd(xblocks);  var_x = var_x(:);
  else
    xblocks=reshape(x,[blkSize,blkSize,prod(nblocks)]);
    mu_x = NaN(prod(nblocks),1);
    var_x=mu_x;
    for k=1:prod(nblocks)
      xblock = xblocks(:,:,k); 
      mu_x(k)  = nanmean(xblock(:));
      var_x(k) = nanstd(xblock(:));
    end;

  end;
  
  var_x = var_x.^2; % now var_x is what it is supposed to be: variance of x
  
  
  % force it through the origin: neither slope nor offset can be negative for
  % image data !
  fo_ = fitoptions('method','LinearLeastSquares','Lower',[0,0]);

  ok_ = isfinite(mu_x) & isfinite(var_x);
  
  if ~all( ok_ )
    warning(generatemsgid('NaN&Inf'),'Ignoring NaNs and Infs in data' );
  end
  
  ft_ = fittype('poly1');

  % Fit this model using new data
  fitObj = fit(mu_x(ok_),var_x(ok_),ft_,fo_);
 
  cf=coeffvalues(fitObj);
  
  alpha=cf(1);
  
  if(isnan(alpha))
    alpha=1;
  end;
                
  delta = nanmean(x(:)); % delta = sample mean of background (i.e. of AWGN)
  
%   beta=cf(2); % beta = sigma2 - delta*alpha; 
  
  sigma2 = cf(2) + alpha.*delta;
  
end % function pawgn_noise_est