利用多个点拟合圆

给定多个点，根据最小二乘原理，可以拟合出一个圆形，拟合的精度很高，

#include
#include
#include 
#include 
int N;
double abr*z=NULL;

double Average（double A[]int n）
{
	int i;
	double tmp=0;
	for （i=0;i	tmp=tmp/（double）n;
	return tmp;
}

void AccuCircleSim（double x[]double y[]int n）
{
    int i=0;
    double X1Y1X2Y2Z1X3Y3Z2Z3;
    double p;
   	double *tp=NULL;
	tp= （double*）malloc（（unsigned）（n*sizeof（double）））;
    X1=Average（xn）;
	Y1=Average（yn）;
	
	for（i=0;i	Y2=Average（tpn）;

    for（i=0;i	X2=Average（tpn）;
        
    for（i=0;i	Z1=Average（tpn）;

	for（i=0;i	Z2=Average（tpn）;
    
    for（i=0;i	X3=Average（tpn）;
    
	for（i=0;i	Y3=Average（tpn）;

	for（i=0;i	Z3=Average（tpn）;

	p=（double）2*（X1*X1-X2）*（Y1*Y1-Y2）-（double）2*（X1*Y1-Z1）*（X1*Y1-Z1）;
	
	a=（（X2*X1+X1*Y2-X3-Z3）*（Y1*Y1-Y2）-（X2*Y1+Y1*Y2-Z2-Y3）*（X1*Y1-Z1））/p;
	
	b=（（X2*Y1+Y1*Y2-Z2-Y3）*（X1*X1-X2）-（X2*X1+X1*Y2-X3-Z3）*（X1*Y1-Z1））/p;	
    r=sqrt（a*a-（double）2*X1*a+b*b-（double）2*Y1*b+X2+Y2）;
    
	free（tp）;
	}

void SndCircleSim（double x[]double y[]int n）
{
    int i=0;
    double u2v2u3u4u5;
    double p;
   	double *tp=NULL*u=NULL*v=NULL*w=NULL;
	u= （double*）malloc（（unsigned）（n*sizeof（double）））;
	v= （double*）malloc（（unsigned）（n*sizeof（double）））;
	w= （double*）malloc（（unsigned）（n*sizeof（double）））;
    	
	for（i=0;i		{
		 u[i]=x[i+1]-x[i];
		 v[i]=y[i+1]-y[i];
		 w[i]=（x[i+1]*x[i+1]-x[i]*x[i]+y[i+1]*y[i+1]-y[i]*y[i]）/（double）2;
		}
	tp= （double*）malloc（（unsigned）（n*sizeof（double）））;
	
	for（i=0;i	u2=Average（tpn-1）;
    for（i=0;i	v2=Average（tpn-1）;       
    for（i=0;i	u3=Average（tpn-1）;
	for（i=0;i	u4=Average（tpn-1）;
	for（i=0;i	u5=Average（tpn-1）;
    
	p=u2*v2-u3*u3;
	a=（v2*u4-u3*u5）/p;
	b=（u2*u5-u3*u4）/p;
	r=0;
	for（i=0;i	r=r/（double）n;

	free（u）;
	free（v）;
	free（w）;
	free（tp）;
	}

//===========求矩阵乘积

double * MatrixProduct（double A[]double B[]int lint mint n）
    {
	 int ijk;
	 double *C=NULL;
     C=（double *）malloc（m*m*sizeof（double））;
     for（i=0;i        for（j=0;j	 for（i=0;i	    for（j=0;j			for（k=0;k     return C;
}

double Surplus（double A[]int mint n） //求矩阵行列式
    {
     int ijkpr;
     double Xtemp=1.0temp1=1.0s=0s1=0;

     if（n==2） X=A[0]*A[3]-A[1]*A[2];
     else
	 {
		for（k=0;k		{
			for（i=0j=k;i         if（m>i）
			{
			 for（p=m-ir=m-1;p>0;p--r--）temp*=A[r*n+p-1];
			}
         s+=temp;
        temp=（double）1;
		}
		for（k=n-1;k

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       5976  2004-02-15 10:46  拟和圆程序\CircleLS.c

     文件       3316  2004-02-13 16:56  拟和圆程序\circles.c

     文件       1450  2004-02-10 14:03  拟和圆程序\CircleSim.c

     文件      10218  2008-01-23 15:56  拟和圆程序\CircleSimD.c

     文件       3451  2007-06-09 11:30  拟和圆程序\CircleSimD.dsp

     文件        543  2007-01-19 16:37  拟和圆程序\CircleSimD.dsw

     文件      50176  2008-04-08 17:39  拟和圆程序\CircleSimD.ncb

     文件     245760  2008-04-08 17:39  拟和圆程序\CircleSimD.opt

     文件       1383  2008-03-01 16:02  拟和圆程序\CircleSimD.plg

     文件     208954  2008-03-01 16:02  拟和圆程序\Debug\CircleSimD.exe

     文件      21204  2008-03-01 16:02  拟和圆程序\Debug\CircleSimD.obj

     文件     492544  2008-03-01 16:02  拟和圆程序\Debug\CircleSimD.pdb

     文件          0  2008-03-01 16:02  拟和圆程序\Debug\CircleSimD.sbr

     文件      45056  2008-03-01 16:02  拟和圆程序\Debug\vc60.pdb

     目录          0  2008-06-14 11:23  拟和圆程序\Debug

     目录          0  2008-06-14 11:23  拟和圆程序

----------- ---------  ---------- -----  ----

              1090031                    16