%GAUSSIANBLOB  Adds Gauss shaped spots to an  image
%
% SYNOPSIS:
%  image_out = gaussianblob(image_in, coordinates, sigma, strength, domain, truncation)
%
%  image_in:    input image, where the spots will be set
%  image_out:   image_in + spots
%  domain:      'spatial' or 'frequency'
%    SPATIAL: the image is in the spatial domain
%      coordinates: array containing the location of the spots 
%                   range:(0,image size), array size must be form [Nx2] for fast 2D placement
%      sigma:       the sigmas of the gaussians
%      strength:    the integrated intensity of the spots
%      spot = strength * (2pi sigma^2)^(-imagedimension/2) * exp(-1/2 (x/sigma)^2)
%    FREQUENCY: the image is in the frequency domain
%      coordinates: array containing the location of the spot in frequency domain
%                   units, which are in the range (-0.5,0.5)
%      sigma:       the sigma of the gauss (in spatial domain units)
%      strength:    the integrated intensity of the spot in the spatial domain.
%      spot = strength / (sqrt(sum(pixels))) * exp(-1/2 (sigma 2 pi f)^2)
%  truncation:  gives the region in which the blob is computed. For the spatial
%    domain, truncation*sigma is the radius of the region being computed. For the
%    frequency domain, truncation/(2*pi*sigma) is the radius.
%
% DEFAULTS:
%  sigma = 2
%  strength = 255
%  domain = 'spatial'
%  truncation = 3

% (C) Copyright 1999-2012               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Bernd Rieger, Apr 2001
% 30 April 2004: Added 'truncation' parameter. (CL)
% 5 May 2004   : Output bug fix (BR)
% November 2005: Added possibility of anisotropic sigma (BR)
% October 2006:  Improved speed for large images and small blobs. (CL)
% 5 April 2007:  Fixed bug for non-integer coordinates.
%                Fixed bug for normalization in spatial coordinates.
%                Fixed bug for truncation in frequency domain. (CL)
% 1 Feb 2013:    The option 'freq' was called 'frequency' in documentation. Fixed. (CL)
% Aug 2013:      Added possibility to add _fast_ multiple Gaussians with one call (BR)

function in = gaussianblob(varargin)

d = struct('menu','Generation',...
           'display','Add Gaussian blob',...
           'inparams',struct('name',       {'in',         'coord',      'sigma','strength','domain','truncation'},...
                             'description',{'Input image','Coordinates','Sigma','Strength','Domain','Truncation'},...
                             'type',       {'image',      'array',      'array','array',   'option','array'},...
                             'dim_check',  {0,            [-1,2],            -1,      -1,         0,       0},...
                             'range_check',{[],           'R',          'R+',   'R',{'spatial','frequency'},'R+'},...
                             'required',   {1,            1,            0,      0,         0,       0},...
                             'default',    {'a',          '[128 128]',  2,      255,      'spatial',3}...
                            ),...
           'outparams',struct('name',{'in'},...
                              'description',{'Output image'},...
                              'type',{'image'}...
                              )...
         );
if nargin == 1
   s = varargin{1};
   if ischar(s) & strcmp(s,'DIP_GetParamList')
      in = d;
      return
   end
end
if nargin>=5
   if isequal(varargin{5},'freq')
      varargin{5} = 'frequency';
   end
end
try
   [in,coord,sigma,strength,domain,truncation] = getparams(d,varargin{:});
catch
   if ~isempty(paramerror)
      error(paramerror)
   else
      error(firsterr)
   end
end

sz = size(coord);
% route via dipimage object: for one blob and Fourier space option
if strcmp(domain,'frequency') || ndims(coord)~=2 || ~all(sz>1)
   di = ndims(in);
   if length(sigma)<di; sigma = repmat(sigma,1,di); end
   mask = newim(in,'bin');
   switch domain
      case 'spatial'
         normalisation = prod(sigma) * (sqrt(2*pi))^di;
         sigma = 1./sigma;
         origin = 'corner';
      case 'frequency'
         normalisation = sqrt(prod(size(in)));
         coord = 2*pi*coord;
         origin = 'radfreq';
         % make sure ift(frequency)=spatial therefore normalisation 1/sqrt(N)
      otherwise
         error('Unkown domain.');
   end
   if isinf(truncation)
      mask(:) = 1;
   else
      s = substruct('()',{});
      for ii = 1:di
         I = coord(ii)+(truncation/sigma(ii)*[-1,1]);
         if strcmp(domain,'frequency')
            I = (I+pi)*(size(in,ii)/(2*pi));
         end
         I = [floor(I(1)),ceil(I(2))];
         I(1) = max(I(1),0);
         I(2) = min(I(2),size(in,ii)-1);
         s.subs{ii} = I(1):I(2);
      end
      mask = subsasgn(mask,s,1);
   end

   exponent = newim(in);
   for ii = 1:di
      ramp = createramp(size(in),ii,origin);
      exponent(mask) = exponent(mask) + ((ramp(mask)-coord(ii)).*sigma(ii))^2;
   end
   in(mask) = in(mask) + (strength/normalisation) * exp(-0.5*exponent(mask));

else
%route via matlab double object for multiple spatial coordinates in 2D (code can be extended to 3D)
   if ndims(coord)~=2
       error('input coordinates must be 2D');
   end
   sz = size(coord);
   if length(sigma)==1
       sigma = repmat(sigma,sz(1),1);
   end
   if length(strength)==1
       strength = repmat(strength,sz(1),1);
   end
   if sz(1)~=length(sigma)
       error('Sigmas and coordinates sizes do not match.')
   end
   if sz(1)~=length(strength)
       error('Strength and coordinates sizes do not match.')
   end

   b = double(in);
   Asz = size(b);
   coT = coord-floor(coord);
   coP = floor(coord);
   %for double the coordiantes range from 1:N, for a dipimage from 0:N-1
   for ii=1:sz(1)
     gsz = ceil(truncation*2*sigma(ii));
     if ~mod(gsz,2); gsz=gsz+1; end %make odd size
     if gsz<3; gsz=3; end %for very small sigma, otherwise not normalized
     [y,x] = meshgrid(-floor(gsz/2):floor(gsz/2), -floor(gsz/2):floor(gsz/2) );
     % not taken into account the pixellation
     %z = strength(ii)/(2*pi*sigma(ii)^2)*exp( -( (x-coT(ii,1)).^2 + (y-coT(ii,2)).^2)/(2*sigma(ii)^2) );
     
     s2 = sqrt(2)*sigma(ii); %only allow isotrop Gaussian
     z =  strength(ii)/4.*((erf((x-coT(ii,1)+.5)./s2) - ...
     erf((x-coT(ii,1)-.5)./s2)).* ...
        (erf((y-coT(ii,2)+.5)./s2)- ...
        erf((y-coT(ii,2)-.5)./s2)));

     id1 = coP(ii,:)-floor(gsz/2)+1;
     id2 = coP(ii,:)+floor(gsz/2)+1;
     if any(id1<=0) || any(id2>=Asz) %remove all spots that are touching the border, could add them cut if needed?
         continue;
     end
     b(id1(1):id2(1) , id1(2):id2(2)) = b(id1(1):id2(1) , id1(2):id2(2)) + z;
   end
   in = mat2im(b');
end
