球经过摄像机投影为椭圆,而由于射影变换不能保证简比不变,因此球心的投影并不在椭圆中心,而是会向主点一侧偏移。运动捕捉实践中,常常将椭圆中心作为球心投影进行定位,这样会产生一定的误差,下面分析球心投影的真实位置和传统方法的误差大小。

公式推导

下面分别通过欧氏几何和射影几何的方法给出球(二次曲面)的成像结果。

欧氏几何方法推导球心投影位置

只要球心不在主轴上,其投影就会变成椭圆。按照对称性,此问题可简化为平面问题。如图 1 所示,设圆椎为球的视椎,视椎的半椎角为 \(\theta\) ,球心视线与像平面的夹角为 \(\alpha\) ,\(O\) 为视点,\(P\) 为球心,\(A\) 为椭圆的近心点,\(B\) 为椭圆的远心点,\(C\) 为球心的真实投影位置,\(D\) 为椭圆中心。其中 \(AB=2a\) 为椭圆长轴(显然 \(AD=DB=a\)),\(OC=h\) 为光心到像点的距离,\(AC=l\) 为椭圆中心到近心点的距离。

Figure 1: 球(线段)在像平面上的投影

Figure 1: 球(线段)在像平面上的投影

因此圆椎的方程可写为 \(x^2+z^2=y^2\tan^2\theta\) ,\(A\) 点坐标为 \((\frac{-h}{\cot\theta + \cot\alpha},\ \frac{h\cot\theta}{\cot\theta + \cot\alpha})\) ,\(B\) 点坐标为 \((\frac{h}{\cot\theta - \cot\alpha},\ \frac{h\cot\theta}{\cot\theta - \cot\alpha})\) ,\(D\) 点坐标为 \((\frac{h\cot\alpha}{\cot^2\theta - \cot^2\alpha},\ \frac{h\cot^2\theta}{\cot^2\theta - \cot^2\alpha})\) ,半长轴 \(a=\frac{h\cot\theta}{\sin\alpha(\cot^2\theta - \cot^2\alpha)}\) 。椭圆短轴过 \(D\) 点并垂直于 \(XY\) 平面,将 \(D\) 点的横纵坐标代入圆椎方程,求得坚坐标满足 \(z^2=\frac{h^2}{\cot^2\theta - \cot^2\alpha}\) 。所以椭圆的半短轴为 \(b=\frac{h}{\sqrt{\cot^2\theta - \cot^2\alpha}}\) ,进而求得椭圆的离心率为 \(e=\sqrt{1-\frac{b^2}{a^2}} = \frac{\cos\alpha}{\cos\theta}\) 。另有几何关系 \(l=\frac{h}{\sin\alpha\cot\theta + \cos\alpha}\) ,因此得到真实投影点 \(C\) 的位置关系:\(\frac{l}{a}=1-\frac{\cot\alpha}{\cot\theta}\) 。

综上所述,球的成像椭圆和球心投影位置满足关系:

\begin{equation} \label{eqn.sphere_projection} \left\{\begin{aligned} a &= \frac{h\cos\theta}{\sin\alpha(\cot^2\theta - \cot^2\alpha)} \\\
b &= \frac{h}{\sqrt{\cot^2\theta - \cot^2\alpha}} \\\
e &= \frac{\cos\alpha}{\cos\theta} \\\
\frac{l}{a} &= 1-\frac{\cot\alpha}{\cot\theta} \end{aligned}\right. \end{equation}

显然,球心投影位置并不在椭圆焦点处。当且仅当 \(\alpha = \theta\) 时,球心投影位置与像平面二次曲线的焦点重合,而这时二次曲线为抛物线(\(e = 1\))。即球与过光心平行于像平面的平面相切且深度为正,这显然不符合物理实际。

射影几何方法推导二次曲面成像

在射影几何中, \(n\) 维空间中的点一般用 \(n+1\) 维矢量的齐次坐标来表示。两个平行的矢量,由于只相差一个因子,在射影空间中代表同一个点。按照摄像机投影方程,三维空间点与二维像满足关系 \(s\pmb{x} = \pmb{PX}\) 。方程左边的因子 \(s\) 代表空间点在摄像机坐标系中的深度,可舍去。在射影平面内点和直线互为对偶,三维射影空间中点与平面互为对偶。直线 \(\pmb{l}\) 上的点 \(\pmb{x}\) 与其存在关系 \(\pmb{l}^T\pmb{x} = \pmb{l}^T\pmb{PX} = \pmb{\pi}^T\pmb{X} = 0\) ,其中 \(\pmb{X}\) 是过直线和投影中心的平面 \(\pmb{\pi}\) 上的点,如此就得到了像平面上直线的反投影平面方程 \(\pmb{\pi} = \pmb{P}^T\pmb{l}\) 。

在空间中的二次曲面 \(\pmb{Q}\) 上的点 \(\pmb{X}\) 满足二次型 \(\pmb{X}^T\pmb{QX} = 0\) ,其中 \(\pmb{Q}\) 为二次型矩阵,因此为实对称阵(\(\pmb{Q}^T = \pmb{Q}\))。空间中的任意点 \(\pmb{X}\) 与平面 \(\pmb{\pi} = \pmb{QX}\) 构成对应的极点和极平面,一般来说极点不在极平面上 (\(\pmb{X}^T\pmb{QX} \neq 0\))。当且仅当 \(\pmb{X}\) 在二次曲面上时,极平面经过极点(\(\pmb{X}^T\pmb{QX} = 0\)),此时平面 \(\pmb{\pi} = \pmb{QX}\) 是二次曲面上 \(\pmb{X}\) 点处的切平面。将 \(\pmb{X} = \pmb{Q}^{-1}\pmb{\pi}\) 代入二次曲面方程可以得到其对偶方程为 \(\pmb{\pi}^T\pmb{Q}^*\pmb{\pi} = 0\) ,其中二次曲面的对偶 \(\pmb{Q}^* = \pmb{Q}^{-1}\) 是由它的无数切平面围成的包络。类似的,在平面上,二次曲线 \(\pmb{C}\) 上的点 \(\pmb{x}\) 满足关系 \(\pmb{x}^T\pmb{Cx} = 0\) ,其对偶 \(\pmb{C}^* = \pmb{C}^{-1}\) 是由它的无数切线 \(\pmb{l}\) 围成的包络,满足关系 \(\pmb{l}^T\pmb{C}^*\pmb{l} = 0\) 。

显然空间二次曲面经过摄像机投影的轮廓为平面二次曲线,则组成对偶二次曲线 \(\pmb{C}^*\) 的包络切线 \(\pmb{l}\) 经过反投影得到的平面一定是空间对偶二次曲面 \(\pmb{Q}^*\) 的包络切平面。因此 \(\pmb{\pi}^T\pmb{Q}^*\pmb{\pi} = \pmb{l}^T\pmb{PQ}^*\pmb{P}^T\pmb{l} = \pmb{l}^T\pmb{C}^*\pmb{l} = 0\) ,所以二次曲面成像的外轮廓曲线可由下式给出:

\begin{equation} \label{eqn.quadric_projection} \pmb{C}^* = \pmb{PQ}^*\pmb{P}^T \end{equation}

此方程适用于包括球在内的所有二次曲面成像过程,更具一般意义 Hartley-2004-1

椭球的标准形式为 \(\frac{X^2}{a^2} + \frac{Y^2}{b^2} + \frac{Z^2}{c^2} = 1\) ,对方程进行齐次变换 \(X = X_1/X_4,\ Y = X_2/X_4,\ Z = X_3/X_4\) 得到 \(\frac{X_1^2}{a^2} + \frac{X_2^2}{b^2} + \frac{X_3^2}{c^2} - X_4^2= 0\) ,其二次型矩阵为对角阵:

\begin{equation} \label{eqn.standard_ellipsoid} \pmb{\Lambda} = \begin{bmatrix} 1/a^2 & & & \\ & 1/b^2 & & \\ & & 1/c^2 & \\ & & & -1 \end{bmatrix} \end{equation}

经过在空间中旋转和平移变换,就能得到椭球的一般形式。例如 \(\pmb{X}’ = \pmb{HX}\) ,其中 \(\pmb{H} = \begin{bmatrix}\pmb{R} & \pmb{T} \\ & 1\end{bmatrix}\) ,则有 \(\pmb{Q} = \pmb{H}^{-T} \pmb{\Lambda H}^{-1}\) 。旋转 \(\pmb{R}\) 为正交矩阵,\(\pmb{T}\) 为平移矢量。要从椭球的一般形式得到标准形式 \(\pmb{\Lambda} = \pmb{H}^T\pmb{QH}\) ,需首先求出欧氏变换 \(\pmb{H}\) 。稍加分析,可知矢量 \(\pmb{T}\) 是椭球的中心坐标。根据二次曲面对偶关系,中心点对应的极平面是无穷远平面,因此 \(\pmb{\pi}_\infty = \pmb{Q}\tilde{\pmb{T}}\) 或者 \(\tilde{\pmb{T}} = \pmb{Q}^*\pmb{\pi}_\infty\) ,其中 \(\tilde{\pmb{T}} = [\pmb{T}^T,\ 1]^T,\ \pmb{\pi}_\infty = [0,,0,,0,,1]^T\) 。变换 \(\pmb{H}\) 可分解为:

\begin{equation} \label{eqn.H-TR} \pmb{H} = \begin{bmatrix} \pmb{I} & \pmb{T} \\ & 1 \end{bmatrix} \begin{bmatrix} \pmb{R} & \\ & 1 \end{bmatrix} \end{equation}

令 \(\pmb{T}_h = \begin{bmatrix}\pmb{I} & \pmb{T} \\ & 1\end{bmatrix},\ \pmb{R}_h = \begin{bmatrix}\pmb{R} & \\ & 1\end{bmatrix}\) ,则有 \(\pmb{\Lambda} = \pmb{R}_h^T\pmb{T}_h^T\pmb{QT}_h\pmb{R}_h\) 。在求出椭球中心后,可通过变换 \(\pmb{Q}’ = \pmb{T}_h^T\pmb{QT}_h\) 将椭球平移到原点。再对二次型矩阵 \(\pmb{Q}'\) 进行特征值分解,即可求得正交矩阵 \(\pmb{R}_h\) 和椭球的标准形式\eqref{eqn.standard_ellipsoid}。 同理,平面上的椭圆拥有类似的标准化过程。例如通过式\eqref{eqn.quadric_projection}求出球的像(椭圆),然后进行二次型标准化即可得到椭圆的半轴长度和方向等信息。

对于球来说,给定球心位置 \((a,\ b,\ c)\) 和半径 \(r\) ,则球面方程可写为 \((X-a)^2+(Y-b)^2+(Z-c)^2 = r^2\) 。对其进行齐次变换得到 \(X_1^2 + X_2^2 + X_3^2 + (a^2+b^2+c^2-r^2)X_4^2 - 2aX_1X_4 - 2bX_2X_4 - 2cX_3X_4 = 0\) ,则球的二次型矩阵可写为:

\begin{equation} \label{eqn.sphere_matrix} \pmb{Q} = \begin{bmatrix} 1 & & & -a \\ & 1 & & -b \\ & & 1 & -c \\ -a & -b & -c & a^2+b^2+c^2-r^2 \end{bmatrix} \end{equation}

其对偶二次曲面为 \(\pmb{Q}^* = \pmb{Q}^{-1}\) 。

MATLAB 仿真分析

给定摄像机参数:视场角为 60~80 度,图像大小 1280x1024,摄像机无畸变,视场范围 1~10 m,另外 Marker 小球的直径 16 mm。下面分别按照式\eqref{eqn.sphere_projection}和\eqref{eqn.quadric_projection}求解椭球的投影,并对比两种方法得到的结果,进行互相印证。

仿真的结果如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
fc =  768.46
err =

   4.7092e-08
   4.1481e-08
   4.1059e-04
   9.4295e-06
   1.7283e-07
   2.1432e-07

stddx =  0.0024044
maxdx =  0.042632
id =  1900
exid =  0.67850
axid =  8.1111
bxid =  5.9873
maxex =  0.71267

经过仿真对比可以验证两种方法殊途同归,结果中变量 err 表示两种方法结果的误差大小,用来衡量结果的接近程度。变量 stddx 为球心投影位置与椭圆中心距离的标准差,为传统方法的平均误差;变量 maxdx 为球心投影位置与椭圆中心的最大距离,即传统方法的最大误差。作为本文最重要的结论,当虚拟 Marker 半径为 8 mm 时,最大误差为 0.02~0.05 个像素;当 Marker 半径为 16 mm 时,最大误差达到 0.12~0.18 个像素。

结论

所以对于球形 Marker 的定位来说,球半径越大则结果的误差越大。当球半径不超过 1 cm 时,最大误差不超过 0.08 个像素,利用椭圆中心代替球心投影是完全可以接受的;当球的半径超过 1 cm,若要取得高定位精度则最好计算这个误差。

Bibliography

[Hartley-2004-1] Hartley & Zisserman, Multiple view geometry in computer vision, Cambridge university press (2004).