Hough變換

 Hough變換可以檢測出直線與圓。

 原理介紹：

       假設直線l的斜截式方程為y=kx+b.(x,y)為圖像中的像素值，如果直線上的所有點對（x，y）都滿足這個式子，即他們有相同的參數（b，k），所以他們在同一條直線上，所以思路很清晰：例如某10個點在同一條直線上，那么這10個點共享一個參數集（b1，k1），又有另外幾個點構成了另一條直線，那么這幾個點又共享另一組參數（b2，k2）。所以說有多少個這樣的參數集，就有多少條直線！
       Hough這樣理解：將斜截式改為：b=-xk+y。將（x,y）空間轉為(b,k)空間坐標系下，將k軸等分i份，將b軸等分j份，那么可以將每一個單元稱為一個累加器單元，其值用A(i,j)表示，初值為零。

                                        

          對於圖像中每個點(x,y)，令參數k依次取值為k軸上的每個細分值，將其帶入b=-xk+y，得到b，通過對b近似將其划分至距離累加器中最近的單元格b中。每得到一對（k，b），將其相應的累加器單元的值進行累加，即A(p,q)=A(p,q)+1.那么很好理解：A的非零值個數為直線個數，A(i,j)得值即為直線上的點個數。

           因為當直線垂直時，斜率無窮大，所以采用直線的標准表達式：ρ=xcosθ+ysinθ.基本沒什么變化，就是將(k,b)變為(ρ,θ).θ一般取值為[-90,90],或[0,180]度。ρ取值為[-D,D]，D為圖像對角線長度。

                                     

 代碼:      

function Hough_LineDetect(filename,degree_range,line_length,dtheta,drho)
%Hough變換檢測直線
%參數：filename：被檢測圖像文件名
%     degree_range：檢測直線的角度范圍，1*2矩陣：默認值為[-90 90]
%     line_length：檢測直線的最小長度；默認值為100
%     dtheta：theta的步長；默認值為1(角度)
%     drho：rho的步長；默認值為1(像素)
%功能：從被檢測圖像中檢測出滿足指定角度和長度的直線。
if nargin < 5
    drho=1;
end
if nargin<4
    dtheta=1;
end
if nargin<3
    line_length=100;
end
if nargin<2
    degree_range=[-90 90]
end

I=imread(filename);
[width,height]=size(I);   %圖像尺寸
BI=edge(I);   %邊緣提取

dtheta=dtheta*pi/180;
radian_upper=max(degree_range*pi/180);
radian_lower=min(degree_range*pi/180);
radian_range=radian_upper-radian_lower;

rho_max=(sqrt(width^2+height^2));
nrho=ceil(2*rho_max/drho);
theta_value=[radian_lower:dtheta:radian_upper];    %theta值
ntheta=length(theta_value);

rho_matrix=zeros(nrho,ntheta);     %對於某一rho和theta值，點的個數
hough_line=zeros(width,height);

%Hough 變換
%將圖像空間(x,y)變換到參數空間(rho,theta)：rho=xcos(theta)+ysin(theta)
[rows,cols]=find(BI);    
pointcount=length(rows);       %邊緣點數個數
rho_value=zeros(pointcount,ntheta);   %rho取值，對於每個點對應一個theta
for i=1:pointcount
    m=rows(i);
    n=cols(i);
    for k=1:ntheta
        rho=(m*cos(theta_value(k)))+(n*sin(theta_value(k)));  %根據x，y，theta求rho
        rho_index = round((rho+rho_max)/drho);  %求索引
        rho_matrix(rho_index,k)=rho_matrix(rho_index,k)+1;
        rho_value(rho_index,k)=rho;  
    end
end

%搜索同一直線上的點
index=find(rho_matrix>line_length);  %求滿足條件的索引值
for k=1:length(index)
    [rho_th,theta_th]=ind2sub(size(rho_matrix),index(k)); %根據索引求參數值
    theta=theta_value(theta_th);
    rho=rho_value(rho_th,theta_th);
    for i=1:pointcount
        x=rows(i);
        y=cols(i);
        rate=(x*cos(theta)+y*sin(theta))/rho;   %求rate來設定條件
        if (rate>1-10^-3 & rate<1+10^-3)   %滿足條件
            hough_line(x,y)=1;  %將滿足條件即直線的點塗黑
        end
    end
end

figure;imshow(I);title('原始圖像');
figure;imshow(BI);title('邊緣檢測后的圖像');
figure;imshow(hough_line);title('hough變換檢測到的直線');

                   

     檢測圓：原理類似，和（ρ,θ）來表示一條直線相似，使用（a,b,r）來確定一個圓心為（a,b）半徑為 r  的圓。

                 某個圓過點（x1,y1），則有：(x1-a1)^2 + (y1-b1)^2 = r1^2 。所以在圓（a1,b1,r1）上的點共用這一套參數集。參數集的個數代表圓的個數，共用某一參數集的點數目即為此圓上的點數。

    代碼：

function hough_CircleDetect(filename,radius_range,step_angle,step_radius)
%Hough變換檢測圓
%參數：filename：被檢測圖像文件名
%     radius_range：檢測圓的半徑范圍，1*2矩陣，默認值為[10 100]
%     step_angle：角度步長；默認值為5（角度）
%     step_radius：半徑步長；默認值為1（像素）
%功能：從被檢測圖像中檢測出滿足指定半徑的圓。
if nargin<4
    step_radius=1;
end
if nargin<3
    step_angle=5;
end
if nargin<2
    radius_range=[10 100];
end
radius_min=min(radius_range);
radius_max=max(radius_range);
step_angle=step_angle*pi/180;

I=imread(filename);
[m,n,l]=size(I);
if l>1
    I=rgb2gray(I);
end
BI=edge(I);
[rows,cols]=find(BI);
PointCount=size(rows);

RadiusCount=ceil((radius_max-radius_min)/step_radius);
AngleCount=ceil(2*pi/step_angle);
hough_space=zeros(m,n,RadiusCount);

%Hought變換
%將圖像空間（x,y）變換到參數空間（a,b,r）
%a=x-r*cos(theta);
%b=y-r*sin(theta);
for i=1:PointCount
    for r=1:RadiusCount
        for k=1:AngleCount
            a=round(rows(i)-(radius_min+(r-1)*step_radius)*cos(k*step_angle));
            b=round(cols(i)-(radius_min+(r-1)*step_radius)*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

%搜索同一圓上的點
thresh = 0.7;
max_PointCount=max(max(max(hough_space)));
index=find(hough_space>=max_PointCount*thresh);
length=size(index);
hough_circle=zeros(m,n);
size_hough_space = size(hough_space);
for i=1:PointCount
    for k=1:length
        [a,b,r]=ind2sub(size_hough_space,index(k));
        rate=((rows(i)-a)^2+(cols(i)-b)^2)/(radius_min+(r-1)*step_radius)^2;
        if(rate<1.1)
            hough_circle(rows(i),cols(i))=1;
        end
    end
end

figure;imshow(I);title('原始圖像');
figure;imshow(BI);title('邊緣檢測后的圖像');
figure;imshow(hough_circle);title('hough變換檢測到的圓');

%hough_CircleDetect('line_circle.bmp')
%hough_CircleDetect('line_circle.bmp',[5,20])
%hough_CircleDetect('line_circle.bmp',[5,100])
%hough_CircleDetect('line_circle.bmp',[5,200])

              

 

其他參考資料：

http://cn.mathworks.com/help/images/ref/imfindcircles.html

http://cn.mathworks.com/matlabcentral/fileexchange/26978-hough-transform-for-circles?focused=6019424&tab=example

https://stackoverflow.com/questions/27171527/circle-detection-from-gray-level-image-in-matlab