matlab 霍夫變換——檢測圓弧

類似的代碼網上一搜全都是，不少代碼運作起來各種錯誤，真正正确的還是要靠自己慢慢嘗試和修改。以下是本人親測，絕對正确可行的代碼。

首先需要建立兩個腳本，功能分别是圖像的前期處理和檢測圓弧，話不多說，代碼實作如下：

function [bw] = image_deal(bw)
 
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
%------------------------------輸入參數-----------------------------  
% BW:處理前圖像；  
%-------------------------------------------------------------------  
  
%------------------------------輸出參數-----------------------------  
% BW:二值後圖像；  
%-------------------------------------------------------------------  
  
% From Nanjing University of Science and Technology，ZLy 
%%建立sobel算子
sobel_xl=[-1,0,1;-2,0,2;-1,0,1];
sobel_xr=[1,0,-1;2,0,-2;1,0,-1];
sobel_yl=[1,2,1;0,0,0;-1,-2,-1];
sobel_yr=[-1,-2,-1;0,0,0;1,2,1];
%%圖像灰階化
bw=rgb2gray(bw);
%%把圖像資料類型轉換為雙精度浮點類型
bw=im2double(bw);
%%圖像與sobel算子進行卷積
A=conv2(bw,sobel_xl);
B=conv2(bw,sobel_xr);
C=conv2(bw,sobel_yl);
D=conv2(bw,sobel_yr);
%%圖像二值化
A=im2bw(A,graythresh(A));
B=im2bw(B,graythresh(B));
C=im2bw(C,graythresh(C));
D=im2bw(D,graythresh(D));
%%圖像相加
bw=A+B+C+D;
imshow(bw);title('卷積後的圖像');
%消除噪點
se=strel('disk',2);
bw=imclose(bw,se);
figure;imshow(bw);title('消除噪點後的圖像');
end
           

function [hough_space,hough_circle,para] = hough_circle(BW,step_r,step_angle,r_min,r_max,p,AW) 
%[HOUGH_SPACE,HOUGH_CIRCLE,PARA] = HOUGH_CIRCLE(BW,STEP_R,STEP_ANGLE,R_MAX,P)  
%------------------------------算法概述-----------------------------  
% 該算法通過a = x-r*cos(angle)，b = y-r*sin(angle)将圓圖像中的邊緣點  
% 映射到參數空間(a,b,r)中，由于是數字圖像且采取極坐标，angle和r都取  
% 一定的範圍和步長，這樣通過兩重循環（angle循環和r循環）即可将原圖像  
% 空間的點映射到參數空間中，再在參數空間（即一個由許多小立方體組成的  
% 大立方體)中尋找圓心，然後求出半徑坐标。  
%-------------------------------------------------------------------  
  
%------------------------------輸入參數-----------------------------  
% BW:二值圖像；  
% step_r:檢測的圓半徑步長  
% step_angle:角度步長，機關為弧度  
% r_min:最小圓半徑  
% r_max:最大圓半徑  
% p:以p*hough_space的最大值為門檻值，p取0，1之間的數  
% AW:原始圖像；
%-------------------------------------------------------------------  
  
%------------------------------輸出參數-----------------------------  
% 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); %傳回BW中非零元素的一個列或行向量v，同時傳回行和列的索引值。
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)));  %求出hough_space的最大值,三個方向上a,b,r
index = find(hough_space>=max_para*p);  %找到大于界定門檻值的聚集點
length = size(index);  
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)) = 1;  
        end  
    end  
end  
% 列印結果  
fprintf('圓的個數 %d\n',length);
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]';
    %%找出對應參數，圓中心加藍
    AW(par1,par2,1)=0;
    AW(par1,par2,2)=232; 
    AW(par1,par2,3)=232; 
    %%圓周加紅
    for angle=0:0.005:2*pi
        x = round(par1+par3*cos(angle));
        y = round(par2+par3*sin(angle));
        if(x>0&y>0)
           AW(x,y,1)=232;
           AW(x,y,2)=0;            
           AW(x,y,3)=0;
        end   
    end
end 
%輸出識别後圖像
imshow(AW);title('圓心加藍,圓周加紅');
           

在指令視窗輸入：

>> bw=imread('p.jpg');
>> aw=bw;
>> bw=image_deal(bw);
>> f=hough_circle(bw,10,0.03,100,150,0.5,aw);
           

原圖：

檢測結果：