三维几何投影


小孔成像

早在公元前四世纪,中国古代思想家墨子在《墨经》中就记录了人类历史上最早的小孔成像实验

“景倒,在午有端,与景长。说在端。”

当光线穿过暗室墙壁上的一个微小针孔时,来自物体顶部的光线沿直线穿过小孔射在暗室底部的墙壁上,来自底部的光线射在顶部,形成一个上下颠倒、左右相反的实像。这一成像原理也被称为针孔相机模型(Pinhole Model),具体叙述如下

设相机光心(针孔位置)位于三维空间的原点 OcO_c,规定方向

  • 向前观察的方向为 ZcZ_c 轴,又称光轴
  • 水平向右的方向为 XcX_c 轴
  • 竖直向下的方向为 YcY_c 轴

这样就构成了三维欧式空间中的相机坐标系,现假定距离光心 ff 处放置一块垂直于光轴的成像传感器平面 z=fz=f,我们称这个距离为焦距,该平面为焦平面。由相似三角形知识:三维空间中任意一个物理点 Pc(xc,yc,zc)P_c(x_c,y_c,z_c),其发出的光线经过光心,在焦平面上的成像点为 P(x,y,f)P(x,y,f),具体表示为

Pc(xc,yc,zc)⇒P(x,y,f)=P(fxczc,fyczc,f)P_c(x_c,y_c,z_c)\Rightarrow P(x,y,f)=P(f\dfrac{x_c}{z_c},f\frac{y_c}{z_c},f)

忽略掉 ZcZ_c 轴方向的 ff,这样就借助针孔相机模型实现了一个映射

C:R3(xc,yc,zc)→R2(x,y)C:\mathbb{R}^3(x_c,y_c,z_c)\rightarrow\mathbb{R}^2(x,y)

上述坐标映射可以通过添加齐次项记为

zc[xy1]=[f0000f000010][xcyczc1]z_c\begin{bmatrix}x\\y\\1\end{bmatrix}=\begin{bmatrix}f&0&0&0\\0&f&0&0\\0&0&1&0\end{bmatrix}\begin{bmatrix}x_c\\y_c\\z_c\\1\end{bmatrix}

称 (x,y,1)(x,y,1) 系为图像系,这种齐次化处理将十分有助于后续的坐标变换。

内参

成像传感器上能够接受的物理尺寸 (x,y)(x,y) 通常为毫米级,但计算机读取数字图像时以像素为单位,且像素坐标系原点 (0,0)(0,0)位于图像的左上角处。现在考虑从相机系坐标 (x,y)(x,y) 转化到像素系坐标 (u,v)(u,v)

设 xx 方向每毫米包含 mxm_x 个像素,yy 方向每毫米包含 mym_y 个像素。光轴与该平面的交点在像素坐标系下的坐标为 (cx,cy)(c_x,c_y),通常称其为主点(Principal Point)。现在来求主点的表达式,根据像素坐标系的定义可知点 (x,y)(x,y)

u=mxx+cx,v=myy+cyu=m_xx+c_x,\quad v=m_yy+c_y

移项后得到

cx=u−mxfxczc,cy=v−myfyczcc_x=u-m_xf\frac{x_c}{z_c},\quad c_y=v-m_yf\frac{y_c}{z_c}

其中定义 fx=mxf,fy=myff_x=m_xf,f_y=m_yf 为以像素为单位的等效水平焦距和垂直焦距。为了简化该表示方法,我们引入齐次坐标,可以得到下述等价矩阵乘法形式

zc[uv1]=[fx0cx0fycy001][xcyczc]=K[xcyczc]z_c\begin{bmatrix}u\\v\\1\end{bmatrix}=\begin{bmatrix}f_x&0&c_x\\0&f_y&c_y\\0&0&1\end{bmatrix}\begin{bmatrix}x_c\\y_c\\z_c\end{bmatrix}=\boldsymbol{K}\begin{bmatrix}x_c\\y_c\\z_c\end{bmatrix}

其中的 K\boldsymbol{K} 称为相机内参矩阵(Intrinsic Matrix),这是针孔相机模型中给的一个重要参数。

外参刚体变换

在自动驾驶或机器人工作站中,物体通常定义在统一的世界坐标系中,而相机会随车身或机械臂运动。因此还需要考虑从世界系到相机系的坐标变换,这可以借助三维刚体变换知识得到

Pc=RPw+tP_c=\boldsymbol{R}P_w+\boldsymbol{t}

其中 R∈SO(3),t∈R3\boldsymbol{R}\in SO(3),\boldsymbol{t}\in\mathbb{R}^3。借助 44 维齐次坐标,可以将上述变换记为

P~c=[Pc1]=[Rt0T1][Pw1]=TcwP~w\tilde{P}_c=\begin{bmatrix}P_c\\1\end{bmatrix}=\begin{bmatrix}\boldsymbol{R}&\boldsymbol{t}\\\boldsymbol{0}^T&1\end{bmatrix}\begin{bmatrix}P_w\\1\end{bmatrix}=\boldsymbol{T}_{cw}\tilde{P}_w

其中的 4×44\times 4 矩阵 Tcw\boldsymbol{T}_{cw} 称为相机外参矩阵(Extrinsic Matrix),通常也被记作 3×43\times 4 的形式 [R∣t][\boldsymbol{R}|\boldsymbol{t}]。此时有

Pc=[R∣t][Pw1]P_c=\begin{bmatrix}\boldsymbol{R}|\boldsymbol{t}\end{bmatrix}\begin{bmatrix}P_w\\1\end{bmatrix}

将相机内参矩阵和外参矩阵级联,就可以将一个三维世界点 Pw(xw,yw,zw)P_w(x_w,y_w,z_w) 直接投影至二维像素系,即有

zc[uv1]=K[R∣t][xcyczc1]=P[xwywzw1]z_c\begin{bmatrix}u\\v\\1\end{bmatrix}=\boldsymbol{K}[\boldsymbol{R}|\boldsymbol{t}]\begin{bmatrix}x_c\\y_c\\z_c\\1\end{bmatrix}=\boldsymbol{P}\begin{bmatrix}x_w\\y_w\\z_w\\1\end{bmatrix}

其中的 P=K[R∣t]∈R3×4\boldsymbol{P}=\boldsymbol{K}[\boldsymbol{R}|\boldsymbol{t}]\in\mathbb{R}^{3\times 4} 称为相机的投影矩阵(Projection Matrix)。

下表直观展示了相机内参和外参的区别。

内参 K\boldsymbol{K} 外参 [R∣t][\boldsymbol{R}|\boldsymbol{t}]
描述 相机自身 相机和世界关系
是否随相机运动改变 否 是
单位 像素 米/角度
作用 相机系到像素系 世界系到相机系
矩阵大小 3×33\times3 3×43\times4 或 4×44\times4

相机位姿

相机外参矩阵记为 4×44\times 4 的形式时,其逆矩阵可以称作相机位姿(Pose)。根据三维刚体运动,其同样可以写作旋转矩阵 Rc\boldsymbol{R}_c 和平移 C\boldsymbol{C} 的组合。根据定义有位姿的 3×43\times 4 形式

Pw=[Rc∣C][Pc1]P_w=\begin{bmatrix}\boldsymbol{R}_c|\boldsymbol{C}\end{bmatrix}\begin{bmatrix}P_c\\1\end{bmatrix}

现在考虑相机位姿和外参矩阵的换算,以 4×44\times 4 形式转换为例有

Tcw=Twc−1=[RcC0T1]−1\boldsymbol{T}_{cw}=\boldsymbol{T}_{wc}^{-1}=\begin{bmatrix}\boldsymbol{R}_c&\boldsymbol{C}\\\boldsymbol{0}^T&1\end{bmatrix}^{-1}

由于 Twc\boldsymbol{T}_{wc} 可以被拆分为

Twc=[RcC0T1]=[IC0T1][Rc00T1]\boldsymbol{T}_{wc}=\begin{bmatrix}\boldsymbol{R}_c&\boldsymbol{C}\\\boldsymbol{0}^T&1\end{bmatrix}=\begin{bmatrix}\boldsymbol{I}&\boldsymbol{C}\\\boldsymbol{0}^T&1\end{bmatrix}\begin{bmatrix}\boldsymbol{R}_c&\boldsymbol{0}\\\boldsymbol{0}^T&1\end{bmatrix}

代入上式有

[RcC0T1]−1=[Rc00T1]−1[IC0T1]−1=[Rc−100T1][I−C0T1]\begin{bmatrix}\boldsymbol{R}_c&\boldsymbol{C}\\\boldsymbol{0}^T&1\end{bmatrix}^{-1}=\begin{bmatrix}\boldsymbol{R}_c&\boldsymbol{0}\\\boldsymbol{0}^T&1\end{bmatrix}^{-1}\begin{bmatrix}\boldsymbol{I}&\boldsymbol{C}\\\boldsymbol{0}^T&1\end{bmatrix}^{-1}=\begin{bmatrix}\boldsymbol{R}_c^{-1}&\boldsymbol{0}\\\boldsymbol{0}^T&1\end{bmatrix}\begin{bmatrix}\boldsymbol{I}&-\boldsymbol{C}\\\boldsymbol{0}^T&1\end{bmatrix}

由于 Rc∈SO(3)\boldsymbol{R}_c\in SO(3),因此有 Rc−1=RcT\boldsymbol{R}_c^{-1}=\boldsymbol{R}_c^T,故整理有

Tcw=[Rt0T1]=[RcT−RcTC0T1]\boldsymbol{T}_{cw}=\begin{bmatrix}\boldsymbol{R}&\boldsymbol{t}\\\boldsymbol{0}^T&1\end{bmatrix}=\begin{bmatrix}\boldsymbol{R}_c^T&-\boldsymbol{R}_c^T\boldsymbol{C}\\\boldsymbol{0}^T&1\end{bmatrix}

因此有互换公式 R=RcT,t=−RcTC\boldsymbol{R}=\boldsymbol{R}_c^T,\boldsymbol{t}=-\boldsymbol{R}_c^T\boldsymbol{C},反之亦然。

一般而言,使用w2c矩阵指代相机外参矩阵,使用c2w矩阵指代位姿。

相机畸变


在上一节中,我们得到了一个从世界系到像素系的投影式

zc[uv1]=K[R∣t][xwywzw1]z_c\begin{bmatrix}u\\v\\1\end{bmatrix}=\boldsymbol{K}[\boldsymbol{R}|\boldsymbol{t}]\begin{bmatrix}x_w\\y_w\\z_w\\1\end{bmatrix}

然而,为了获得好的成像效果,人们通常在相机的前方添加透镜。透镜的加入对成像过程会产生新的影响

  1. 透镜自身的形状对光线传播的影响
  2. 在机械组装过程中,透镜和成像平面不可能完全平行

这些都会使得光线穿过透镜投影到成像面时的位置发生变化,从而产生畸变。

  • 径向畸变:在实际拍摄的照片中,摄像机的透镜往往使得真实环境中的一条直线在图片中变成了曲线。越靠近图像的边缘,这种现象越明显。由于实际加工制作的透镜往往是中心对称的,这使得不规则的畸变通常径向对称。主要分为两大类,桶形畸变和枕形畸变。
  • 切向畸变:在相机的组装过程中由于不能使得透镜和成像面严格平行,而引入切向畸变。

我们一般使用下述数学模型描述相机畸变

xd=x(1+k1r2+k2r4+k3r6)+2p1xy+p2(r2+2x2)x_d=x(1+k_1r^2+k_2r^4+k_3r^6)+2p_1xy+p_2(r^2+2x^2)

yd=y(1+k1r2+k2r4+k3r6)+2p1(r2+2y2)+2p2xyy_d=y(1+k_1r^2+k_2r^4+k_3r^6)+2p_1(r^2+2y^2)+2p_2xy

其中 x,yx,y 表示畸变发生时实际观察到的图像系坐标,xd,ydx_d,y_d 表示修正畸变后的图像系坐标,r=x2+y2r=x^2+y^2。上式包含了径向畸变与切向畸变项,具体有

  • k1,k2,k3k_1,k_2,k_3 描述了相机的径向畸变
  • p1,p2p_1,p_2 描述了相机的切向畸变

此外,还有一种不属于相机畸变的相机误差项——像素网格倾斜,即在内参矩阵 K\boldsymbol{K} 内添加倾斜系数 ss。

K=[fxscx0fycy001]\boldsymbol{K}=\begin{bmatrix}f_x&s&c_x\\0&f_y&c_y\\0&0&1\end{bmatrix}

在去除畸变时,通常采用上述 xd,ydx_d,y_d 的表达式进行修正,前提是我们知晓畸变系数。

理想相机标定


标定

大多数条件下,一个相机的内外参数等必须通过对实际物理点坐标的实验与计算才能得到,这个求解参数的过程就称之为相机标定(Camera Calibration)。对于一个理想相机模型有

zc[uv1]=[fx0cx0fycy001][r11r12r13t1r21r22r23t2r31r32r33t3][xwywzw1]z_c\begin{bmatrix}u\\v\\1\end{bmatrix}=\begin{bmatrix}f_x&0&c_x\\0&f_y&c_y\\0&0&1\end{bmatrix}\begin{bmatrix}r_{11}&r_{12}&r_{13}&t_1\\r_{21}&r_{22}&r_{23}&t_2\\r_{31}&r_{32}&r_{33}&t_3\end{bmatrix}\begin{bmatrix}x_w\\y_w\\z_w\\1\end{bmatrix}

但由于 R\boldsymbol{R} 正交,即 RTR=I\boldsymbol{R}^T\boldsymbol{R}=\boldsymbol{I},故旋转矩阵可以概括为三个方向的旋转量 (θx,y,z)(\theta_{x,y,z})。通过观察可以发现,上述方程中含有十个未知量,这些就是相机标定中的求解目标量。

  • fx,fy,cx,cyf_x,f_y,c_x,c_y:相机内参
  • θx,θy,θz\theta_{x},\theta_y,\theta_z:旋转外参
  • t1,t2,t3t_1,t_2,t_3:平移外参

三维标靶标定

可以发现,每个物理点 (xw(i),yw(i),zw(i))(x_w^{(i)},y_w^{(i)},z_w^{(i)}) 都能够提供一个二元约束 (u(i),v(i))(u^{(i)},v^{(i)}),因此求解该十元一次方程需要至少五个三维物理点。此时通常使用三维标靶,其每一个立体面上的靶点空间坐标已知,只需测量相机内像素系坐标并带入方程求解即可。

即使不考虑相机畸变,限于各类误差,上述方程的准确解总是无法求得,因此还需要利用最小二乘法等优化算法进行修正。更多内容限于篇幅不在此讲述。

张正友标定


单应矩阵

三维标靶通常较难制作,张正友博士于1999年提出了一种在二维棋盘形平面上实现相机标定的方法,现称张正友标定法。因其高效、便捷的特性,且能够对相机的径向畸变系数 k1,k2,k3k_1,k_2,k_3、内参倾斜系数 ss 等进行标定,现已成为大多数相机标定入门方法的标准。

考虑一个 zw=0z_w=0 的二维棋盘形平面,使用非理想的针孔相机模型,那么棋盘上的点坐标满足

zc[uv1]=[fxscx0fycy001][r11r12r13t1r21r22r23t2r31r32r33t3][xwyw01]z_c\begin{bmatrix}u\\v\\1\end{bmatrix}=\begin{bmatrix}f_x&s&c_x\\0&f_y&c_y\\0&0&1\end{bmatrix}\begin{bmatrix}r_{11}&r_{12}&r_{13}&t_1\\r_{21}&r_{22}&r_{23}&t_2\\r_{31}&r_{32}&r_{33}&t_3\end{bmatrix}\begin{bmatrix}x_w\\y_w\\0\\1\end{bmatrix}

此时可以将 R\boldsymbol{R} 按列表示为

zc[uv1]=K[r1r2t][xwyw1]z_c\begin{bmatrix}u\\v\\1\end{bmatrix}=\boldsymbol{K}\begin{bmatrix}\boldsymbol{r}_1&\boldsymbol{r}_2&\boldsymbol{t}\end{bmatrix}\begin{bmatrix}x_w\\y_w\\1\end{bmatrix}

因此定义单应矩阵(Homography Matrix)

H=K[r1r2t]=[h11h12h13h21h22h23h31h32h33]\boldsymbol{H}=\boldsymbol{K}\begin{bmatrix}\boldsymbol{r}_1&\boldsymbol{r}_2&\boldsymbol{t}\end{bmatrix}=\begin{bmatrix}h_{11}&h_{12}&h_{13}\\h_{21}&h_{22}&h_{23}\\h_{31}&h_{32}&h_{33}\end{bmatrix}

单应矩阵描述了两个平面之间的射影变换,在此例子中被看做二维棋盘到像素系的变换。引入单应矩阵 H\boldsymbol{H} 后,由于其作用于齐次坐标上,因此有

[uv1]∝H[xwyw1],H∼λH\begin{bmatrix}u\\v\\1\end{bmatrix}\propto\boldsymbol{H}\begin{bmatrix}x_w\\y_w\\1\end{bmatrix},\quad\boldsymbol{H}\sim\lambda\boldsymbol{H}

因此可以固定某一元素 hij=1h_{ij}=1,此时的单应矩阵相当于只有 88 个自由度。又因为每个棋盘点能够提供两个约束,即有二元一次方程组

u=h11xw+hh12yw+h13h31xw+hh32yw+h33,v=h21xw+hh22yw+h23h31xw+hh32yw+h33u=\frac{h_{11}x_w+h_{h12}y_w+h_{13}}{h_{31}x_w+h_{h32}y_w+h_{33}},\quad v=\frac{h_{21}x_w+h_{h22}y_w+h_{23}}{h_{31}x_w+h_{h32}y_w+h_{33}}

因此,当一张图片上的标定板角点数量等于 44 时,即可求得对应的 H\boldsymbol{H} 。当一张图片上的标定板角点数量大于 44时,利用最小二乘法回归最佳的矩阵 H\boldsymbol{H}。

内参矩阵求解

现在已经求得单应矩阵 H=[h1h2h3]=K[r1r2t]\boldsymbol{H}=\begin{bmatrix}\boldsymbol{h}_1&\boldsymbol{h}_2&\boldsymbol{h}_3\end{bmatrix}=\boldsymbol{K}\begin{bmatrix}\boldsymbol{r}_1&\boldsymbol{r}_2&\boldsymbol{t}\end{bmatrix},因此不难有

r1=K−1h1,r2=K−1h2\boldsymbol{r}_1=\boldsymbol{K}^{-1}\boldsymbol{h}_1,\quad \boldsymbol{r}_2=\boldsymbol{K}^{-1}\boldsymbol{h}_2

由于旋转矩阵 R\boldsymbol{R} 的正交性,有

r1Tr2=0,r1Tr1=r2Tr2=1\boldsymbol{r}_1^T\boldsymbol{r}_2=0,\quad \boldsymbol{r}_1^T\boldsymbol{r}_1=\boldsymbol{r}_2^T\boldsymbol{r}_2=1

代入可以得到

h1T(KT)−1K−1h2=0,h1T(KT)−1K−1h1=h2T(KT)−1K−1h2=1\boldsymbol{h}_1^T(\boldsymbol{K}^T)^{-1}\boldsymbol{K}^{-1}\boldsymbol{h}_2=0,\quad \boldsymbol{h}_1^T(\boldsymbol{K}^T)^{-1}\boldsymbol{K}^{-1}\boldsymbol{h}_1=\boldsymbol{h}_2^T(\boldsymbol{K}^T)^{-1}\boldsymbol{K}^{-1}\boldsymbol{h}_2=1

上式约束中均出现了由内参矩阵组成的矩阵 (KT)−1K−1(\boldsymbol{K}^T)^{-1}\boldsymbol{K}^{-1},记 B=(KT)−1K−1=(KKT)−1\boldsymbol{B}=(\boldsymbol{K}^T)^{-1}\boldsymbol{K}^{-1}=(\boldsymbol{KK}^T)^{-1} 是一个对称矩阵。因此我们可以尝试先求出矩阵 B\boldsymbol{B},这样就顺便求得了相机的内参矩阵 K\boldsymbol{K}。

将新矩阵 B\boldsymbol{B} 代入约束关系式可以得到

h1TBh2=0,h1TBh1=h2TBh2=1\boldsymbol{h}_1^T\boldsymbol{B}\boldsymbol{h}_2=0,\quad \boldsymbol{h}_1^T\boldsymbol{B}\boldsymbol{h}_1=\boldsymbol{h}_2^T\boldsymbol{B}\boldsymbol{h}_2=1

利用二次型的展开,设有

vij=[h1ih1jh1ih2j+h2ih1jh2ih2jh1ih3j+h3ih1jh2ih3jh3ih2j+h3ih3j],b=[b11b12b22b13b23b33]\boldsymbol{v}_{ij}=\begin{bmatrix}h_{1i}h_{1j}\\h_{1i}h_{2j}+h_{2i}h_{1j}\\h_{2i}h_{2j}\\h_{1i}h_{3j}+h_{3i}h_{1j}\\h_{2i}h_{3j}\\h_{3i}h_{2j}+h_{3i}h_{3j}\end{bmatrix},\quad \boldsymbol{b}=\begin{bmatrix}b_{11}\\b_{12}\\b_{22}\\b_{13}\\b_{23}\\b_{33}\end{bmatrix}

上述方程可表示为

{vijTb=0v11Tb=v22Tb=1⇒Vb=0\begin{cases}\boldsymbol{v}_{ij}^T\boldsymbol{b}=0\\\boldsymbol{v}_{11}^T\boldsymbol{b}=\boldsymbol{v}_{22}^T\boldsymbol{b}=1\end{cases}\Rightarrow \boldsymbol{Vb}=0

其中的向量 vij\boldsymbol{v}_{ij} 由已知矩阵 H\boldsymbol{H} 的元素给出,因此只需解出向量 b\boldsymbol{b} 就可以得到矩阵 B\boldsymbol{B}。这里的 b\boldsymbol{b} 只有 66 个元素,但由于 B\boldsymbol{B} 是一个对称矩阵,且每张标定棋盘图能够提供一个二元约束关系。因此,我们只要取 33 张标定板照片即可求解向量 b\boldsymbol{b} 和矩阵 B\boldsymbol{B}。当照片更多时,可以采用最小二乘法拟合。

利用 B\boldsymbol{B} 和 K\boldsymbol{K} 的关系,可以求得相机内参矩阵。

外参矩阵求解

当求得 K,H\boldsymbol{K},\boldsymbol{H} 时,利用单应矩阵的性质有

[r1r2t]=K−1[h1h2h3]\begin{bmatrix}\boldsymbol{r}_1&\boldsymbol{r}_2&\boldsymbol{t}\end{bmatrix}=\boldsymbol{K}^{-1}\begin{bmatrix}\boldsymbol{h}_1&\boldsymbol{h}_2&\boldsymbol{h}_3\end{bmatrix}

由于 R\boldsymbol{R} 需要满足旋转矩阵的性质,因此可以设

r3=r1×r2\boldsymbol{r}_3=\boldsymbol{r}_1\times \boldsymbol{r}_2

由于棋盘标定板 zw=0z_w=0,因此旋转矩阵的第三列没有起到坐标转换作用,这也解释了为什么我们只能显式的求得 r1,r2\boldsymbol{r}_1,\boldsymbol{r}_2。

畸变参数

最后考虑一个由径向畸变参数引入的相机畸变模型。根据像素系的定义有

u=mxx+cx,v=myy+cyu=m_xx+c_x,\quad v=m_yy+c_y

其中的 (u,v)(u,v) 为理想的无畸变像素坐标,并把倾斜系数的影响忽略不计。畸变后的像素坐标 (u^,v^)(\hat{u},\hat{v}) 表示为

u^=mxx^+cx,v^=myy^+cy\hat{u}=m_x\hat{x}+c_x,\quad \hat{v}=m_y\hat{y}+c_y

代入径向畸变公式有

u^−cx=(u−cx)(1+kr2+k2r4+k3r6)\hat{u}-c_x=(u-c_x)(1+kr^2+k_2r^4+k_3r^6)

v^−cy=(v−cy)(1+kr2+k2r4+k3r6)\hat{v}-c_y=(v-c_y)(1+kr^2+k_2r^4+k_3r^6)

变形得到

[(u−cx)r2(u−cx)r4(u−cx)r6(v−cy)r2(v−cy)r4(v−cy)r6][k1k2k3]=[u^−uv^−v]\begin{bmatrix}(u-c_x)r^2&(u-c_x)r^4&(u-c_x)r^6\\(v-c_y)r^2&(v-c_y)r^4&(v-c_y)r^6\\\end{bmatrix}\begin{bmatrix}k_1\\k_2\\k_3\end{bmatrix}=\begin{bmatrix}\hat{u}-u\\\hat{v}-v\end{bmatrix}

使用最小二乘法即可求得 k1,2,3k_{1,2,3} 的近似拟合值。