100分，求思路！由点得到拟合球

[jianxinlee]

现有离散空间点若干，如何得到拟合意义上的球，即要求算法给出球心坐标和半径！

[jianxinlee]

关于提到“拟合意义”，是我自己感觉这个问题就像常见的直线拟合或者多项式拟合。

当然一个很显然的思路是:（x-x0）^2 +（y-y0）^2 +（z-z0）^2 = R^2, 然后对R求和S,构造S达到最小时的(x0, y0, z0)。但是这个构造过程如何进行呢？

或者那位牛牛有更好的办法，感激不尽！

[saint001]

如果进行最小二乘拟合，就是求这些点的坐标平均值

[saint001]

问题不是你说的对R求和，使和最小，如果是这个意思，平均值就差不多了(对R^2求和最小)
而是要求Ri的偏差(方差)比较小

[jianxinlee]

saint001，应该你是对的，我的描述有问题！

但是这个Ri怎么构造呢？球心半径都是待求量

[saint001]

很简单的问题
昨天把它想复杂了

直接对x^2+y^2+z^2+ax+by+cz+d=0拟合就行了
拟合的四个系数是a,b,c,d，最小二乘法
本来很简单，可是昨天局限在了把球心和半径作为自变量上

有了a,b,c,d球心坐标和半径也就知道了

[saint001]

前些天编写的最小二乘拟合的程序，这几天就用了好几次，从d:\CircleFit.txt中读取数据
#include "math.h"
#include "stdio.h"
#define ABS(x) (x)>0?(x):-(x)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

void solve(double **a,double *b,double *x,int n)
{
int i,j,k,ik;
double mik,temp;
for(k=0;k<n;k++)
{
mik=-1;
for(i=k;i<n;i++)
if(ABS(a[i][k])>mik)
{
mik=ABS(a[i][k]);
ik=i;
}
for(j=k;j<n;j++)
SWAP(a[ik][j],a[k][j]);
SWAP(b[k],b[ik]);
b[k]/=a[k][k];
for(i=n-1;i>=k;i--)
a[k][i]/=a[k][k];
for(i=k+1;i<n;i++)
{
b[i]-=a[i][k]*b[k];
for(j=n-1;j>=k;j--)
a[i][j]-=a[i][k]*a[k][j];
}
}
for(i=n-1;i>=0;i--)
{
x[i]=b[i];
for(j=i+1;j<n;j++)
x[i]-=a[i][j]*x[j];
}
}

void linear(double **x,double *y,double *beta,int n,int p)
{
double **a,*b;
int i,j,k;
a=new double*[p];
for(i=0;i<p;i++)
a[i]=new double[p];
for(i=0;i<p;i++)
for(j=0;j<p;j++)
{
a[i][j]=0;
for(k=0;k<n;k++)
a[i][j]+=x[k][i]*x[k][j];
}
b=new double[p];
for(i=0;i<p;i++)
{
b[i]=0;
for(j=0;j<n;j++)
b[i]+=x[j][i]*y[j];
}
solve(a,b,beta,p);
for(i=0;i<p;i++)
delete a[i];
delete a,b;
}

void main()
{
int i,n,p=4;
double x0,y0,z0,r;
double **x,*y,beta[4];

FILE *fp=fopen("d:\\CircleFit.txt","r");
fscanf(fp,"%d",&n);
x=new double*[n];
y=new double[n];
for(i=0;i<n;i++)
{
x[i]=new double[4];
fscanf(fp,"%lf",x[i]);
fscanf(fp,"%lf",x[i]+1);
fscanf(fp,"%lf",x[i]+2);
x[i][3]=1;
y[i]=-x[i][0]*x[i][0]-x[i][1]*x[i][1]-x[i][2]*x[i][2];
}
fclose(fp);
linear(x,y,beta,n,p);
x0=-beta[0]/2;
y0=-beta[1]/2;
z0=-beta[2]/2;
r=sqrt(x0*x0+y0*y0+z0*z0-beta[3]);
printf("%f\t%f\t%f\t%f\t\n",x0,y0,z0,r);
}

[saint001]

刚才一段时间上不去论坛
d:\CircleFit.txt中的数据
20
2.8708 1.3399 3.2936
1.3901 0.9456 4.7586
-0.2118 1.4578 4.4892
-0.7987 1.6279 3.5585
2.3814 0.9131 3.9640
2.4616 0.8509 2.2743
-0.7898 2.6383 3.6532
2.5070 3.3155 3.3500
-0.6203 0.9506 3.3615
-0.8986 2.4616 2.6251
2.5346 3.2571 2.5281
2.5465 1.8771 4.4302
3.0246 2.2083 3.3029
3.0213 2.4154 2.6920
0.1899 2.2169 4.8677
2.6968 1.4048 2.1883
-0.7191 3.0063 2.8426
2.7501 1.5074 2.1484
2.8091 2.7637 2.3159
-0.9991 2.1089 3.0588
数据是用matlab产生的球心为1，2，3，半径为2的球面上的随机点，并有随机误差(0.1以内)，所以拟合出来球心近似是[1，2，3]，半径应该近似为2
程序输出结果是
1.039589 2.051142 3.058152 2.009236
前三项为球心坐标，后一项为半径

[jianxinlee]

强人，呵呵，谢了，这就结贴！