本文将介绍张正友于 2000 年提出的一种被广泛使用的相机标定方法。

使用该方法时，标定物可视为一个平面，而世界坐标系到图像坐标系的映射可以简化为单应矩阵（Homography）。

相机投影矩阵 P 的简化

在上一篇中，我们已经知道如何将世界坐标系中的三维点映射为图像中的像素：

X~i∼K[R∣t]X~w=K[r11r12r13txr21r22r23tyr31r32r33tz]X~w=K[r1r2r3t][XwYwZw1]\begin{equation} \begin{aligned} \bm{\tilde{X}_{\text{i}}} & \sim K \begin{bmatrix} R | \bm{t} \end{bmatrix} \bm{\tilde{X}_{\text{w}}} \\ & = K \begin{bmatrix} r_{11} & r_{12} & r_{13} & t_x \\ r_{21} & r_{22} & r_{23} & t_y \\ r_{31} & r_{32} & r_{33} & t_z \\ \end{bmatrix} \bm{\tilde{X}_{\text{w}}} \\ & = K \begin{bmatrix} \bm{r_{1}} & \bm{r_{2}} & \bm{r_{3}} & \bm{t} \end{bmatrix} \begin{bmatrix} X_{\text{w}} \\ Y_{\text{w}} \\ Z_{\text{w}} \\ 1 \end{bmatrix} \end{aligned} \end{equation}

其中 X~i=(x,y,1)T\bm{\tilde{X}_{\text{i}}} = (x, y, 1)^T，X~w=(Xw,Yw,Zw,1)T\bm{\tilde{X}_{\text{w}}} = (X_{\text{w}}, Y_{\text{w}}, Z_{\text{w}}, 1)^T。符号 ∼\sim 表示两个齐次坐标只相差一个非零尺度。

在张正友标定法中，标定物可视为平面。我们把世界坐标系建立在标定板上，使标定平面满足 Zw=0Z_{\text{w}}=0。只有在这个坐标系选择下，才能把 X~w\bm{\tilde{X}_{\text{w}}} 中的 ZwZ_{\text{w}} 在计算中省略；相应地，矩阵中的 r3\bm{r_3} 也可去除。于是有

X~i∼K[r1r2r3t][XwYwZw1]=K[r1r2t][XwYw1]\begin{equation} \begin{aligned} \bm{\tilde{X}_{\text{i}}} & \sim K \begin{bmatrix} \bm{r_{1}} & \bm{r_{2}} & \red{\sout{\bm{r_{3}}}} & \bm{t} \end{bmatrix} \begin{bmatrix} X_{\text{w}} \\ Y_{\text{w}} \\ \red{\sout{Z_{\text{w}}}} \\ 1 \end{bmatrix} & = K \begin{bmatrix} \bm{r_{1}} & \bm{r_{2}} & \bm{t} \end{bmatrix} \begin{bmatrix} X_{\text{w}} \\ Y_{\text{w}} \\ 1 \end{bmatrix} \end{aligned} \end{equation}

用 MM 表示 (Xw,Yw)T(X_{\text{w}}, Y_{\text{w}})^T，M~\tilde{M} 表示 (Xw,Yw,1)T(X_{\text{w}}, Y_{\text{w}}, 1)^T，则

H∼K[r1r2t]\begin{equation} H \sim K \begin{bmatrix} \bm{r_{1}} & \bm{r_{2}} & \bm{t} \end{bmatrix} \end{equation} X~i∼HM~\begin{equation} \bm{\tilde{X}_{\text{i}}} \sim H\tilde{M} \end{equation}

这里的矩阵 HH 是一个单应矩阵，求解 (4)(4) 即可得到 HH 的估计。

通过调用 cv::findHomography() 之类的函数即可轻松求得单应矩阵 HH。

因为 HH 有 8 个自由度，而每个点提供 2 个观测（xx 和 yy），所以至少需要 4 个点，并且其中任意 3 个点不能共线。实际标定通常使用更多角点，并通过归一化 DLT 或带鲁棒估计的方法计算 HH。

得到 HH 后，我们便可以进一步通过 HH 来得到 KK。

内参的约束

根据 (3)(3)，有

H=[h1h2h3]∼K[r1r2t]\begin{equation} H = \begin{bmatrix} \bm{h_1} & \bm{h_2} & \bm{h_3} \end{bmatrix} \sim K \begin{bmatrix} \bm{r_{1}} & \bm{r_{2}} & \bm{t} \end{bmatrix} \end{equation}

估计得到的 HH 只确定到尺度。选定 H=[h1 h2 h3]H=[\bm{h_1}\ \bm{h_2}\ \bm{h_3}] 的一个代表后，存在非零尺度 λ\lambda，使

[r1r2t]=λK−1H.\begin{equation} \begin{bmatrix}\bm{r_1}&\bm{r_2}&\bm{t}\end{bmatrix} =\lambda K^{-1}H. \end{equation}

因为 r1\bm{r_1} 与 r2\bm{r_2} 正交且各自都已单位化，即

r1⋅r2=r1Tr2=0\begin{equation} \bm{r_1}\cdot \bm{r_2} = \bm{r_1}^T\bm{r_2} = 0 \end{equation} ∥r1∥=∥r2∥=1\begin{equation} \lVert \bm{r_1} \rVert = \lVert \bm{r_2} \rVert = 1 \end{equation}

由上面的尺度关系和正交关系，有

λK−1[h1h2h3]=[r1r2t](λK−1h1)T(λK−1h2)=0\begin{equation} \begin{gather*} \lambda K^{-1} \begin{bmatrix} \bm{h_1} & \bm{h_2} & \bm{h_3} \end{bmatrix} = \begin{bmatrix} \bm{r_{1}} & \bm{r_{2}} & \bm{t} \end{bmatrix} \\ (\lambda K^{-1}\bm{h_1})^T (\lambda K^{-1}\bm{h_2}) = 0 \end{gather*} \end{equation} h1TK−TK−1h2=0\begin{equation} \begin{gather*} \blue{\bm{h_1}^T K^{-T} K^{-1}\bm{h_2} = 0} \end{gather*} \end{equation}

由上面的尺度关系和单位范数条件可知，两个向量的范数相等：

∥K−1h1∥2=∥K−1h2∥2=1λ2(K−1h1)T(K−1h1)=h1TK−TK−1h1(K−1h2)T(K−1h2)=h2TK−TK−1h2\begin{equation} \begin{gather*} \lVert K^{-1}\bm{h_1} \rVert ^2 = \lVert K^{-1}\bm{h_2} \rVert ^2 = \frac{1}{\lambda^2} \\ (K^{-1}\bm{h_1})^T (K^{-1}\bm{h_1}) = \bm{h_1}^T K^{-T} K^{-1}\bm{h_1} \\ (K^{-1}\bm{h_2})^T (K^{-1}\bm{h_2}) = \bm{h_2}^T K^{-T} K^{-1}\bm{h_2} \end{gather*} \end{equation} h1TK−TK−1h1−h2TK−TK−1h2=0\begin{equation} \blue{\bm{h_1}^T K^{-T} K^{-1}\bm{h_1} - \bm{h_2}^T K^{-T} K^{-1}\bm{h_2} = 0} \end{equation}

蓝色高亮的两个式子，便是给定一个单应矩阵时，内参所满足的两个基本约束。

定义 BB 为

B=K−TK−1\begin{equation} B = K^{-T} K^{-1} \end{equation}

上述约束可简写为

h1TBh2=0\begin{equation} \bm{h_1}^T \blue{B} \bm{h_2} = 0 \end{equation}

以及

h1TBh1−h2TBh2=0\begin{equation} \bm{h_1}^T \blue{B} \bm{h_1} - \bm{h_2}^T \blue{B} \bm{h_2} = 0 \end{equation}

估计相机标定矩阵 K

B=K−TK−1B=K^{-T}K^{-1} 本身总是对称正定矩阵，因此，只要从数据中估计出的 BB 仍然正定，便可以通过 Cholesky 分解恢复矩阵 KK：

B=AAT\begin{equation} B = AA^T \end{equation}

其中 A=K−TA=K^{-T}。换句话说，如果已知 BB 的一个正定代表，我们便可通过矩阵分解求出 KK，再消除其未知整体尺度。

因此我们先尝试通过 HH 得到 BB。

由于 BB 对称（矩阵与其转置之积总是对称矩阵），

B=[b11b12b13b12b22b23b13b23b33]\begin{equation} B = \begin{bmatrix} \blue{b_{11}} & \blue{b_{12}} & \blue{b_{13}} \\ b_{12} & \blue{b_{22}} & \blue{b_{23}} \\ b_{13} & b_{23} & \blue{b_{33}} \\ \end{bmatrix} \end{equation}

可写为 6 维向量：

b=(b11,b12,b13,b22,b23,b33)T\begin{equation} \bm{b} = (\blue{b_{11}, b_{12}, b_{13}, b_{22}, b_{23}, b_{33}})^T \end{equation}

展开正交约束，有

h1TBh2=[h11h21h31][b11b12b13b12b22b23b13b23b33][h12h22h32]=[h11h21h31][h12b11+h22b12+h32b13h12b12+h22b22+h32b23h12b13+h22b23+h32b33]=h11(h12b11+h22b12+h32b13)+h21(h12b12+h22b22+h32b23)+h31(h12b13+h22b23+h32b33)=h11h12b11+(h11h22+h21h12)b12+(h31h12+h11h32)b13+h21h22b22+(h31h22+h21h32)b23+h31h32b33=0\begin{equation} \begin{split} \bm{h_1}^T B \bm{h_2} &= \begin{bmatrix} h_{11} & h_{21} & h_{31} \end{bmatrix} \begin{bmatrix} b_{11} & b_{12} & b_{13} \\ b_{12} & b_{22} & b_{23} \\ b_{13} & b_{23} & b_{33} \\ \end{bmatrix} \begin{bmatrix} h_{12} \\ h_{22} \\ h_{32} \end{bmatrix} \\ &= \begin{bmatrix} h_{11} & h_{21} & h_{31} \end{bmatrix} \begin{bmatrix} h_{12}b_{11} + h_{22}b_{12} + h_{32}b_{13} \\ h_{12}b_{12} + h_{22}b_{22} + h_{32}b_{23} \\ h_{12}b_{13} + h_{22}b_{23} + h_{32}b_{33} \end{bmatrix} \\ &= h_{11}(h_{12}b_{11} + h_{22}b_{12} + h_{32}b_{13}) \\ &\quad + h_{21}(h_{12}b_{12} + h_{22}b_{22} + h_{32}b_{23}) \\ &\quad + h_{31}(h_{12}b_{13} + h_{22}b_{23} + h_{32}b_{33}) \\ &= h_{11}h_{12} \blue{b_{11}} + (h_{11}h_{22} + h_{21}h_{12}) \blue{b_{12}} \\ &\quad + (h_{31}h_{12} + h_{11}h_{32}) \blue{b_{13}} + h_{21}h_{22} \blue{b_{22}} \\ &\quad + (h_{31}h_{22} + h_{21}h_{32}) \blue{b_{23}} + h_{31}h_{32} \blue{b_{33}} \\ &= 0 \end{split} \end{equation}

我们可以将上式改写为更简洁的向量点积的形式：

[h11h12h11h22+h21h12h31h12+h11h32h21h22h31h22+h21h32h31h32]Tb=0\begin{equation} \begin{bmatrix} h_{11}h_{12} \\ h_{11}h_{22} + h_{21}h_{12} \\ h_{31}h_{12} + h_{11}h_{32} \\ h_{21}h_{22} \\ h_{31}h_{22} + h_{21}h_{32} \\ h_{31}h_{32} \\ \end{bmatrix} ^T \bm{b} = 0 \end{equation}

我们甚至还可以更简洁一点，记

vij=[h1ih1jh1ih2j+h2ih1jh3ih1j+h1ih3jh2ih2jh3ih2j+h2ih3jh3ih3j]\begin{equation} \bm{v_{ij}} = \begin{bmatrix} h_{1i}h_{1j} \\ h_{1i}h_{2j} + h_{2i}h_{1j} \\ h_{3i}h_{1j} + h_{1i}h_{3j} \\ h_{2i}h_{2j} \\ h_{3i}h_{2j} + h_{2i}h_{3j} \\ h_{3i}h_{3j} \end{bmatrix} \end{equation}

因此上面的正交约束可写为：

v12Tb=0\begin{equation} \bm{v_{12}}^T \bm{b} = 0 \end{equation}

类似地，根据等范数约束，有

v11Tb−v22Tb=0\begin{equation} \begin{split} \bm{v_{11}}^T \bm{b} - \bm{v_{22}}^T \bm{b} = 0 \end{split} \end{equation}

将上述两式合并在一个矩阵中，便得到了

Vi=[v12Tv11T−v22T]\begin{equation} V_i = \begin{bmatrix} \bm{v_{12}}^T \\ \bm{v_{11}}^T - \bm{v_{22}}^T \end{bmatrix} \\ \end{equation}

此时我们就得到单张图像的方程：

Vib=0\begin{equation} V_i \bm{b} = \bm{0} \end{equation}

假设有 nn 张图像，将这些矩阵堆叠成一个 2n×62n\times6 矩阵：

V=[V1V2⋯Vn]\begin{equation} V = \begin{bmatrix} V_1 \\ V_2 \\ \cdots \\ V_n \end{bmatrix} \end{equation}

现在需要求解如下线性方程组：

Vb=0\begin{equation} V \bm{b} = \bm{0} \end{equation}

显然，该线性方程组有平凡解 b=0\bm{b}=\bm{0}，但这样的平凡解并没有意义。因为 BB 本来就只能从这些齐次约束中确定到一个非零尺度，我们可施加归一化约束：

∥b∥=1\begin{equation} \lVert \bm{b} \rVert = 1 \end{equation}

在有噪声的情况下，实际求解的是

min⁡∥b∥=1∥Vb∥2.\begin{equation} \min_{\lVert\bm{b}\rVert=1}\lVert V\bm{b}\rVert^2. \end{equation}

对 VV 做 SVD，取最小奇异值对应的右奇异向量即可得到 b\bm{b}。由于 b\bm{b} 有 5 个尺度无关的自由度，而每张图像只提供 2 个约束，所以我们通常至少需要 3 张具有不同且非退化姿态的图像。标定的过程中，如果标定板姿态变化不足，方程可能退化或严重病态，这就是我们应该采集更多不同视角图像的原因。

得到 KK 后，可以从每张图像的单应矩阵 HH 中恢复该视角所对应的外参。

相机作为测量设备

求出相机标定矩阵 KK 和畸变参数后，我们可将图像坐标系中的一个像素映射为相机坐标系中的一条射线。如果输入是原始畸变图像上的像素，必须先用标定得到的畸变模型将其反畸变；以下用 X~i\bm{\tilde{X}_{\text{i}}} 表示反畸变后的像素齐次坐标。

若在相机坐标系中指定一个测量平面，还可进一步把该像素映射为平面上的一个点。

设测量平面为：

AXc+BYc+CZc+D=0AX_{\text{c}} + BY_{\text{c}} + CZ_{\text{c}} + D = 0

然后计算射线与平面的交点，即可得到该点在相机坐标系中的三维坐标。

投影射线的方向可以表示为：

d=K−1X~i\begin{equation} \bm{d} = K^{-1} \bm{\tilde{X}_{\text{i}}} \end{equation}

因此射线可以参数化为：

Xc(λ)=λd,λ>0\begin{equation} \bm{X_{\text{c}}}(\lambda) = \lambda\bm{d}, \quad \lambda > 0 \end{equation}

其中 λ\lambda 为标量。

将射线方程代入平面方程即可求出 λ\lambda。这里要求射线不与平面平行，即 Adx+Bdy+Cdz≠0Ad_x+Bd_y+Cd_z\neq0：

A(λdx)+B(λdy)+C(λdz)+D=0\begin{equation} A(\lambda d_x) + B(\lambda d_y) + C(\lambda d_z) + D = 0 \end{equation} λ∗=−DAdx+Bdy+Cdz\begin{equation} \lambda^* = -\frac{D}{Ad_x + Bd_y + Cd_z} \end{equation}

只有当 λ∗>0\lambda^*>0 时，交点才位于相机前方的射线上。此时交点为：

Xc∗=Xc(λ∗)=λ∗d\begin{equation} \bm{X_{\text{c}}}^* = \bm{X_{\text{c}}}(\lambda^*) = \lambda^*\bm{d} \end{equation}

如果已知把世界坐标变换到相机坐标的同一组外参 R,tR,\bm{t}，还可将该点从相机坐标系变换回世界坐标系：

Xw∗=RT(Xc∗−t)\begin{equation} \bm{X_{\text{w}}}^* = R^T (\bm{X_{\text{c}}}^* - \bm{t}) \end{equation}

因此，在测量平面及相机外参已知的前提下，可以把二维像素映射为该平面上的世界坐标点。单个像素本身只确定一条射线；如果没有已知平面、深度或其他视角，它不能唯一确定任意三维点。

参考文献

  1. 张正友标定法原文