鱼眼参数的数值计算优化方法

        在下面两篇文章中，我们大概介绍了开源图像处理软件 Hugin 中鱼眼图像的矫正方法，即如何将目标全景图上的坐标映射到鱼眼图像上，从而获取相应的像素信息。

1. https://blog.csdn.net/qq_33552519/article/details/107362517
2. https://blog.csdn.net/qq_33552519/article/details/121039023

对于以上前面两篇文章，我们忽略了一个问题，即以上矫正所用到的参数大部分都是不准确或者未知的，这样会导致不同的鱼眼图像经过矫正后不能很好地拼接在一起，最直接的后果是在拼接缝处出现断层，从而使得全景图像的画面不连续，影响观感。因此我们需要对一些参数进行优化。这里主要还是参考 Hugin 中的方法，一种基于数值计算迭代优化的方法。

图1 鱼眼图像矫正流程

        首先，我们对鱼眼矫正的流程做个总结，如图 1 所示。那么，我们可优化的一些参数如表 1 所示。一般来说，为了覆盖整张全景图像即达到 $360°\times 180°$ 的视场，需要至少两张鱼眼图像，例如典型的用两个视场大于 $180°$ 的鱼眼镜头分别在偏航为 $0°$ 和 $180°$ 的进行拍摄，这样就可以保证覆盖整张全景图，同时保留一定的重叠区域以方便参数的优化以及保证图像融合的效果。当然使用更多的鱼眼镜头可以使得画面更加精细，毕竟鱼眼镜头越远离光轴的景象越被压缩，而且重叠区面积增大也可以提供更多的控制点对，提升参数优化和图像融合的质量，但也会带来成本增加，机械设备误差或者故障发生的风险增大的问题。另外，增加一个相机意味着增加一组参数，那么参数优化的复杂度也会相应提升。

表1 可优化参数列表

	参数类型	参数列表	备注
镜头参数	镜头有效视场	FOV	Field of View，圆形
	镜头径向畸变参数	a, b, c	径向畸变矫正参数
	传感器倾斜参数	tilt, slant, rotate	较少使用
	传感器中心偏移	offset_x, offset_y	–
	剪切效应参数	shear_x, shear_y	较少使用
方位参数	镜头姿态	yaw, pitch, roll	偏航，俯仰，旋转
	光心偏移	X, Y, Z	详见文章 2
	平面偏航、俯仰	plane_y, plane_p	详见文章 2

        在系统参数优化问题中，参数的数量越多，所需要的约束条件就越多，同时优化问题收敛的条件就越苛刻。为了减少需要优化的参数的数量，通常假设这些相机的固有镜头参数是相同的，即包括视场、径向畸变参数、传感器倾斜和中心偏移等这些外部一般不能改变的参数，因此只需要优化一组镜头参数。另外，通常也会选定一个相机作为全景图像的基准，称为主相机，这样就可以减少一组需要优化的方位参数，其他的相机则可以主相机为基准进行优化。注意，这里是参数优化，而不是纯粹的参数从无到有的求解，所以每个相机都应该首先提供一个合理的参数初始值，否则会大大地提升参数优化的难度，导致无法收敛。

        为了进行参数优化，我们首先得确定一个优化的目标。我们知道，物理空间的一个点以不同的方位参数拍摄会落在不同的鱼眼图像的不同位置，把这些对应着物理空间同一个点的不同位置称为控制点，一般来说以两个作为一对。因为目标全景图像只有一张，那么该物理空间点可以通过球面透视表示为唯一的球面坐标。因此，利用这些鱼眼图像各自的参数将控制点对根据图 1 从下至上的过程逆推回全景图像对应的球面时，理论上这些坐标应该是重合的。所以，我们的目标就是在不同的鱼眼图像上找出一系列的控制点对，然后优化出最佳的参数，使得这些控制点对逆映射到球面上时，彼此之间的距离最小。

        假设图 1 从下至上的逆映射过程可定义为一个非线性函数，即

$${\mathbf{x}}_p=\mathbf{g}\left({\mathbf{x}}_f;{\theta }_f\right),{\mathbf{x}}_f\in {\mathbb{R}}^2,{\mathbf{x}}_p\in {\mathbb{R}}^3,\Vert {\mathbf{x}}_p\Vert =1.\text{\hspace{1em}(1)}$$

其中 ${\mathbf{x}}_f$ 代表鱼眼图像上的坐标，${\mathbf{x}}_p$ 代表全景图像对应的球面坐标，${\theta }_f$ 代表表 1 中的参数，不同鱼眼图像对应着不同的参数。我们假设这个函数是连续可导的。虽然在开头两篇文章的分析中，我们主要讨论了图 1 从上至下的推导过程，但我们也能很容易推导出其逆向映射过程，除了镜头径向畸变校正这一步。因为我们采用的畸变校正模型中，实际落点距离 $r$ 是关于入射角 $\theta$ 的多项式函数，那么在逆映射的过程中，我们需要从实际落点距离 $r$ 逆推回入射角 $\theta$，即使用校正模型的反函数。然而，该反函数我们很难推导，因此可以采用数值近似的方法，一种常用的方法是牛顿迭代法，这里不细说。

        如果我们有一对控制点 $({\mathbf{x}}_{f1},{\mathbf{x}}_{f2})$，可以得到它们的球面坐标对 $({\mathbf{x}}_{p1},{\mathbf{x}}_{p2})$，那么定义这一对控制点的代价函数为

$$l\left(\mathbf{x};\theta \right)=\psi \left({\mathbf{x}}_{p1},{\mathbf{x}}_{p2}\right),\mathbf{x}={\left[{\mathbf{x}}_{f1}^T,{\mathbf{x}}_{f2}^T\right]}^T,\theta ={\left[{\theta }_{f1}^T,{\theta }_{f2}^T\right]}^T.\text{\hspace{1em}(2)}$$

图2 球面距离示意图

其中 $\psi$ 是两点的距离函数，可根据实际情况选择。通常来说，因为这两个点是在球面上的，因此一般使用球面距离，而不使用欧几里德距离。如图 2 所示，假如球面上有两点 $A({\phi }_1,{\theta }_1)$, $B({\phi }_2,{\theta }_2)$，其中 $\phi$, $\theta$ 分别是纬度和经度。那么我们可定义其球面距离为 A, B 之间的大圆劣弧长度，即

$$\psi \left(\mathbf{a},\mathbf{b}\right)=R\cdot \angle AOB=R\cdot \arccos \left(\frac{{\mathbf{a}}^T\mathbf{b}}{\Vert \mathbf{a}\Vert \Vert \mathbf{b}\Vert }\right)=R\cdot \arcsin \left(\frac{\Vert \mathbf{a}\times \mathbf{b}\Vert }{\Vert \mathbf{a}\Vert \Vert \mathbf{b}\Vert }\right).\text{\hspace{1em}(3)}$$

那么，根据式（3）我们就得到了一条约束方程。另外，我们也可以分别计算 A, B 两点的南北和东西方向距离，其中东西方向距离可定义为平均纬度上的小圆劣弧长度，南北方向距离则可定义为纬度之间沿着经线的弧线长度，即

$${\psi }_{\phi }\left(\mathbf{a},\mathbf{b}\right)=R\cdot \left(\left|{\phi }_1-{\phi }_2\right|\right),\text{\hspace{1em}(4.1)}$$

$${\psi }_{\theta }\left(\mathbf{a},\mathbf{b}\right)=R\cdot \cos \left({\phi }_1+{\phi }_2\right)\cdot \left(\left|{\theta }_1-{\theta }_2\right|\right).\text{\hspace{1em}(4.2)}$$

这样，我们实际上就得到了两个约束方程。一般来说，在优化过程中，我们不会同时使用式（3）和式（4）。

        由此，我们可以定义参数优化所用到的代价函数为

$$F\left(\mathbf{\Theta }\right)=\frac{1}{2}\sum _{i=1}^m{\left(f_i\left(\mathbf{\Theta }\right)\right)}^2=\frac{1}{2}{\Vert \mathbf{f}\left(\mathbf{\Theta }\right)\Vert }^2,,\text{\hspace{1em}(5.1)}$$

$$f_i\left(\mathbf{\Theta }\right)=l_i\left({\mathbf{x}}_i;{\theta }_i\right),\mathbf{\Theta }={\theta }_1\cup {\theta }_2\cup \cdots \cup {\theta }_m\in {\mathbb{R}}^n.\text{\hspace{1em}(5.2)}$$

其中 $m$ 代表总共有 $m$ 个约束方程，注意其不一定等于控制点的对数，因为在式（4）中一个控制点对可以定义两个约束方程。$n$ 代表总共有 $n$ 个需要优化的参数，主要每个控制点对用到的鱼眼图像可能不同，所以其参数也就不一定相同，但也有一些公共的参数，例如镜头参数，因此这里使用了它们的并集。尽管在式（3）和（4）中我们使用的是球面距离，但在求所有控制点对的整体距离时还是使用了平方和，这相当于把球面距离当作一种样本误差，我们的目标是要最小化整体的误差平方和。对于 $F\left(\mathbf{\Theta }\right)$ 来说我们其实并不知道其真正的函数形式，但为了方便一般假设其是一个高度非线性的连续函数，这时可将问题形式化为

Θ ∗ = arg ⁡ min ⁡ Θ { F ( Θ ) } . (6) {{\bf{\Theta }}^ * } = \arg {\min _{\bf{\Theta }}}\left\{ {F\left( {\bf{\Theta }} \right)} \right\}.\tag{6} Θ∗=argΘmin​{F(Θ)}.(6)

        由于 $F\left(\mathbf{\Theta }\right)$ 为平方和形式，因此也将式（6）称为最小二乘问题，相关的优化算法也有很多，比较常用的是 LM 算法，其原理可参考我的另一篇博客。注意，在 LM 算法优化过程中需要用到 $F\left(\mathbf{\Theta }\right)$ 的 Jacobian 矩阵，也就是各 $f_i\left(\mathbf{\Theta }\right)$ 关于各个参数的一阶偏导组成的矩阵，由于我们并不知道各 $f_i\left(\mathbf{\Theta }\right)$ 的具体形式，因此可通过一些数值近似方法，例如比较简单的有限差分法，即给某一个参数增加非常小的步长 $\delta$，其他的参数保持不变，求解出新的函数值，那么函数值的差值除以 $\delta$ 即为关于该参数的近似导数，即

$$f(x+\delta )=f(x)+f^{\prime }(x)\delta +o(\delta ),\text{\hspace{1em}(7.1)}$$

$$f^{\prime }(x)\approx \frac{f(x+\delta )-f(x)}{\delta },\delta \text{ is sufficient samll.}\text{\hspace{1em}(7.2)}$$

不过这种方法计算量比较大，而且 $\delta$ 应该足够小，否则不一定能满足有限差分的收敛性条件。

        在 LM 算法中，我们需要保证 $F\left(\mathbf{\Theta }\right)$ 的 Jacobian 矩阵是列满秩的，一般来说也就要求约束方程要比需要优化的参数多，所以我们应该提供足够多的质量好的控制点对。PanoTools 全景工具库中使用了两种优化策略，也就是式（3）和（4）的区别。式（3）中一个控制点对只能提供一个约束方程，由于约束性相对不那么强，因此收敛速度会慢一些，但是稳定性相对较好。而式（4）中一个控制点对提供了两个约束方程，约束性要强一些，因此在初始参数设置较好时会有较快的收敛速度，但如果初始参数较差时稳定性会差一些。因此 PanoTools 中也采用了混合的优化策略，即前期使用式（3），只要每次迭代 $F\left(\mathbf{\Theta }\right)$ 的下降比例都达到一定要求，否则切换到式（4）继续优化，这样就可以提高收敛速度和稳定性。