球面拟合算法

最新推荐文章于 2022-10-11 11:52:16 发布

原创最新推荐文章于 2022-10-11 11:52:16 发布 · 3.4k 阅读

17 ·

CC 4.0 BY-SA版权

数据结构与算法专栏收录该内容

10 篇文章

订阅专栏

博客围绕一组分散在球面上且混杂噪声的离散数据点，探讨如何拟合球面参数。通过构造函数求使函数取得最小值的参数，经过一系列化简得到方程并求解。还给出了使用Eigen库及封装结构体的代码，最后展示了实验计算出的球心和半径。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

一组离散数据点分布在一个球的球面上的，该如何拟合这个球的参数?

简单描述一下问题：我们有一组数据：, 这些数据是分散在一个球面上的，当然数据中混杂着一些噪声。球面方程:

$(x-x_0)^2+(y-y_0)^2+(z-z_0)^2=R^2$

其中 $(x_0, y_0, z_0, R)$ 是我们需要求出的参数。

我们构造一个函数：

$H(x_0, y_0, z_0 , R)=\sum_{i=1}^{N}\left((x_i-x_0)^2+(y_i-y_0)^2+(z_i-z_0)^2-R^2\right)^2$

我们求的 $(x_0, y_0, z_0, R)$ ,是使 H 取得最小值的那组参数。也就是满足下面四个方程的.

$\frac{\partial H}{\partial x_0} =-4\sum _{i=1}^N \left(x_i-x_0\right) \left(\left(x_i-x_0\right){}^2+\left(y_i-y_0\right){}^2+\left(z_i-z_0\right){}^2-R^2\right)= 0\\ \frac{\partial H}{\partial y_0} = -4\sum _{i=1}^N \left(y_i-z_0\right) \left(\left(x_i-x_0\right){}^2+\left(y_i-y_0\right){}^2+\left(z_i-z_0\right){}^2-R^2\right)=0 \\ \frac{\partial H}{\partial z_0} = -4\sum _{i=1}^N \left(z_i-z_0\right) \left(\left(x_i-x_0\right){}^2+\left(y_i-y_0\right){}^2+\left(z_i-z_0\right){}^2-R^2\right)= 0\\ \frac{\partial H}{\partial R} = -4 \sum _{i=1}^N R \left(\left(x_i-x_0\right){}^2+\left(y_i-y_0\right){}^2+\left(z_i-z_0\right){}^2-R^2\right) = 0$
最后一个方程可以化简为： $\sum _{i=1}^N \left(\left(x_i-x_0\right){}^2+\left(y_i-y_0\right){}^2+\left(z_i-z_0\right){}^2-R^2\right) = 0$

利用这个式子可以把前三个方程化简为：

$\sum _{i=1}^N x_i \left(\left(x_i-{x_0}\right){}^2+\left(y_i-{y_0}\right){}^2+\left(z_i-{z_0}\right){}^2-R^2\right)= 0\\ \sum _{i=1}^N y_i \left(\left(x_i-{x_0}\right){}^2+\left(y_i-{y_0}\right){}^2+\left(z_i-{z_0}\right){}^2-R^2\right)=0 \\ \sum _{i=1}^N z_i \left(\left(x_i-{x_0}\right){}^2+\left(y_i-{y_0}\right){}^2+\left(z_i-{z_0}\right){}^2-R^2\right)= 0\\$

为了进一步化简，再做一次变换：

$u_i = x_i - \bar x \\ v_i = y_i - \bar y \\ w_i = z_i - \bar z$

这样替换之后有：

$\sum _{i=1}^N u_i \left(\left(u_i-{u_0}\right){}^2+\left(v_i-{v_0}\right){}^2+\left(w_i-{w_0}\right){}^2-R^2\right)= 0\\ \sum _{i=1}^N v_i \left(\left(u_i-{u_0}\right){}^2+\left(v_i-{v_0}\right){}^2+\left(w_i-{w_0}\right){}^2-R^2\right)=0 \\ \sum _{i=1}^N w_i \left(\left(u_i-{u_0}\right){}^2+\left(v_i-{v_0}\right){}^2+\left(w_i-{w_0}\right){}^2-R^2\right)= 0\\ \sum _{i=1}^N \left(\left(u_i-{u_0}\right){}^2+\left(v_i-{v_0}\right){}^2+\left(w_i-{w_0}\right){}^2-R^2\right) = 0$

再次化简后有：

$(\sum u_i^2 )u_0 +(\sum{u_i v_i})v_0+(\sum {u_i w_i}) w_0= \frac{\sum{(u_i^3+u_i v_i^2+u_i w_i^2)}}{2} \\ (\sum u_i v_i )u_0 +(\sum{ v_i^2})v_0+(\sum {v_i w_i}) w_0= \frac{\sum{(u_i^2 v_i+ v_i^3+v_i w_i^2)}}{2} \\ (\sum u_i w_i )u_0 +(\sum{v_i w_i})v_0+(\sum {w_i^2}) w_0= \frac{\sum{(u_i^2 w_i+v_i^2 w_i+ w_i^3)}}{2}$

解出这个方程就可以得到 $(u_0, v_0, w_0)$ , 然后就可以得到 $(x_0, y_0, z_0)$

了。之后带入第四个式子，就可以得到 $R$ 的值：

$R =\sqrt{ \frac{\sum _{i=1}^N\left(\left(u_i-{u_0}\right){}^2+\left(v_i-{v_0}\right){}^2+\left(w_i-{w_0}\right){}^2\right)} {N}}$

下面是代码：
1. 其中用到了Eigen这个第三方开源库。

2. 自己封装了三维点和标准球的结构体

#include <Eigen/Dense>

bool FittingSphere(const vector<XYZPoint> &vecPoint, CStandardSphere &StandardSphere)
{
	bool bRes = true;
	Eigen::Vector3d Y;
	Eigen::Matrix3d X;
	Y.setZero();
	X.setZero();
	int nCount = vecPoint.size();

	for (int i = 0; i < nCount; i++)
	{
		Y[0] += pow(vecPoint[i].dX, 3) + pow(vecPoint[i].dY, 2) * vecPoint[i].dX + pow(vecPoint[i].dZ, 2) * vecPoint[i].dX;
		Y[1] += pow(vecPoint[i].dX, 2) * vecPoint[i].dY + pow(vecPoint[i].dY, 3) + pow(vecPoint[i].dZ, 2) * vecPoint[i].dY;
		Y[2] += pow(vecPoint[i].dX, 2) * vecPoint[i].dZ + pow(vecPoint[i].dY, 2) * vecPoint[i].dZ + pow(vecPoint[i].dZ, 3);

		X(0, 0) += pow(vecPoint[i].dX, 2);
		X(0, 1) += vecPoint[i].dX * vecPoint[i].dY;
		X(0, 2) += vecPoint[i].dX * vecPoint[i].dZ;
		X(1, 1) += pow(vecPoint[i].dY, 2);
		X(1, 2) += vecPoint[i].dY * vecPoint[i].dZ;
		X(2, 2) += pow(vecPoint[i].dZ, 2);
	}
	//对称矩阵
	X(1, 0) = X(0, 1);
	X(2, 0) = X(0, 2);
	X(2, 1) = X(1, 2);
	//判断是否可逆
	double dX = 0.0, dY = 0.0, dZ = 0.0;
	double dRadius = 0.0;
	if (X.determinant() == 0)
	{
		bRes = false;
	}
	else
	{
		Eigen::Matrix3d XI = X.inverse();
		for (int ni = 0; ni < XI.RowsAtCompileTime; ni++)
		{
			dX += XI(0, ni) * Y[ni];
			dY += XI(1, ni) * Y[ni];
			dZ += XI(2, ni) * Y[ni];
		}
		dX = dX / 2.0;
		dY = dY / 2.0;
		dZ = dZ / 2.0;
		for (int i = 0; i < nCount; i++)
		{
			dRadius += (pow(vecPoint[i].dX - dX, 2) + pow(vecPoint[i].dY - dY, 2) + pow(vecPoint[i].dZ - dZ, 2));
		}
		dRadius = pow(dRadius / nCount, 0.5); 		
	}
	StandardSphere.SetX(dX);
	StandardSphere.SetY(dY);
	StandardSphere.SetZ(dZ);

	StandardSphere.SetRadius(dRadius);
	return bRes;
}

实验设定的球心是（-60，25，0），半径是12.5
计算出来的球心是：（-60.003079，24.999393，0.002547），半径12.500950