edge_detection.rar_matlab例程_matlab_

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+');