hough.txt

在Matlab上基于Hough变换开发的快速圆检测方法

 function   [x,y,r,count]=hough(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;