%GAUSSIANBLOB  Adds a Gauss shaped spot to an  image
%
% SYNOPSIS:
%  image_out = gaussianblob(image_in, coordinates, sigma, strength, domain, truncation)
%
%  image_in:    input image, where the spot will be set
%  image_out    = image_in + spot
%  domain:      'spatial' or 'frequency'
%    SPATIAL: the image is in the spatial domain
%      coordinates: array containing the location of the spot (0,image size)
%      sigma:       the sigma of the gauss
%      strength:    the integrated intensity of the spot
%      spot = strength * (2pi sigma^2)^(-imagedimension/2) * exp(-1/2 (x/sigma)^2)
%    FREQ: 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-2007               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)

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,            1,      0,         0,       0},...
                             'range_check',{[],           'R',          'R+',   'R',{'spatial','freq'},'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
try
   [in,coord,sigma,strength,domain,truncation] = getparams(d,varargin{:});
catch
   if ~isempty(paramerror)
      error(paramerror)
   else
      error(firsterr)
   end
end

di = ndims(in);
mask = newim(in,'bin');
switch domain
   case 'spatial'
      normalisation = prod(sigma) * (sqrt(2*pi))^di;
      sigma = 1./sigma;
      origin = 'corner';
   case 'freq'
      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,'freq')
         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));
