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;