【鱼眼相机模型】鱼眼相机投影模型理解

一、从普通镜头到鱼眼镜头

如图1所示，普通镜头下的光线依据针孔相机模型进行成像（该部分可参考相机投影关系）。但该模型存在一个缺陷：相机视野范围越大，所需的成像平面也越大，当相机视野范围要求在180°时，所需的成像平面要求为无限大。

图1.针孔相机模型

在一些需要大角度视野的场景下，为解决相机视野需求和成像平面之间的矛盾，人们通过将一系列透镜进行组合，使得光线出射角小于入射角度，将大角度视野中的空间投影到有限尺寸的成像平面上。

图2.简化版本的针孔模型与鱼眼模型光路图对比

二、鱼眼相机模型详述

图3.鱼眼相机

鱼眼相机的一般结构： 如图3.a所示，鱼眼相机一般由若干不同的透镜组合而成，在成像过程中，入射光线经过不同程度的折射，投影到尺度有限的成像平面上。
鱼眼相机模型： 由于鱼眼相机的多元件结构使对鱼眼相机的折射关系分析变得相当复杂，如图3.b所示，在文¹中提出单位球面投影模型，用以简化该折射关系，该模型将鱼眼相机的成像过程分解为两步：

1. 将三维空间点线性投影到球心与相机坐标系原点重合的虚拟单位球面上
2. 将单位球面的点投影到图像平面上，该过程是非线性的。根据投影函数的不同，可将投影模型进一步划分为以下表，投影模型的图示见图4。

投影模型	投影函数	特征
i. 体视投影 (stereographic projection)	$r=2ftan\frac{\theta }{2}$	任何直线相交的角度，在交换后保持不变
ii.等距投影 (equidistance projection)	$r=f\theta$	物体成像面上距离画面中心的距离与入射角成正比
iii. 等立体角投影 (equisolid angle projection)	$r=2fsin\frac{\theta }{2}$	在变换前后，物体所占的立体角大小不变
iv. 正交投影 (orthogonal projection)	$r=fsin\theta$	投影畸变最大，而且最大视场角不能大于180°

图4.投影模型

图3.b与图4的变量释义：

• $X_c\_Y_c\_Z_c$:相机坐标系
• $x\_y$:图像坐标系
• P:物点，其在相机坐标系下的坐标$P(X_c,Y_c,Z_c)$
• $\theta :$入射角
• $p$:鱼眼模型的像点$p$(x,y),即畸变的像点
• $p^{{}^{\prime }}$:针孔模型的像点$p^{{}^{\prime }}$$(x^{{}^{\prime }},y^{{}^{\prime }})$，即未畸变的像点，（针孔模型的入射角等于出射角）
• $r$:畸变像点p在极坐标系下的弧长
• $r^{{}^{\prime }}$:非畸变像点$p^{{}^{\prime }}$在极坐标系下的弧长

Kannala-Brandt 模型：
为了利于标定，文¹提出Kannala-Brandt 模型，将上述四个投影模型的入射角通过泰勒展开，取前5项，表示为如下形式，标定时即标定$k_1,k_2,k_3,k_4$参数：
$$2tan\frac{\theta }{2}=\theta +k_1{\theta }^3+k_2{\theta }^5+k_3{\theta }^7+k_4{\theta }^9={\theta }_d$$
$$2sin\frac{\theta }{2}=\theta +k_1{\theta }^3+k_2{\theta }^5+k_3{\theta }^7+k_4{\theta }^9={\theta }_d$$
$$sin\theta =\theta +k_1{\theta }^3+k_2{\theta }^5+k_3{\theta }^7+k_4{\theta }^9={\theta }_d$$

则上述四个投影函数可统一用下式表示：
$$r=f{\theta }_d$$

A:物点P的成像过程：

1. 世界坐标系到相机坐标系的坐标转换：$P=R{P}_w+T$，其中$P_w$为世界坐标系坐标，$P(X_c,Y_c,Z_c)$为相机坐标系
2. 利用针孔模型，求解非畸变点图像坐标系下坐标：$x_{p^{{}^{\prime }}}=\frac{X_c}{Z_c}f$，$y_{p^{{}^{\prime }}}=\frac{Y_c}{Z_c}f$
3. 求解非畸变像点的弧长：$r^{{}^{\prime }2}=x_{p^{{}^{\prime }}}^2+y_{p^{{}^{\prime }}}^2$
4. 求解入射角$\theta :$：由于非畸变像点的出射角等于入射角，$\theta =atan(r^{{}^{\prime }}/f)$;
5. 求解${\theta }_d$：${\theta }_d=\theta +k_1{\theta }^3+k_2{\theta }^5+k_3{\theta }^7+k_4{\theta }^9$
6. 求解畸变点弧长:$r=f{\theta }_d$
7. 求解畸变像点图像坐标系坐标：利用相似三角形$x_p=\frac{r}{r^{{}^{\prime }}}{x}_{p^{{}^{\prime }}}$，$y_p=\frac{r}{r^{{}^{\prime }}}{y}_{p^{{}^{\prime }}}$
8. 利用相机内参将畸变像点：从图像坐标系转至像素坐标系： $u=\frac{x_p}{{\Delta }_x}+c_x$,$v=\frac{y_p}{{\Delta }_y}+c_y$;$c_x、c_y$是相机内参，${\Delta }_x、{\Delta }_y$为像元大小（图像像素大小）

B:上述投影在opencv中简化如下（文末提供证明过程）：

1. 世界坐标系到相机坐标系的坐标转换：$P=R{P}_w+T$，其中$P_w$为世界坐标系坐标，$P(X_c,Y_c,Z_c)$为相机坐标系
2. 利用针孔模型，非畸变点图像坐标系下坐标：${x_{p^{{}^{\prime }}}}_{new}=\frac{X_c}{Z_c}$，${y_{p^{{}^{\prime }}}}_{new}=\frac{Y_c}{Z_c}$
3. 求解非畸变像点的弧长：${r^{{}^{\prime }2}}_{new}={x_{p^{{}^{\prime }}}}_{new}^2+{y_{p^{{}^{\prime }}}}_{new}^2$
4. 求解入射角$\theta :$由于非畸变像点的出射角等于入射角，${\theta }_{new}=atan({r^{{}^{\prime }}}_{new})$ (ps:结合B.2、B3、B4，A.2，A.3、A4可知${\theta }_{new}=\theta$)
5. 求解畸变点弧长：$r_{new}={\theta }_{new}+k_1{\theta }_{new}^3+k_2{\theta }_{new}^5+k_3{\theta }_{new}^7+k_4{\theta }_{new}^9$
6. 求解畸变像点图像坐标系坐标：${x_p}_{new}=\frac{r_{new}}{r_{new}^{{}^{\prime }}}{x}_{p_{new}^{{}^{\prime }}}$，$y_p=\frac{r_{new}}{r_{new}^{{}^{\prime }}}{y}_{p_{new}^{{}^{\prime }}}$
7. 利用相机内参将畸变像点：从图像坐标系转至像素坐标系$u=f_x\ast {x_p}_{new}+c_x$,$v=f_y\ast {y_p}_{new}+c_y$

（11月16日）再补充：上述B中前面的物理含义注释只是为了和A中做一个对应，读者不必在意其代表的真实物理含义，opencv做这个简化的目的是因为真实使用时我们得到的相机内参为fx整体，无法独立获得焦距f, 而通过B的过程，可以不依赖f，仅依赖fx计算得到三维点的二维投影像素坐标。读者在理解的时候，着重体会A.7 = B.7的最终结果即可！！！！！！！！！ 换句话说B可以看做以下过程（但直接像下面这么写，容易让人看得云里雾里，故原文加了前面的对应注释）：

1. 世界坐标系到相机坐标系的坐标转换：$P=R{P}_w+T$，其中$P_w$为世界坐标系坐标，$P(X_c,Y_c,Z_c)$为相机坐标系
2. 计算新变量：${x_{p^{{}^{\prime }}}}_{new}=\frac{X_c}{Z_c}$，${y_{p^{{}^{\prime }}}}_{new}=\frac{Y_c}{Z_c}$
3. 计算新变量：${r^{{}^{\prime }2}}_{new}={x_{p^{{}^{\prime }}}}_{new}^2+{y_{p^{{}^{\prime }}}}_{new}^2$
4. 计算新变量，${\theta }_{new}=atan({r^{{}^{\prime }}}_{new})$ (ps:结合B.2、B3、B4，A.2，A.3、A4可知${\theta }_{new}=\theta$)
5. 计算新变量：$r_{new}={\theta }_{new}+k_1{\theta }_{new}^3+k_2{\theta }_{new}^5+k_3{\theta }_{new}^7+k_4{\theta }_{new}^9$
6. 计算新变量：${x_p}_{new}=\frac{r_{new}}{r_{new}^{{}^{\prime }}}{x}_{p_{new}^{{}^{\prime }}}$，$y_p=\frac{r_{new}}{r_{new}^{{}^{\prime }}}{y}_{p_{new}^{{}^{\prime }}}$
7. 得到像素坐标：$u=f_x\ast {x_p}_{new}+c_x$,$v=f_y\ast {y_p}_{new}+c_y$

等价性证明：通过简单的等式代入可得A中 $\frac{r}{r^{{}^{\prime }}}$等于B中$\frac{r_{new}}{r_{new}^{{}^{\prime }}}$，而A.8展开如下（@JasonGao1991，A.8式确实无法简单抵消掉f,现提供完整过程）

$$u=\frac{x_p}{{\Delta }_x}+c_x=\frac{f}{{\Delta }_x}\ast \frac{r}{r^{{}^{\prime }}}\ast \frac{X_c}{Z_c}+c_x==\frac{f}{{\Delta }_x}\ast \frac{r_{new}}{r_{new}^{{}^{\prime }}}\ast \frac{X_c}{Z_c}+c_x=f_x\ast {x_p}_{new}+c_x=B.7$$
$$v=\frac{y_p}{{\Delta }_y}+c_y=\frac{f}{{\Delta }_y}\ast \frac{r}{r^{{}^{\prime }}}\ast \frac{Y_c}{Z_c}+c_x=\frac{f}{{\Delta }_y}\ast \frac{r_{new}}{r_{new}^{{}^{\prime }}}\ast \frac{Y_c}{Z_c}+c_x=f_y\ast {y_p}_{new}+c_y=B.7$$

opencv原文：

    Eigen::Vector2d KannalaBrandt8::project(const Eigen::Vector3d &v3D) {
        const double x2_plus_y2 = v3D[0] * v3D[0] + v3D[1] * v3D[1];
        const double theta = atan2f(sqrtf(x2_plus_y2), v3D[2]);
        const double psi = atan2f(v3D[1], v3D[0]);

        const double theta2 = theta * theta;
        const double theta3 = theta * theta2;
        const double theta5 = theta3 * theta2;
        const double theta7 = theta5 * theta2;
        const double theta9 = theta7 * theta2;
        const double r = theta + mvParameters[4] * theta3 + mvParameters[5] * theta5
                        + mvParameters[6] * theta7 + mvParameters[7] * theta9;

        Eigen::Vector2d res;
        res[0] = mvParameters[0] * r * cos(psi) + mvParameters[2];
        res[1] = mvParameters[1] * r * sin(psi) + mvParameters[3];

        return res;

    }

文章部分参考²：

---

1. A Generic Camera Model and Calibration Method for Conventional, Wide-Angle, and Fish-Eye Lenses ↩︎ ↩︎

2. 鱼眼镜头的成像原理到畸变矫正 ↩︎