hough变换提取椭圆

1、hough 变换提取椭圆任意椭圆方程表达式xa coscosb sinsinpya cossinb sincosq步骤： 1、读入图片（图片大小为256256 ），将图片二值化并提取边缘（用 edge）；2、设置一个 5维的参数空间并使初值为0；3、对边缘点集中的每一点进行hough 变换，如果 p、q在图像范围内，则 hough_space(p, q, a, b,)hough_space(p, q, a, b,)1；4、在参数空间中寻找超过阈值的p, q, a, b,，就是椭圆的参数；5、再对椭圆参数进行求精；（直接讲下面的程序复制到.m 文件里，再改 a、b、theta参数 运行）程序：椭

2、圆提取（图片大小256*256）f1=imread(' 你的图片 .JPG');%读入要提取椭圆的图片imshow(f1)f2=im2bw(f1,0.72);%二值化f3=im2uint8(f2);figureimshow(f3)BW=edge (f3,'canny',0,0.7,0.9);%提取边缘figureimshow(BW)n,m = size(BW);a_max=80;%a为椭圆长轴，范围需要自己根据图片中的椭圆形状确定a_min=70;b_max=50;%b为椭圆短轴b_min=40;step_a=2;%步长step_b=2;theta_min=-p

3、i/18;%theta倾斜角theta_max=pi/9;step_theta=pi/60;step_angle=1;yz =0.82;%阈值（自己定）size_a = round(a_max-a_min)/step_a)+1;size_b = round(b_max-b_min)/step_b)+1;size_theta=round(theta_max-theta_min)/step_theta)+1;size_angle = round(2*pi/step_angle);hough_space = zeros(m,n,size_a,size_b,size_theta); %设定 5维的参数

4、空间并使初值为 0cols,rows = find(BW);%find找出非零元素的索引ecount = size(rows);% Hough%图像空间 (x,y)对应到参数空间(a,b,p、 q 、 theta)% p = x-a*cos(angle)*cos(theta)+b*sin(angle)*sin(theta)% q = y-a*cos(angle)*sin(theta)-b*sin(angle)*cos(theta)%rows(i)行坐标fori=1:ecountfora=1:size_aforb=1:size_bfortheta=1:size_thetafork=1:size_

5、angle%hough 变换p =round(rows(i)-(a_min+(a-1)*step_a)*cos(k*step_angle)*cos(theta_min+(theta-1)*step_theta)+(b_min+(b-1)*step_b)*sin(k*step_angle)*sin(theta_min+(theta-1)*step_theta);q =round(cols(i)-(a_min+(a-1)*step_a)*cos(k*step_angle)*sin(theta_min+(theta-1)*step_theta)-(b_min+(b-1)*step_b)*sin(k*

6、step_angle)*cos(theta_min+(theta-1)*step_theta);if(p>0&p<=m&q>0&q<=n)hough_space(p,q,a,b,theta) =hough_space(p,q,a,b,theta)+1;endendendendendend% 搜索超过阈值的聚焦点max_para = max(max(max(max(max(hough_space);index = find(hough_space>max_para*yz); %find 找出 hough_space 中大于阈值的缩引并存入 i

7、ndexlength = size(index);hough_circle1=zeros(m,n); %确定为椭圆上的点的坐标 hough_circle2=zeros(m,n);%找出峰值对应的参数空间坐标fork=1:lengthpar5 =floor(index(k)-1)/(m*n*size_a*size_b)+1;%theta增量par4 =floor(index(k)-(par5-1)*(m*n*size_a*size_b)/(m*n*size_a)+1;%b增量par3 =floor(index(k)-(par5-1)*(m*n*size_a*size_b)-(par4-1)*(m

8、*n*size_a)/(m*n)+1;%a 增量par2 =floor(index(k)-(par5-1)*(m*n*size_a*size_b)-(par4-1)*(m*n*size_a)-(par3-1)*(m*n)/m)+1;%p增量par1 =index(k)-(par5-1)*(m*n*size_a*size_b)-(par4-1)*(m*n*size_a)-(par3-1)*(m*n)-(par2-1)*m;%q增量par5=theta_min+(par5-1)*step_theta;par4 = b_min+(par4-1)*step_b;par3 = a_min+(par3-1

9、)*step_a;theta(k)=par5;b(k)=par4;a(k)=par3;q(k)=par2;p(k)=par1;end%求出两圆参数平均值row1 col1=size(p);count=1;theta=sort(theta);p=sort(p);q=sort(q);a=sort(a);b=sort(b);THETA(count)=theta(1);P(count)=p(1);A(count)=a(1);B(count)=b(1);Q(count)=q(1);fort=1:1:col1ifabs(P(count)-p(t)<=10THETA(count)=(theta(t)+

10、THETA(count)/2;A(count)=(a(t)+A(count)/2;B(count)=(b(t)+B(count)/2;P(count)=(p(t)+P(count)/2;Q(count)=(q(t)+Q(count)/2;elsecount=count+1;THETA(count)=theta(t);A(count)=a(t);B(count)=b(t);P(count)=p(t);Q(count)=q(t);endendTHETAABPQ%绘制椭圆TYZ=zeros(1,2);TYY=zeros(1,2);ct_z=1;ct_y=1;fori=1:ecountifround(

11、rows(i)-P(1)*cos(THETA(1)+(cols(i)-Q(1)*sin(THETA(1)2/(A(1)2)+(-(rows(i)-P(1)*sin(THETA(1)+(cols(i)-Q(1)*cos(THETA(1)2/(B(1)2)<1.5.&round(rows(i)-P(1)*cos(THETA(1)+(cols(i)-Q(1)*sin(THETA(1)2/ (A(1)2)+(-(rows(i)-P(1)*sin(THETA(1)+(cols(i)-Q(1)*cos(THETA(1) 2/(B(1)2)>0.5TYY(ct_y,1)=rows(i);

12、TYY(ct_y,2)=cols(i);hough_circle1(cols(i),rows(i)=1;ct_y=ct_y+1;endifround(rows(i)-P(2)*cos(THETA(2)+(cols(i)-Q(2)*sin(THETA(2)2/( A(2)2)+(-(rows(i)-P(2)*sin(THETA(2)+(cols(i)-Q(2)*cos(THETA(2)2/(B(2)2)<1.5.&round(rows(i)-P(2)*cos(THETA(2)+(cols(i)-Q(2)*sin(THETA(2)2/ (A(2)2)+(-(rows(i)-P(2)*

13、sin(THETA(2)+(cols(i)-Q(2)*cos(THETA(2) 2/(B(2)2)>0.5TYZ(ct_z,1)=rows(i);TYZ(ct_z,2)=cols(i);hough_circle2(cols(i),rows(i)=1;ct_z=ct_z+1;endendfigureimshow(hough_circle1),title('?ì 2a?á 1?' )figureimshow(hough_circle2),title('?ì 2a?á 1?' )%分别计算两个椭圆的参数row_TYZ,co

14、l_TYZ=size(TYZ);fori1=1:1:row_TYZWTZ(i1,:)=TYZ(i1,1)2 TYZ(i1,1)*TYZ(i1,2) TYZ(i1,2)2 TYZ(i1,1)TYZ(i1,2) 1; ;endv_z1,d_z1=svd(WTZ'*WTZ);v_z1=vpa(v_z1,8)d_z1=double(d_z1);row_TYY,col_TYY=size(TYY);forj1=1:1:row_TYYWTY(j1,:)=TYY(j1,1)2 TYY(j1,1)*TYY(j1,2) TYY(j1,2)2 TYY(j1,1)TYY(j1,2) 1; ;endv_y1,d

15、_y1=svd(WTY'*WTY);v_y1=vpa(v_y1,8)d_y1=double(d_y1);在原图上绘制拟合椭圆，以查看拟合精度f1=imread(' 你的图片 .JPG');f2=im2bw(f1,0.72);f3=im2uint8(f2);BW=edge (f3,'canny',0,0.7,0.9);%绘制左边椭圆fori2=1:1:256forj2=1:1:256if0<double(.25527789e-4*i22-.62429860e-5*i2*j2+.52621684e-4*j22-.86887797e-2*i2-.78888347e-2*j2+.99993113).&double(.25527789e-4*i22-.62429860e-5*i2*j2+.52621684e-4*j22-.86887 797e-2*i2-.78888347e-2*j2+.99993113)<0.01f1(j2,i2)=0;f2(j2,i2)=0;BW(j2,i2)=1;endendend%绘制右边椭圆fori3=1:1:256forj3=1:1:256if0<double(.46540784e-4*i32+.19567656e-4*i3*j3+.10878161e-3*j32-.9180