gen_susan(I,w,OPS)

function J = gen_susan(I,w,OPS)
% J = gen_susan(I,regsize,sigma,tau)
% 
% GEN_SUSAN - Generalized SUSAN 2-D filtering
% SUSAN filtering with filter kernel W scaled with generalized
% Gaussian of intensity difference. Different prefiltering
% functions can be selected as well as width and exponent of the
% intensity scaling. GEN_SUSAN can produce filtering with
% caracteristics similar to wiener2 and medfilt2 and "everything
% inbetween"
% 
% INPUT:
% I - 2-D matrix,
% W - Filter kernel, 
%       Default: ones(3)/9 
% OPS - options struct with fields:
%   TAU - width of intensity weighting. 1/eps -> filtering similar to
%         linear filter, std(I(:)) and I give results similar to
%         sigma-filter/wiener2. 
%         Default: 1, more useful values have to be set manually,
%         ex: I.^0.5
%   GAMMA - shape factor, 2 <-> Gaussian, 1 <-> double sided
%         exponential, Wi(i,j) = W(i,j)*exp(-abs((Ip-I(i,j))/TAU).^GAMMA)
%         Default: 2
%   PRE_FILTER - [{'n'}|'f'|'w'|'m'|'s'], 'n' for none, 'f' for
%         FILTER2 (using same filter kernel), 'w' for
%         WIENER2 (using size(W) as neighbourhood region),
%         'm' for MEDFILT2, 's' for GEN_SUSAN, using same
%         filter kernel but others options default.
%   NO_CENTER - [{0}|1] to include the center point or not to
%         include the center point, default is to include it.
% 
% OUTPUT:
%   J - Filtered image, always of type DOUBLE; ro
%       with no input arguments GEN_SUSAN returns the default
%       OPS-struct. 
% 
% GEN_SUSAN is a generalization of "Smoothing over Univalue Segment
% Assimilating Nucleus" S.M. Smith and J.M. Brady. SUSAN - a new
% approach to low level image processing.  Int. Journal of Computer
% Vision, 23(1):45--78, May 1997.
% 
% No error checking.
%
% See also: WIENER2, MEDFILT2, FILTER2,

% Copyright Bjorn Gustavsson 20050129


% If no input arguments give the OPS-struct with default options.
if nargin == 0
  
  J.gamma = 2;
  J.pre_filter = 'n';
  J.tau = 1;
  J.no_center = 0;
  return
  
end

% Default filter kernel
if nargin == 1 | isempty(w)
  w = ones(3)/9;
end

% If no OPS given creat a default one.
if nargin <= 2
  OPS = gen_susan;
else  % If there are an OPS-struct complete it with defaults for
      % missing fields.
  if ~isfield(OPS,'gamma')
    OPS.gamma = 2;
  end
  if ~isfield(OPS,'tau')
    OPS.tau = 1;  % Pretty useless value but there has to be something...
  end
  if ~isfield(OPS,'pre_filter')
    OPS.pre_filter = 'n';
  end
  if ~isfield(OPS,'no_center')
    OPS.no_center = 0;
  end
end

if ~isa(I, 'double')
  I = double(I);
end

regsize = size(w);

% Constant extrapolation at the edges. Higher order extrap too
% complicated and since there is a need for filtering there should
% be some noise and that makes higher order extrapolation even worse?
I_internal = I([ones(1,ceil((regsize(1)-1)/2)) 1:end size(I,1)*ones(1,(regsize(1)-1)/2)],...
               [ones(1,ceil((regsize(2)-1)/2)) 1:end size(I,1)*ones(1,(regsize(2)-1)/2)]);

% Pre-filtering
switch OPS.pre_filter
 case 'f'
  I_internal = filter2(w,I_internal);
 case 'w'
  try 
    I_internal = wiener2(I_internal,size(w));
  catch
    warning('susan:nonexistinprefilter',...
            'Could not find WIENER2 for SUSAN prefiltering')
  end
 case 'm'
  try 
    I_internal = medfilt2(I_internal,size(w));
  catch
    warning('susan:nonexistinprefilter',...
            'Could not find MEDFILT2 for SUSAN prefiltering')
  end
 case 's'
    I_internal = gen_susan(I_internal,w);
 otherwise
  % sit quiet and happy
end

[Iy,Ix] = size(I);

% Initialization
J = 0;
S_W = 0;


%% Filtering
%          sy,sx
%          ___
%         \
% J(y,x) = >  I(y+i-sy/2,x+j-sx/2)*W(i,j)*exp(-|Ip(y,x)-I(y+i-sy/2,x+j-sx/2)|^gamma/tau^gamma)
%         /___
%
%         i,j = 1,1
%        -------------------------------------------------------------------------------------
%          sy,sx
%          ___
%         \
%          > W(i,j)*exp(-|Ip(y,x)-I(y+i-sy/2,x+j-sx/2)|^gamma/tau^gamma)
%         /___
%
%         i,j = 1,1
%
%
for i = 1:regsize(1),
  
  for j = 1:regsize(2),
    
    if OPS.no_center & abs(i-(regsize(1)+1)/2)<eps  & abs(j-(regsize(2)+1)/2)<eps 
      % then do not include that point (center point)
    else
      W = w(1+end-i,1+end-j)*exp( - abs((I-I_internal(i:(i-1+Iy),j:(j-1+Ix)))./OPS.tau).^OPS.gamma );
      J = J+I_internal(i:(i-1+Iy),j:(j-1+Ix)).*W;
      S_W = S_W + W;
    end
  end
  
end

J = J./S_W;