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