hough.m

hough2005用于检测圆的半径和圆的位置

function res=hough100(im,RHO_MAX,THETA_MAX)
if (nargin~=3 | nargout~=1)
	fprintf(1,'Correct use: res=hough2(im,RHO_MAX,THETA_MAX).\n');
	error('Exiting...\n');
end

[X,Y]=size(im);
if X~=Y
	fprintf(1,'Input image is not square.\n');
	error('Exiting...\n');
elseif rem(X,2)==1
	fprintf(1,'Input image size has to be even in pixels.\n');
	error('Exiting...\n');
end

%hack to make it work when there is one feature point
if sum(im(:)) == 1
   [i,j]=find(im);
   
   is=i-1:i+1;
   js=j-1:j+1;
   
   is=is(find(is>0&is<size(im,1)));
   js=js(find(js>0&js<size(im,2)));
   
   im_new=im;
   im_new(is,js)=1;
   
   im2=im_new;
   im2(i,j)=0;
   
   h_new=hough3(im_new,RHO_MAX,THETA_MAX);
   h2=hough3(im2,RHO_MAX,THETA_MAX);
   
   res=h_new-h2;
   
   return;
elseif sum(im(:)) == 0
   
   res=zeros(RHO_MAX,THETA_MAX);
   
   return;
end

tic

d_rho=X/RHO_MAX;
d_theta=pi/THETA_MAX;

theta=0:d_theta:pi-d_theta;

smat=sin(theta);
cmat=cos(theta);

fprintf('Finding feature points.\n');
[x,y]=find(im);

fprintf('Translating so the origin is in the middle of the image.\n');
fprintf('Doing the Hough Transform.\n');
h1=((y-Y/2-1) * smat + (x-X/2-1) * cmat )/d_rho;

fprintf('Rounding.\n');
h3=round(h1+(RHO_MAX+1)/2);

clear h1 im x y
temp=flipud(sort(h3));
 
difference=(diff([max(temp)+1;temp])~=0);

fprintf('calc sorting-differences time:\n');

% step 3
K=find(difference);

[i,j]=find(difference);
rho=temp(K);

votes=zeros(size(i));
for counter=1:length(j)-1
   
   if j(counter)==j(counter+1)
      votes(counter)=i(counter+1)-i(counter);
   else
      votes(counter)=size(h3,1)+1-i(counter);
   end
   
end
votes(length(j))=size(h3,1)+1-i(length(j));

% step 4
res=zeros(RHO_MAX,THETA_MAX);

for counter=1:length(j)
   
   if rho(counter)>-1 & rho(counter)<RHO_MAX
      res(rho(counter)+1,j(counter))=votes(counter);
   end
end
toc