【论文翻译】Unbiased Estimator for Distorted Conics in Camera Calibration-相机标定中畸变圆锥曲线的无偏估计器

---

Unbiased Estimator for Distorted Conics in Camera Calibration-相机标定中畸变圆锥曲线的无偏估计器

阅读论文，顺带把翻译贴上来 ，方便后面再看。

---

摘要

在相关文献中，点和圆锥曲线一直是相机几何标定的主要特征。尽管圆锥曲线比点包含更多的信息，但在畸变情况下，圆锥曲线的属性损失严重限制了其在相机标定中的实用性。现有的许多方法都试图通过忽略畸变或引入3D球形目标来解决基于圆锥曲线的标定问题，以此规避这一限制。在本文中，我们提出了一种利用矩进行基于圆锥曲线的标定新方法。我们的方法建立在数学发现的基础之上，即即使在畸变条件下，也可以无偏地估计一阶矩。这使得我们能够追踪投影和畸变过程中的矩变化，确保畸变后的圆锥曲线的一阶矩能够得以保留。通过无偏估计器，我们可以在亚像素级别上精确地检测到圆形图案，并可以充分利用整个标定流程，从而显著提高标定效果。整个代码可从以下网址获取：https://github.com/ChaehyeonSong/discocal。

---

1. 引言

相机标定在3D计算机视觉和基于视觉的感知中至关重要。显然，从图像中理解几何结构的过程依赖于准确的相机标定，包括3D密集重建[17, 18]、视觉同步定位与地图构建（SLAM）[12, 19, 22]以及深度估计[9]。这个基本问题已通过多种方式得到解决[4, 5, 14, 26, 27]，并且大多数标定方法都利用覆盖有特殊图案的平面目标，如方格网格结构[25, 28]或圆形图案[10, 13]。这种基于平面目标的方法需要利用控制点（即正方形的角或圆的质心）的无偏投影模型进行精确测量，以获得准确的标定结果。许多文献研究了这两种图案的明显优缺点。

棋盘格图案可以确保在投影变换和畸变中进行精确（无偏）估计，但其检测精度仅限于像素级别。另一方面，圆形图案[11]在实现亚像素级别检测精度方面表现出色，然而有偏的投影模型导致了较差的标定结果。仅仅使用圆锥曲线特征来补偿透视偏差是不够的，更主要的因素是镜头的非线性畸变[16]。它会降低圆锥的几何特性，使得难以估计投影圆的质心，如图1所示。

图1. 畸变图像中的质心估计。
本图说明了相机投影和镜头畸变对图像中圆形目标的影响。由于这些畸变，圆形会失去其圆锥曲线属性并变形为畸变的椭圆，使得中心点跟踪变得具有挑战性。（红色）表示传统控制点估计方法的失败，如错误跟踪的中心点所示。（绿色）所提出的无偏估计器可以从闭合形式计算中准确识别出变换后椭圆的中心点。

在本文中，我们通过提出一个闭合形式的无偏估计器来处理畸变偏差，从而推动了圆形图案的发展。据我们所知，这是第一个描述径向多项式畸变下投影圆锥曲线特征的分析解决方案。受到矩表示可以在任何多项式映射下描述一般分布变换的事实的启发，我们证明了畸变圆锥的特性，如质心，可以表示为未畸变圆锥的n阶矩的线性组合。我们的工作不仅限于针孔相机标定的畸变模型。此外，我们提出了一种通用且可微分的方法来跟踪任意非线性变换下的圆锥曲线，这些变换可以近似为多项式函数。

通过进行合成实验和真实实验，我们验证了我们的估计器在各种畸变和已知相机参数的圆半径下的重投影误差是无偏的。我们的估计器还应用于RGB和热红外（TIR）相机的标定，这些相机由于边界模糊效应而难以检测控制点（见附录9.1）。因此，我们的方法在重投影误差和6D姿态方式上优于现有的基于圆形图案的方法和棋盘格方法。我们工作的主要贡献总结如下：

• 我们率先为畸变圆锥曲线提供了无偏估计器，以充分利用圆形图案的优点，继承了一个简单、鲁棒且准确的检测器。我们的工作结合了数学的优雅和实用的巧思，完成了基于圆锥的标定流程中缺失的部分。
• 我们利用了矩的概率概念，这是以前在标定中未曾尝试过的，并提供了圆锥曲线的一般矩作为解析形式，并进行了彻底的数学推导和证明。这种方法使我们能够为圆形图案设计无偏估计器。
• 当在合成图像和真实图像上进行测试时，我们的无偏估计器提高了整体标定性能。特别是，我们展示了我们的方法为TIR图像产生了显著改善的标定结果，这些图像通常包括高水平的模糊、噪声和显著的畸变。我们公开了我们的算法，以支持社区在其标定中使用圆锥特征。

2. 相关工作

平面图案和控制点。
Zhang [28]以及Sturm和Maybank [25]介绍了一种利用印有特定图案的平面目标进行标定的方法。平面目标应包含可以从特定图案中轻松且准确提取的控制点。最初考虑的是由黑白方块组成的棋盘格图案，然后由Heikkila [10]引入了带有网格结构的圆形图案作为替代。棋盘格图案使用方块的角作为控制点，而圆形图案则使用图像中投影圆的质心作为控制点，具有亚像素精度 [11]。由于控制点的精确位置至关重要，因此也考虑了通过交互式方式细化控制点的位置 [7]。

无偏估计器。
为了获得准确的标定结果，控制点的测量质量和估计值都至关重要。与在大多数情况下被证明具有无偏估计器的棋盘格图案不同，对圆形图案应用相同的估计器会导致偏差 [16]。具体来说，目标中的圆心并非投影到图像中投影圆的中心点上。有几项研究 [10, 14] 通过采用基于圆锥的变换来挑战这个问题。这些方法在线性变换（如透视变换）下实现了无偏估计器。然而，由于圆锥特性在非线性变换下不保持不变，因此畸变引起的偏差仍未解决。畸变引起的偏差在大多数相机中是一个主要因素 [16]。为了解决这个问题，Kannala和Brandt [13]为圆形图案引入了无偏估计器的广义概念，但无法为给定的积分方程推导出解析解。

无控制点方法。
同时，也开发了直接使用圆锥特性的方法。通过利用同心圆，早期的研究 [4, 14] 可以从圆锥方程中获得绝对圆锥的图像。有些研究 [26, 27] 专注于球体，其投影形状仅取决于中心点和相机之间的配置。不幸的是，这些方法也受到了畸变偏差的影响，这种偏差会破坏圆锥或球体的所有关键几何形状。为了处理非线性畸变函数，Devernay和Faugeras [8]利用了在未失真图像中始终保持直线的线特征。尽管这项工作没有明确估计内在参数，但最近的研究 [5] 引入了使用未失真图像中的线特征来求解内在参数的闭合形式解。

通用模型方法。
现有的相机模型假设了一个理想的镜头，并用少数参数来近似非线性映射。为了缓解这一假设，引入了更一般的射线追踪模型。Schops等人 [23]提出了一种使用B样条插值和最近点进行一般反投影的模型。另一项研究 [20] 提出了像素级焦距，以考虑径向方向的一般非线性映射。这些替代方法具有潜力；然而，现有的3D视觉应用仍然需要几何相机模型。因此，在本文中，我们专注于基于几何模型的标定，尤其是针孔模型。

3. 预备知识

3.1 符号表示

向量$p\in {\Re }^n$和矩阵$Q\in {\Re }^{(n\times n)}$分别用小写和大写的粗体表示，下标为坐标。我们将目标平面设置为与世界坐标的xy平面相同。因此，$p_w=(x_w,y_w{)}^{\top }$ 是目标平面上以目标坐标表示的二维点，$p_i$是图像中对应的点。该目标点$p_w=(x_w,y_w{)}^{\top }$通过透视投影投影到归一化平面上的点$p_n$，然后$p_n$在畸变和内部矩阵的作用下映射到图像平面上的点$p_i$。另外，$Q_w$是目标平面上椭圆的矩阵表示，$Q_n$是归一化平面上的椭圆。这些符号和坐标之间的几何关系的更多详细信息如图2所示。我们使用∼来表示齐次向量。同时，$\tilde{p}\simeq \tilde{q}$意味着两个向量在尺度上是相同的。

图2. 图像投影几何。
由于投影和镜头畸变，目标平面上的圆投影中心（圆圈点）、归一化平面上的椭圆投影中心（三角形标记）和图像平面上的实际形状质心（交叉标记）之间存在一些不匹配。由归一化平面中的非线性畸变引起的变换形状不能通过分析描述为原始圆锥曲线。然而，尽管存在这些畸变，我们的算法还是通过矩跟踪成功地定位了真实的质心。

3.2 圆锥曲线的矩阵表示

一个圆锥曲线是由一系列满足以下方程的点 $(x,y)$ 构成的集合：

$$a{x}^2+2bxy+c{y}^2+2dx+2ey+f=0(1)$$

这个方程也可以以矩阵的形式来表达：

$$x^{\top }Qx=0$$

其中，

$$Q=\left[\begin{array}{ccc}a& b& d\\ b& c& e\\ d& e& f\end{array}\right],x=\left[\begin{array}{c}x\\ y\\ 1\end{array}\right](2)$$

对称矩阵 $Q$ 是该圆锥曲线的特征矩阵。一般而言，圆锥曲线可能包括椭圆、抛物线和双曲线，但本文我们主要关注椭圆，因为在大多数情况下（除了极端情形），圆的图像都可以视为椭圆。给定特征矩阵 $Q$，椭圆的中心可以通过以下公式计算：

$$\tilde{p}=Q^{-1}\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right](3)$$

关于椭圆几何特征的详细推导，请参阅附录 8.1。

3.3 相机模型

采用针孔相机模型，目标平面上的二维点 ${\tilde{p}}_w$ 通过以下方式投影到图像上：

$$K=\left[\begin{array}{ccc}{f}_x& \eta & c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right],(4)$$

$$T_{cw}=\left[\begin{array}{cc}R& t\\ 0& 1\end{array}\right]=\left[\begin{array}{cccc}{r}_1& r_2& r_3& t\\ 0& 0& 0& 1\end{array}\right],(5)$$

$${\tilde{p}}_i=\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]\simeq K\left[\begin{array}{cccc}{r}_1& r_2& r_3& t\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ z_w\\ 1\end{array}\right](6)$$

$$\simeq K\left[\begin{array}{ccc}{r}_1& r_2& t\end{array}\right]\left[\begin{array}{c}{x}_w\\ y_w\\ 1\end{array}\right](\because z_w=0)(7)$$

$$\simeq KE{\tilde{p}}_w({\tilde{p}}_w\triangleq [x_w,y_w,1{]}^{\top })(8)$$

$$\simeq K{\tilde{p}}_n\simeq H{\tilde{p}}_w.(9)$$

这里，$K$ 是相机的内参矩阵，$E$ 是目标坐标系与相机坐标系之间的外参，而 $H$ 是单应性矩阵。由于 $p_w$ 是位于 $xy$ 平面上的一个点，所以 $z_w$ 为零。上述方程完全确定了目标点与归一化平面上的投影点（${\tilde{p}}_n$）之间的关系。不幸的是，当这个投影点（${\tilde{p}}_n$）被投影到图像平面上时，会进一步经历非线性变换（镜头畸变）。在径向畸变和切向畸变这两种类型中，我们只考虑径向畸变，这在大多数相机中已知是足够的。径向畸变通常使用多项式函数进行建模：

$$s_n=x_n^2+y_n^2,(10)$$

$$k=\sum _{i=0}^{n_d}{d}_i{s}_n^i(d_0=1),(11)$$

$${\tilde{p}}_d=D({\tilde{p}}_n)=[x_d,y_d,1{]}^{\top }=[k{x}_n,k{y}_n,1{]}^{\top },(12)$$

$${\tilde{p}}_i\simeq K{\tilde{p}}_d=KD({\tilde{p}}_n)=KD(E{\tilde{p}}_w),(13)$$

其中，$d_i$ 是畸变参数，$D$ 是畸变函数。$n_d$ 是畸变参数的最大阶数，在现有的标定方法中通常小于四。

3.4 测量模型（控制点）

图像中黑点的质心可以通过以下方式获得：

$${\tilde{p}}_i=\frac{\int w{p}_{i}d{A}_i}{\int wd{A}_i}\approx \frac{{\sum }_j{w}_j{p}_i(j)}{{\sum }_j{w}_j}(14)$$

在理想情况下，权重 $w$ 是 (255 - 强度)。然而在实际应用中，其他光源经常会扭曲投影圆的颜色值，因此为了鲁棒性，通常将 $w$ 设置为1。

4. 圆形图案的无偏估计器

4.1 矩方法

本节的主要目标是提供控制点的无偏估计器的详尽数学推导，该控制点是来自目标图案在图像中投影圆的中心点。在同态变换下跟踪控制点是很直观的。一个椭圆在同态变换 $H$ 下被投影到另一个椭圆，如下所示：
$${\tilde{p}}_f\simeq H{\tilde{p}}_i\text{\hspace{1em}(15)}$$
$$Q_f\simeq H^{-\top }{Q}_i{H}^{-1}\text{\hspace{1em}(16)}$$
下标 $i$ 表示原始帧，$f$ 表示变换帧。通过结合等式（3）和（16），变换后的椭圆中心计算为 $H{Q}_i^{-1}{H}^{\top }(0,0,1{)}^{\top }$，它与原始椭圆的变换中心点 $H{Q}_i^{-1}(0,0,1{)}^{\top }$ 不同。在我们采用这种方法之前，我们先定义一些几何学的初步概念。

定义1（形状）。形状 $A_{xy}$ 是由 $xy$ 平面上的封闭曲线包围的点集。$|A_{xy}|$ 是封闭曲线的内部区域。$A_k$ 表示 $A_{x_k{y}_k}$。

圆锥特征在非线性变换（如畸变）下会失去其特征。为了在畸变下构建一个无偏估计器，我们利用矩理论。在概率论中，随机变量 $X$ 和 $Y$ 的 $(i+j)$ 阶矩定义为 $E[X^i{Y}^j]$。假设给定形状的点服从均匀分布，则空间平均值对应于随机变量的期望值。我们在 $xy$ 坐标中定义一个二维形状的 $(i+j)$ 阶矩。

定义2 (矩)。对于任何 $m,n\in {\mathbb{Z}}^{\ast }$ (非负整数集)，$M_{m,n}^{xy}$ 是 $A_{xy}$ 的 $(m+n)$ 阶矩，定义为：

$$M_{xy}^{m,n}\triangleq \frac{1}{|A_{xy}|}\int x^m{y}^{n}d{A}_{xy}\text{\hspace{1em}(17)}$$

一阶矩因子（$M_{1,0}$ 和 $M_{0,1}$）代表形状的中心点，而二阶矩因子（$M_{2,0}$, $M_{1,1}$, 和 $M_{0,2}$）与形状的协方差有关。因此，一个椭圆可以仅通过一阶和二阶矩因子来描述。椭圆的矩阵形式是 $x^{\top }Qx\le 0$，而矩形式为 $[M_{1,0},M_{0,1},M_{2,0},M_{1,1},M_{0,2}]$。每种表示都具有相同的五个自由度（DoF）。

定理1。给定一个在域 $X$ 中可逆的多项式函数 $D:(x,y)\to (x^{\prime },y^{\prime })$，设 $A_{xy}\subset X$ 通过 $D$ 变换为 $A_{x^{\prime }{y}^{\prime }}$。对于任何 $m,n\in {\mathbb{Z}}^{\ast }$，存在 $c_{ij}(D)\in \mathbb{R}$ 和 $p,q\in {\mathbb{Z}}^{\ast }$，使得：

$$|A_{x^{\prime }{y}^{\prime }}|M_{x^{\prime }{y}^{\prime }}^{m,n}=|A_{xy}|\sum _{i=0}^p\sum _{j=0}^q{c}_{ij}{M}_{xy}^{i,j}$$

证明：见附录7.1。

定理1表明，对于任何多项式映射，变换后的形状的矩总是可以通过原始形状的矩的线性组合来表示。所需的 $p$ 和 $q$ 的数量取决于多项式函数的阶数。

4.2. 追踪畸变下的控制点

尽管无法以任何解析形式描述畸变的椭圆，但根据定理1，它的矩表示总是存在的。尽管可以使用高阶矩，但我们只使用一阶矩，即形状的质心。其原因将在下一节4.3中解释。

为了检查畸变映射下的形状，设$A_n$为归一化平面上的一个形状，由点$(x_n,y_n)$组成，$A_d$为畸变后的形状，由点$(x_d,y_d)$组成。

使用定理1，可以在归一化平面上计算畸变椭圆的中心，如下所示：

$$w_{0r}=\sum _i(2i+1)d_i{d}_{r-i}\text{\hspace{1em}(18)}$$

$$w_{1r}=\sum _i\sum _j(2i+1)d_i{d}_j{d}_{r-i-j}\text{\hspace{1em}(19)}$$

$$1=M_d^{0,0}=\frac{|A_n|}{|A_d|}\sum _{r=0}^{2n_d}{w}_{0r}\left[\frac{1}{|A_n|}\int s_n^{r}d{A}_n\right]\text{\hspace{1em}(20)}$$

$${\overline{x}}_d=M_d^{1,0}=\frac{|A_n|}{|A_d|}\sum _{r=0}^{3n_d}{w}_{1r}\left[\frac{1}{|A_n|}\int x_n{s}_n^{r}d{A}_n\right]\text{\hspace{1em}(21)}$$

$${\overline{y}}_d=M_d^{0,1}=\frac{|A_n|}{|A_d|}\sum _{r=0}^{3n_d}{w}_{1r}\left[\frac{1}{|A_n|}\int y_n{s}_n^{r}d{A}_n\right]\text{\hspace{1em}(22)}$$

参考等式(10)，其中符号$s_n$定义为 $(x_n^2+y_n^2)$。等式(20)用于计算 $\frac{|A_n|}{|A_d|}$。$2n_d$和$3n_d$来自于计算黎曼度量，详细信息见附录8.2。

计算图像平面上投影形状的质心可以如下简化计算，详细过程见附录8.3。

质心的x坐标计算公式为：
$${\overline{x}}_i=f_x{\overline{x}}_d+\eta {\overline{y}}_d+c_x\text{\hspace{1em}(23)}$$

质心的y坐标计算公式为：
$${\overline{y}}_i=f_y{\overline{y}}_d+c_y.\text{\hspace{1em}(24)}$$

因此，为了构建$[{\overline{x}}_i,{\overline{y}}_i]$的无偏估计量，我们只需要对每个整数$r$从0到$3n_d$计算以下三个积分：
$$\frac{1}{|A_n|}\int s^{r}d{A}_n,\frac{1}{|A_n|}\int x_n{s}^{r}d{A}_n,\text{和}\frac{1}{|A_n|}\int y_n{s}^{r}d{A}_n,$$
其中$(x_n,y_n)$是集合 A n = { p n ∣ p ~ n T Q n p ~ n ≤ 0 } A_n = \{p_n | \tilde{p}_n^T Q_n \tilde{p}_n \leq 0\} An​={pn​∣p~​nT​Qn​p~​n​≤0}中的一个点。

我们将这些值定义为一个向量$v_{rn}$：
$$v_n^r\triangleq \left[\begin{array}{c}\frac{1}{|A_n|}\int x_n{s}^{r}d{A}_n\\ \frac{1}{|A_n|}\int y_n{s}^{r}d{A}_n\\ \frac{1}{|A_n|}\int s^{r}d{A}_n\end{array}\right]\text{\hspace{1em}(25)}$$

直接从归一化平面上的椭圆计算$v_{rn}$并不简单；因此，我们将这个过程分为两步，如图3所示。

图3 任意椭圆$Q_n$的矩计算策略。
利用定理2，可以从未旋转的标准椭圆$Q_s$的矩得到旋转后的椭圆$Q_n$的矩。而标准椭圆$Q_s$的矩，则可以通过位于原点的另一个标准椭圆$Q_0$的矩，利用定理5来得到。

定理2（旋转等变性）。对于任何二维旋转变换$R:(x_s,y_s)\to (x_n,y_n)$，都存在一个$\alpha \in [0,2\pi ]$，使得：

$$x_n=\cos \alpha x_s-\sin \alpha y_s$$

$$y_n=\sin \alpha x_s+\cos \alpha y_s$$

$$v_n^r=\left(\begin{array}{ccc}\cos \alpha & -\sin \alpha & 0\\ \sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right)v_s^r\text{\hspace{1em}(26)}$$

证明：见附录7.2。

定理2意味着，如果我们能计算出任意未旋转椭圆$A_s$（其长轴和短轴平行于坐标系轴）的向量$v_s^r$，那么就有可能获得任意椭圆$A_n$的$v_n^r$。

引理3。考虑$I^{m,n}=\frac{1}{2\pi }{\int }_0^{2\pi }{\cos }^m\theta {\sin }^n\theta d\theta$，其中$m,n,i,j\in {\mathbb{Z}}^{\ast }$，那么：

$$I^{m,n}=\{\begin{array}{ll}\frac{\left(\genfrac{}{}{0ex}{}{2i+2j}{i+j}\right)\left(\genfrac{}{}{0ex}{}{i+j}{i}\right)}{\left(\genfrac{}{}{0ex}{}{2i+2j}{2i}\right)2^{2i+2j}}& \text{if }m=2i,n=2j\\ 0& \text{otherwise}\end{array}$$

证明：见附录7.3。
（注意，这里和后面的形如$\left(\genfrac{}{}{0ex}{}{n}{k}\right)$的符号，代表组合运算$C(n,k)=\frac{n!}{k!(n-k)!}$。）

引理4。$A_0=\{(x_0,y_0)|(\frac{x_0}{a}{)}^2+(\frac{y_0}{b}{)}^2\le 1\}$的$(m+n)$阶矩为：

$$M_0^{m,n}=\frac{a^m{b}^n}{1+(m+n)/2}{I}^{m,n}$$

证明：见附录7.4。

我们可以利用引理3和引理4以及定理5来获得$v_s^r$。$M_0^{m,n}$表示位于原点的未旋转椭圆的$(m+n)$阶矩的解析解。

定理5（$v_s^r$的解）。当$x_s$和$y_s$满足$((\frac{x_s-t_x}{a}{)}^2+(\frac{y_s-t_y}{b}{)}^2\le 1)$时，

$$v_s^r[0]=\sum _{i=0}^r\sum _{j=0}^{r-i}{M}_0^{2i,2j}\sum _{k=i}^{r-j}\left(\genfrac{}{}{0ex}{}{r}{k}\right)\left(\genfrac{}{}{0ex}{}{2k+1}{2i}\right)\left(\genfrac{}{}{0ex}{}{2r-2k}{2j}\right)t_x^{2k-2i+1}{t}_y^{2r-2k-2j}$$

$$v_s^r[1]=\sum _{i=0}^r\sum _{j=0}^{r-i}{M}_0^{2i,2j}\sum _{k=i}^{r-j}\left(\genfrac{}{}{0ex}{}{r}{k}\right)\left(\genfrac{}{}{0ex}{}{2k}{2i}\right)\left(\genfrac{}{}{0ex}{}{2r-2k+1}{2j}\right)t_x^{2k-2i}{t}_y^{2r-2k-2j+1}$$

$$v_s^r[2]=\sum _{i=0}^r\sum _{j=0}^{r-i}{M}_0^{2i,2j}\sum _{k=i}^{r-j}\left(\genfrac{}{}{0ex}{}{r}{k}\right)\left(\genfrac{}{}{0ex}{}{2k}{2i}\right)\left(\genfrac{}{}{0ex}{}{2r-2k}{2j}\right)t_x^{2k-2i}{t}_y^{2r-2k-2j}$$

证明：见附录7.5。

我们已经通过使用动态规划降低了计算$v_s^r$的计算成本。系数包括二项式的组合，而$I^{m,n}$对于圆锥曲线形状是不变的；因此，可以提前计算并存储这些值。因此，可以在$O(1)$时间内获得$v_s^r$。

$v_n^r$可以通过将旋转矩阵乘以$v_s^r$来简单获得，如定理2所述。使用前面的推导来计算相关等式，无偏估计器的整个过程可以被描述为算法1。

4.3. 第一矩的稳健性

如第4.1节所述，圆锥曲线是用二阶矩定义的。利用上述结果，可以计算出图像平面上畸变圆锥曲线的二阶矩。然而，边界模糊效应很容易污染高阶矩。例如，如果椭圆存在一些膨胀或侵蚀，其长轴和短轴的长度会变短或变长，而形状的质心则保持不变。对于标定来说，无偏估计和精确测量都是必不可少的；因此，仅利用第一矩更有利于精确标定。第一矩的另一个优点是对图像噪声的稳健性。

假设形状边界点存在一些噪声，且噪声服从均值为零、方差为${\sigma }^2$的正态分布，那么形状第一矩的方差将减少$\frac{1}{n}$。对于遵循$X_i\sim N(\mu ,{\sigma }^2)$的边界点，第一矩$M_1^X$及其方差的计算公式如下：

$$M_1^X=\frac{1}{n}\sum _{i=1}^n{X}_i（第一矩）\text{\hspace{1em}(27)}$$

$$Var(M_1^X)=Var\left(\frac{1}{n}\sum _i{X}_i\right)\text{\hspace{1em}(28)}$$

$$=\frac{1}{n^2}Var\left(\sum _i{X}_i\right)=\frac{1}{n^2}\cdot n\cdot Var(X_i)=\frac{{\sigma }^2}{n}$$

这是圆形图案比棋盘格对边界模糊效应更稳健的原因之一，因为棋盘格的控制点是直接从单点获得的。这一发现将通过第5.1节中的一系列实验得到验证。

4.4. 使用无偏估计器进行标定

标定过程分为两个阶段。首先，我们使用Zhang[28]的方法，在封闭形式中建立无畸变假设下的内参初始值。由于初始值没有考虑畸变参数，因此质量较差，需要通过优化进行细化。

在优化过程中，我们最小化了重投影误差，该误差是观测到的控制点位置与使用无偏估计器估计的控制点位置之间的平方差，如下所示：

$$K,D=\underset{K,D{,}^{1,2,\dots n}}{\mathrm{argmin}}\sum _{j=1}^n\sum _{k=1}^m{\Vert {\overline{p}}_i^{jk}-{\hat{p}}_i^{jk}\Vert }^2,\text{\hspace{1em}(29)}$$

$${\hat{p}}_i^{jk}=\text{UnbiasedEstimator}(K,E^j,Q_w^{jk}),\text{\hspace{1em}(30)}$$

其中，${\overline{p}}_i^{jk}$ 是第$j$张图像中的第$k$个控制点，而${\hat{p}}_i^{jk}$ 是与$p_i^{jk}$ 对应的估计控制点。与目标平面中对应于$p_i^{jk}$ 的圆是$Q_w^{jk}$。上述优化问题使用Ceres Solver[1]解决。

5. 结果

5.1. 估计器的比较

我们将我们的无偏估计器与棋盘格方法以及其他现有的估计器进行了比较，如基于点的、基于圆锥曲线的和数值方法[13]，定义如下：

$$基于点的：\widetilde{{\hat{p}}_i}=KD(E{Q}_w^{-1}[0,0,1{]}^{\top }),\text{\hspace{1em}(31)}$$

$$基于圆锥曲线的：\widetilde{{\hat{p}}_i}=KD(E{Q}_w^{-1}{E}^{\top }[0,0,1{]}^{\top }),\text{\hspace{1em}(32)}$$

$$数值方法：\widetilde{{\hat{p}}_i}=\text{Eq.}(19)\text{ in }[13],\text{\hspace{1em}(33)}$$

$$我们的方法：算法1$$

Numerical(n)表示迭代步数为n。

大多数现有的标定算法使用基于点的估计器，它忽略了圆锥曲线几何，并直接将扭曲图像中的形状中心与目标中的圆心相匹配。基于圆锥曲线的估计器补偿了由透视变换引入的偏差。然而，基于圆锥曲线的估计器仍然获得畸变偏差。数值方法在理论上是无偏的，但也显示出由数值积分引起的相当大的误差。这些限制在图4中得到了很好的描述。每个数据点表示在固定的24个场景中的重投影误差的平均值，我们用标准偏差绘制了误差图。

图4. 合成图像中的重投影误差比较。
我们检查了每种方法的重投影误差的一致性和大小，同时改变了圆的半径和表示为$d_1$的畸变量。在图表中，每条曲线代表平均重投影误差，包络线代表相应的标准差边界。第一行使用原始图像，第二行将高斯模糊应用于图像以检查鲁棒性。与其他圆形图案方法不同，我们的方法是无偏的，无论畸变和半径如何变化，误差都接近零。

我们基于矩的无偏估计器无论半径大小和畸变如何，都保持了非常小的重投影误差。然而，其他方法的误差较大，而像圆锥曲线和点这样的圆形图案方法的误差会随着半径或畸变的增加而逐渐增加。尽管畸变量对棋盘格方法的结果影响不大，但与圆形图案相比，它由于不准确的控制点检测器而产生的固有误差。棋盘格的控制点是正方形的角，由于像素的不连续性，角位置应该通过图案边界附近的插值来近似。相比之下，圆形图案的控制点是通过计算圆内数千个内点的平均值获得的，从而得到小数点后精确的位置。

由于基于圆锥曲线和点的估计器具有由圆锥几何形状引起的偏差，因此随着圆的大小和畸变的增加，误差也会增加。这两个估计器有略微不同的趋势。基于圆锥曲线的估计器由于畸变而非透视变换而产生偏差。因此，重投影误差仅取决于畸变系数的绝对值。相比之下，基于点的估计器同时具有透视和畸变偏差。透视偏差导致误差图中标准差较高且不对称。

5.2. 合成图像的标定精度

为了评估我们的无偏估计器对标定性能的提升，我们在合成图像上开展了一系列实验。我们设定了焦距$f_x,f_y$和主点坐标$c_x,c_y$的真实值（Ground Truth, GT），以确保视场（Field of View, FOV）大约为90度。我们为畸变系数准备了两种场景，基于公式（11）：一种是低畸变场景（$d_1=-0.2$），另一种是高畸变场景（$d_1=-0.4$）。在高畸变情况下，畸变函数在给定的FOV范围内可能不可逆，因此我们稍微调整了$d_2$以确保其可逆性。在我们的分析中，由于$d_2$的值极小，几乎可以忽略不计，因此我们主要评估$d_1$。我们生成了100张图像，并随机选择30张进行标定。为了确保结果的可靠性，我们重复了这一过程30次，并计算出平均值和标准差，如表1所示。

表1. 合成图像的标定结果。
此表显示了标定结果的平均值和标准差。我们生成了100张合成图像，并使用每种方法从30张随机选择的图像中获得了内部参数。通过重复此操作30次，我们能够计算出与每种方法对应的标定结果的平均值和标准差。我们的方法显示的结果最接近真实值，且方差远低于其他方法。同时，其他圆形图案方法，如基于点和基于圆锥曲线的方法，尽管使用了与我们相同的控制点，但却收敛到了不合适的值。更多的迭代次数可以使数值方法表现得无偏，但计算时间会显著增加。棋盘格方法也是无偏的，收敛到真实值；然而，由于测量噪声，其方差很大。值得注意的是，当图像质量由于高斯模糊而降低时，测量精度会进一步恶化，并且几乎有一半的图像无法被检测到。

实验结果显示，在所有测试场景中，我们的方法均表现最佳。我们的方法不仅提供了准确的平均值，还展现出了极低的方差，这是其特别的价值所在。由于棋盘格图案的估计器是无偏的，在无噪声干扰的情况下，其标定结果能够非常接近真实值。然而，由于控制点测量的不精确性，标定结果呈现出显著的方差，这一点在第5.1节的实验中已得到验证。当图像中引入噪声时，这种情况会进一步恶化，控制点检测在许多情况下会失败，从而导致性能下降。

对于使用圆形图案的基于点和基于圆锥曲线的方法，虽然其测量结果是准确且稳健的，但由于估计器存在偏差，其性能相较于棋盘格方法仍有所不足。随着畸变程度的增加，这种偏差会变得更加明显，从而导致标定结果的不理想。数值方法需要在精度和运行时间之间做出权衡，而且为了获得有意义的结果，往往需要花费大量时间。相比之下，我们的方法有效地解决了这些限制，展现出了比棋盘格方法和其他方法更高的精度、效率和稳健性。值得强调的是，本文的一个重要贡献是证明了圆形图案的价值。尽管圆形图案具有更多的信息特征，但由于算法上的限制，其应用一直受到忽视。

5.3. 真实图像的标定精度

我们进一步使用从RGB和TIR摄像机捕获的真实图像来评估我们的方法。由于无法找到真实摄像机的GT（真值）内部参数，因此需要将实验扩展到现实世界，这需要另一种评估指标。我们没有直接比较内部参数，而是评估了不同场景中每个目标之间的重投影误差分布以及相对平移和旋转值。为了突出各方法之间的明显差异，我们使用了两种具有高畸变的摄像机。一种是分辨率为1200×930的RGB摄像机Trition 5.4 MP。另一种是分辨率为640×512的TIR摄像机FLIR A65。这两个摄像机的样本图像如图5所示。由于TIR在识别彩色图案方面受到限制，只能通过红外能量来区分物体，因此我们使用了一种由具有不同导热性的正方形或圆形组成的印刷电路板（PCB），这是在[24]中引入的。

重投影误差分布。当总重投影误差较小且误差分布均匀时，我们期望标定结果更好。我们从三个不同的距离（即近、中、远）收集了24张图像，以研究这两个方面。如图5所示，基于点和基于圆锥曲线的方法导致高重投影误差，以及不同距离之间误差值的显著差异，这意味着计算出的内部参数不能在所有距离上一致地解释投影模型。相反，棋盘格方法和我们的方法，其无偏估计器在所有距离上都显示出较低的重投影误差。棋盘格方法和我们的方法在RGB摄像机上表现出相似的性能。同时，由于我们的方法能够准确检测控制点，因此在TIR摄像机上优于棋盘格方法。如第5.2节所示，随着迭代步骤的增加，数值方法的准确性收敛于我们的方法。Numerical(1600)显示出与我们的方法相似的结果，但收敛需要近五分钟。更多信息，如2D图像上误差向量的分布，见附录9.2。

图5. 真实世界实验：重投影误差。
第一行显示了来自（a）RGB和（b）TIR摄像机的样本图像。第二行显示了重投影误差分布除以摄像机到目标的距离。有偏方法的误差分布，如基于点和圆锥曲线的方法，随距离的变化而变化很大。无偏方法，如棋盘格和我们的方法，显示出较低且均匀的误差分布。在TIR图像中，由于圆形图案的鲁棒性控制点，我们的方法明显优于棋盘格方法。

6D位姿误差。为了验证每种方法的反投影性能，我们研究了相对的6-DoF位姿误差。每个摄像机捕获了20个不同的场景，我们用这些图像进行了标定。在标定过程中，还优化了每个场景中目标的位姿，这些位姿用于评估标定精度。相对位姿的真值是通过OptiTrack的运动捕捉系统获得的。设$T_{t_i}^c$是从第$i$个目标坐标到摄像机坐标的$SE(3)$矩阵，$T_{o_i}^m$是从第$i$个物体坐标到运动捕捉系统坐标的$SE(3)$矩阵。请注意，对象框架是附加到目标平面上的框架。但是，由于运动捕捉系统为目标定义了另一个坐标系，如图6所示，因此这两个框架并不相同。

如果没有坐标变换矩阵$T_t^o$和$T_c^m$，就不可能直接评估相对位姿。因此，我们首先计算变换矩阵，最终误差定义为：

$$\text{error}=\sum _{i=1}^n||T_{o_i}^m\ominus {\hat{T}}_c^m{T}_{t_i}^c({\hat{T}}_t^o{)}^{-1}||(34)$$

其中，${\hat{T}}_c^m$和${\hat{T}}_t^o$是估计值。估计$T_t^o$和$T_c^m$的详细信息见附录9.3。根据表2，我们的方法在所有情况下都表现最佳，而棋盘格方法在反投影任务中性能较低，尽管与基于点和基于圆锥曲线的方法相比，其重投影误差较低。

表2. 真实世界实验：6D姿态误差。
从校准结果和运动捕捉系统获得的目标姿态之间的旋转和平移误差。我们的方法在旋转部分优于其他方法。在RGB摄像机上，棋盘格方法在平移部分显示出相当的结果；然而，它在TIR摄像机上出现了故障。

图6. 实验装置。
我们利用动作捕捉系统来确定目标的真实6D姿态。满足$T_c^m{T}_t^c=T_o^m{T}_t^o$的关系。

6. 结论

在本文中，我们介绍了圆锥曲线的第n阶封闭形式矩，并证明了扭曲椭圆的质心可以通过原始矩的线性组合来表示。基于这些结果，我们开发了一种用于摄像机校准的扭曲下圆形图案的无偏估计器。使用这个估计器，基于圆形图案的校准克服了其算法限制，并且表现优于基于棋盘格的校准。

附录（略）