鱼眼镜头的成像原理到畸变矫正（完整版）

  最近刚接触鱼眼相机，发现网上资料还是比较零散的，于是把搜罗到的资料汇总梳理了一下。这里推荐大家直接看链接6的论文，从成像模型到畸变矫正整个过程讲的比较清楚，网上很多版本其实都是根据这两篇论文来的，而且绝大部分博客只写到了成像畸变过程，至于如何进行畸变矫正完全没有提。

1. 摄影入门之相机镜头的分类
2. 鱼眼镜头是怎么「鱼眼」的？
3. 鱼眼相机成像模型
4. 鱼眼相机成像模型学习
5. 2004 - A generic camera calibration method for fish-eye lenses
6. 2006 - A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses
7. OpenCV 4.5.1 – Fisheye camera model
8. 鱼眼摄像头标定与畸变校正（双OPENCV版本）
9. {畸变矫正}图像去畸变opencv各种方法与matlab各种方法-联合分析（一统江湖）
10. 【数值分析】一元方程的迭代解法
11. 数值分析笔记整理（1）——求解方程的方法

1. 什么是鱼眼镜头

1.1 镜头分类

  相机镜头大致上可以分为变焦镜头和定焦镜头两种。顾名思义，变焦镜头可以在一定范围内变换焦距，随之得到不同大小的视野；而定焦镜头只有一个固定的焦距，视野大小是固定的。鱼眼镜头是定焦镜头中的一种视野范围很大的镜头，视角通常大于180°。 如下图所示，在获取更大视野范围的同时，鱼眼镜头成像的畸变也更大。

1.2 鱼眼镜头和普通镜头的区别

  那么鱼眼镜头为什么可以获得比普通镜头更大的视野范围呢，他们两者之间的差别在哪？
  其实，我们平常接触的大多数镜头都可以近似看做针孔相机模型，该模型下，光线沿直线传播，像与物之间是相似的，或者更严格地用数学语言来说，像与物之间是经过了透视变换（Perspective Transform）。在透视变换下，直线经过变换仍是直线，曲线经过变换仍是曲线，两直线交点经过变换仍是两直线相交的点等。正因为投影变换保持了很多几何性质不变，所以我们看照片是能够与原场景联系起来的，照片与原场景之间存在某些相似的特性。

  从某种意义上来说，相机镜头所起的作用，就是做了一个数学变换，将物空间变换为像空间，成像平面就是在像空间内切了一刀，截取了一个平面，成为拍下的照片。
  但基于针孔相机模型的镜头存在一个缺陷——光线始终沿直线传播使得镜头难以捕捉位于边缘的物体。如下图所示，对于同样长度的红色箭头，越靠近边缘的经过镜头成像后就变得越长，而实际上我们底片的尺寸是有限的，所以极端接近边缘的物体普通的镜头就无法成像记录了。

  于是人们想到的水下的鱼。由于水的折射率比空气大，光线从空气进入水中，折射角比入射角更小，并且入射角越大，这个变小的程度也越大。由于这个特性，使得在水中向上看时，能一眼看到整个水面上的这个半球形空间，整个空间的影像都背扭曲、压缩到了一个半顶角约为48°的锥形内。

在这个锥形空间内部，是来自水面上的空间的光线，在这个锥形外部，是来自水面下景色的反射。也就是说，在水下向上看，在一个圈之外，只能看到水底的景色；所有水面上的景色，都被压缩在一个圈内，如下图所示。鱼眼镜头也是人们根据这种特性发明的，另外，鱼眼镜头的前镜片直径很短，且呈抛物状像前部凸出，与鱼的眼睛十分相似，“鱼眼镜头”因此而得名。

2. 鱼眼镜头的投影模型

  鱼眼镜头一般是由十几个不同的透镜组合而成的，如下图所示，在成像的过程中，入射光线经过不同程度的折射，投影到尺寸有限的成像平面上，使得鱼眼镜头与普通镜头相比起来拥有了更大的视野范围。

  在研究鱼眼相机成像时，可以将上面的镜头组简化为一个球面，如下图 (b) 所示， O 1 − X c Y c Z c O_1-X_cY_cZ_c O1​−Xc​Yc​Zc​是相机坐标系， O 2 − x y O_2-xy O2​−xy是成像平面。现实世界有一点 P P P，入射角为 θ \theta θ，如果按照普通相机的针孔相机模型，入射光线 P O 1 PO_1 PO1​经过镜头后不改变路线， P 、 O 1 、 p ′ P、O_1、p' P、O1​、p′三点共线， p ′ p' p′为 P P P的像点；但对于鱼眼相机，入射光线 P O 1 PO_1 PO1​经过镜头后会发生折射，因此 P P P的像点为 p p p点，极坐标表示为 ( r , φ ) (r, \varphi) (r,φ)。

  为了将尽可能大的场景投影到有限的图像平面内，鱼眼相机会按照一定的投影函数来设计，如上图 (a) 所示。根据投影函数的不同，鱼眼相机的设计模型大致能被分为五种：透视投影（即针孔相机模型）、等积投影、等距投影、体视投影、正交投影。

3. 鱼眼镜头的成像过程

  实际的镜头因为各种原因并不会精确的符合投影模型，为了方便鱼眼相机的标定，一般取 r r r关于 θ \theta θ泰勒展开式的前5项来近似鱼眼镜头的实际投影函数： r ( θ ) ≈ k 0 θ + k 1 θ 3 + k 2 θ 5 + k 3 θ 7 + k 4 θ 9 (1) r(\theta)\approx k_0\theta+k_1\theta^3+k_2\theta^5+k_3\theta^7+k_4\theta^9\tag{1} r(θ)≈k0​θ+k1​θ3+k2​θ5+k3​θ7+k4​θ9(1)

  假设相机坐标系下有一点 P ( x , y , z ) P(x,y,z) P(x,y,z)，点 P ( x , y , z ) P(x,y,z) P(x,y,z)如果按照针孔相机模型投影，则不存在畸变，像点为 P 0 ( a , b ) P_0(a,b) P0​(a,b)。不妨假设 f = 1 f=1 f=1（最终可以求得 r d r_d rd​和 r r r的比值与 f f f无关）,可求得 P 0 P_0 P0​点坐标y以及入射角 θ \theta θ： { a = x / z b = y / z r 2 = a 2 + b 2 θ = a t a n ( r ) (2) \left \{ \begin{aligned}a&=x/z \\ b&=y/z \\ r^2&=a^2+b^2 \\ \theta&=atan(r)\end{aligned}\right.\tag{2} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​abr2θ​=x/z=y/z=a2+b2=atan(r)​(2)
  由于畸变的存在，像点到图像中心的距离 r r r被压缩成 r d r_d rd​，实际的像点位置为 p ′ ( x ′ , y ′ ) p'(x',y') p′(x′,y′)，有 ∣ O p ′ ∣ = r d , ∣ O P 0 ∣ = r |Op'|=r_d,|OP_0|=r ∣Op′∣=rd​,∣OP0​∣=r。结合等距投影函数和式（1）有： r d = f θ ≈ k 0 θ + k 1 θ 3 + k 2 θ 5 + k 3 θ 7 + k 4 θ 9 (3) r_d=f\theta \approx k_0\theta+k_1\theta^3+k_2\theta^5+k_3\theta^7+k_4\theta^9\tag{3} rd​=fθ≈k0​θ+k1​θ3+k2​θ5+k3​θ7+k4​θ9(3)因为 f = 1 f=1 f=1，且 θ d \theta_d θd​的一次项系数 k 0 k_0 k0​可以为1，最终可以得到OpenCV中使用的鱼眼相机模型：
r d = f ⋅ t a n θ d = θ d = θ ( 1 + k 1 θ 2 + k 2 θ 4 + k 3 θ 6 + k 4 θ 8 ) (4) r_d=f \cdot tan \theta_d = \theta_d=\theta(1+k_1\theta^2+k_2\theta^4+k_3\theta^6+k_4\theta^8)\tag{4} rd​=f⋅tanθd​=θd​=θ(1+k1​θ2+k2​θ4+k3​θ6+k4​θ8)(4)由相似三角形原理： r d r = θ d r = x ′ a = y ′ b (5) \frac{r_d}{r}=\frac{\theta_d}{r}=\frac{x'}{a}=\frac{y'}{b}\tag{5} rrd​​=rθd​​=ax′​=by′​(5)所以就求得了畸变后点 p ′ p' p′的坐标为： { x ′ = θ d r a y ′ = θ d r b (6) \left \{ \begin{aligned}x'&=\frac{\theta_d}{r}a \\ y'&=\frac{\theta_d}{r}b \end{aligned}\right.\tag{6} ⎩⎪⎪⎨⎪⎪⎧​x′y′​=rθd​​a=rθd​​b​(6)最后利用相机内参将像平面的点转换到像素坐标系就得到了最终图像上的点： { u = f x x ′ + c x v = f y y ′ + c y (7) \left \{ \begin{aligned}u&=f_xx'+c_x \\ v&=f_yy'+cy\end{aligned}\right.\tag{7} {uv​=fx​x′+cx​=fy​y′+cy​(7)

对于公式(3)的解释

  首先感谢评论区几位认真的小伙伴指出这个问题，这里我确实表示的有错误，导致不少人对公式3有误解，这里单独再拎出来解释一下。
  这里对比下面两张图一起看，左边的是论文的截图，右边是我按照自己理解画的，两者的符号表示会不太一样，不要混淆在一起。

  论文指出（左图），根据投影函数的不同，鱼眼相机的设计模型大致能被分为五种：透视投影、等积投影、等距投影、体视投影、正交投影，最常用的是等距投影模型，即 r = f θ (eq1) r=f\theta\tag{eq1} r=fθ(eq1)，但投影模型这么多是不利于自动标定的，于是将上述几种投影模型统一用 r r r关于 θ \theta θ的泰勒展开式来近似表示，这么做会带来一些精度损失，但基本可以忽略： r ( θ ) = f θ ≈ k 0 θ + k 1 θ 3 + k 2 θ 5 + k 3 θ 7 + k 4 θ 9 (eq2) r(\theta) = f\theta \approx k_0\theta+k_1\theta^3+k_2\theta^5+k_3\theta^7+k_4\theta^9\tag{eq2} r(θ)=fθ≈k0​θ+k1​θ3+k2​θ5+k3​θ7+k4​θ9(eq2)这里 θ \theta θ是入射光线 ∣ P O 1 ∣ |PO_1| ∣PO1​∣和光轴的夹角，即入射角， r r r表示相机空间任意点 P P P在相机成像平面的像点 p p p距离光心的距离 ∣ O 2 p ∣ |O_2p| ∣O2​p∣。根据 e q 3 eq3 eq3我们能够知道，成像点到光心的距离 r r r是关于入射角 θ \theta θ的函数，但光线入射后以什么角度射出我们是难以计算的，因为鱼眼镜头是由一组透镜组成的，光线入射后的光路非常复杂，会在不同透镜间反复折射，可以看到论文的截图中入射光线经过 O 1 O_1 O1​射到点 p p p的光路画的也是一条曲线。

  然后再来看看我自己画的图（右图），畸变前的像点 P 0 ( a , b ) P_0(a,b) P0​(a,b)和畸变后的像点 p ′ ( x ′ , y ′ ) p'(x',y') p′(x′,y′)到光心 O O O的距离分别为 r r r和 r d r_d rd​。实际中我们一般知道相机空间的点 P P P坐标，如果不考虑畸变，知道相机焦距 f f f的情况下，根据出射角等于入射角，畸变点的像点 P 0 P_0 P0​也很容易可以得到，所以 r r r的值也知道，根据相似三角形原理： r d r = x ′ a = y ′ b (eq3) \frac{r_d}{r}=\frac{x'}{a}=\frac{y'}{b}\tag{eq3} rrd​​=ax′​=by′​(eq3)从而可以求得，一般记 s c a l e = r d / r scale=r_d/r scale=rd​/r: { x ′ = r d r a = s c a l e ⋅ a y ′ = r d r b = s c a l e ⋅ b (eq4) \left \{ \begin{aligned}x'&=\frac{r_d}{r}a=scale\cdot a \\ y'&=\frac{r_d}{r}b=scale\cdot b \end{aligned}\right.\tag{eq4} ⎩⎪⎨⎪⎧​x′y′​=rrd​​a=scale⋅a=rrd​​b=scale⋅b​(eq4)取 k 0 = 1 k_0=1 k0​=1，将 r d r_d rd​代入 e q 3 eq3 eq3有： r d ≈ θ + k 1 θ 3 + k 2 θ 5 + k 3 θ 7 + k 4 θ 9 (eq5) r_d \approx \theta+k_1\theta^3+k_2\theta^5+k_3\theta^7+k_4\theta^9\tag{eq5} rd​≈θ+k1​θ3+k2​θ5+k3​θ7+k4​θ9(eq5)到这一步为止应该没有什么歧义，至于 k 0 k_0 k0​为什么取1我暂时没找到比较官方的解释，个人感觉是为了少标定一个参数，实际上一次项系数对于整个表达式的值影响也没那么大， k 0 = 1 k_0=1 k0​=1带来的精度损失靠后面高次项系数补回来就好了。然后重点来了，这也是我之前理解错的地方，如果不是评论区的小伙伴指出我就一直想当然地错下去了。
  翻阅过opencv文档的小伙伴们会发现opencv文档中是这么写的： θ d = θ ( 1 + k 1 θ 2 + k 2 θ 4 + k 3 θ 6 + k 4 θ 8 ) (eq6) \theta_d = \theta(1+k_1\theta^2+k_2\theta^4+k_3\theta^6+k_4\theta^8)\tag{eq6} θd​=θ(1+k1​θ2+k2​θ4+k3​θ6+k4​θ8)(eq6)如果和上面我们推导出来的 e q 5 eq5 eq5比较，会得到 r d = θ d (eq7) r_d=\theta_d\tag{eq7} rd​=θd​(eq7)的结论，之前我的解释是错误的，正确理解应该是：

   θ \theta θ是入射角， θ d \theta_d θd​畸变后的等效折射角（不是实际的折射角）， r d r_d rd​和 θ d \theta_d θd​满足 r d = f ⋅ t a n θ d r_d=f \cdot tan\theta_d rd​=f⋅tanθd​，如果取等效焦距 f = 1 f=1 f=1（由eq4可以得到，畸变和相机焦距无关，仅和 s c a l e scale scale值有关，所以 f f f取多少都无所谓， f f f变了， r d r_d rd​相应的也会跟着变化），就有 r d = t a n θ d r_d= tan\theta_d rd​=tanθd​，考虑到相机的成像CCD平面尺寸一般都是几毫米，焦距在几百毫米左右，所以相机实际成像过程中 θ d \theta_d θd​是比较小的， t a n θ d tan \theta_d tanθd​可以近似用 θ d \theta_d θd​表示( x → 0 , t a n x = x x\rightarrow0, tanx=x x→0,tanx=x)，所以才有了这么一个等式 r d = θ d = θ ( 1 + k 1 θ 2 + k 2 θ 4 + k 3 θ 6 + k 4 θ 8 ) (eq7) r_d=\theta_d=\theta(1+k_1\theta^2+k_2\theta^4+k_3\theta^6+k_4\theta^8)\tag{eq7} rd​=θd​=θ(1+k1​θ2+k2​θ4+k3​θ6+k4​θ8)(eq7)  所以不管鱼眼镜头使用的什么模型，推导过程都是一样的，至于 t a n θ d tan\theta_d tanθd​是否要近似成 θ d \theta_d θd​，我个人觉得只要畸变参数的标定过程和解畸变过程保持一致就不会有太大问题。

  另外，本文讨论的只有径向畸变，如果考虑切向畸变的话 e q 6 eq6 eq6会有所不同。

4. 鱼眼镜头的畸变矫正

  谈畸变矫正之前先简单回顾一下畸变成像过程。
  相机坐标系存在一点 P ( x , y , z ) P(x,y,z) P(x,y,z)，现在要获得该点在鱼眼相机像平面的投影，需要经过如下步骤：

Step1：根据针孔相机模型成像原理，可以求得未发生畸变时，点 P P P的像点 P 0 ( a , b ) P_0(a,b) P0​(a,b)，极坐标形式表示为 ( r , φ ) (r,\varphi) (r,φ)，以及点 P P P的投影入射角 θ \theta θ
Step2：实际上由于畸变的存在，光线出射角 θ d ≠ θ \theta_d \neq \theta θd​​=θ，实际的像点为 p ′ ( x ′ , y ′ ) p'(x',y') p′(x′,y′)
Step3：利用等距投影公式以及泰勒展开式，可以近似求得 θ d = θ ( 1 + k 1 θ 2 + k 2 θ 4 + k 3 θ 6 + k 4 θ 8 ) \theta_d=\theta(1+k_1\theta^2+k_2\theta^4+k_3\theta^6+k_4\theta^8) θd​=θ(1+k1​θ2+k2​θ4+k3​θ6+k4​θ8)
Step4：因为 r d = θ d r_d=\theta_d rd​=θd​，所以点 p ′ p' p′的极坐标为 ( θ d , φ ) (\theta_d,\varphi) (θd​,φ)，从而求得笛卡尔坐标值 x ′ = ( θ d / r ) a , y ′ = ( θ d / r ) b x'=(\theta_d/r)a,y'=(\theta_d/r)b x′=(θd​/r)a,y′=(θd​/r)b
Step5：最后根据相机内参将 p ′ p' p′转换到像素坐标系： u = f x x ′ + c x , v = f y y ′ + c y u=f_xx'+c_x,v=f_yy'+c_y u=fx​x′+cx​,v=fy​y′+cy​

  鱼眼相机的成像过程是已知入射角 θ \theta θ求出射角 θ d \theta_d θd​，而鱼眼相机的畸变矫正则是已知畸变后的像点位置 ( x ′ , y ′ ) (x',y') (x′,y′)求实际的入射角 θ \theta θ。由于相机参数已知，可以根据 ( x ′ , y ′ ) (x',y') (x′,y′)以及相机焦距求得 θ d \theta_d θd​的值，所以畸变矫正的本质问题是求解关于 θ \theta θ的一元高次方程： θ d = θ ( 1 + k 1 θ 2 + k 2 θ 4 + k 3 θ 6 + k 4 θ 8 ) \theta_d=\theta(1+k_1\theta^2+k_2\theta^4+k_3\theta^6+k_4\theta^8) θd​=θ(1+k1​θ2+k2​θ4+k3​θ6+k4​θ8)其中 k 1 , k 2 , k 3 , k 4 k_1,k_2,k_3,k_4 k1​,k2​,k3​,k4​是畸变参数，由相机标定结果提供。
  常用的求解一元高次方程的方法有二分法、不动点迭代、牛顿迭代法。这里使用牛顿迭代法求解： f ( θ ) = θ ( 1 + k 1 θ 2 + k 2 θ 4 + k 3 θ 6 + k 4 θ 8 ) − θ d θ 0 = θ d θ n + 1 = θ n − f ( θ n ) f ′ ( θ n ) f(\theta)=\theta(1+k_1\theta^2+k_2\theta^4+k_3\theta^6+k_4\theta^8) -\theta_d \\ \theta_0=\theta_d \\ \theta_{n+1}=\theta_n-\frac{f(\theta_n)}{f'(\theta_n)} f(θ)=θ(1+k1​θ2+k2​θ4+k3​θ6+k4​θ8)−θd​θ0​=θd​θn+1​=θn​−f′(θn​)f(θn​)​循环迭代直到 f ( θ ) ≈ 0 f(\theta)\approx0 f(θ)≈0（具体精度根据需要自行设置），或达到迭代次数上限。求得 θ \theta θ之后，未畸变像点 P 0 P_0 P0​到像平面中心的距离也可以求得： r = t a n ( θ ) r=tan(\theta) r=tan(θ)根据式(6)可以求得 P 0 P_0 P0​的坐标： { a = r θ d x ′ b = r θ d y ′ (8) \left \{ \begin{aligned}a&=\frac{r}{\theta_d}x' \\ b &=\frac{r}{\theta_d}y' \end{aligned}\right.\tag{8} ⎩⎪⎨⎪⎧​ab​=θd​r​x′=θd​r​y′​(8)最后利用相机内参将 P 0 ( a , b ) P_0(a, b) P0​(a,b)转换到像素坐标系即可得到未畸变的像素坐标。