如何进行圆拟合?

[wrcluomo]

我有近似的圆上一批点,想把他拟合成圆?如何做?找了一下,有一个matlab中的,可惜matlab不会,

===================================
7)Matlab中如何作圆回归？
:#Peter Boettcher ([e-mail]),2002/5/16, comp.soft-sys.matlab#

Q5.5: How can I fit a circle to a set of XY data?
=================================================

An elegant chunk of code to perform least-squares circle fitting
was written by Bucher Izhak and has been floating around the
newgroup for some time. The first reference to it that I can
find is in:

function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.A is an
% optional output describing the circle's equation:
%
% x^2+y^2+a(1)*x+a(2)*y+a(3)=0

% by Bucher izhak 25/oct/1991

n=length(x); xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy)...
sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R = sqrt((a(1)^2+a(2)^2)/4-a(3));

Tom Davis provided a more sophisticated approach that works
for more cases in and Code included.

[syy64]

我提供的思路请参考：
1、首先根据这批点寻找一个闭合多边形；
2、计算多边形的面积，将这个面积作为圆的面积，这样可以计算出圆的半径，半径是固定的，改变的是圆心的位置，
3、找出多边形中距离最长的两点，以这两点的中点为圆心；
4、计算该多边形与圆周的交点，这样圆周和多边形共同可能形成多个多边形；
5、计算这些多边形的面积；
6、在圆外的多边形面积记为负，在内为正；
7、分别累加正负面积，计算正负面积的绝对值之差，沿着最长距离按域值改变圆心的位置，差值最小者就是要找的圆心。

[wrcluomo]

TO syy64:我觉着你这个方法有一定的可行性,我以前也写过一个类似的算法,但这都不是拟合,只能说是改善.因为你这个算法是把半径固定找圆心.其实在拟合的过程呀 x0,y0,r三个都是变量.固定任何一个都可能达不到目的.
我觉得这个可行,可惜数学学得不好,都忘了:
三参数逆合, 是圆心坐标x,y以及半径r
求下面公式
N
E | (xi - x)^2 + (yi-y)^2 - r*r |
i=0
就是求到预计的圆的距离的和
让这个数值最小
就是把上述表达式对x,y和r分别求偏导, 使得其数值为0, 从而求解获得x,y,r的数值
我求的大家看看对不对,之后如何做?
对X:2*x-2*x(i)+x(i)^2+y^2-2*y(i)*y+y(i)^2-r^2

[syy64]

用数学方法来解决空间问题，主要是看最终效果，和要达到的目的，楼主提出的距离最小法也是一种效果，如果楼主追求的是这种效果，那你就达到目的了。

[wrcluomo]

唉,不是我提出的,是人家提出的啊.

[syy64]

我提供的思路请参考：
1、首先根据这批点寻找一个闭合多边形；
2、找出闭合多边形的重心，
3、以某点为半径；
4、计算该多边形与圆周的交点，这样圆周和多边形共同可能形成多个多边形；
5、计算这些多边形的面积；
6、在圆外的多边形面积记为负，在里为正；
7、累加这些面积，如累加后面积为正，就按一定域值扩大半径，否则缩小半径；
8、重复1--7步，知道累加值为最小。

[wrcluomo]

谢谢SYY64,但是你的方法感觉只是找最佳半径,这个多边形的重心可不一定是所找的圆心呀.

[折腾_苏州]

VC做可能比较麻烦，Matlab很方便啊，以前模糊数学的时候用过，把上面的function函数保存成
m文件，进行编译

或者Matlab与VC混合编程的方式，网上也有很多
http://www.vckbase.com/document/viewdoc/?id=1435
http://www.simwe.com/jour/prog/p001008.htm

[tiaoci]

这个不是很难吧

[booklove]

连接距离最远的点，作中垂线，形成许多交点，圆心在交点最密集的地方。
求圆心到点的距离，取中间的80%，求平均差不多是半径吧。

[jerry]

以前做过多项式拟合. 应该差不多吧.
就是求点到圆的均方差最小. 然后解线性议程组

看看数值计算方法-> 最小二乘拟合算法

[syy64]

我上面的提法不够严密，应该是这样，先计算原始多边形的面积，将这个面积作为圆的面积，这样可以计算出圆的半径，半径是固定的，改变的是圆心的位置；
所以上面第7步应该这样：
7、分别累加正负面积，计算正负面积绝对值之差，....
吃饭了下午再想想；

[wrcluomo]

关于这个问题,我现在还不想用MATLAB解决.因为我不只是拟合圆,还要得到他的一些参数.我查了许多关于拟合方面的,但圆拟合的少呀.我贴两个大家讨论一下:
BBS水木清华站∶精华区
发信人: FangQ (诗经·小雅·鹤鸣), 信区: MathTools
标 题: Re: 求用圆拟合数据点的方法 (转载)
发信站: BBS 水木清华站 (Mon Mar 15 22:53:56 1999)

/***********************************************/
/*以前我写的一个小程序中的一段，你可以看看。*/
/*圆回归是在西电梁昌洪的计算微波一书中有公式*/
/***********************************************/

double Sum(double (*D)[][2],int num,int flag)
{
int i;
double SUM=0;
double ID;

for(i=0;i<num;i++)
{
switch(flag)
{
case 1:ID= (*D)[i][0]; break; /*X*/
case 2:ID= (*D)[i][1]; break; /*Y*/
case 3:ID= (*D)[i][0]*(*D)[i][0]; break; /*X*X*/
case 4:ID= (*D)[i][1]*(*D)[i][1]; break; /*Y^2*/
case 5:ID= (*D)[i][0]*(*D)[i][1]; break; /*XY*/
case 6:ID= (*D)[i][0]*(*D)[i][1]*(*D)[i][1];break; /*X*Y^2*/
case 7:ID= (*D)[i][1]*(*D)[i][0]*(*D)[i][0];break; /*Y*X^2*/
case 8:ID= (*D)[i][0]*(*D)[i][0]*(*D)[i][0];break; /*X^3*/
case 9:ID= (*D)[i][1]*(*D)[i][1]*(*D)[i][1];break; /*Y^3*/
case 10:ID= (*D)[i][0]*(*D)[i][0]; break; /*X^2*/
}
SUM+=ID;
}
return SUM;
}
/*以下是在主函数中,DATA为数据点集,Xc,Yc,R为回归圆参数*/
/************************************************************/
Xc=(((Sum(DATA,MAXNUM,8)+Sum(DATA,MAXNUM,6))*Sum(DATA,MAXNUM,4))-\
((Sum(DATA,MAXNUM,7)+Sum(DATA,MAXNUM,9))*Sum(DATA,MAXNUM,MAXNUM)))\
/(2*((Sum(DATA,MAXNUM,10)*Sum(DATA,MAXNUM,4)+Sum(DATA,MAXNUM,3))*\
Sum(DATA,MAXNUM,3)));
Yc=(((Sum(DATA,MAXNUM,7)+Sum(DATA,MAXNUM,9))*Sum(DATA,MAXNUM,10))-\
((Sum(DATA,MAXNUM,8)+Sum(DATA,MAXNUM,7))*Sum(DATA,MAXNUM,3)))\
/(2*((Sum(DATA,MAXNUM,10)*Sum(DATA,MAXNUM,4)+Sum(DATA,MAXNUM,3))*\
Sum(DATA,MAXNUM,3)));

for(i=0;i<MAXNUM;i++)
R+=(pow((DATA[i][0]-Xc),2)+pow((DATA[i][1]-Yc),2));
/************************************************************/

[wrcluomo]

这一个他这里我想是有错误的,我调试了一下没调试出来.其码这个就不对Sum(DATA,MAXNUM,MAXNUM))),还有这个: case 3:ID= (*D)[i][0]*(*D)[i][0]; break; /*X*X*/
和 case 10:ID= (*D)[i][0]*(*D)[i][0]; break; /*X^2*/ 是一样的.

[wrcluomo]

还有这个:没调试过,是拟合曲线的:
//最小二乘法曲线拟合typedef CArray<double,double>CDoubleArray;BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,CDoubleArray *A){ //X,Y -- X,Y两轴的坐标 //M -- 结果变量组数 //N -- 采样数目 //A -- 结果参数 register long i,j,k; double Z,D1,D2,C,P,G,Q; CDoubleArray B,T,S; B.SetSize(N); T.SetSize(N); S.SetSize(N); if(M>N)M=N; for(i=0;i<M;i++) (*A)[i]=0; Z=0; B[0]=1; D1=N; P=0; C=0; for(i=0;i<N;i++) { P=P+(*X)[i]-Z; C=C+(*Y)[i]; } C=C/D1; P=P/D1; (*A)[0]=C*B[0]; if(M>1) { T[1]=1; T[0]=-P; D2=0; C=0; G=0; for(i=0;i<N;i++) { Q=(*X)[i]-Z-P; D2=D2+Q*Q; C=(*Y)[i]*Q+C; G=((*X)[i]-Z)*Q*Q+G; } C=C/D2; P=G/D2; Q=D2/D1; D1=D2; (*A)[1]=C*T[1]; (*A)[0]=C*T[0]+(*A)[0]; } for(j=2;j<M;j++) { S[j]=T[j-1]; S[j-1]=-P*T[j-1]+T[j-2]; if(j>=3) { for(k=j-2;k>=1;k--) S[k]=-P*T[k]+T[k-1]-Q*B[k]; } S[0]=-P*T[0]-Q*B[0]; D2=0; C=0; G=0; for(i=0;i<N;i++) { Q=S[j]; for(k=j-1;k>=0;k--) Q=Q*((*X)[i]-Z)+S[k]; D2=D2+Q*Q; C=(*Y)[i]*Q+C; G=((*X)[i]-Z)*Q*Q+G; } C=C/D2; P=G/D2; Q=D2/D1; D1=D2; (*A)[j]=C*S[j]; T[j]=S[j]; for(k=j-1;k>=0;k--) { (*A)[k]=C*S[k]+(*A)[k]; B[k]=T[k]; T[k]=S[k]; } } return TRUE;}

[tiaoci]

设已知的数据点为 P[i] (i = 1 - N)

圆心坐标为 O(x, y) 半径 R

于是圆方程为 (X - xo)^2 + (Y - yo)^2 - R^2 = 0

或者用向量表示为 |P[i] - O|^2 - R^2 = 0

建立每个数据点的误差函数为 f = |P[i] - O|^2 - R^2

总误差函数 F = Sigma[(|P[i] - O|^2 - R^2)^2] (i = 1 - N)

对总误差函数的O点的 x, y 和 R 求导数得到三个高次方程

然后去接这个方程组得数值解，估计可行，就是数值解不打好求 :)

[wrcluomo]

http://topic.csdn.net/t/20050113/17/3723480.html
还有这个讨论贴

[wrcluomo]

UP

[ReverseEngineering]

好难的问题。