数学建模(9)——霍夫(Hough)变换圆检测

来源:互联网 发布:乐视手机移动网络设置 编辑:程序博客网 时间:2024/06/07 22:49

本篇文章主要是针对霍夫变换在数学建模2014A题嫦娥那道题的应用,借鉴网友的代码http://blog.csdn.net/huangquanming/article/details/7621557?locationNum=2&fps=1
1

针对这一阶段,从最小的高程值往上检测,获取某一高程的平面图,然后进行圆检测,再找到一个合适的位置作为预定着陆的中心位置

%tif文件处理clear,clcg=imread('fujian3.tif');     gg=double(g);             % 转为数值矩阵%gg=gg-1/255;              % 将彩色值转为 0-1 的渐变值[x,y]=size(gg);             % 取原图大小%[X,Y]=meshgrid(1:y,1:x); % 以原图大小构建网格%mesh(X,Y,gg);            % 高程图像bb=gg;k=80;    for i=1:x        for j=1:y            if bb(i,j)==k                bb(i,j)=1;            else                bb(i,j)=0;            end        end    endcircleParaXYR=[];  I = bb;  %figure(1),imshow(I),title('原图')  [m,n,l] = size(I);  if l>1      I = rgb2gray(I);  end  BW = edge(I,'sobel');  step_r = 1;  step_angle = 0.1;  minr = 3;  maxr = 30;  thresh = 0.51;  [hough_space,hough_circle,para] = hough_circle(BW,step_r,step_angle,minr,maxr,thresh);  figure(1),imshow(I),title('原图')  figure(2),imshow(BW),title('边缘')  figure(3),imshow(hough_circle),title('检测结果')  circleParaXYR=para;  %输出  fprintf(1,'\n---------------圆统计----------------\n');  [r,c]=size(circleParaXYR);%r=size(circleParaXYR,1);  fprintf(1,'  检测出%d个圆\n',r);%圆的个数  fprintf(1,'  圆心     半径\n');%圆的个数  for n=1:r  fprintf(1,'%d (%d,%d)  %d\n',n,floor(circleParaXYR(n,1)),floor(circleParaXYR(n,2)),floor(circleParaXYR(n,3)));  end  %标出圆  figure(4),imshow(I),title('检测出图中的圆')  hold on;  [row,~]=size(circleParaXYR); %plot(circleParaXYR(:,2), circleParaXYR(:,1), 'r+');   for k = 1 : size(circleParaXYR, 1)    t=0:0.01*pi:2*pi;    x=cos(t).*circleParaXYR(k,3)+circleParaXYR(k,2);y=sin(t).*circleParaXYR(k,3)+circleParaXYR(k,1);    plot(x,y,'r-');   end  picmap=zeros(2300,2300);for k = 1 : size(circleParaXYR, 1)    picmap(circleParaXYR(k,2),circleParaXYR(k,1))=1;end  % length(find(picmap==1))for k =1 : size(circleParaXYR, 1)    for t=0:0.01*pi:2*pi;    x=ceil(cos(t).*circleParaXYR(k,3)+circleParaXYR(k,2));  y=ceil(sin(t).*circleParaXYR(k,3)+circleParaXYR(k,1));    if x<1||x>2300      x=1;  end  if y<1||y>2300      y=1;  end  picmap(x,y)=1;  endend  R_r=150;%寻找中心着陆位置的圆的半径可以根据需要修改for i = R_r : size(picmap, 1)-R_r      for j=R_r : size(picmap, 1)-R_r          temp=0;        for m = i-R_r+1 : i+R_r              for n=j-R_r +1: j+R_r                  if sqrt((m-i).^2+(n-j).^2)<R_r                    if picmap(m,n)==1                    temp=1;                    break;                    end                end            end        end        if temp==0            mm=[i,j]%中心着陆位置        end    endend
function [hough_space,hough_circle,para] = hough_circle(BW,step_r,step_angle,r_min,r_max,p);  %[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之间的数  %-------------------------------------------------------------------  %------------------------------输出参数-----------------------------  % hough_space:参数空间,h(a,b,r)表示圆心在(a,b)半径为r的圆上的点数  % hough_circl:二值图像,检测到的圆  % para:检测到的圆的圆心、半径  %-------------------------------------------------------------------  % From Internet,Modified by mhjerry,2011-12-11  circleParaXYR=[];  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);%查找非零元素的行列坐标  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;%h(a,b,r)的坐标,圆心和半径              end          end      end  end  % 搜索超过阈值的聚集点。对于多个圆的检测,阈值要设的小一点!通过调此值,可以求出所有圆的圆心和半径  max_para = max(max(max(hough_space)));%返回值就是这个矩阵的最大值  index = find(hough_space>=max_para*p);%一个矩阵中,想找到其中大于max_para*p数的位置  length = size(index);%符合阈值的个数  hough_circle = false(m,n);  %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)) = true;   %检测的圆          end      end  end                 % 从超过峰值阈值中得到  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;      circleParaXYR = [circleParaXYR;par1,par2,par3];      hough_circle(par1,par2)= true; %这时得到好多圆心和半径,不同的圆的圆心处聚集好多点,这是因为所给的圆不是标准的圆      %fprintf(1,'test1:Center %d %d \n',par1,par2);  end  %集中在各个圆的圆心处的点取平均,得到针对每个圆的精确圆心和半径!  while size(circleParaXYR,1) >= 1      num=1;      XYR=[];      temp1=circleParaXYR(1,1);      temp2=circleParaXYR(1,2);      temp3=circleParaXYR(1,3);      c1=temp1;      c2=temp2;      c3=temp3;      temp3= r_min+(temp3-1)*step_r;     if size(circleParaXYR,1)>1            for k=2:size(circleParaXYR,1)        if (circleParaXYR(k,1)-temp1)^2+(circleParaXYR(k,2)-temp2)^2 > temp3^2           XYR=[XYR;circleParaXYR(k,1),circleParaXYR(k,2),circleParaXYR(k,3)];  %保存剩下圆的圆心和半径位置        else          c1=c1+circleParaXYR(k,1);        c2=c2+circleParaXYR(k,2);        c3=c3+circleParaXYR(k,3);        num=num+1;        end       end     end         %fprintf(1,'sum %d %d radius %d\n',c1,c2,r_min+(c3-1)*step_r);        c1=round(c1/num);        c2=round(c2/num);        c3=round(c3/num);        c3=r_min+(c3-1)*step_r;        %fprintf(1,'num=%d\n',num)        %fprintf(1,'Center %d %d radius %d\n',c1,c2,c3);           para=[para;c1,c2,c3]; %保存各个圆的圆心和半径的值        circleParaXYR=XYR;  end    
原创粉丝点击