亚像素边缘提取MATLAB码

matlab亚像素边缘提取代码，供学习！

clear
clc
II=imread（‘Sample_BW.jpg‘）;
figure（1）imshow（II）
I=double（II）;
[mn]=size（I）;
%用Canny算子进行像素边缘检测，判别阈值为0.1
ff=edge（I‘canny‘0.1）;
figure（2）imshow（ff）
f=double（ff）;
%求原函数梯度函数
R=zeros（mn）;
for xi=2:1:m-1
    for yi=2:1:n-1
        R（xiyi）=abs（I（xi+1yi+1）+2*I（xiyi+1）+I（xi-1yi+1）-I（xi-1yi-1）-2*I（xiyi-1）-I（xi+1yi-1））+abs（I（xi-1yi-1）+2*I（xi-1yi）+I（xi-1yi+1）-I（xi+1yi+1）-2*I（xi+1yi）-I（xi+1yi-1））;
    end
end
%开始求坐标
nn=1;
W=zeros（mn）;
for x=2:1:m-1
    for y=2:1:n-1
        if f（xy）==0
           W（xy）=0;
        else
            i=[R（x-1y）R（xy）R（x+1y）];
            j=[R（xy-1）R（xy）R（xy+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）;%求亚像素边缘坐标[XeYe]
                Ye=y+（j（1）-j（3））/（（j（1）-2*j（2）+j（3））+eps）;
                if （Xe0）&（Ye0）
                    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）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       1474  2009-05-18 09:24  yaxiangsu.m

     文件       7378  2005-11-22 10:20  Sample_BW.jpg

----------- ---------  ---------- -----  ----

                 8852                    2