% runEZD2.m        heavily modified by:    Dion Monstavicius
%                  originally created by:  Anat Caspi
% This function runs the Elder and Zucker multi-scale
% edge detector and returns the edgel image.
%
% INPUTS:  im, an image.
%          sigma--an estimate of the sensor noise of the image
%
% RETURNS: zc:         an edge map, not necessarily optimized.
%          lap_of_gau: the 2nd derivative estimate along
%                      the gradient direction
%          scale2:     the minimum reliable scale of the
%                      second derivative estimator
%          blur:       the blur estimate at every pixel

function [zc, lap_of_gau, scale2, blur] = runEZD2(im, sigma)

DEBUG = 1;
DISP  = 1;

if(DISP)
  figure(1);
  imagesc(im);
  colormap gray;
  print -deps origim.eps
end

if (nargin <2)
  sigma = 0.05;
end

global sn;
global alphaP;

sn = sigma;
%alphaP = 1 - (1 - sigma).^(1./length(im))
alphaP = 2e-7;
% creating the different scales. These represent 
% the different scale sizes supported here, expressed
% as the standard deviation of the filter.
sd1 = [16 8 4 2 1 0.5];
szim = size(im);
map = zeros(szim);
mag = zeros(szim);
ang = zeros(szim);

% we iterate. For each standard deviation, 
% we compute the 'gradient'- which in this
% case uses the steering functions defined
% by Elder and Zucker.

for h = 1:length(sd1)
  [mg ag map crit] = gradient(sd1(h), im, map);
  mag(find(mg ~= 0)) = 0;
  mag = mag + mg;
  ang(find(ag ~= 0)) = 0;
  ang = ang + ag;
end
%scaind = log2(map) + repmat(2, size(map));


sd2 = [4 2 1 0.5];
szim = size(im);
scale2 = repmat(0, szim);
lap_of_gau = repmat(0, szim);

for h = 1:length(sd2)
  [l_o_g scale2] = laplace(sd2(h), im, ang, scale2);
  lap_of_gau(find(l_o_g ~= 0)) = 0;
  lap_of_gau = lap_of_gau + l_o_g;
end
%scaind2 = log2(scale2) + repmat(2, size(scale2));

zc = edge(lap_of_gau, 'canny');

% now compute the distance between extrema

[A,B] = find(zc ~= 0);
[Alim, Blim] = size(im);
window = 20;
w = window;
D = zeros(size(im));
blur = zeros(size(im));
for i= 1:length(A),

%define search window
  Alower = A(i) - window;
  Blower = B(i) - window;
  Aupper = A(i) + window;
  Bupper = B(i) + window;
%clip to edge of image
  if Alower < 1
   Alower = 1;
  end
  if Blower < 1
   Blower = 1;
  end
  if Aupper>Alim
   Aupper = Alim;
  end
  if Bupper > Blim
   Bupper = Blim;
  end

currentmax = lap_of_gau(A(i),B(i));
currentmin = currentmax;
nextmin = currentmin;
nextmax = currentmax;

currentpixelx = A(i);
currentpixely = B(i);
currentangle = 45*round(rad2deg(ang(currentpixelx, currentpixely))/45);
if (currentangle == 0 )| (currentangle == 360 )| (currentangle == -360),
 pixelincx = 1;
 pixelincy = 0;
end
if (currentangle == 45 )| (currentangle == -315),
 pixelincx = 1;
 pixelincy = 1;
end
if (currentangle == 90 )| (currentangle == -270),
 pixelincx = 0;
 pixelincy = 1;
end
if (currentangle == 135 )| (currentangle == -225),
 pixelincx = -1;
 pixelincy = 1;
end
if (currentangle == 180 )| (currentangle == -180),
 pixelincx = -1;
 pixelincy = 0;
end
if (currentangle == 225 )| (currentangle == -135),
 pixelincx = -1;
 pixelincy = -1;
end
if (currentangle == 270 )| (currentangle == -90),
 pixelincx = 0;
 pixelincy = -1;
end
if (currentangle == 315 )| (currentangle == -45),
 pixelincx = 1;
 pixelincy = -1;
end
counter = 0;
while(nextmax >= currentmax),
  nextmaxx = currentpixelx + pixelincx;
  nextmaxy = currentpixely + pixelincy;
  if nextmaxx < 1 | (nextmaxx > size(im, 1)) | (nextmaxy < 1) | (nextmaxy > size(im,2))
    break;
  end
  nextmax = lap_of_gau(nextmaxx, nextmaxy);
  counter = counter +1;
  if counter > w,
   break;
  end
  currentpixelx = nextmaxx;
  currentpixely = nextmaxy;
end
maxx = currentpixelx;
maxy = currentpixely;

currentpixelx = A(i);
currentpixely = B(i);
counter = 0;
while(nextmin <= currentmin),
  nextminx = currentpixelx - pixelincx;
  nextminy = currentpixely - pixelincy;
  if nextminx < 1 | (nextminx > size(im, 1)) | (nextminy < 1) | (nextminy > size(im,2))
    break;
  end
  nextmin = lap_of_gau(nextminx, nextminy);
  counter = counter +1;
  if counter > w,
   break
  end
  currentpixelx = nextminx;
  currentpixely = nextminy;

end
minx = currentpixelx;
miny = currentpixely;

d = sqrt((maxx-minx)^2 + (maxy-miny)^2);
  temp = sqrt((d/2)^2 - scale2(A(i),B(i))^2);
  if isreal(temp) 
   blur(A(i),B(i)) = temp;
  else
   blur(A(i),B(i)) = 0;
  end
end




