霍夫直线变换，利用梯度方向加快计算速度.
%确定accumulate大小。初始化
%计算hough变化，增加相应accumulator的vote.
%通过阀值，过滤，之后再计算得到局部最大值。
%利用梯度方向加快计算速度.
function lines=houghline2(img,ratio)
lines.x1=0;lines.y1=0;lines.x2=0;lines.y2=0;
BW = edge(img,'canny');
gradient=gradientdirection(img);
[ry,cx]=size(img);
%theta=[1,180];
rmin=-cx;%负数
rmax=max(ry,cx)*sqrt(2);
accumulator=zeros(round(rmax-rmin+1),180);%确定accumulate大小。初始化
for y=1:ry
for x=1:cx
if(BW(y,x)==1)
theta=0;
if(gradient(y,x)>=0)
theta=gradient(y,x);
else if (abs(gradient(y,x))==90)
theta=0;
else
theta=180-abs(gradient(y,x));
end
end
for k=theta-20+1:theta+20
if(k>=1 && k<=180)
r=round(x*cos(pi/180*k)+y*sin(pi/180*k)-rmin+1); %r存在负数,归一化，从1索引开始
accumulator(r,k)=accumulator(r,k)+1;
end
end
end
end
end
[acc_h, acc_w]=size(accumulator);
% average=sum(accumulator(:))/numel(accumulator);
% threshold=round(average*1.3);
[maxaa,indexmax]=max(accumulator(:));
threshold=round(max(accumulator(:))*ratio);
localmax=0;
dx=3;dy=3;
n=1;%the point index
for y=1:acc_h
for x=1:acc_w
if (accumulator(y,x)>threshold)
localmax=accumulator(y,x);
for ly=-dy:dy
for lx=-dx:dx
if( (y+ly)>=1 && (y+ly)<=acc_h && (x+lx)>=1 && (x+lx)<=acc_w)
if(accumulator(y+ly,x+lx)>localmax)
localmax=accumulator(y+ly,x+lx);
ly=dy+1;%记号，退出循环
lx=dx+1;
end
end
end
end
if(localmax>accumulator(y,x))
continue;
end
if x==0 || x==180
lines(n).x1=y;
lines(n).y1=1;
lines(n).x2=y;
lines(n).y2=ry;
n=n+1;
else if x==90
lines(n).x1=1;
lines(n).y1=y;
lines(n).x2=cx;
lines(n).y2=y;
n=n+1;
else
lines(n).y1=1;
lines(n).x1=(y+rmin-1)/cos(pi/180*x);
lines(n).y2=ry;
lines(n).x2=(y+rmin-1-(ry)*sin(pi/180*x))/cos(pi/180*x);
n=n+1;
end
end
end
end
end
梯度方向计算：
function gradientimage=gradientdirection360(img)
[height,width]=size(img);
gradientimage=zeros(height,width);
gradient_temp=zeros(height+2,width+2);
gradient_temp(2:height+1,2:width+1)=img;
gradient_temp(1,2:width+1)=img(1,:);
gradient_temp(height+2,2:width+1)=img(height,:);
gradient_temp(2:height+1,1)=img(:,1);
gradient_temp(2:height+1,2+width)=img(:,width);
theta=0;
for y=2:height+1
for x=2:width+1
gx=gradient_temp(y,x+1)-gradient_temp(y,x-1);
gy=gradient_temp(y+1,x)-gradient_temp(y-1,x);
theta=floor(atan(gy/(gx+1e-6))*(180/pi));
if (gx>0 && gy>0)
theta=theta;
else if(gx<0 && gy>0)
theta=theta+180;
else if(gx<0 && gy<0)
theta=theta+180;
else if(gx>0 && gy<0)
theta=theta+360;
else if(gx==0 && gy>0)
theta=theta; %90
else if(gx==0 && gy<0)
theta=theta+360; %270
else if(gx>0 && gy==0)
theta=theta; % 0
else if(gx<0 && gy==0)
theta=theta+180; %180
else if(gx==0 && gy==0)
theta=-1; %-1±íʾûÓÐÌݶȡ£ÊÇƽ»¬µÄÇøÓò
end
end
end
end
end
end
end
end
end
gradientimage(y-1,x-1)=round(theta);
end
end运行脚本文件：
close all;
clc;
img=imread('example.png');
img=rgb2gray(img);
%img=imresize(img,0.5);
imshow(img);
tic
lines=houghline2(img,0.4);
toc
hold on;
m=length(lines);
for i=1:m
plot([lines(i).x1,lines(i).x2],[lines(i).y1,lines(i).y2],'b');
end
霍夫圆变换：
%threshhold:votes的阀值 0.8: 0.8*maxvote
%dp:a,b采样的频率参数 1:0-360 2; 0:2:360
%dr:半径的分辨率 1：minR:maxR,2: minR:2:maxR
%minR,最小圆半径
%maxR,最大圆半径
%similarity:圆的相似度 c1=(x1 y1 r1), c2=(x2 y2 r2), smilarity = |x1-x2|+|y1-y2|+|r1-r2|
function circles=houghcircle2(img,minR,maxR,threshhold)
circles(1).x=0;circles(1).y=0;circles(1).radii=0;
[height,width]=size(img);
%gradient=gradientdirection360(img);
BW = edge(img,'canny');
figure;
imshow(BW);
%初始化类accumulator
accumulator=zeros(2*maxR+height,2*maxR+width,maxR-minR+1);
%得到满足在[minR,maxR]之间的点，(realx,realy)；对图像进行canny变化
%得到边缘点，（ex,ey）,然后对(realx,realy)进行平移得到：（realx+ex,realy+ey）
[X Y] = meshgrid(-maxR:maxR, -maxR:maxR);
Rmap=round((sqrt((X).^2+(Y).^2)));
realX=X((Rmap>=minR) & (Rmap<=maxR)); %满足条件的x和y value,向量
realY=Y((Rmap>=minR) & (Rmap<=maxR));
r=Rmap((Rmap>=minR) & (Rmap<=maxR));
%得到边缘的点
[Ey,Ex]=find(BW==1);
mEdge=length(Ex);
for i=1:mEdge
%圆心可能所在的位置
centerX=Ex(i)+realX;%vector
centerY=Ey(i)+realY;
%归一化,索引从1开始
index=sub2ind(size(accumulator),centerY+maxR,centerX+maxR,r-minR+1);
%vote加1
accumulator(index)=accumulator(index)+1;
end
[acc_rows,acc_cols,acc_3dim]=size(accumulator);
localmax=0;
%求局部极大值的领域尺寸
dx=3;dy=3;
n=1;
for r=minR:maxR
%选出满足条件的candidate点。条件：周长的大小
slice=accumulator(:,:,r-minR+1);
slice(slice<round(threshhold*2*pi*r))=0;
[y,x,votes]=find(slice);
%求局部极大值
m=length(y);
for i=1:m
y2=y(i);
x2=x(i);
localmax=votes(i);
for ly=-dy:dy
for lx=-dx:dx
if( (y2+ly)>=1 && (y2+ly)<=acc_rows && (x2+lx)>=1 && (x2+lx)<=acc_cols)
index=sub2ind(size(accumulator),y2+ly,x2+lx,r-minR+1);
if(accumulator(index)>localmax)
localmax=accumulator(index);
ly=dy+1;%记号trick，退出循环
lx=dx+1;
end
end
end
end
if(localmax>votes(i))
continue;
end
end
if m==0
continue;
end
circles(n).x=x2-maxR;
circles(n).y=y2-maxR;
circles(n).radii=r;
n=n+1;
end
figure;
imshow(img);
hold on;
for i = 1:length(circles)
x = circles(i).x-circles(i).radii;
y = circles(i).y-circles(i).radii;
if((x<=0) || (y<=0))
continue;
end
w = 2*circles(i).radii;
rectangle('Position', [x y w w], 'EdgeColor', 'red', 'Curvature', [1 1]);
end
运行脚本文件：
close all;
clc;
img=imread('test2.jpg');
img=imresize(img,0.6);
img=rgb2gray(img);
figure;
imshow(img);
houghcircle2(img,0.8);
pause;
采用梯度方向来定位可能的圆心，这样就不用找一圈（360度）所有可能的点，节省计算时间。
但是1.程序却不能找到圆。。不知道什么原因。而且2.运行速度比之前的还慢。。写代码能力不行啊。。
把这段找不到圆的代码贴上去先。以后自己发现错误了再看。。
%注意：：：：：！！
% 这个程序用梯度方向来定位可能的圆心，这样就不用找一圈（360度）所有可能的点，节省计算时间。
% 但是
% 1.程序却不能找到圆。。不知道什么原因。而且2.运行速度比之前的还慢。。写代码能力不行啊。。
%threshhold:votes的阀值 0.8: 0.8*maxvote
%dp:a,b采样的频率参数 1:0-360 2; 0:2:360
%dr:半径的分辨率 1：minR:maxR,2: minR:2:maxR
%minR,最小圆半径
%maxR,最大圆半径
%similarity:圆的相似度 c1=(x1 y1 r1), c2=(x2 y2 r2), smilarity = |x1-x2|+|y1-y2|+|r1-r2|
function circles=houghcircle3(img,minR,maxR,threshhold)
circles(1).x=0;circles(1).y=0;circles(1).radii=0;
[height,width]=size(img);
gradient=gradientdirection360(img);
BW = edge(img,'canny');
figure;
imshow(BW);
%初始化类accumulator
accumulator=zeros(2*maxR+height,2*maxR+width,maxR-minR+1);
%得到满足在[minR,maxR]之间的点，(realx,realy)；对图像进行canny变化
%得到边缘点，（ex,ey）,然后对(realx,realy)进行平移得到：（realx+ex,realy+ey）
[X Y] = meshgrid(-maxR:maxR, -maxR:maxR);
Rmap=round((sqrt((X).^2+(Y).^2)));
realX=X((Rmap>=minR) & (Rmap<=maxR)); %满足条件的x和y value,向量
realY=Y((Rmap>=minR) & (Rmap<=maxR));
r=Rmap((Rmap>=minR) & (Rmap<=maxR));
%
ration=zeros(360,1);
count=zeros(360,1);
alpha=[0:359];
aa=length(realX(realX>minR*cos(alpha(2)) & realX<maxR*cos(alpha(2))));
%找到每个角度满足条件的realX(realY)坐标的索引，并存储在ration中。
for i=1:360
% BB=realX>minR*cos(alpha(i)) && realX<maxR*cos(alpha(i));
if i==1 || i==181
index_y=find(realY==0);
else
index_y= find(realY>minR*sin(alpha(i)) & realY<maxR*sin(alpha(i)));
end
% vectorx=realX(index_y);
if i==91 || i==271
index_x=find(realX==0);
else
index_x=find(realX>minR*cos(alpha(i)) & realX<maxR*cos(alpha(i)));
end
k=1;
if length(index_y)<length(index_x)
for j=1:length(index_y)
if(find(index_x==index_y(j)))
ration(i,k)=index_y(j);
k=k+1;
end
end
else
for j=1:length(index_x)
if(find(index_y==index_x(j)))
ration(i,k)=index_x(j);
k=k+1;
end
end
end
count(i)=k-1;
end
%得到边缘的点
[Ey,Ex]=find(BW==1);
mEdge=length(Ex);
for i=1:mEdge
%圆心可能所在的位置
degree=gradient(Ey(i),Ex(i))+1;
index_zb=ration(degree,:);
index_zb=index_zb(1:count(degree));
if isempty(index_zb)
continue;
end
realX2=realX(index_zb);
realY2=realY(index_zb);
rr=r(index_zb);
centerX=Ex(i)+realX2;%vector
centerY=Ey(i)+realY2;
%归一化,索引从1开始
index=sub2ind(size(accumulator),centerY+maxR,centerX+maxR,rr-minR+1);
%vote加1
accumulator(index)=accumulator(index)+1;
end
[acc_rows,acc_cols,acc_3dim]=size(accumulator);
localmax=0;
%求局部极大值的领域尺寸
dx=3;dy=3;
n=1;
for r=minR:maxR
%选出满足条件的candidate点。条件：周长的大小
slice=accumulator(:,:,r-minR+1);
slice(slice<round(threshhold*2*pi*r))=0;
[y,x,votes]=find(slice);
%求局部极大值
m=length(y);
for i=1:m
y2=y(i);
x2=x(i);
localmax=votes(i);
for ly=-dy:dy
for lx=-dx:dx
if( (y2+ly)>=1 && (y2+ly)<=acc_rows && (x2+lx)>=1 && (x2+lx)<=acc_cols)
index=sub2ind(size(accumulator),y2+ly,x2+lx,r-minR+1);
if(accumulator(index)>localmax)
localmax=accumulator(index);
ly=dy+1;%记号trick，退出循环
lx=dx+1;
end
end
end
end
if(localmax>votes(i))
continue;
end
end
if m==0
continue;
end
circles(n).x=x2-maxR;
circles(n).y=y2-maxR;
circles(n).radii=r;
n=n+1;
end
aaa=circles;
figure;
imshow(img);
hold on;
for i = 1:length(circles)
x = circles(i).x-circles(i).radii;
y = circles(i).y-circles(i).radii;
if((x<=0) || (y<=0))
continue;
end
w = 2*circles(i).radii;
rectangle('Position', [x y w w], 'EdgeColor', 'red', 'Curvature', [1 1]);
end
参考资料：
http://en.wikipedia.org/wiki/Hough_transform
http://homepages.inf.ed.ac.uk/rbf/HIPR2/hough.htm
http://matlabtricks.com/post-39/understanding-the-hough-transform
5295
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
