X=imread('im1.bmp');
X=im2double(X);
[H W] = size(X);
sigma=1;
kernelsz=sigma.*6+1;
maxthresh=0.9;
minthresh=maxthresh.*(0.04);
PSF = fspecial('gaussian',kernelsz,sigma);
Blurred1 = imfilter(X,PSF,'conv');
PSF = fspecial('gaussian',25,4);
Blurred4 = imfilter(X,PSF,'conv');
for r=2:H-1
for c=2:W-1
dx(r,c)=0.25*Blurred1(r, c+1)+0.125*Blurred1(r-1, c+1)+0.125*Blurred1(r+1, c+1)-0.25*Blurred1(r, c-1)-0.125*Blurred1(r-1, c-1)-0.125*Blurred1(r+1, c-1);
dy(r,c)=-0.25*Blurred1(r+1, c)-0.125*Blurred1(r+1, c-1)-0.125*Blurred1(r+1, c+1)+0.25*Blurred1(r-1, c)+0.125*Blurred1(r-1, c-1)+0.125*Blurred1(r-1, c+1);
end
end
for r=2:H-1
for c=2:W-1
GRADIENT(r,c) = sqrt (dx(r,c)^2 + dy(r,c)^2 );
end
end
%%Perform Non maximum suppression:
non_max = GRADIENT;
for r=1+ceil(kernelsz/2):H-ceil(kernelsz/2)
for c=1+ceil(kernelsz/2):W-ceil(kernelsz/2)
if (dx(r,c) == 0) tangent = 10;
else tangent = (dy(r,c)/dx(r,c));
end
if (-0.4142<tangent && tangent<=0.4142)
if(GRADIENT(r,c)<GRADIENT(r,c+1) || GRADIENT(r,c)<GRADIENT(r,c-1))
non_max(r,c)=0;
end
end
if (0.4142<tangent && tangent<=2.4142)
if(GRADIENT(r,c)<GRADIENT(r-1,c+1) || GRADIENT(r,c)<GRADIENT(r+1,c-1))
non_max(r,c)=0;
end
end
if ( abs(tangent) >2.4142)
if(GRADIENT(r,c)<GRADIENT(r-1,c) || GRADIENT(r,c)<GRADIENT(r+1,c))
non_max(r,c)=0;
end
end
if (-2.4142<tangent && tangent<= -0.4142)
if(GRADIENT(r,c)<GRADIENT(r-1,c-1) || GRADIENT(r,c)<GRADIENT(r+1,c+1))
non_max(r,c)=0;
end
end
end
end
subplot(2,4,1);
imshow(X);
title ('Original Image');
subplot(2,4,2);
imshow(Blurred1);
title ('sigma=1')
subplot(2,4,3);
imshow(Blurred4);
title ('sigma=4')
subplot(2,4,5);
imshow(dx);
title ('Sobel X');
subplot(2,4,6);
imshow(dy);
title ('Sobel Y')
subplot(2,4,7);
imshow(GRADIENT);
title ('Gradient')
subplot(2,4,8);
imshow(non_max);
title ('Non-max')
hysteresis = non_max;
for r=1+ceil(kernelsz/2):H-ceil(kernelsz/2)
for c=1+ceil(kernelsz/2):W-ceil(kernelsz/2)
if(hysteresis(r,c)>=maxthresh)
hysteresis(r,c)=1;
end
if(hysteresis(r,c)<maxthresh && hysteresis(r,c)>=minthresh)
hysteresis(r,c)=2;
end
if(hysteresis(r,c)<minthresh)
hysteresis(r,c)=0;
end
end
end
for r=1+ceil(kernelsz/2):H-ceil(kernelsz/2)
for c=1+ceil(kernelsz/2):W-ceil(kernelsz/2)
if (hysteresis(r,c)>0)
if(hysteresis(r,c)==2)
if( hysteresis(r-1,c-1)==1 || hysteresis(r-1,c)==1 || hysteresis(r-1,c+1)==1 || hysteresis(r,c-1)==1 || hysteresis(r,c+1)==1 || hysteresis(r+1,c-1)==1 || hysteresis(r+1,c)==1 || hysteresis(r+1,c+1)==1 )
hysteresis(r,c)=1;
end
end
end
end
end
figure;imshow(hysteresis);
EIGIMAGE = zeros(H,W);
result = zeros(H,W);
R = zeros(H,W);
winsz=5;
Rmax = 0;
for r = 2:H-1
for c = 2:W-1
A = zeros(2,2);
for rr=-(floor(winsz/2)):floor(winsz/2)
for cc=(-floor(winsz/2)):floor(winsz/2)
if(r+rr>0 && c+cc>0 && r+rr<H && c+cc<W)
A=A+[dx(r+rr,c+cc)*dx(r+rr,c+cc) dx(r+rr,c+cc)*dy(r+rr,c+cc); dy(r+rr,c+cc)*dx(r+rr,c+cc) dy(r+rr,c+cc)*dy(r+rr,c+cc)];
[V,D] = eig(A);
EIGV1(r,c)=D(1,1);
EIGV2(r,c)=D(2,2);
clear V;
clear D;
if EIGV1(r,c) < EIGV2(r,c)
EIGIMAGE(r,c)=EIGV1(r,c);
else
EIGIMAGE(r,c)=EIGV2(r,c);
end;
R(r,c) = det(A)-0.06*(trace(A))^2;
if R(r,c) > Rmax
Rmax = R(r,c);
end;
end;
end;
end;
end;
end;
figure,imshow(EIGIMAGE)
cnt = 0;
for r = 2:H-1
for c = 2:W-1
if R(r,c) > 0.01*Rmax && R(r,c) > R(r-1,c-1) && R(r,c) > R(r-1,c) && R(r,c) > R(r-1,c+1) && R(r,c) > R(r,c-1) && R(r,c) > R(r,c+1) && R(r,c) > R(r+1,c-1) && R(r,c) > R(r+1,c) && R(r,c) > R(r+1,c+1)
result(r,c) = 1;
cnt = cnt+1;
end;
end;
end;
[pc, pr] = find(result == 1);
figure,imshow(X)
hold on;
plot(pr,pc,'r+');