hough变换检测直线和圆

close all
filename='hough/data/circle.bmp';
degree_range = [0 90];  %检测0到90度范围的直线
line_length = 100;      %检测长度大于100的直线
dtheta=1;               %theta的步长
drho=1;                 %rho的步长
%图像预处理
I=rgb2gray(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];
ntheta  =  length(theta_value);

rho_matrix = zeros(nrho,ntheta);
hough_line = zeros(width,height);

%霍夫变换，将空间域变换到参数域  rho=x*cos(theta)+y*sin(theta),rho和theta是参数
[rows,cols] = find(BI);
pointcount = length(rows);
rho_value  = zeros(pointcount,ntheta);
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)));
        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_rho,index_theta] = find(rho_matrix>line_length);
for k = 1:length(index_rho)
    
    theta = theta_value(index_theta(k));
    rho   = rho_value(index_rho(k),index_theta(k));
   
    for i = 1:pointcount
        x=rows(i);
        y=cols(i);
        rate = (x*cos(theta)+y*sin(theta))/rho;
        if(rate>1-10^-3&rate<1+10^-3)
            hough_line(x,y)=1;
        end
    end
end
figure(1);imshow(I);title('src image');
figure(2);imshow(BI);title('edge image');
figure(3);imshow(rho_matrix,[]);title('t');
figure(4);imshow(hough_line);title('hough image');

close all
filename     = 'hough\data\circle.bmp';                           %文件名
radius_range = [5 80];                                            %检测半径范围
step_angle   = 5;                                                 %角度步长
step_radius  = 1;                                                 %半径步长
radius_min   = min(radius_range);
radius_max   = max(radius_range);
step_angle   = step_angle*pi/180;

I            = rgb2gray(imread(filename));
[m,n]        = size(I);
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);

%将原图转换到参数域，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.48;
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(1),imshow(I); title('原图');
figure(2),imshow(BI);title('边缘图');
figure(3),imshow(hough_circle);title('结果图');