希望可以加分，这能鼓励我们对生活工作中遇到的问题应用matlab来解决。

在工作中遇到这样的一个问题：当用三坐标测量圆特征，得到平面上离散分布的若干点，这些点近似分布在圆周上。如何评价该圆特征的圆度？

圆度误差的评定方法有4种：最小包容区域法，最小外接圆法，最大内切圆法，最小二乘法。最小包容区域法是符合圆度定义的方法，其他方法实际上是不符合圆度定义的，但由于最小包容区域法很难实现，故只好用其他方法来近似代替。

假设平面上有若干个测量点，大概在它们的中间随便找一点，以该点为中心创建一个矩形区域，该矩形区域内创建出均匀分布的很多点，充满矩形区域。然后计算这些构造的“中心点”到周围测量点之间的距离。对某一个构造的“中心点”来说，总可以找到距离它最远和最近的两个测量点，记下这它们的距离差，认为该距离差是圆度的一个备选值。对于矩形框内所有的中心点来说，必然有一个中心点(假设为o1)对应着最小的一个距离差。理论上，只要矩形框内构造的“中心点”足够密，则一定可以找到符合精度的圆度值。

但如果圆特征较大，按照上面所说，可能计算就要花费很多时间。于是对该方案进行了修改，先创建一个较疏的矩形点阵，找到o1后，再以o1为中心，创建另一个矩形框，该矩形框尺寸比第一个的减小一半，但里面点数目不变，这样点的密度就增加到原来的4倍。由于点数目没变，比较少，所以计算速度应该比较快。而且由于点的密度增加了，所找到的“中心点”也更精确。

根据需要设定循环次数。可以找到符合精度要求的中心点，以及它对应的圆度误差值。

以下是matlab程序：经过测试，所得到的圆度值随着循环次数n的增大而减小，而且下降速度越来越慢。可以认为，最终收敛于真正的圆度值。附件将圆度值与循环次数n的关系用图表显示出来，比较直观。

% 求圆特征的圆度及相应的中心点

xi=input('输入测量点的横坐标:');

yi=input('输入测量点的纵坐标:');  %yi与xi长度相等

xmin=min(xi);

xmax=max(xi);

ymin=min(yi);

ymax=max(yi);

xoc=0.5;  %可以修正

yoc=0.5;  %可以修正

xco=xmin+xoc*(xmax-xmin);  %设定的圆环中心点群的中心点的横坐标

yco=ymin+yoc*(ymax-ymin);  %设定的圆环中心点群的中心点的纵坐标

x_range=(xmax-xmin)/2;

y_range=(ymax-ymin)/2;

incr_x=x_range/10;  %x轴上有20个刻度

incr_y=y_range/10;  %y轴上有20个刻度

xo=[xco-x_range:incr_x:xco+x_range];  %设定的中心搜索点的横坐标组成的向量

yo=[yco-y_range:incr_y:yco+y_range];  %设定的中心搜索点的纵坐标组成的向量

R=[];

for n=1:10  %循环次数，可以修正

for j=1:length(xo)  %某个搜索点在x轴上的索引号

for k=1:length(yo)  %该搜索点在y轴上的索引号

for i=1:length(xi)  %测量点索引号

R(j,k,i)=sqrt((xi(i)-xo(j)).^2+(yi(i)-yo(k)).^2);

end

rmax(j,k)=max(R(j,k,:));

rmin(j,k)=min(R(j,k,:));

delta_r(j,k)=rmax(j,k)-rmin(j,k);

end

end;

[delta_rmi,row]=min(delta_r);

[delta_rmin,column]=min(delta_rmi);

x_range=x_range*0.5;  %搜索区域减小到原来的25%,可以修正

y_range=y_range*0.5;  %搜索区域减小到原来的25%,可以修正

n=n+1;

incr_x=incr_x/2;  %x方向密度增加到原来的2倍

incr_y=incr_y/2;  %y方向密度增加到原来的2倍

xco=xo(row(column));

yco=yo(column);

xo=[xco-x_range:incr_x:xco+x_range];

yo=[yco-y_range:incr_y:yco+y_range];

end;

x=xco

y=yco

for i=1:length(xi)

r(i)=sqrt((xi(i)-x)^2+(yi(i)-y)^2);

end

r_o=max(r)

r_i=min(r)

delta_radi=r_o-r_i