亚像素边缘检测源代码，亲测可用

clear clc II=imread('Sample_BW.jpg'); figure(1),imshow(II) I=double(II); [m,n]=size(I); %用Canny算子进行像素边缘检测，判别阈值为0.1 ff=edge(I,'canny',0.1); figure(2),imshow(ff) f=double(ff); %求原函数梯度函数 R=zeros(m,n); for xi=2:1:m-1 for yi=2:1:n-1 R(xi,yi)=abs(I(xi+1,yi+1)+2*I(xi,yi+1)+I(xi-1,yi+1)-I(xi-1,yi-1)-2*I(xi,yi-1)-I(xi+1,yi-1))+abs(I(xi-1,yi-1)+2*I(xi-1,yi)+I(xi-1,yi+1)-I(xi+1,yi+1)-2*I(xi+1,yi)-I(xi+1,yi-1)); end end %开始求坐标 nn=1; W=zeros(m,n); for x=2:1:m-1 for y=2:1:n-1 if f(x,y)==0 W(x,y)=0; else i=[R(x-1,y),R(x,y),R(x+1,y)]; j=[R(x,y-1),R(x,y),R(x,y+1)]; if (i(2)>i(1))&(i(2)>i(3))&(j(2)>j(1))&(j(2)>j(3)) Xe=x+(i(1)-i(3))/((i(1)-2*i(2)+i(3))+eps);%求亚像素边缘坐标[Xe,Ye] Ye=y+(j(1)-j(3))/((j(1)-2*j(2)+j(3))+eps); if (Xe<m)&(Xe>0)&(Ye<n)&(Ye>0) Xe1(nn)=Xe; %得到的亚像素边缘点 Ye1(nn)=Ye; Xe2(nn)=round(Xe); %得到的逼近的边缘点 Ye2(nn)=round(Ye); W(Xe2(nn),Ye2(nn))=1; %把检测到的边缘点直观的表示出来 nn=nn+1; %设置循环，用数组的形势表示Xe2，Ye2，方便后期查看处理亚像素边缘点坐标 end end end end end figure(3),imshow(W)