椭球与椭圆方程拟合

1 背景
有时，需要将一堆数据拟合成椭圆(椭球)方程，有时需要将一般椭圆(椭球)方程标准化。这时就需要使用专门的算法进行处理，本文介绍算法进行处理。

2 椭球
椭球是3D空间方程求解。
2.1 椭球拟合
对于椭球的一般方程，可表示为如下形式  
(2.1-1)
$$
s_1x^2+s_2y^2+s_3z^2+2s_4yz+2s_5xz+2s_6xy+2s_7x+2s_8y+2s_9z+s_{10}=0
$$
上面的方程要表示为椭球，必须满足如下条件 
(2.1-2)
$$
\begin{cases}
W&=\begin{bmatrix}s_1&s_6&s_5
\\
s_6&s_2&s_4
\\
s_5&s_4&s_3\end{bmatrix}
\\
\\
|W|(s_1+s_2+s_3)&>0\\\\s_1s_2+s_2s_3+s_3s_1&>s_4^2+s_5^2+s_6^2
\end{cases}
$$
现在已知有n组\((x,y,z)\)数据，为了拟合椭球，先设置如下中间变量  
(2.1-3)
$$
\begin{cases}
\mathbf{v}_1 &=\begin{bmatrix}s_1 &s_2&s_3&s_4&s_5&s_6\end{bmatrix}^T
\\
\mathbf{v}_2 &=\begin{bmatrix}s_7 &s_8&s_9&s_{10}\end{bmatrix}^T
\\
\\
\eta &=\dfrac{\alpha}{2}-1
\\
B_1 &=\dfrac{1}{3\eta^2+2\eta^3-1}\begin{bmatrix}
1-\eta^2 &\eta+\eta^2&\eta+\eta^2
\\
\eta+\eta^2 &1-\eta^2&\eta+\eta^2
\\
\eta+\eta^2&\eta+\eta^2&1-\eta^2
\end{bmatrix}
\\
\\
B_2&=-\dfrac{1}{\alpha}\begin{bmatrix}
1&0&0
\\
0&1&0
\\
0&0&1
\end{bmatrix}
\\
\\
B&=\begin{bmatrix}
B_1 &\mathbf{0}
\\
\mathbf{0} &B_2
\end{bmatrix}
\\
\\
D&=\begin{bmatrix}
x_1^2&y_1^2&z_1^2&2y_1z_1&2x_1z_1&2x_1y_1&2x_1&2y_1&2z_1&1
\\
x_2^2&y_2^2&z_2^2&2y_2z_2&2x_2z_2&2x_2y_2&2x_2&2y_2&2z_2&1
\\
…
\\
x_n^2&y_n^2&z_n^2&2y_nz_n&2x_nz_n&2x_ny_n&2x_n&2y_n&2z_n&1
\end{bmatrix}
\end{cases}
$$
然后根据文献[1]，设置一个\(\alpha \geq 4\)，然后如下算法求解\(\mathbf{v}_1,\mathbf{v}_2\).  
(2.1-4)
$$
\begin{cases}
D^TD &=\begin{bmatrix}
S_{11} & S_{12}
\\
S_{12}^T &S_{22}
\end{bmatrix}
\\
M&=-S_{22}^{-1}S_{12}^T
\\
P&=B(S_{11}+S_{12}M)
\\
\mathbf{v}_2&=M\mathbf{v}_1
\end{cases}
$$
上面公式里\(S_{11},S_{12},S_{22}\)的维度分别为\(6\times6,6\times4,4\times4\)。在实际使用时，取\(P\)的特征向量作为\(\mathbf{v}_1\)，因为特征向量有多个，因此选择满足条件(2.1-2)的特征向量作为解。注意，(2.1-2)不能完全排除非椭球方程，从(2.2-7)可知道，椭球方程还与其它系数相关，因此严格判断结果是否为椭球，需要确定最终求解的系数能否通过椭球标准化，通过则为椭球系数。

注意，在实际使用过程中，需要设置\(\alpha\),一般地，设置为4即可，如果4不满足，则设置其它值。

2.2 椭球标准化
对于标准化椭球方程，可表示为如下形式  
(2.2-1)
$$
\dfrac{x'^2}{a^2}+\dfrac{y'^2}{b^2}+\dfrac{z'^2}{c^2}=1
$$
对于参数化标准方程有如下表示方式  
(2.2-2)
$$
\begin{cases}
x'&=a\sin(\theta)\cos(\phi)
\\
y'&=b\sin(\theta)\sin(\phi)
\\z'&=c\cos(\theta)\end{cases}\ \ \ ,0\leq \theta\leq\pi,0\leq\phi<2\pi
$$
现在考虑，(2.1-1)的椭球方程，当经过\([R,\mathbf{t}]\)的刚体变换后，可以在另一个坐标系下表示成(2.2-1)的标准椭球方程。则有如下表示  
(2.2-3)
$$
\begin{cases}
\mathbf{x}' &=R(\mathbf{x}+\mathbf{t})
\\
\mathbf{x} &=\begin{bmatrix}x&y&z\end{bmatrix}^T
\\
\mathbf{x}' &=\begin{bmatrix}x'&y'&z'\end{bmatrix}^T
\end{cases}
$$
因此，标准化过程需要求的参数为\(a,b,c,R,\mathbf{t}\)。 
这里设置中间矩阵  
(2.2-4)
$$
\begin{aligned}
S&=\begin{bmatrix}
\dfrac{1}{a^2} &0&0
\\
0&\dfrac{1}{b^2}&0
\\
0&0&\dfrac{1}{c^2}
\end{bmatrix}
\\
Q&=R^TSR
\end{aligned}
$$
因此，求得\(S\)矩阵，即可求得\(a,b,c\)。这里将(2.2-3)代入(2.2-1)有  
(2.2-5)
$$
\begin{aligned}
1 &=\dfrac{x'^2}{a^2}+\dfrac{y'^2}{b^2}+\dfrac{z'^2}{c^2}
\\
&=(R(\mathbf{x}+\mathbf{t}))^TS(R(\mathbf{x}+\mathbf{t}))
\\
&=(\mathbf{x}+\mathbf{t})^TR^TSR(\mathbf{x}+\mathbf{t})
\\
&=(\mathbf{x}+\mathbf{t})^TQ(\mathbf{x}+\mathbf{t})
\\
&=\mathbf{x}^TQ\mathbf{x}+\mathbf{t}^TQ\mathbf{t}+2\mathbf{t}^TQ\mathbf{x}
\end{aligned}
$$
由于(2.2-5)与(2.1-1)为同一个椭球，因此总存在一个非0系数\(k\)使得\(W=kQ\)情况下，下面算式\(x,y,z\)以及常数项系数与(2.1-1)对应  
(2.2-6)
$$
\mathbf{x}^TW\mathbf{x}+2\mathbf{t}^TW\mathbf{x}+\mathbf{t}^TW\mathbf{t}-k=0
$$
根据系数对应关系，可以得到  
(2.2-7)
$$
\begin{cases}
W&=R^T(kS)R
\\
\\
\mathbf{t} &=W^{-1}\begin{bmatrix}
s_7
\\
s_8
\\
s_9
\end{bmatrix}=R^T(kS)^{-1}R\begin{bmatrix}
s_7
\\
s_8
\\
s_9
\end{bmatrix}
\\
\\
k &=\mathbf{t}^TW\mathbf{t}-s_{10}
\end{cases}
$$
因此，在上面公式里先奇异值分解\(W\)获得\(R,kS\)(注意，这里按第4部分进行分解，只有第4部分对称矩阵分解提到的\(m=0\)或者\(m = n\)时，才有解，否则此参数不满足椭球方程。也可以通过特征分解来获得)，然后计算\(\mathbf{t}\)，最后计算\(k\)，进而得到\(S\)(求解得到的\(S\)对角线元素必须都大于0，否则椭球方程无解)。需要注意的是，奇异值分解不唯一，可以想象，在标准椭球体\(x',y',z'\)坐标系下，坐标轴对调后的方程依然为标准椭球体方程，所以需要注意选择。

最终，在求得标准椭球方程\(a,b,c\)系数后，就可以求椭球体积，以及表面积。为了方便求椭球体积\(V\)与表面积\(A\)，这里约定系数\(a\geq b \geq c >0\)，且设置如下中间变量  
(2.2-8)
$$
\begin{cases}
e_1 &=\sqrt{\dfrac{a^2-c^2}{a^2}}
\\
\\
e_2&=\sqrt{\dfrac{b^2-c^2}{b^2}}
\\
\\
k &=\dfrac{e_2}{e_1}
\\
\\
\phi &=\arcsin(\sqrt{1-\dfrac{c^2}{a^2}})
\end{cases}
$$
最终，使用如下公式计算体积与面积[3]  
(2.2-9)
$$
\begin{cases}
V&=\dfrac{4\pi}{3}abc
\\
A&=2\pi(c^2+\dfrac{bc^2}{\sqrt{a^2-c^2}}\mathbf{EllipticF}(\phi,k)+b\sqrt{a^2-c^2}\ \mathbf{EllipticE}(\phi,k))
\end{cases}
$$
公式里 \(\mathbf{EllipticF}、\mathbf{EllipticE}\) 分别为第一类不完全椭圆积分函数，第二类不完全椭圆积分函数。

注意，如果椭球有两个轴相等，这里设\(a=b\)，则表面积\(A\)存在如下精确解

(2.2-10)
$$
A=\begin{cases}
2\pi(a^2+c^2\dfrac{\mathbf{arctanh}(\sin(\phi))}{\sin(\phi)})&,\phi=\arccos(\dfrac{c}{a}),a>c
\\
2\pi(a^2+c^2\dfrac{\phi}{\tan(\phi)})&,\phi=\arccos(\dfrac{a}{c}),a\leq c
\end{cases}
$$

3 椭圆
椭圆是2D空间求解。
3.1 椭圆拟合
对于椭圆的一般方程，可表示为如下形式  
(3.1-1)
$$
s_1x^2+s_2y^2+2s_3xy+2s_4x+2s_5y+s_6=0
$$
上面的方程要表示为椭圆，要求如下条件成立  
(3.1-2)
$$
\begin{cases}
W&=\begin{bmatrix}s_1&s_3
\\
s_3&s_2\end{bmatrix}\\\\
G&=\begin{bmatrix}
s_1 &s_3&s_4
\\
s_3&s_2&s_5
\\
s_4&s_5&s_6
\end{bmatrix}\\\\
|W|&>0\\\\|G|(s_1+s_2)&<0
\end{cases}
$$
现在已知有n组\((x,y)\)数据，为了拟合椭圆，先设置如下中间变量  
(3.1-3)
$$
\begin{cases}
\mathbf{v}_1 &=\begin{bmatrix}s_1 &2s_3&s_2\end{bmatrix}^T
\\
\mathbf{v}_2 &=\begin{bmatrix}2s_4 &2s_5&s_6\end{bmatrix}^T
\\
\\
B&=\dfrac{1}{2}\begin{bmatrix}
0&0&1
\\
0&-2&0
\\
1&0&0
\end{bmatrix}
\\
\\
D_1&=\begin{bmatrix}
x_1^2&x_1y_1&y_1^2
\\
x_2^2&x_2y_2&y_2^2
\\
…
\\
x_n^2&x_ny_n&y_n^2
\end{bmatrix}
\\
\\
D_2&=\begin{bmatrix}
x_1&y_1&1
\\
x_2&y_2&1
\\
…
\\
x_n&y_n&1
\end{bmatrix}
\end{cases}
$$
然后根据文献[2]，使用如下算法求解\(\mathbf{v}_1,\mathbf{v}_2\).  
(3.1-4)
$$
\begin{cases}
S_1 &=D_1^TD_1
\\
S_2 &=D_1^TD_2
\\
S_3&=D_2^TD_2
\\
M&=-S_3^{-1}S_2^T
\\
P&=B(S_1+S_2M)
\\
\mathbf{v}_2&=M\mathbf{v}_1
\end{cases}
$$
在实际使用时，取\(P\)的特征向量作为\(\mathbf{v}_1\)，因为特征向量有多个，因此选择满足条件(3.1-2)的特征向量作为解。

3.2 椭圆标准化
对于标准化椭圆方程，可表示为如下形式  
(3.2-1)
$$
\dfrac{x'^2}{a^2}+\dfrac{y'^2}{b^2}=1
$$

而对于参数化标准椭圆可表示为如下方式

(3.2-2)
$$
\begin{cases}
x'&=a\cos \theta
\\
y'&=b\sin\theta
\end{cases}\ \ \ ,0\leq \theta<2\pi
$$
现在考虑，(3.1-1)的椭圆方程，当经过\([R,\mathbf{t}]\)的刚体变换后，可以在另一个坐标系下表示成(3.2-1)的标准椭圆方程。则有如下表示  
(3.2-3)
$$
\begin{cases}
\mathbf{x}' &=R(\mathbf{x}+\mathbf{t})
\\
\mathbf{x} &=\begin{bmatrix}x&y\end{bmatrix}^T
\\
\mathbf{x}' &=\begin{bmatrix}x'&y'\end{bmatrix}^T
\end{cases}
$$
因此，标准化过程需要求的参数为\(a,b,R,\mathbf{t}\)。 
这里设置中间矩阵  
(3.2-4)
$$
\begin{aligned}
S&=\begin{bmatrix}
\dfrac{1}{a^2} &0
\\
0&\dfrac{1}{b^2}
\end{bmatrix}
\\
Q&=R^TSR
\end{aligned}
$$
因此，求得\(S\)矩阵，即可求得\(a,b\)。这里将(3.2-3)代入(3.2-1)有  
(3.2-5)
$$
\begin{aligned}
1 &=\dfrac{x'^2}{a^2}+\dfrac{y'^2}{b^2}
\\
&=(R(\mathbf{x}+\mathbf{t}))^TS(R(\mathbf{x}+\mathbf{t}))
\\
&=(\mathbf{x}+\mathbf{t})^TR^TSR(\mathbf{x}+\mathbf{t})
\\
&=(\mathbf{x}+\mathbf{t})^TQ(\mathbf{x}+\mathbf{t})
\\
&=\mathbf{x}^TQ\mathbf{x}+\mathbf{t}^TQ\mathbf{t}+2\mathbf{t}^TQ\mathbf{x}
\end{aligned}
$$
由于(3.2-5)与(3.1-1)为同一个椭圆，因此总存在一个非0系数\(k\)使得\(W=kQ\)情况下，下面算式\(x,y\)以及常数项系数与(3.1-1)对应  
(3.2-6)
$$
\mathbf{x}^TW\mathbf{x}+2\mathbf{t}^TW\mathbf{x}+\mathbf{t}^TW\mathbf{t}-k=0
$$
根据系数对应关系，可以得到  
(3.2-7)
$$
\begin{cases}
W&=R^T(kS)R
\\
\\
\mathbf{t} &=W^{-1}\begin{bmatrix}
s_4
\\
s_5
\end{bmatrix}=R^T(kS)^{-1}R\begin{bmatrix}
s_4
\\
s_5
\end{bmatrix}
\\
\\
k &=\mathbf{t}^TW\mathbf{t}-s_6
\end{cases}
$$
因此，在上面公式里先奇异值分解\(W\)获得\(R,kS\)(注意，这里按第4部分进行分解，只有第4部分对称矩阵分解提到的\(m=0\)或者\(m = n\)时，才有解，否则此参数不满足椭圆方程。也可以通过特征分解来获得)，然后计算\(\mathbf{t}\)，最后计算\(k\)，进而得到\(S\)(求解得到的\(S\)对角线元素必须都大于0，否则椭圆方程无解)。需要注意的是，奇异值分解不唯一，可以想象，在标准椭圆\(x',y'\)坐标系下，坐标轴对调后的方程依然为标准椭圆方程，所以需要注意选择。

有的情况下，只想快速确认椭圆长短半轴，这时可直接使用如下公式计算[4]  
(3.2-8)
$$
\begin{cases}
p_1&=\sqrt{(s_1-s_2)^2+4s_3^2}
\\
p_2&=(s_1+s_2)
\\
p_3&=s_3^2-s_1s_2
\\
p_4&=2(s_1s_5^2+s_2s_4^2+s_6s_3^2-2s_3s_4s_5-s_1s_2s_6)
\\
a&=\sqrt{\dfrac{p_4}{p_3(p_1-p_2)}}
\\
b&=\sqrt{\dfrac{p_4}{p_3(-p_1-p_2)}}
\end{cases}
$$

最终，在求得标准椭圆方程\(a,b\)系数后，就可以求椭圆面积，以及周长。为了方便求椭圆面积\(A\)与周长\(C\)，这里约定系数\(a\geq b\)，可使用如下公式计算面积与周长[4]  
(3.2-9)
$$
\begin{cases}
A &=\pi a b
\\
C &=4a\ \mathbf{EllipticE}(\sqrt{\dfrac{a^2-b^2}{a^2}}) &\approx \pi(a+b)(1+\dfrac{3\gamma^2}{10+\sqrt{4-3\gamma^2}}) \\ \gamma &=\dfrac{a-b}{a+b}
\end{cases}
$$
公式里 \(\mathbf{EllipticE}\) 为第二类完全椭圆积分函数，周长后面的近似公式来自拉马努金公式，这个公式精度较高。

4 对称矩阵分解
对于\(n\)阶对称矩阵\(A\),总可以表示为如下形式  
(4-1)
$$
A = UDWU^T
$$
其中\(U\)为正交矩阵,\(D\)为对角线上元素非负的对角矩阵,\(W\)为对角线元素只能为1或-1的对角矩阵。 
对矩阵\(A\)进行奇异值分解有  
(4-2)
$$
A = UDV^T=UDWU^T \Rightarrow V=UW
$$
这里设\(\mathbf{u}_i、\mathbf{v}_i\)分别为\(U、V\)的第\(i\)列，且\(W=\mathbf{Diag}(w_1,w_2,…)\)则可以得到  
(4-3)
$$
\mathbf{v}_i=w_i\mathbf{u}_i
\\
\Downarrow
\\
||\mathbf{v}_i-\mathbf{u}_i||^2=\begin{cases}
0 &,w_i=1
\\
||2\mathbf{v}_i||^2=4||\mathbf{v}_i||^2=4&,w_i=-1
\end{cases}
$$
因此有  
(4-4)
$$
m=\dfrac{1}{4}\sum_{i=1}^n||\mathbf{v}_i-\mathbf{u}_i||^2=\dfrac{1}{4}||U-V||^2
$$
从上面可以发现\(m\)反应了\(W\)对角线元素为-1的个数，并且可以发现，对称矩阵进行奇异值分解后得到的\(U、V\)每一列有对应关系，特别地\(m=0\)时，它们完全一致。

另一方面，如果对称矩阵\(A\)存在如下(4-5)分解(\(U\)为正交矩阵，\(S\)为对角矩阵，且对角线元素同号)，则\(A\)的奇异值分解必须满足\(||U-V||^2=0\text{ or }||U+V||^2=0\)。

(4-5)
$$
A=USU^T
$$
参考
[1] Li Q , Griffiths J G .Least squares ellipsoid specific fitting[C]//Geometric Modeling & Processing.IEEE, 2004.DOI:10.1109/GMAP.2004.1290055.

[2] Radim H R .Numerically Stable Direct Least Squares Fitting Of Ellipses[J].  1999.DOI:doi:http://dx.doi.org/.

[3] 椭球函数

[4] 椭圆函数