hough1.m

国外matlab研究站点下载的hough变换直线及圆检测代码

function [x,y,r,count]=hough1(BW1,x0,y0,r0,rr)
% BW1为输入的需检测的圆图象的矩阵,[x0,y0]为估算的定位中心的位置,r0为估算的定位的半径,
%rr是检测范围的外扩值


[m,n]=size(BW1);
%确定搜索范围，搜索的矩形的四个定点为[e3,e1],[e3,e2],[e4,e1],[e4,e2]
e1=floor(x0-r0)-rr;
if e1<1
    e1=1;
end
    e2=ceil(x0+r0)+rr;
if e2>n
    e2=n;
end
    e3=floor(y0-r0)-rr;
if e3<1
    e3=1;
end
    e4=ceil(y0+r0)+rr;
if e4>m
    e4=m;
end
%Hough变换求瞳孔中心
flag=1;count=0; %flag为循环控制变量，count记录落在所得圆上的点
while(flag)
    flag=0;
    BW2=zeros(m,n); %设0矩阵BW2，以BW2为变换空间
for j=e1:e2
    for i=e3:e4
        if BW1(i,j)==1 %在搜索范围内找到一个值为1点
            k=(j-x0)^2+(i-y0)^2-r0^2;
        if k>-1600&k<1600 %这个点在有效区内
            real=j-x0;imag=y0-i;
            theta=atan2(imag,real); %theta为BW1(i,j)相对于(x0,y0)的圆心角
            ii=round(i+r0*sin(theta));jj=round(j-r0*cos(theta));
%BW2(ii,jj)为以BW1(i,j)的确定的圆心的位置
            BW2(ii,jj)=BW2(ii,jj)+1;
                end
        end
    end
end
A=max(BW2); %BW2中的各列最大值组成向量A
count1=max(A); %count1为A中的最大值
if count1>=count %当求出的值较大时，用新值代替原来的count值,并将r减1再求一次
count=count1;
r0=r0-1;
BW3=BW2;
flag=1;
end
end
[y,x]=find(BW3==count); %最大值的位置坐标
n=length(y);
if n>1 %最大值的位置坐标不唯一时取所有位置的算术平均值
    yy=sum(y)/n;
    xx=sum(x)/n;
    y=yy;
    x=xx;
end
    r=r0+1; %定位的半径
return;