% [mg ag marker] = gradient(sd1, im, marker)
% This function returns two matrices the size of 'im' that
% contain the magnitude and the angle of the intensity
% gradient of 'im' using a gaussian directional derivative
% filter of standard deviation 'sd1'.
% If given, the optional 'marker' matrix will be set to sd1
% everywhere the gradient magnitude exceeds the critical
% value function of sd1 and left alone elsewhere.  This
% matrix will then be returned.  If no marker matrix is
% given as an argument, the function won't return the third
% value.
function [mg, ag, marker, crit] = gradient(scale, im, marker)

if (nargin > 3)
  error('Wrong number of arguments to gradient.');
end
if (nargin > 2)
  if (size(marker) ~= size(im))
    error('Marker image not the same size as im.');
  end
else
  if (nargout > 2)
    error('No marker matrix can be returned unless you supply one.');
  end
end

% We want 2 sd on either side of the filter.
% this means that in all, our filter is 
%  4*scale by 4*scale
tail = ceil(2 * scale);

% construct the filter.
[X Y] = meshgrid([-tail:tail]);
gx = g1x(X, Y, scale);
gy = g1y(X, Y, scale);
%convolve with the filter.
% this becomes very expensive with large sized filters.
gimx = conv2(im, gx, 'same');
gimy = conv2(im, gy, 'same');
ag = angle(gimx + i*gimy);
mg = steer1(ag, gimx, gimy);

abmag = abs(mg);
%compute the critical threshold for the gaussian of this scale.
% then threshold the magnitude by this value.
crit = c1(scale);
if (nargin > 2)
  list = find(abmag >= crit);
  marker(list) = scale;
end

% output magnitude and angle of only those points 
%   with absolute magnitude of gaussian greater than
%   the critical threshold.
list = find(abmag < crit);
mg(list) = 0;
ag(list) = 0;





