传统Hough算法检测圆的MATLAB实现

clear;tic%计时开始 pic1 = imread('coin.png'); pic2 = rgb2gray(pic1); sigma=2;%标准差大小 window=double(uint8(3*sigma)*2+1);%窗口大小一半为3*sigma H=fspecial('gaussian', window, sigma);%fspecial('gaussian', hsize, sigma)产生滤波模板 imag=imfilter(pic2,H,'replicate'); subplot(221);imshow(imag);title('高斯滤波后'); rmin = 10; rmax = 50; value_r(rmax) = 0; flag_r(rmax) = 0; area_num = zeros(1,rmax); xmax = zeros(1,rmax); ymax = zeros(1,rmax); [high,width] = size(imag); % 获得图像的高度和宽度 vote = zeros(high,width,rmax);%投票结果初始化 flag = zeros(high,width,rmax); two_d_vote = zeros(high,width); two_d_center = zeros(high,width); flag_num = 0; sum_d = 0; U = double(imag); uSobel = imag; for i = 2:high - 1 %sobel边缘检测 for j = 2:width - 1 Gy(i,j) = (U(i+1,j-1) + 2*U(i+1,j) + U(i+1,j+1)) - (U(i-1,j-1) + 2*U(i-1,j) + U(i-1,j+1)); Gx(i,j) = -(U(i-1,j+1) + 2*U(i,j+1) + U(i+1,j+1)) + (U(i-1,j-1) + 2*U(i,j-1) + U(i+1,j-1)); uSobel(i,j) = sqrt(Gx(i,j)^2 + Gy(i,j)^2); if Gx(i,j) == 0; a(i,j) = 0; else if Gx(i,j)>0 if Gy(i,j)>0 a(i,j) = atan(Gy(i,j)/Gx(i,j));%1 else a(i,j) = 2*pi + atan(Gy(i,j)/Gx(i,j));%4 end else if Gy(i,j)>0 a(i,j) = pi + atan(Gy(i,j)/Gx(i,j));%2 else a(i,j) = pi + atan(Gy(i,j)/Gx(i,j));%3 end end end end end uSobel(1,:) = 0;uSobel(:,1) = 0;%去黑边 subplot(222);imshow(im2uint8(uSobel));title('边缘检测后'); %画出边缘检测后的图像 for rr = rmin:rmax for i = 1:high for j = 1:width if uSobel(i,j) > 80 %二值化的阈值 Sobel2(i,j) = 1; %=====================%一般方向的vote for ii = 1:high; for jj = 1:width; if 0.95*rr^2<=((ii-i)^2+(jj-j)^2)&&((ii-i)^2+(jj-j)^2)<= 1.05*rr^2;%1 if rand(1)<10/rr %不同半径的投票权重相同 vote(ii,jj,rr) = vote(ii,jj,rr) + 1; end end end end else Sobel2(i,j) = 0; end end end end subplot(223);imshow(im2uint8(Sobel2));title('二值化后'); %画出二值化后的图像 %subplot(333);mesh(vote);title('投票结果'); %============================%找出几个可能的半径 vote_maxall = max(max(max(vote))); for nn = rmin:rmax vote_max(nn)=max(max(vote(:,:,nn))); if vote_max(nn)>=0.75*vote_maxall;%可能存在圆阈值~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ flag_num = flag_num +1; flag_r(nn) = 1; value_r(flag_num) = nn;%圆的半径 end end %===========================% %=====================%找出圆心区域 for nn = rmin:rmax if flag_r(nn) == 1 for i = 1:high-1; for j = 1:width-1; if vote(i,j,nn) >= 0.95*vote_max(nn) %圆心区域阈值~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ center(i,j,nn) = 1; else center(i,j,nn) = 0; end end end end end %subplot(334);mesh(center);title('圆心区域'); %=====================%圆心区域分区 for nn = rmin:rmax if flag_r(nn) == 1 two_d_vote = vote(:,:,nn); two_d_center = center(:,:,nn); while max(max(two_d_center)) == 1 [xmax,ymax] = find(two_d_vote==max(max(two_d_vote))); [width1 hight1]=size(xmax);%不允许多个极值同时存在 [width2 hight2]=size(ymax); if width1*hight1>1 xmax = xmax(1); end if width2*hight2>1 ymax = ymax(1); end area_num(nn) = area_num(nn) + 1; for i = 1:high-1; for j = 1:width-1; if two_d_center(i,j) == 1 if ((i-xmax).^2 + (j-ymax).^2) < 1000 flag(i,j,nn) = area_num(nn); two_d_vote(i,j) = 0; two_d_center(i,j) = 0; end end end end end end end %subplot(335);mesh(flag); title('圆心区域分区'); %=====================%圆心区域极值 area_num_max = max(area_num); area_max_x = zeros(area_num_max,flag_num); area_max_y = zeros(area_num_max,flag_num); area_min_x(:,:) = zeros(area_num_max,flag_num); area_min_x(:,:) = width; area_min_y(:,:) = zeros(area_num_max,flag_num); area_min_y(:,:) = high; n = 0; for nn = rmin:rmax if flag_r(nn) == 1 n = n + 1; %for n = 1:flag_num for m = 1:area_num(nn) for i = 1:high-1; for j = 1:width-1; if(flag(i,j,nn) == m)%如果区域标号相等 if area_max_x(m,n)<j area_max_x(m,n) = j; end if area_max_y(m,n)<i area_max_y(m,n) = i; end if area_min_x(m,n)>j area_min_x(m,n) = j; end if area_min_y(m,n)>i area_min_y(m,n) = i; end end end end %end end end end %=============================%圆心坐标 for n = 1:flag_num nn = value_r(n); for m = 1:area_num(nn) center_x(m,n) = (area_max_x(m,n) + area_min_x(m,n))/2; center_y(m,n) = (area_max_y(m,n) + area_min_y(m,n))/2; end end %=============================% %======================%画出检测的圆 subplot(224); imshow(pic1);hold on;title('检测结果'); for n = 1:flag_num nn = value_r(n); for m = 1:area_num(nn) flag_re = 0; for ii = 1: n mm = value_r(ii); for jj = 1:m %排除重复画圆的情况 if m ==1 && n == 1 draw_circle(center_x(m,n),center_y(m,n),value_r(n));axis([0 width 0 high]);%circle函数是调用的 else if (0.95*mm<nn && nn<1.05*mm) && (0.95*center_x(m,n)<center_x(jj,ii) && center_x(m,n)<center_x(jj,ii)*1.05) && (0.95*center_y(m,n)<center_y(jj,ii) && center_y(m,n)<center_y(jj,ii)*1.05) flag_re = flag_re + 1; end if flag_re < 2 && jj == m && ii == n draw_circle(center_x(m,n),center_y(m,n),value_r(n));axis([0 width 0 high]);%circle函数是调用的 end end end end hold on; end end %subplot(334);mesh(vote(:,:,value_r(1))); toc;%计时结束