⭐ Zeal's Blog
🛠 Zeal's 知乎专栏
🌀 Zeal's Github
本文内容
本文旨在较为直观地介绍相机成像背后的数学模型，主要的章节组织如下：
第一章用最简单的针孔投影模型为例讲解一个三维点是如何映射到图像中的一个像素
第二章介绍除了针孔投影模型外其他一些经典投影模型，旨在让读者建立不同投影模型之间的建模过程
第三章介绍如何把不同的投影模型用一个统一的投影过程表达
第四章进一步补充第三章的统一投影模型，并介绍畸变的定义和去畸变的原理
第五章针对全景相机的基本概念和两种应用广泛的全景相机模型做出介绍
第六章用代码示例介绍如何使用OpenCV的接口对图像去畸变
一、相机是如何成像的？
祭上一张经典图
针孔相机模型
一个三维空间中的点怎样映射到图片上的一个像素？
我们从一个最简单的问题开始：一个三维空间中的点\(P=[X,Y,Z]\)是如何经过相机成像变为图像上的一个像素\(\boldsymbol{p}=[u,v]\)的？
我们最常见的投影模型Perspective Projection Model 描述的就是针孔相机 的成像原理。从上面的图根据相似三角形可以得出
\[
\frac{X_c}{Z_c} = \frac{u - c_x}{f}\tag{1-1}
\]
其中\(c_x\)为光轴在图像中的\(x\)坐标，如果相机的光轴与感光元器件完全对齐的话，\(c_x = \frac{1}{2}W\), \(W\)是图像的宽度（单位是像素）
从上面的关系可以得出，将一个三维点投影到像素坐标系的时候，可以直接使用下面的公式。对应的就是针孔相机模型。下标\((\cdot)_c\) 表示这个点是在相机坐标系下的点
\[
\begin{split}
\lambda
\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}
\end{split}\tag{1-2}
\]
公式（1-2）在代码中也经常对\(u,v\)分开进行计算
\[
\begin{split}
u = f_x \frac{X_c}{Z_c} + c_x
\\
\\
v = f_y \frac{Y_c}{Z_c} + c_y
\end{split}\tag{1-3}
\]
当然，为了简化记号，公式（1-2）也通常记为矩阵相乘的形式
\[
\lambda\boldsymbol{p}=KP_c\tag{1-4}
\]
如果我们提前对\(P_c\)做归一化处理，也就是除以\(Z_c\)（假设点位于相机前，即\(Z_c>0\)），则可以去掉\(\lambda\)系数，即如下形式：
\[
\begin{split}
&\boldsymbol{p}_n=\frac{P_c}{Z_c}=
\begin{bmatrix}
\frac{X_c}{Z_c}&\frac{Y_c}{Z_c}&1
\end{bmatrix}^T
\\
&\boldsymbol{p}=K\boldsymbol{p}_n
\end{split}\tag{1-5}
\]
公式（1-4）中\(\boldsymbol{p}_n\)位于\(Z=1\)的平面，又叫归一化平面（nomalized plane），后文会再次讲到这个平面。
公式（1-2）到公式（1-5）其实都是等价的
如果给的三维点是在世界坐标系下，也就是\(P_w=[X_w,Y_w,Z_w]\)，那么我们只需要先把该点用相机的外参转换到相机坐标系下（刚性变换）即可：
\[
\begin{split}
\lambda\boldsymbol{p} = KT_{cw}P_w\\
\\
T_{cw} = [R|t]
\end{split}\tag{1-6}
\]
由于刚性变换过程不影响对相机投影模型的讨论，因此后面都假设三维点是处于相机坐标系
相机坐标系
在OpenCV中以及大部分视觉处理库中，相机坐标系的规定都与上述的图一致，就是相机光轴往前为\(Z\)，图像水平往右为\(X\)，图像垂直往下为\(Y\)。不过要注意的是在一些仿真渲染器或者特定任务的数据集中可能会规定图像垂直往上为\(Z\)，前为\(X\)，朝左为\(Y\)，但是这一点是无关紧要的，这一点差别可以反映在相机的外参里，也就是公式（1-6）中的\(T_{cw}\)，只要按照OpenCV的方式规定相机坐标系，总是可以找到一个外参矩阵\(T_{cw}\)将世界坐标系下的点变换到相机坐标系（前\(Z\)右\(X\)下\(Y\)）
大白话总结
问：相机是如何成像的？ 答：光束从物体表面反射，经过相机镜头，到达感光原件，这一系列物理过程可以通过数学公式表达，最终变成一个简单的矩阵操作将三维空间中的点对应到图片的一个像素。
二、不同的相机投影模型
第一节介绍的是针孔投影模型，但是事实上相机镜头都是多种多样的，不可能都是符合针孔投影模型。本节会介绍经典的相机投影模型，并从直观感受和形式化定义上介绍不同的投影模型是如何联系在一起的
什么是相机投影模型
相机投影模型用数学的方式描述了一个真实世界中的三维点到图像上像素坐标的映射关系
相机投影模型实际上就是对相机成像过程（物理）的数学建模。建模的目的是为了能够尽量符合真实的成像过程。不同的建模方式就对应不同的相机投影模型
经典的相机投影模型
我们回头看看公式（1-3），并暂时只关注\(X\)轴的映射关系
\[
\begin{split}
u = f_x \frac{X_c}{Z_c} + c_x
\end{split}\tag{1-3}
\]
上式中\(f_x\)称为相机焦距，反应了一个单位长度应该映射为几个像素，单位是\(pixels/m\)，\(c_x\)是相机坐标系（以光轴点为原点）到图像坐标系（以左上角为原点），这两个参数都是相机的固有参数。 而上式轴\(\frac{X_c}{Z_c}\)刚好是点\(P_c\)到相机光心连线与光轴角度的正切值，我们记光束与光轴的夹角为\(\theta\)，并将图片原点移动到图片中心，则公式（1-3）可以写为
\[
u^{\prime}=f_x\tan(\theta)\tag{2-1}
\]
上式的示意图如下（图中的\(r\)在只考虑\(X\)轴的时候就是\(u^{\prime}\)）：
光束与光轴夹角
在公式（2-1）中，\(X\)轴的投影坐标是\(\theta\)的函数，于是，我们是否可以用不同的函数表达这个过程？答案是肯定的！不同的函数就对应了不同的投影模型。下图就给出了在经典投影模型中对\(\theta\)的不同映射方式
经典投影模型
事实上，\(f\)作为相机的焦距，在上图中的不同投影模型都统一出现，于是我们可以舍弃焦距符号。于是，上图中不同的函数关系与投影模型的对应关系如下：
\(r(\theta)=\tan(\theta)\)，perspective projection/针孔投影模型/rectilinear model
\(r(\theta)=2\tan{\big(\frac{\theta}{2}\big)}\), stereographic projection
\(r(\theta)=\theta\), equidistant projection
\(r(\theta)=\sin{(\theta)}\), sine-law projection
\(r(\theta)=2\sin{\big(\frac{\theta}{2}\big)}\), equi-solid angle projection
上面两幅图出自:Models for the various classical lens projections，这篇文章比较形象地介绍各种经典的相机投影模型，并给出他们的函数曲线分析。不过总体偏形象化，没有引入更形式化的描述。
三、相机投影模型的统一表达形式
上一节我们将投影关系限制在\(X\)轴，并且给出了较为直观的图示。目的在于两个：1）给读者建立更深刻的相机投影过程;2）让读者对几种经典的投影模型有初步的直观了解。在这一小节中，我们给出更为统一的相机投影表达方式，同时为后文讨论相机的畸变 建立必要的理论基础。
我们将相机的投影过程拆分为三个过程：1）将空间中的点投影到单位球表面;2）将单位球上的点投影到归一化平面;3）将归一化平面上的点利用针孔模型投影到图像坐标系。下面详细介绍这几个过程
单位球投影
想象一下一束光束从相机光心射出，经过图像中的一个像素然后往外无限延伸，可以想象到，这个光束经过的任何点到图像的投影都是经过的那个像素。这个简单的事实告诉我们，我们可以对一个三维点进行任意的放缩，其在图像上的成像点都不会改变。于是，我们将三维点\(P_c\)除以它自身的膜，将其投影到一个单位球，其投影坐标为\(P_s\)
\[
\begin{split}
P_s = \frac{P_c}{||P_c||}=
\begin{bmatrix}
\frac{X_c}{||P_c||}&\frac{Y_c}{||P_c||}&\frac{Z_c}{||P_c||}
\end{bmatrix}^T
\end{split}\tag{3-1}
\]
示意图如下：
单位球投影示意图
在上图中，\(\theta\)是光束与光轴的夹角，\(\alpha\)为光束与水平轴的夹角。
两个角度有如下关系：
\[
\begin{split}
&\rho = \sqrt{(X_s^2+Y_s^2)}\\
&\theta=\operatorname{atan}(\frac{\rho}{Z_s})=\operatorname{asin}(\rho)\\
&\alpha=\operatorname{atan}(\frac{Y_s}{X_s})
\end{split}\tag{3-2}
\]
之所以构建\(\theta\)和\(\alpha\)是因为，我们后续可以将投影过程建立为这两个角度的函数，也就是只与光束的角度有关，而与具体的点坐标无关，这也是符合直观的。
另外，从前面的叙述以及常识，我们知道针孔成像结果是一个倒立的像，为了方便叙述，我们可以将相机做一个镜像，如下图所示：
对相机做镜像
这个单位球有时候又叫视球（Viewing Sphere）
归一化平面&&模型平面
将世界坐标系的点投影到单位球后，我们进一步将其映射到\(Z=1\)的平面上，这个平面又叫归一化平面（Normalized Plane）。 此时不同的投影模型会对在归一化平面上的点到原点的连线做放缩\(r(\theta)\)，为了后文叙述的统一性，我们再拆分出一个模型平面。在归一化平面上的点只是\(P_c\)与光心点的连线和平面的交点，即：
\[
p_n = \frac{P_c}{||Z_c||}\tag{3-3}
\]
在模型平面上（Model Plane），对归一化平面上的点做半径放缩，即
\[
\begin{split}
\boldsymbol{p}_d=
\begin{bmatrix}
r(\theta)\cos(\alpha)\\
r(\theta)\sin(\alpha)\\
1
\end{bmatrix}=
\begin{bmatrix}
u_d\\
v_d\\
1
\end{bmatrix}
\end{split}\tag{3-4}
\]
这两个平面的变换过程如下图所示
归一化平面投影&&模型平面
从第二节中经典的投影模型我们发现，其实不同的投影模型都没有对\(\alpha\)产生影响，而是对投影点到原点的距离\(r(\theta)\)建立不同的函数形式
稍后我们会看到，为了能够将所有模型都统一为针孔投影模型，会将模型平面变为畸变平面，用更具有一般性的多项式代替前面\(r(\theta)\)的表达形式。（因此，现在模型平面上的点用下标\(d\)标明，表示distortion）
透视变换
得到模型平面上的坐标\(\boldsymbol{p}_d=[u_d,v_d,1]^T\)后，我们可以用相机内参\(K\)将其变换到图像平面，这个步骤实际上就只是坐标的变换了。
\[
\begin{split}
u = f_xu_d+c_x\\
v = f_yv_d+c_y
\end{split}\tag{3-5}
\]
最后，我们可以得到一张完整地表达相机投影模型的示意图：
相机投影模型统一表达形式
将针孔投影带入上述统一模型
我们可以将针孔模型，即\(r(\theta)=\tan(\theta)\)代入上面的模型，联立公式（3-1）-（3-5）,最后可以发现，公式（3-5）的结果其实就是公式（1-3），也就是针孔投影模型的公式。这个过程比较简单，就不再展开公式。读者自己推导一下这个过程非常有利于理解上述的过程。
不同投影模型的函数图
我们可以将不同的经典投影模型的函数画出来（横轴为\(\theta\)，纵轴为\(r(\theta)\)），结果如下
经典投影模型的函数曲线
从上图我们应该至少要观察到一个重要的事实：针孔投影模型无法对\(>=90\degree\)的视野成像。因为\(\tan(\cdot)\)在90度的时候会趋于无穷大。实际上，从上面的函数图可以看到，针孔投影模型只能在大约水平140度以内的视野成像
相机投影模型总结
到目前位置，我们应该建立至少以下几个方面的认识：
统一化的投影模型经过：单位球投影->归一化平面->透视变换几个过程，将一个三维点投影到图像上的像素
不同的经典投影模型在投影过程中不会改变光束与\(X\)轴的夹角，只是对像素到图像原点的距离建立不同的方程
四、相机的畸变
针孔模型的优越性
首先描述一下我们人类直观上对于“标准的图像”这个词的一个感性认识，是不是我们会觉得横平竖直，真实中是直线则图像中也是直线，这样的图片会比较“标准”？事实上，针孔投影模型就刚好具有这样的性质。这个性质也可以从其投影方程看出来。经过针孔成像的物体，好像就是把整个物体缩小放在图片上，因此圆是圆，直线是直线。而其他投影模型就可能会呈现膨胀、紧缩的成像效果。如下图
标准的针孔成像与其他投影模型的图像
定义
“畸变”这个词从词语上应该理解成由于镜头加工等因素造成镜头与投影模型的差异。但实际上，相机畸变现在描述的是相机成像过程与针孔投影模型的差异（a deviation from the pinhole model） ，也就是：针孔投影+畸变模型=实际成像
而“去畸变”则是使用畸变模型对图像进行逆操作，使得图像就像用针孔投影模型成像出来的一样。
在归一化平面插入畸变模型
为了能够利用针孔模型的性质，我们在前文给出的统一相机投影模型中将模型平面用畸变平面代替。这里的核心思想是，用一个更具一般性的多项式替代各个投影模型中的模型函数，以此达到用一个方程表达多个投影模型的目的。
替换后的示意图如图所示：
加入畸变平面的投影模型
注意上图与上一节最后的统一投影模型其实是一样的，不过Model Plane名字换成Distortion Plane
我们这一节采用OpenCV实现的畸变模型来讲解畸变过程
我们还是先将\(P_c\)点投影到归一化平面得到\(\boldsymbol{p}_n=[u_n,v_n]\)，并令其到原点的半径为\(r=\sqrt{(u_n^2+v_n^2)}\)
径向畸变（Radial Distortion）
所谓的径向畸变（Radial Distortion）就是指只对点\(p_n\)做半径上的伸缩变化，这一点其实就跟经典投影模型一样用一个函数建立半径的变化，在OpenCV中，标准镜头的径向畸变可以用如下方程表达：
\[
\begin{split}
&u^{\prime}=u_n(1+k_1r^2+k_2r^4+k_3r^6)\\
&v^{\prime}=v_n(1+k_1r^2+k_2r^4+k_3r^6)
\end{split}\tag{4-1}
\]
径向畸变对于环绕光轴一周的改变是一致的，因此也叫做半径对称畸变（Radial Symmetric Distortion）或者畸变的对称部分（Symmetric Part Of Distortion Model）
切向畸变（Tangential Distortion）
真实的镜头由于加工误差和安装误差，会导致镜头与感光原件的中心不是完全对齐的，因此在平行的方向上会与标准的针孔成像模型有差异，这种差异对于光轴不是旋转对称的，也叫做切向畸变（Tangential Distortion）
OpenCV中镜头的切向畸变方程如下：
\[
\begin{split}
u_{d}=u^{\prime}+2p_1u_nv_n+p_2(r^2+2u_n)\\
v_{d}=v^{\prime}+p_1(r^2+2{v_n}^2) + 2p_2u_nv_n
\end{split}\tag{4-2}
\]
公式（4-1）和公式（4-2）中的\(k_1,k_2,k_3,p_1,p_2\)叫做畸变参数（Distortion Coefficients）
去畸变
有了畸变模型之后，我们就可以将一个三维点经过：单位球投影->归一化平面投影->畸变模型->透视变换，得到该点在图像中的像素位置。我们可以通过一些标定物，经过求解得到上述的畸变系数，进而得到畸变过程的反过程。我们将图像上的点先经过透视变换的反变换得到在畸变平面上的点，再经过畸变过程的反过程，再投回图像中，这样我们就得到一副没有畸变的图像，其看起来就像是用完全标准的针孔投影模型成像的照片。
下面的图片展示了从畸变到去畸变的图像变化
去畸变图像变化
做完去畸变后，整个的成像过程就像直接用针孔投影模型（公式1-3）成像一般，就像下图一样。
去畸变的等效投影过程
提出畸变模型的文章需要给出如何标定出畸变参数，同时如何计算从畸变点到非畸变点，这里面主要是一些数学求解，这里就不展开讨论，有兴趣的可以看看wiki:Distortion
五、全景相机（Omidirectional Camera）
所谓的全景相机广泛概念上指成像角度能够大于等于180度的相机，他们看起来大概像是下面这样
全景相机类型
从第二节（不同的相机投影模型）我们可以知道，根据针孔投影模型设计出的镜头无法对大于等于90度（左右）的视野成像，通常由于进光量等问题，这类相机一般只有140度左右的成像视野。
一般的畸变模型的设计以针孔相机模型为基础，通过参数模型修正真实的成像与针孔成像的差异，因为在很多的应用中我们希望能够通过“去畸变”的方式，将图片变成“直线还是直线（Straight lines are straight），在KannalaBrandt论文中是这样描述的：
It is impossible to project the hemispherical field of view on a finite image plane by a perspective projection so fish-eye lenses are designed to obey some other projection model.This is the reason why the inherent distortion of a fish-eye lens should not be considered only as a deviation from the pinhole model
于是，就有很多专门针对全景相机（Omidirection Camera/Fisheye Camera/Wide-Angle camera）的建模研究出现。
重新理清一下我们的目的：
拥有一种统一的表达方式能够尽量拟合真实的全景相机的成像过程
这种表达方式应该简洁有效
能够对模型的参数求解，并将图像通过“去畸变”变成像是由针孔相机拍摄出来的横平竖直的图像
OpenCV中针对全景相机的标定和去畸变给出了两种实现
KannalaBrandt模型.对应实现OpenCV::Fisheye
CMei模型. 对应实现OpenCV::Omnidir
下面简单介绍这两种模型
KannalaBrandt模型
第一步还是先将点投影到单位球模型，这样我们就得出了\(\theta,\alpha\)两个角度，后续的畸变模型就是关于这两个角度的函数
这里这里摘录OpenCV的描述方式OpenCV_Fisheye，相比于论文，在模型参数上简化了很多
KannalaBrandt模型使用一个多项式描述径向畸变（畸变的对称部分）
\[
\begin{split}
&\boldsymbol{p}_n=\frac{P_c}{P_z}=[u_n,v_n,1]^T\\
\\
&r = \sqrt{u_n^2+v_n^2}\\
\\
&\theta_d = k_1\theta + k_2\theta^3 + k_3\theta^5+k_4\theta^7+k_5\theta^9+....
\end{split}\tag{5-1}
\]
径向畸变后坐标变为：
\[
\begin{split}
u^{\prime} = \frac{\theta_d}{r}u_n\\
v^{\prime} = \frac{\theta_d}{r}v_n
\end{split}\tag{5-2}
\]
最后再进行非对称畸变
\[
\begin{split}
&u_d=u^{\prime}+\alpha v^{\prime}\\
&v_d=v^{\prime}
\end{split}\tag{5-3}
\]
最后再经过公式（3-5）变换到图像坐标系
CMei模型
CMei模型相比于其他模型有一个较大不同之处在于在从单位球投影到归一化平面时将相机光心往后移动了距离\(\xi\)，总体的投影过程如下图：
CMei投影模型
上图出自引用CMei的一篇论文：Design and Calibration of an Omni-RGB+D Camera
六、代码实例
实例一：针孔投影模型去畸变
123456789101112# undistort imageh, w = img.shape[:2]newcameramtx, roi = cv2.getOptimalNewCameraMatrix(K, distortions, (w,h), alpha, (w,h))undistorted_img = cv2.undistort(img, K, distortions, None, newcameramtx)x, y, w, h = roiundistorted_img = undistorted_img[y:y+h, x:x+w]# undistort image pointsif points2d.ndim == 2: points2d = points2d[:, None, :]undistorted_points = cv2.undistortPoints(points2d, K, distortions, P=K)undistorted_points = undistorted_points.reshape(-1, 2)
实例二：CMei模型去畸变
123456789101112131415161718192021222324252627# undistort imageif newK: K_new = np.zeros((3,3), np.float64) h,w = img.shape[:2] K_new[0, 0] = w/4 K_new[0, 2] = w/2 K_new[1, 1] = h/4 K_new[1, 2] = h/2 K_new[2, 2] = 1.0else: K_new = Noneundistorted_img = cv2.omnidir.undistortImage(img, K, distortions, Xi, cv2.omnidir.RECTIFY_PERSPECTIVE, Knew=K_new)# undistort image pointsif points2d.ndim == 2: points2d = points2d[:, None, :]undistorted_points = cv2.omnidir.undistortPoints(points2d, K, distortions, Xi, None)undistorted_points = undistorted_points.squeeze()f0 = K_new[0,0]f1 = K_new[1,1]c0 = K_new[0,2]c1 = K_new[1,2]undistorted_points[:, 0] = f0* undistorted_points[:, 0] + c0undistorted_points[:, 1] = f1 * undistorted_points[:, 1] + c1
PS:cv2::omnidir::undisortPoints的旧版本有bug，参考这个issue。omnidir空间目前还没有成为opencv的正式稳定接口，因此维护在opencv-contrib-python包中，最新的包（4.6.x）已经修复了bug。一定要检查一下是不是有那个bug！！！！
实例三：KannalaBrandt模型去畸变
123456789101112131415# You should replace these 3 lines with the output in calibration stepDIM=XXXK=np.array(YYY)D=np.array(ZZZ)def undistort(img_path): img = cv2.imread(img_path) h,w = img.shape[:2] map1, map2 = cv2.fisheye.initUndistortRectifyMap(K, D, np.eye(3), K, DIM, cv2.CV_16SC2) undistorted_img = cv2.remap(img, map1, map2, interpolation=cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT) cv2.imshow("undistorted", undistorted_img) cv2.waitKey(0) cv2.destroyAllWindows()if __name__ == '__main__': for p in sys.argv[1:]: undistort(p)
相机标定方法
最经典的是使用张正友标定法。原文;中文讲解
参考资料
OpenCV:Camera Calibration
OpenCV2.4:Camera Calibration
Models for the various classical lens projections
Wiki:Distortion
KannalaBrandt:A Generic Camera Model and Calibration Method for Conventional, Wide-Angle, and Fish-Eye Lenses
CMei: Single View Point Omnidirectional Camera Calibration from Planar Grids
OpenCV::Fisheye
OpenCV::Omnidir