球心投影位置

球心投影位置

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

公式推导

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

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

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

图 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{array}{rl}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{array}\text{\hspace{1em}(1)}$$

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

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

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

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

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

$${\mathbf{C}}^{\ast }={\mathbf{P}\mathbf{Q}}^{\ast }{\mathbf{P}}^T\text{\hspace{1em}(2)}$$

此方程适用于包括球在内的所有二次曲面成像过程，更具一般意义 ^{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$ ，其二次型矩阵为对角阵：

$$\mathbf{\Lambda }=\left[\begin{array}{cccc}1/a^2& & & \\ & 1/b^2& & \\ & & 1/c^2& \\ & & & -1\end{array}\right]\text{\hspace{1em}(3)}$$

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

$$\mathbf{H}=\left[\begin{array}{cc}\mathbf{I}& \mathbf{T}\\ & 1\end{array}\right]\left[\begin{array}{cc}\mathbf{R}& \\ & 1\end{array}\right]$$

令 ${\mathbf{T}}_h=\left[\begin{array}{cc}\mathbf{I}& \mathbf{T}\\ & 1\end{array}\right],{\mathbf{R}}_h=\left[\begin{array}{cc}\mathbf{R}& \\ & 1\end{array}\right]$ ，则有 $\mathbf{\Lambda }={\mathbf{R}}_h^T{\mathbf{T}}_h^T{\mathbf{Q}\mathbf{T}}_h{\mathbf{R}}_h$ 。在求出椭球中心后，可通过变换 ${\mathbf{Q}}^{\prime }={\mathbf{T}}_h^T{\mathbf{Q}\mathbf{T}}_h$ 将椭球平移到原点。再对二次型矩阵 ${\mathbf{Q}}^{\prime }$ 进行特征值分解，即可求得正交矩阵 ${\mathbf{R}}_h$ 和椭球的标准形式 (3)。 同理，平面上的椭圆拥有类似的标准化过程。例如通过式 (2) 求出球的像（椭圆），然后进行二次型标准化即可得到椭圆的半轴长度和方向等信息。

对于球来说，给定球心位置 $(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-2a{X}_1{X}_4-2b{X}_2{X}_4-2c{X}_3{X}_4=0$ ，则球的二次型矩阵可写为：

$$\mathbf{Q}=\left[\begin{array}{cccc}1& & & -a\\ & 1& & -b\\ & & 1& -c\\ -a& -b& -c& a^2+b^2+c^2-r^2\end{array}\right]$$

其对偶二次曲面为 ${\mathbf{Q}}^{\ast }={\mathbf{Q}}^{-1}$ 。

MATLAB 仿真分析

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

仿真的结果如下：

fc =  961.51
err =

   3.5116e-08
   3.1730e-08
   3.2576e-04
   8.9195e-06
   2.8262e-08
   1.3897e-07

stddx =  0.0018535
maxdx =  0.039864
id =  1284
exid =  0.56819
axid =  9.0695
bxid =  7.4859
maxex =  0.61848

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

结论

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

参考文献

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