如何进行圆拟合?

Question

我有近似的圆上一批点,想把他拟合成圆?如何做?找了一下,有一个matlab中的,可惜matlab不会,===================================
7)Matlab中如何作圆回归？
:#Peter Boettcher ([email protected]),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/1991n=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.

Answer 1

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

Answer 2

VC做可能比较麻烦，Matlab很方便啊，以前模糊数学的时候用过，把上面的function函数保存成
m文件，进行编译或者Matlab与VC混合编程的方式，网上也有很多
http://www.vckbase.com/document/viewdoc/?id=1435
http://www.simwe.com/jour/prog/p001008.htm

Answer 3

关于这个问题,我现在还不想用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)); 
/************************************************************/ 

Answer 4

这一个他这里我想是有错误的,我调试了一下没调试出来.其码这个就不对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*/ 是一样的.

Answer 5

还有这个:没调试过,是拟合曲线的:
//最小二乘法曲线拟合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;}

Answer 6

设已知的数据点为 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 求导数得到三个高次方程然后去接这个方程组得数值解，估计可行，就是数值解不打好求 :)

Answer 7

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

Answer 8

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