function [lap_of_gau, marker] = laplace(sd2, im, gau_ang, marker)
% [lap_of_gau marker] = laplace(sd2, im, gau_ang, marker)
% This function returns a matrix representing the laplacian of
% 'im' at the angles given by gau_ang.
% If given, the optional 'marker' matrix will be set to sd2
% everywhere the gradient magnitude exceeds the critical
% value function of sd2 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.

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

% We want 2 sd on either side.
tail = ceil(2 * sd2);

[X Y] = meshgrid([-tail:tail]);
gx = g2x(X, Y, sd2);
gy = g2y(X, Y, sd2);
gxy= g2xy(X, Y, sd2);
gimx = conv2(im, gx, 'same');
gimy = conv2(im, gy, 'same');
gimxy= conv2(im, gxy, 'same');
lap_of_gau = steer2(gau_ang, gimx, gimy, gimxy);

ablog = abs(lap_of_gau);
crit = c2(sd2);
if (nargin > 3)
  list = find(ablog >= crit);
  marker(list) = sd2;
end

list = find(ablog < crit);
lap_of_gau(list) = 0;


