OpenCV Stitching图像拼接算法详解

OpenCV-4.x Stitching

流程图

模块详解

核心模块：

1. 特征点检测与图像匹配（matchers.hpp、matchers.cpp）
2. 估计相机运动，计算相机内外参（motion_estimators.hpp、motion_estimators.cpp）
3. 图像变形（warpers.hpp、warpers.cpp）
4. 计算接缝（seam_finders.hpp、seam_finders.cpp）
5. 补偿曝光（exposure_compensate.hpp、exposure_compensate.cpp）
6. 图像融合（blenders.hpp、blenders.hpp）

图像拼接模式：

1. PANORAMA（全景模式），此模式假设各个图像对应的相机光心位置一致，各相机光轴之间存在旋转关系。

2. SCANS（扫描模式），此模式假设各个图像对应的相机光心位于同一平面，各相机光轴与此平面垂直，各图像之间符合仿射变换的约束，各相机拍摄的是与光心所在平面平行的平面上的场景。

   **注：**官方文档要求图像之间满足仿射变换，假设同一空间点$P$在图像A和在图像B中的像素坐标分别为$(u_1,v_1)$和$(u_2,v_2)$，内参分别为$\mathbf{K_1}$和$\mathbf{K_2}$，在图像A对应的相机坐标系下的坐标和在图像B对应的相机坐标系下的坐标分别为$(X_1,Y_1,Z_1)$和$(Z_2,Y_2,Z_2)$，矩阵$\mathbf{R}$和平移$\mathbf{t}$描述图像B对应的相机相对于图像A对应的相机的位姿。则有

   $$\begin{bmatrix} u_1\\ v_1\\ 1 \end{bmatrix} = \frac{Z_2}{Z_1}\mathbf{K_1}\mathbf{R}\mathbf{K_2}^{-1} \begin{bmatrix} u_2\\ v_2\\ 1 \end{bmatrix} + \frac{1}{Z_1}\mathbf{K_1}\mathbf{t}$$

   根据官方要求，矩阵$\frac{Z_2}{Z_1}\mathbf{K_1}\mathbf{R}\mathbf{K_2}^{-1}$第三行应为$[0,0,1]$，平移$\frac{1}{Z_1}\mathbf{K_1}\mathbf{t}$第三维应为$0$. 由此可知$Z_1=Z_2$，且对于任意点$P$，$Z_1$为常数，因此可知SCANS模式下要求：各相机光轴与此平面垂直，各相机拍摄的是与光心所在平面平行的平面上。

特征点检测与图像匹配

模块过程：

1. 计算图像集合中每个图像的特征点集合及其描述子集合；
2. 根据对应图像的特征点和描述子，对图像集合进行两两配对；
3. 对匹配上的两张图像，根据对应的特征点集合，计算其变换矩阵$\mathbf{H}$。

PANORAMA模式：

图像匹配的class：BestOf2NearestMatcher
H是两个特征点集之间的单应变换矩阵，利用RANSAC算法，拟合出H，计算矩阵H的函数为：
cv::Mat cv::findHomography( InputArray _points1, InputArray _points2,
                           OutputArray _mask, int method, double ransacReprojThreshold )

SCANS模式：

图像匹配的class：AffineBestOf2NearestMatcher
H是两个特征点集之间的仿射变换矩阵，opencv的具体实现中，只包含旋转、缩放和平移，利用RANSAC算法，拟合出H，计算矩阵H的函数为：
CV_EXPORTS_W cv::Mat estimateAffinePartial2D(InputArray from, InputArray to, OutputArray inliers = noArray(),
                                  int method = RANSAC, double ransacReprojThreshold = 3,
                                  size_t maxIters = 2000, double confidence = 0.99,
                                  size_t refineIters = 10);

估计相机运动，计算相机内外参

模块过程：

1. 利用图像匹配过程中计算得到的$\mathbf{H}$，估计相机外参$\mathbf{R}$，并初始化相机内参$\mathbf{K}$。其中$\mathbf{H}$和$\mathbf{R}$是相对值，首先由多个图像为节点，各个图像之间的匹配关系为边，构造Spanning Tree，然后树中心节点对应的图像外参为单位矩阵，其他图像外参是基于此图像的位姿变化。
2. 利用Bundle Adjustment方法优化相机内外参；
3. 波形矫正，进一步优化相机外参。

PANORAMA模式：

估计相机外参的class：HomographyBasedEstimator
根据H计算R的过程：CalcRotation结构体及其operator()成员函数
优化相机参数的clas：BundleAdjusterRay
波形矫正优化策略：WAVE_CORRECT_HORIZ

PANORAMA模式下，根据两个图像像素之间的单应变换矩阵$\mathbf{H}$，计算相机旋转矩阵$\mathbf{R}$的过程如下：

假设同一空间点$P$在图像A和在图像B中的像素坐标分别为$(u_1,v_1)$和$(u_2,v_2)$，在图像A对应的相机坐标系下的坐标和在图像B对应的相机坐标系下的坐标分别为$(X_1,Y_1,Z_1)$和$(Z_2,Y_2,Z_2)$，则根据opencv stitching模块中的定义，$\mathbf{H}$满足

而根据相机成像模型，有

其中$\mathbf{K_1}$和$\mathbf{K_2}$分别为图像A对应的相机的内参和图像B对应的相机内参。而$\mathbf{R}$为图像B对应的相机相对于图像A对应的相机的旋转，因此有

由上述等式，可以求出$\mathbf{R}$满足

图像B对应的相机相对于图像A对应的相机的旋转矩阵$\mathbf{R}$由上式给出。

（注：opencv stitching模块的具体实现中旋转矩阵用$\mathbf{K_1}^{-1}\cdot\mathbf{H}^{-1}\mathbf{K_2}$来表示，相当于在实际上的旋转矩阵$\mathbf{R}$基础上多一个系数$\frac{Z_1}{Z_2}$. 此系数在利用Bundle Adjustment方法优化相机参数的过程中，通过cv::Rodrigues()函数将旋转矩阵转换为旋转向量的时候，去除了尺度信息）

SCANS模式：

估计相机外参的class：AffineBasedEstimator
根据H计算R的过程：CalcAffineTransform结构体及其operator()成员函数
优化相机参数的class：BundleAdjusterAffinePartial
波形矫正策略：不进行波形矫正

SCANS模式下，根据两个图像像素之间的仿射变换矩阵$\mathbf{H}$，估计相机旋转矩阵$\mathbf{R}$和平移$\mathbf{t}$的过程如下：

假设同一空间点$P$在图像A和在图像B中的像素坐标分别为$(u_1,v_1)$和$(u_2,v_2)$，则根据opencv stitching模块中的定义，$\mathbf{H}$是仿射矩阵，满足

根据仿射变换矩阵的性质，$\mathbf{H}$可表示为

根据仿射变换的性质，可将矩阵$\mathbf{H}$拆分为旋转缩放矩阵$\mathbf{R'}$和平移$\mathbf{t'}$，即

根据opencv中的定义，旋转矩阵$\mathbf{R}$和平移$\mathbf{t}$描述的是图像B对应的相机相对于图像A对应的相机的位姿，是与$\mathbf{H}$相反的过程，因此可得

此时得到的$\mathbf{R}$和$\mathbf{t}$都是基于像素坐标，即对于点$P$在图像A和在图像B中的像素坐标$(u_1,v_1)$和$(u_2,v_2)$，有

**经过bundle adjustment过程优化之后，由$\mathbf{H}$得到与真实的相机之间的位姿$\mathbf{R}$和$\mathbf{t}$。**对于$P$点在图像A对应的相机坐标系下的坐标和在图像B对应的相机坐标系下的坐标$(X_1,Y_1,Z_1)$和$(Z_2,Y_2,Z_2)$，因此有

图像变形

模块过程：

1. 建立原图像和其对应的投影面之间的像素对应关系；
2. 将图像本身的像素范围变换为对应到投影面上的ROI框；
3. 对投影面上的ROI框中每个像素坐标，在原图像中找到对应的像素点，将其值赋给投影面上ROI框中对应的坐标。

PANORAMA模式：

图像变形class：SphericalWarper
原图与投影面映射变换class：SphericalProjector

PANORAMA模式下，图像投影方式为球面投影，即将图像投影到以相机光心为球心的球面上。

注1： 以下涉及的齐次坐标，是指相对于相机$\mathbf{z}$轴方向上的齐次坐标。

注2： 球面坐标系即为半径为常数的球坐标系，具体定义可参考球坐标系。

1. 正向投影：

   首先，将图像像素坐标转换到相机坐标系下的齐次坐标。选取原图像中像素位置$(u,v)$，其对应的空间点$P$在当前相机坐标系下的坐标为$(X,Y,Z)$，则根据相机模型有

   $$\begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix} =\mathbf{K}^{-1} \begin{bmatrix} u\\ v\\ 1 \end{bmatrix}$$

   则$(\frac{X}{Z},\frac{Y}{Z},1)$即为$P$点当前相机坐标系下，投影到单位平面的齐次坐标，其中$\mathbf{K}$是当前相机内参。

   然后，将$P$点在当前相机坐标系下的齐次坐标，变换为中心相机坐标系下的齐次坐标$(X',Y',Z')$，即

   $$\begin{bmatrix} X'\\ Y'\\ Z' \end{bmatrix} =\mathbf{R} \begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix}$$

   其中，$\mathbf{R}$为当前相机相对于中心相机的旋转。

   最后，根据利用直角坐标系和球面坐标系的关系，将$P$点在中心相机坐标系下的齐次坐标$(X',Y',Z')$投影到，以中心相机光心为球心以$s$为半径的球面上得到$P$点在球面坐标系下的坐标$(u',v')$

   $$\begin{aligned} &u'=s\cdot \textbf{arctan2}(\frac{X'}{Z'})\\ &v'=s\cdot (\pi - \arccos\frac{Y'}{\sqrt{X'^2+Y'^2+Z'^2}}) \end{aligned}$$

   其中，函数$\textbf{arctan2}()$与c++中的函数定义一致。

2. 反向投影：

   反向投影即是正向投影的逆向操作，我们沿用上述符号定义。

   首先，将$P$点球面坐标$(u',v')$，转为中心相机坐标系下的齐次坐标

   $$\begin{bmatrix} \frac{X'}{\sqrt{X'^2+Y'^2+Z'^2}}\\ \frac{Y'}{\sqrt{X'^2+Y'^2+Z'^2}}\\ \frac{Z'}{\sqrt{X'^2+Y'^2+Z'^2}} \end{bmatrix} = \begin{bmatrix} \sin(\pi-\frac{v'}{s})\cdot\sin(\frac{u'}{s})\\ \cos(\pi-\frac{v'}{s})\\ \sin(\pi-\frac{v'}{s})\cdot\cos(\frac{u'}{s}) \end{bmatrix}$$

   然后，将中心相机坐标系下的齐次坐标转为原图像相机坐标系下的齐次坐标

   $$\begin{bmatrix} \frac{X}{\sqrt{X^2+Y^2+Z^2}}\\ \frac{Y}{\sqrt{X^2+Y^2+Z^2}}\\ \frac{Z}{\sqrt{X^2+Y^2+Z^2}} \end{bmatrix} = \mathbf{R}^{-1} \begin{bmatrix} \frac{X'}{\sqrt{X'^2+Y'^2+Z'^2}}\\ \frac{Y'}{\sqrt{X'^2+Y'^2+Z'^2}}\\ \frac{Z'}{\sqrt{X'^2+Y'^2+Z'^2}} \end{bmatrix}$$

   将原图像相机下的齐次坐标，除去$\mathbf{z}$轴方向的长度，得到

   $$\begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix} = \frac{\sqrt{X^2+Y^2+Z^2}}{Z} \begin{bmatrix} \frac{X}{\sqrt{X^2+Y^2+Z^2}}\\ \frac{Y}{\sqrt{X^2+Y^2+Z^2}}\\ \frac{Z}{\sqrt{X^2+Y^2+Z^2}} \end{bmatrix}$$

   最后，根据原图像相机的内参，计算出齐次坐标投影到成像平面上的像素坐标$(u,v)$，即

   $$\begin{bmatrix} u\\ v\\ 1 \end{bmatrix} = \mathbf{K} \begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix}$$

SCANS模式：

图像变形class：AffineWarper
原图与投影面映射变换class：PlaneProjector

SCANS模式下，图像投影方式仿射变换，即将图像投影到与其平行的平面，并利用相对于中心相机图像的仿射变换关系对图像进行变换。

1. 正向投影：

   首先，将图像像素坐标转换到相机坐标系下的齐次坐标。选取原图像中像素位置$(u,v)$，其对应的空间点$P$在当前相机坐标系下的坐标为$(X,Y,Z)$，则根据相机模型有

   $$\begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix} =\mathbf{K}^{-1} \begin{bmatrix} u\\ v\\ 1 \end{bmatrix}$$

   则$(\frac{X}{Z},\frac{Y}{Z},1)$即为$P$点当前相机坐标系下，投影到单位平面的齐次坐标，其中$\mathbf{K}$是当前相机内参。

   然后，根据$\mathbf{R}$和$\mathbf{t}$​，将$P$点在当前相机坐标系下的齐次坐标，变换为中心相机坐标系下的齐次坐标$(X',Y',Z')$，即

   $$\begin{bmatrix} \frac{X'}{Z'}\\ \frac{Y'}{Z'}\\ 1 \end{bmatrix} =\mathbf{R} \begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix} + \mathbf{t}$$

   其中$Z'=1$.

   最后，将投影面进行尺度$s$的缩放，得到$P$点在投影面上的像素坐标$(u',v')$：$u'=s\frac{X'}{Z'},v'=s\frac{Y'}{Z'}$. 因此正向投影，相当于把所有图像投影到距离相机光心所在平面为$1$的平面上，再将图像缩放$s$倍。

2. 反向投影：

   首先，$P$点将投影面图像坐标进行$\frac{1}{s}$尺度缩放，得到在中心相机下的齐次坐标$(\frac{X'}{Z'},\frac{Y'}{Z'},1)$，然后进行正向投影的逆操作

   $$\begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix} =\mathbf{R}^{-1} \begin{bmatrix} \frac{X'}{Z'}\\ \frac{Y'}{Z'}\\ 1 \end{bmatrix} + \mathbf{R}^{-1}\mathbf{t}$$

   最后，根据相机内参$\mathbf{K}$，得到原图像像素坐标$(u,v)$，满足

   $$\begin{bmatrix} u\\ v\\ 1 \end{bmatrix} =\mathbf{K} \begin{bmatrix} \frac{X}{Z}\\ \frac{Y}{Z}\\ 1 \end{bmatrix}$$

计算接缝、补偿曝光

PANORAMA模式：

计算接缝class：GraphCutSeamFinder，使用图形切割算法，可以查找图像中最佳的拼接缝隙
补偿曝光class：SphericalWarper，对图像中不同区域进行增益补偿

SCANS模式：

计算接缝class：GraphCutSeamFinder，使用图形切割算法，可以查找图像中最佳的拼接缝隙
补偿曝光class：BlocksGainCompensator，一种不使用曝光补偿的方法来校正拍摄照片中的一些失真

图像融合

模块过程：

1. 根据投影之后各个的图像尺寸和位置，确定最终拼接完成之后的图像尺寸；
2. 利用对应模式下的图像变形方式，将图像投影到最终图像上；
3. 将拼接好的图像中重叠的部分，利用多频带混合算法进行融合。

图像融合class：MultiBandBlender，图像拼接的多频带混合算法

多频带混合算法：

多频带混合的基本思想是图像可以分解为不同频率的图像的叠加，在不同的频率上，使用不同的权重来进行融合，在低频部分应该使用波长较宽的加权信号，在高频部分使用较窄的加权信号，其算法过程如下：

1. 计算输入图像的高斯金字塔；
2. 计算输入图像的拉普拉斯金字塔；
3. 将处于同一级的拉普拉斯金字塔进行融合；
4. 将高层的拉普拉斯金字塔依次扩展直至和相同分辨率；
5. 将上一步中得到的图像依次叠加，则得到最终的输出图像。