实验科学离不开拟合，考虑到二次曲面拟合中球拟合和椭球拟合较常使用，并且球拟合有简单的最小二乘方法，因此这里稍微整理一下二者拟合相关的数学工具与处理思路。

球拟合

处理球的拟合时不会遇到处理椭球问题中的非线性问题，可以直接使用最小二乘法来得出结果。

三维空间中的任意一个球面可以表示成：

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

对该式展开并移相得到：

$$ x^2+y^2+z^2=2x_0x+2y_0y+2z_0z+(r^2-x_0^2-y_0^2-z_0^2) $$

不妨记：

$$ A = \begin{pmatrix} 2x_1 & 2y_1 & 2z_1 & 1 \\ 2x_2 & 2y_2 & 2z_2 & 1 \\ \vdots & \vdots & \vdots & \vdots \\ 2x_n & 2y_n & 2z_n & 1 \\ \end{pmatrix}, \vec{b} = \begin{pmatrix} x_0 \\ y_0 \\ z_0 \\ r^2 - x_{0}^2 - y_{0}^2 - z_{0}^2 \end{pmatrix}, \vec{Y} = \begin{pmatrix} x_1^2 + y_1^2 + z_1^2 \\ x_2^2 + y_2^2 + z_2^2 \\ \vdots \\ x_{n}^2 + y_{n}^2 + z_{n}^2 \end{pmatrix} $$

其中$(x_1,y_1,z_1),(x_2,y_2,z_2),...,(x_n,y_n,z_n)$分别是$n$个散点的坐标。这样问题就转化为了求下述超定方程的解：

$$ A\vec{b}=\vec{Y} $$

由于这是一个线性问题，我们可以直接写出该超定方程的解：

$$ \vec{b_0}=(A^\intercal A)^{-1}A^\intercal Y $$

椭球拟合

三维坐标变换

现在我们考虑一个空间中的任意椭球，为了方便计算误差函数我们需要对它进行坐标变换，使其中心与坐标原点对齐；假设该变换形式为：先平移一个向量$\vec{s}=-(x_0,y_0,z_0)^\intercal$，后分别绕$z,y,x$轴旋转$\theta_z,\theta_y,\theta_x$。将其写成矩阵形式为：

$$ \vec{x'}=R_xR_yR_z(\vec{x}+\vec{s})\tag{2.1} $$

此时椭球方程可以转换为标准形式：

$$ \frac{x'^2}{a^2}+\frac{y'^2}{b^2}+\frac{z'^2}{c^2}=1 $$

在这种情况下，我们可以简单定义最小距离平方和为：

$$ S_{err}^2=\sum_i\frac{1}{\frac{x_i'^2}{a^2}+\frac{y_i'^2}{b^2}+\frac{z_i'^2}{c^2}} $$

椭球曲面拟合

求解使$S_{err}^2(\theta_x,\theta_y,\theta_z,a,b,c,x_0,y_0,z_0)$最小的拟合曲面解$\vec{s_0}=(\theta_x,\theta_y,\theta_z,a,b,c,x_0,y_0,z_0)^\intercal$是一个非线性问题，无法使用常规线性代数方法解决，不过由于函数值在任意处数值可求，因此在计算机上可以使用高斯-牛顿方法、梯度下降等方法进行求解得到区域极小值；

例如我们可以用一个试解用$\vec{s}$代入误差函数，然后在试解$\vec{s}$附近对误差函数求方向导数$\vec{p}$，然后用新的试解$\vec{p}+\vec{s}$作为下一次尝试，直到两次试解的误差函数相差小于某个特定的值退出计算。

这种方法会遇到两种困难：

1、参数$\vec{s}$和对应坐标的变换表示形式非常复杂，并且目前的数学工具下暂时找不到更加简洁、有效的表达方式；

2、该方法可能会遇到数值困难(在椭球向球体退化的时候)，点偏差会造成角度的巨大变化(虽然我们可以把这种情况单独判断，并把对应最小二乘法的Jacobi矩阵秩从9改成6)

我们考虑另一种二次曲面的一般表示形式：

$$ Ax^2+By^2+Cz^2+Dxy+Exz+Fyz+Gx+Hy+Kz+L=0\tag{2.4} $$

该方程中L为非零时，我们可以在两边同时除以L，实际上是一个九元方程。

我们将式(2.1)做适当改写：

$$ (R_1R_2R_3)^{-1}\begin{pmatrix} x'\\y'\\z' \end{pmatrix}= \begin{pmatrix} x\\y\\z \end{pmatrix}+ \begin{pmatrix} -x_0\\-y_0\\-z_0 \end{pmatrix} $$

将$(x',y',z')^\intercal$写成极坐标形式：

$$ \begin{pmatrix} x\\y\\z \end{pmatrix}= (R_1R_2R_3)^{-1}\begin{pmatrix} a\cos u\cos v\\b\sin u\cos v\\c\sin v \end{pmatrix}+ \begin{pmatrix} x_0\\y_0\\z_0 \end{pmatrix} $$

则标准椭球方程可以写成：

$$ \begin{aligned} \frac{1}{a^2}&[(x-x_0)\cos\theta_z\cos\theta_y+(y-y_0)(\sin\theta_z\cos\theta_x+\cos\theta_z\sin\theta_y\sin\theta_x)+\\ &(z-z_0)(\sin\theta_z\sin\theta_x-\cos\theta_z\sin\theta_y\cos\theta_z)]^2+\\ \frac{1}{b^2}&[-(x-x_0)\sin\theta_z\cos\theta_y+(y-y_0)(\cos\theta_x\cos\theta_z-\sin\theta_x\sin\theta_y\sin\theta_z)+\\ &(z-z_0)(\cos\theta_z\sin\theta_x+\sin\theta_y\sin\theta_z\cos\theta_x)]^2+\\ \frac{1}{c^2}&[(x-x_0)\sin\theta_y-(y-y_0)\cos\theta_y\sin\theta_x+(z-z_0)\cos\theta_y\cos\theta_x]^2=1 \end{aligned} $$

对比(2.4)式，我们可以得到：

$$ \begin{aligned} A=&\frac{\cos^2\theta_y\cos^2\theta_z}{a^2}+\frac{\cos^2\theta_y\sin^2\theta_z}{b^2}+\frac{\sin^2\theta_y}{c^2}\\ B=&\frac{(\cos\theta_x\sin\theta_z+\sin\theta_x\sin\theta_y\cos\theta_z)^2}{a^2}+ \frac{(\cos\theta_x\cos\theta_z-\sin\theta_x\sin\theta_y\sin\theta_z)^2}{b^2}+ \frac{(\sin\theta_x\cos\theta_y)^2}{c^2}\\ C=&\frac{(\sin\theta_x\sin\theta_y\cos\theta_z-\sin\theta_x\sin\theta_z)^2}{a^2}+ \frac{(\sin\theta_x\cos\theta_z+\cos\theta_x\sin\theta_y\sin\theta_z)^2}{b^2}+ \frac{(\sin\theta_x\cos\theta_y)^2}{c^2}\\ \end{aligned} $$

显然，$A,B,C$是严格正的；相应的我们也可以写出$D,E,...,L$的表达式，不过他们没有什么特殊的性质因此不在此列出。需要注意的是，我们在计算时需要将式子展开，因此$x_0,y_0,z_0$的信息将会包含在$D,E,...,L$中。

由于$A,B,C$的性质，我们做如下变形：

$$ \begin{aligned} \frac{1}{3}&[(A+B+C)(x^2+y^2+z^2)+(A-C)(x^2+y^2-2z^2)+(A-B)(x^2-2y^2+z^2)]+\\ &Dxy+Exz+Fyz+Gx+Hy+Kz+L=0 \end{aligned} $$

在上式两边同乘$\frac{3}{A+B+C}$可以得到：

$$ x^2+y^2+z^2-U(x^2+y^2-2z^2)-V(x^2-2y^2+z^2)-4Mxy-2Nxz-2Pyz-Qx-Ry-Sz-T=0 $$

稍微整理一下，我们可以得到形如以下的式子：

$$ (x^2+y^2-2z^2)X_1+(x^2-2y^2+z^2)X_2+(4xy)X_3+(2xz)X_4+(2yz)X_5+(x)X_6+(y)X_7+(z)X_8+X_9=x^2+y^2+z^2 $$

这里记

$$ \vec{X}=\begin{pmatrix} X_1,X_2,X_3,X_4,X_5,X_6,X_7,X_8,X_9 \end{pmatrix}^\intercal $$

这样，我们可以将上式转换成$A\vec{X}=b$的形式，也就是说，我们可以通过最小二乘法解出与$x_0,y_0,z_0,a,b,c,\theta_x,\theta_y,\theta_z$相关的$X_1,X_2,X_3,X_4,X_5,X_6,X_7,X_8,X_9$，再通过对应关系求出椭球的九个参数信息。

特别感谢Green L对原博文的斧正，由于我最近没什么空写相关的代码，因此这里也贴一下Green L给出的Sample Code：
MATLAB ellipsoid fit(by Yury Petrov)