霍夫變換(Hough)是一個非常重要的檢測間斷點邊界形狀的方法。它通過将圖像坐标空間變換到參數空間，來實作直線與曲線的拟合。
參數空間可以表示為(a,b,r)，圖像坐标空間中的一個圓對應參數空間中的一個點。
A(a,b,r)。
計算過程是讓a,b在取值範圍内增加，解出滿足上式的r值，每計算出一個(a,b,r)值，就對相應的數組元素A(a,b,r)加1。計算結束後，找到的超過門檻值的A(a,b,r)所對應的a,b,r就是所求的圓的參數。
函數程式：
function [Hough_space,Hough_circle_result,Para] = Hough_circle(BW,Step_r,Step_angle,r_min,r_max,p)
%---------------------------------------------------------------------------------------------------------------------------
% input：
% BW:二值圖像;
% Step_r:檢測的圓半徑步長;
% Step_angle:角度步長，機關為弧度;
% r_min:最小圓半徑;
% r_max:最大圓半徑;
% p:以p*Hough_space的最大值為門檻值，p取0，1之間的數.
% a = x-r*cos(angle); b = y-r*sin(angle);
%---------------------------------------------------------------------------------------------------------------------------
% output：
% Hough_space:參數空間，h(a,b,r)表示圓心在(a,b)半徑為r的圓上的點數;
% Hough_circle:二值圖像，檢測到的圓;
% 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)));
%一個矩陣中，想找到其中大于max_para*p數的位置
index = find(Hough_space>=max_para*p);
length = size(index);%符合門檻值的個數
Hough_circle_result=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_result(rows(i),cols(i)) = 1;%檢測的圓
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_result(par1,par2)= 1; %這時得到好多圓心和半徑，不同的圓的圓心處聚集好多點，這是因為所給的圓不是标準的圓
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
c1=round(c1/num);
c2=round(c2/num);
c3=round(c3/num);
c3=r_min+(c3-1)*Step_r;
Para=[Para;c1,c2,c3]; %儲存各個圓的圓心和半徑的值
circleParaXYR=XYR;
end
主程式：
clc
clear all
%------------------------------------------------------------------------------
%模拟圖像
I=ones(256,256);
I(20:20:220,20:20:220)=0;
for i=20:20:220
for j=20:20:220
circle_size=floor(rand*4);
for ci=i-circle_size:i+circle_size
for cj=j-circle_size:j+circle_size
if (ci-i)^2+(cj-j)^2<circle_size^2
I(ci,cj)=0;
end
end
end
end
end
I= imgaussfilt(I, 4);
I(I>0.9)=1;
I(I<0.9)=0;
%---------------------------------------------------------------------------------
%加入直線幹擾
I(64,:)=0;I(:,64)=0;I(128,:)=0;I(:,128)=0;I(192,:)=0;I(:,192)=0;
figure(1);
imshow(I,[]),title('原圖');
%---------------------------------------------------------------------------------
% 用sobel進行邊緣檢測
BW = edge(I,'sobel');
%---------------------------------------------------------------------------------
%設定參數：
%檢測的圓半徑步長為1
Step_r = 0.5;
%角度步長0.1，機關為弧度
Step_angle = 0.1;
%最小圓半徑2
minr =5;
%最大圓半徑30
maxr = 20;
%以thresh*hough_space的最大值為門檻值，thresh取0-1之間的數
thresh = 0.8;
circleParaXYR=[];
%---------------------------------------------------------------------------------
%開始檢測
[Hough_space,Hough_circle_result,Para] = Hough_circle(BW,Step_r,Step_angle,minr,maxr,thresh);
circleParaXYR=Para;
axis equal
figure(2);
imshow(BW,[]),title('邊緣');
axis equal
figure(3);
imshow(Hough_circle_result,[]),title('檢測結果');
axis equal
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(3,k)+circleParaXYR(2,k);
y=sin(t).*circleParaXYR(3,k)+circleParaXYR(1,k);
plot(x,y,'r-');
end
- 效果圖：
如有問題請留言。