如何利用MATLAB实现哈夫变换检测圆

clc,clear all

I = imread('yuan.tif')

[m,n,l] = size(I)

if l>1

I = rgb2gray(I)

end

BW = edge(I,'sobel')

step_r = 1

step_angle = 0.1

r_min = 20

r_max = 30

thresh = 0.7

% %%%%%%%%%%%%%%%%%%%%%%%%%%

% 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:检测到的圆的圆心、半径

[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

end

end

end

end

% 搜索超过阈值的聚集点

max_para = max(max(max(hough_space)))

index = find(hough_space>=max_para*thresh )

length = size(index)

hough_circle = false(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

par3 = r_min+(par3-1)*step_r

fprintf(1,'Center %d %d radius %d\n',par1,par2,par3)

para(:,k) = [par1,par2,par3]

end

subplot(221),imshow(I),title('原图')

subplot(222),imshow(BW),title('边缘')

subplot(223),imshow(hough_circle),title('检测结果')

用Hough变换来检测圆，可以标出圆心，并画出轮廓！！

http://apps.hi.baidu.com/share/detail/2196138这上面有代码，已经试过，可以运行检测。

这个不难的：

- 图像预处理，自动阀值方法二值化，然后滤掉噪声点，得到比较干净的圆形光斑离散点集；

- 用以下这个程序拟合出离散点的圆，并找出圆心。

其中第一步的自动阀值可以用otsu函数(otsu method,大津法)，其余都很基础；第二步的程序如果看不懂，可以进一步看看参考资料连接。

function [xc,yc,R,a] = circfit(x,y)

%

% [xc yx R] = circfit(x,y)

%

% fits a circle in x,y plane in a more accurate

% (less prone to ill condition )

% procedure than circfit2 but using more memory

% x,y are column vector where (x(i),y(i)) is a measured point

%

% result is center point (yc,xc) and radius R

% an optional output is the vector of coeficient a

% describing the circle's equation

%

% x^2+y^2+a(1)*x+a(2)*y+a(3)=0

%

% By: Izhak bucher 25/oct /1991,

x=x(:)y=y(:)

a=[x y ones(size(x))]\[-(x.^2+y.^2)]

xc = -.5*a(1)

yc = -.5*a(2)

R = sqrt((a(1)^2+a(2)^2)/4-a(3))

---