文件1---hough_circle.m
function [hough_space,hough_circle,para] = hough_circle(BW,step_r,step_angle,r_min,r_max,p)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input
% BW:二值图像;
% step_r:检测的圆半径步长
% step_angle:角度步长,单位为弧度
% r_min:最小圆半径
% r_max:最大圆半径
% p:阈值,0,1之间的数 通过调此值可以得到图中圆的圆心和半径
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output
% hough_space:参数空间,h(a,b,r)表示圆心在(a,b)半径为r的圆上的点数
% hough_circl:二值图像,检测到的圆
% para:检测到的圆的圆心、半径
circleParaXYR=[];
para=[];
[m,n] = size(BW);
size_r = round((r_max-r_min)/step_r)+1;%四舍五入
size_angle = round(2*pi/step_angle);
hough_space = zeros(m,n,size_r);
[rows,cols] = find(BW);%查找非零元素的行列坐标
ecount = size(rows);%非零坐标的个数
% Hough变换
% 将图像空间(x,y)对应到参数空间(a,b,r)
% a = x-r*cos(angle)
% b = y-r*sin(angle)
for i=1:ecount
for r=1:size_r %半径步长数
for k=1:size_angle %按一定弧度把圆几等分
a = round(rows(i)-(r_min+(r-1)*step_r)*cos(k*step_angle));
b = round(cols(i)-(r_min+(r-1)*step_r)*sin(k*step_angle));
if(a>0&a<=m&b>0&b<=n)
hough_space(a,b,r) = hough_space(a,b,r)+1;%h(a,b,r)的坐标,圆心和半径
end
end
end
end
% 搜索超过阈值的聚集点。对于多个圆的检测,阈值要设的小一点!通过调此值,可以求出所有圆的圆心和半径
max_para = max(max(max(hough_space)));%返回值就是这个矩阵的最大值
index = find(hough_space>=max_para*p);%一个矩阵中,想找到其中大于max_para*p数的位置
length = size(index);%符合阈值的个数
hough_circle = false(m,n);
%hough_circle = zeros(m,n);
%通过位置求半径和圆心。
for i=1:ecount
for k=1:length
par3 = floor(index(k)/(m*n))+1;
par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;
par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m;
if((rows(i)-par1)^2+(cols(i)-par2)^2<(r_min+(par3-1)*step_r)^2+5&...
(rows(i)-par1)^2+(cols(i)-par2)^2>(r_min+(par3-1)*step_r)^2-5)
hough_circle(rows(i),cols(i)) = true; %检测的圆
end
end
end
% 从超过峰值阈值中得到
for k=1:length
par3 = floor(index(k)/(m*n))+1;%取整
par2 = floor((index(k)-(par3-1)*(m*n))/m)+1;
par1 = index(k)-(par3-1)*(m*n)-(par2-1)*m;
circleParaXYR = [circleParaXYR;par1,par2,par3];
hough_circle(par1,par2)= true; %这时得到好多圆心和半径,不同的圆的圆心处聚集好多点,这是因为所给的圆不是标准的圆
%fprintf(1,'test1:Center %d %d \n',par1,par2);
end
%集中在各个圆的圆心处的点取平均,得到针对每个圆的精确圆心和半径!
while size(circleParaXYR,1) >= 1
num=1;
XYR=[];
temp1=circleParaXYR(1,1);
temp2=circleParaXYR(1,2);
temp3=circleParaXYR(1,3);
c1=temp1;
c2=temp2;
c3=temp3;
temp3= r_min+(temp3-1)*step_r;
if size(circleParaXYR,1)>1
for k=2:size(circleParaXYR,1)
if (circleParaXYR(k,1)-temp1)^2+(circleParaXYR(k,2)-temp2)^2 > temp3^2
XYR=[XYR;circleParaXYR(k,1),circleParaXYR(k,2),circleParaXYR(k,3)]; %保存剩下圆的圆心和半径位置
else
c1=c1+circleParaXYR(k,1);
c2=c2+circleParaXYR(k,2);
c3=c3+circleParaXYR(k,3);
num=num+1;
end
end
end
%fprintf(1,'sum %d %d radius %d\n',c1,c2,r_min+(c3-1)*step_r);
c1=round(c1/num);
c2=round(c2/num);
c3=round(c3/num);
c3=r_min+(c3-1)*step_r;
%fprintf(1,'num=%d\n',num)
%fprintf(1,'Center %d %d radius %d\n',c1,c2,c3);
para=[para;c1,c2,c3]; %保存各个圆的圆心和半径的值
circleParaXYR=XYR;
end
文件2 main.m
clc,clear all
circleParaXYR=[];
I = imread('h6.jpg');
[m,n,l] = size(I);
if l>1
I = rgb2gray(I);
end
BW = edge(I,'sobel');
step_r = 1;
step_angle = 0.1;
minr = 3
maxr = 30;
thresh = 0.51;
[hough_space,hough_circle,para] = hough_circlefu(BW,step_r,step_angle,minr,maxr,thresh);
figure(1),imshow(I),title('原图')
figure(2),imshow(BW),title('边缘')
figure(3),imshow(hough_circle),title('检测结果')
circleParaXYR=para;
%输出
fprintf(1,'\n---------------圆统计----------------\n');
[r,c]=size(circleParaXYR);%r=size(circleParaXYR,1);
fprintf(1,' 检测出%d个圆\n',r);%圆的个数
fprintf(1,' 圆心 半径\n');%圆的个数
for n=1:r
fprintf(1,'%d (%d,%d) %d\n',n,floor(circleParaXYR(n,1)),floor(circleParaXYR(n,2)),floor(circleParaXYR(n,3)));
end
%标出圆
figure(4),imshow(I),title('检测出图中的圆')
hold on;
plot(circleParaXYR(:,2), circleParaXYR(:,1), 'r+');
for k = 1 : size(circleParaXYR, 1)
t=0:0.01*pi:2*pi;
x=cos(t).*circleParaXYR(k,3)+circleParaXYR(k,2);y=sin(t).*circleParaXYR(k,3)+circleParaXYR(k,1);
plot(x,y,'r-');
end
R_max=maxr;
acu=zeros(R_max);
stor =[];
for j=1:R_max
for n=1:r
if j == floor(circleParaXYR(n,3))
acu(j)= acu(j)+1;
end
end
stor=[stor;j,acu(j)];
%fprintf(1,'%d,%d\n',j,acu(j));
end
fprintf(1,'\n------------粒子大小,数目统计---------\n');
fprintf(1,'粒子半径,粒子个数\n');
for j=1:R_max
if acu(j) > 0
fprintf(1,'%4d %8d\n',stor(j,1),stor(j,2));
end
end
fprintf(1,'----------------------------------------\n');
figure(5),plot(stor(:,1),stor(:,2),'-k','LineWidth',2),title('粒径谱');
xlabel('粒子大小');
ylabel('粒子个数');
grid on;
z=[0,10,20,30,40,50,60,70,80,90,11,35,25,42,48,40,20,75,88,94,23,10,20,30,40,78,60,76,84,95,58,10,20,30,40,50,60,70,80,90,100];%给出z的坐标
Z=z(:);
S=floor(abs(Z)*1);
C=floor(abs(Z)*0.5);
figure(6),scatter3(circleParaXYR(:,1),circleParaXYR(:,2),Z,circleParaXYR(:,3)*7,'filled'),title('构建三维粒子场');
注:对网上的代码做了改动,使它更能精准检测每个圆的圆心和半径
2012-5-31