%SCALESPACE   Gaussian scale-space
%
% SYNOPSIS:
%  [sp,bp,pp] = scalespace(image_in,nscales,base)
%
% DEFAULTS:
%  nscales = 7
%  base = sqrt(2)    images are smoothed with base^i, i=0:nscales-1
%
% OUTPUT PARAMETERS:
%  sp: Scale pyramide  
%  bp: Difference between scales
%  pp: Variance between scales
%
% EXAMPLES:
%  [x,y,z] = scalespace(readim) %2D image
%  [x1,y1] = scalespace(readim('chromo3d')) %3D image 
%
% LITERATURE:
%  J.J. Koenderink, The Structure of Images, Biological Cybernetics, 50:363-370, 1984.
%  T. Lindeberg, Scale-Space for Discrete Signals, IEEE Transactions PAMI, 12(3):234-254, 1990.


% (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
%
% Peter Bakker, Oct 2000.
% Jan 2004: make nD, extended help and make faster (BR)
% March 2007: changed default base to sqrt(2)  (BR)

function [varargout] = scalespace(varargin)

d = struct('menu','Analysis',...
   'display','Gaussian scale-space',...
   'inparams',struct('name',       {'in',   'scales','base'},...
                     'description',{'Input image','Scales','Base'},...
                     'type',       {'image',      'array', 'array'},...
                     'dim_check',  {0,            0,       1},...
                     'range_check',{[],           'N+',    'R+'},...
                     'required',   {1,            0,       0},...
                     'default',    {'a',          7,       sqrt(2)}...
                    ),...
   'outparams',struct('name',       {'sp','bp','pp'},...
                      'description',{'Scale','Difference','Variance'},...
                      'type',       {'image','image','image'}...
                      )...
  );
if nargin == 1
   s = varargin{1};
   if ischar(s) & strcmp(s,'DIP_GetParamList')
      varargout{1} = d;
      return
   end
end
try
   [in,scales,base] = getparams(d,varargin{:});
catch
   if ~isempty(paramerror)
      error(paramerror)
   else
      error(firsterr)
   end
end
% Total rewrite of PB code to make nD and 
% make variable number of output arguments
% for faster processing (BR)

%make process strings
di = ndims(in);
str_sp = ['varargout{1}(' repmat(':,',1,di) 'ii)'];
str_sp2 =['varargout{1}(' repmat(':,',1,di) 'ii-1)'];
str_bp = ['varargout{2}(' repmat(':,',1,di) 'ii-1)'];
str_pp = ['varargout{3}(' repmat(':,',1,di) 'ii-1)'];

N = nargout; if N<1, N=1; end
varargout = cell(N,1);

varargout{1} = newim([size(in),scales]);
ii=0;
eval([str_sp '=in;']);

switch N
   case {0, 1}
      for ii=1 : scales-1
         sc = base.^(ii-1);
         %#function gaussf
         eval([str_sp '=gaussf(in,sc);']);
      end
   case 2
      varargout{2} = newim([size(in),(scales-1)]);
      for ii=1 : scales-1
         sc = base.^(ii-1);
         %#function gaussf
         eval([str_sp '=gaussf(in,sc);']);
         eval([str_bp '=' str_sp '-' str_sp2 ';']);
      end
   case 3
      varargout{2} = newim([size(in),(scales-1)]);
      varargout{3} = newim([size(in),(scales-1)]);
      for ii=1 : scales-1
         sc = base.^(ii-1);
         %#function gaussf squeeze
         eval([str_sp '=gaussf(in,sc);']);
         eval([str_bp '=' str_sp '-' str_sp2 ';']);
         eval([str_pp '=gaussf(squeeze(' str_bp '*' str_bp '),sc);']);
      end
   otherwise
      error('Too many output parameters.');
end
