A=imread('Lena.gif','gif');
n=size(A,1);
e_high=80;
e_low=20;
A=double(A);
DX=[-1 0 1;-2 0 2;-1 0 1];
DY=[-1 -2 -1;0 0 0;1 2 1];
gauss3=[1 2 1;2 3 2;1 2 1]/15;
gauss5=[1 4 7 4 1;4 20 33 20 4;7 33 55 33 7;
        4 20 33 20 4;1 4 7 4 1]/331;
gauss7=[1 12 55 90 55 12 1;12 148 665 1097 665 148 12;
        55 665 2981 4915 2981 665 55;90 1097 4915 8103 4915 1097 90;
        55 665 2981 4915 2981 665 55;12 148 665 1097 665 148 12;
        1 12 55 90 55 12 1]/50887;
B=conv2(A,gauss3,'same');
GX=conv2(B,DX,'same');
GY=conv2(B,DY,'same');
G=(GX.*GX+GY.*GY).^(0.5);
T=zeros(n);
for i=1:n
  for j=1:n
    if GX(i,j)==0
      T(i,j)=1; % 1 - vertical direcion
    else
      p=atan(GY(i,j)/GX(i,j));
      if abs(p)>=3*pi/8
        T(i,j)=1; % 1 - vertical direction
      elseif p>=pi/8
        T(i,j)=2; % 2 - SW/NE direction
      elseif p<=-pi/8
        T(i,j)=4; % 4 - NW/SE direction
      else
        T(i,j)=3; % 3 - horizontal direction
      end
    end
  end
end
E=zeros(n);
for i=2:n-1
  for j=2:n-1
    if T(i,j)==1
      u=1;
      v=0;
    elseif T(i,j)==4
      u=-1;
      v=1;
    elseif T(i,j)==3
      u=0;
      v=1;
    else
      u=1;
      v=1;
    end
    if (G(i,j)>=G(i+u,j+v))&(G(i,j)>=G(i-u,j-v))
      E(i,j)=G(i,j);
    end
  end
end
for i=1:n
  for j=1:n
    if E(i,j)<e_high
      E(i,j)=0;
    end
  end
end
TEMP=255-uint8(E);
imwrite(TEMP,'Lena_edges_first_3.bmp','bmp');
