【计算机视觉】：（3）全景图像拼接

1.全景拼接基本流程

① 针对某个场景拍摄多张/序列图像

② 计算第二张图像与第一张图像之间的变换关系

③ 将第二张图像叠加到第一张图像的坐标系中

④ 变换后的融合/合成

⑤ 在多图场景中，重复上述过程

2. 图像映射流程

2.1. 提取特征点，构建关键点描述符（SIFT特征提取）

2.2. 特征匹配（KNN或kd树匹配）

上述用到SIFT特征匹配算法，详情可见https://blog.csdn.net/qq_45802659/article/details/123769669

2.3. 根据图像变换特点，选取合适的变换类型

图像的映射类型有：平移、旋转、仿射、透视映射、尺度变换

刚体变换：平移+旋转，只改变物体位置，不改变物体形状
仿射变换：改变物体位置和形状，但是保持“平直性”（原来平行的边依然平行）
投影变换：彻底改变物体位置和形状

2.4. 根据单应性变换等方法计算变换结构

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。实现单应变换的关键是求出单应性矩阵，而针对不同的变换类型需要相应数量的对应点对来求出相应的单应性矩阵（求解矩阵中的未知数）。
平移变换：2个自由度（需要一个点对）
相似变换：4个自由度（需要两个点对）
仿射变换：6个自由度（需要三个点对）
投影变换：8个自由度（需要四个点对）


进一步转换为：

写成矩阵Ah=0形式：

H单应性矩阵，8个自由度，矩阵可以缩放，使得$h_{22}$为1，所以实际上是8个自由度，至少要4个点对（几个未知数几个方程），在真实的应用场景中，我们计算的点对中都会包含噪声。比如点的位置偏差几个像素，甚至出现特征点对误匹配的现象，如果只使用4个点对来计算单应矩阵，那会出现很大的误差。因此，为了使得计算更精确，一般都会使用远大于4个点对来计算单应矩阵。

可通过最小二乘法求解，但容易受到噪声影响，后面会用到RANSAC随机抽样一致算法计算单应性矩阵。

2.5. 采用正向或逆向映射，利用插值方式实现图像映射变换

给定变换模型 x’= h(x) ，以及输入图像f(x), 如何根据f(x)计算变换后的图像 g(x’) = f(h(x))?即有了单应性矩阵，我们应该就行像素的映射f(x)变为g(x’) 。两种方法：前向映射和逆向映射
前向映射：扫描原图的每个像素(x,y)，变换后的值填入对应的(x’,y’)处。
对于多个像素点映射到同一个(x’,y’)的情况，不太好处理，所以反向映射更有效。
当像素落在两个像素之间时：近邻插值

逆向映射：扫描变换后图像的每个像素（x‘,y’），根据逆变换找到原图中对应的（x,y）
大概率找到的不是整点，所以需要使用图像插值

当像素落在两个像素之间时：线性插值/双线性插值

主要的插值方法:
最近邻插值；
线性插值；
双线性插值；
三线性插值。

3. 图像拼接的几何原理

• 全景融合的 3D 几何解释
– 图像被投影到共同的拼接平面上（同一坐标系）
– 在拼接平面上实现全景融合
– 在拼接的应用中，其实可以简化理解为 2D图像的变换，叠加过程

3. 全景拼接流程

1. 根据给定图像/集，实现特征匹配
2. 通过匹配特征计算图像之间的变换结构(RANSAC算法得出单应性矩阵)
3. 利用图像变换结构，实现图像映射
4. 针对叠加后的图像，采用APAP之类的算法，对齐特征点
5. 通过图割方法，自动选取拼接缝
6. 根据multi-band blending策略实现融合

3.1. SIFT特征匹配

此处不赘述，详情可见https://blog.csdn.net/qq_45802659/article/details/123769669

3.2. RANSAC算法

RANSAC算法的基本假设是样本中包含正确数据(inliers，可以被模型描述的数据)，也包含异常数据(outliers，偏离正常范围很远、无法适应数学模型的数据)，即数据集中含有噪声（SIFT特征匹配虽然效果好但难免因为不同地方的类似图案出现匹配误差）。这些异常数据可能是由于错误的测量、错误的假设、错误的计算等产生的。同时RANSAC也假设，给定一组正确的数据，存在可以计算出符合这些数据的模型参数的方法。简单来说RANSAC反复取样找到误差最小的模型。

RANSAC算法，随机选取k个样本点拟合模型再设定阈值，计算其他点到该拟合模型的距离，小于阈值作为内点，统计内点个数，迭代计算出内点个数最多的模型。可以看出噪声影响小。而最小二乘容易受到噪声干扰。（比如RANSAC和最小二乘拟合直线、曲线等）

RANSAC算法步骤：

1. 随机选择四对匹配特征
2. 根据直接线性变换解法DLT（图像映射中的推导过程）计算单应性矩阵 H (唯一解)
3. 对所有匹配点，计算映射误差ε= ||pi’, H pi|| （可算欧式距离）
4. 根据误差阈值，确定inliers（例如3-5像素）
5. 针对最大inliers集合，重新计算单应矩阵 H（迭代计算，达到迭代次数就停止，求出符合最多匹配点对的H）

3.4. APAP算法

APAP参考：参考1；参考2
论文原文

由于图像是将三维景物显示在2D平面上，图像中的景物的三维信息丢失，所以拼接时经常会有“鬼影”的出现，对目前已有方法提出的单一映射矩阵，只能应用于输入图像在一个平面上，或者相机沿着同一个焦点旋转拍摄的图像。如果违背了这些约束条件，图像拼接就会产生鬼影现象。APAP的算法主旨是得到global projective，但是能容忍变换矩阵不够准确带来的局部偏差。

下图来自论文As-Projective-As-Possible Image Stitching with Moving DLT，来自CVPR 2013。
图中的点表示匹配上的特征点，横纵坐标表示匹配上的特征点的坐标位置，红线是拟合的模型即单应性矩阵。（a）用一个2x2的单一映射矩阵，从图可以看出局部偏差是不可避免的。（c）用APAP的方法，局部偏差更少，并且正确的推断出全局变换的趋势。

APAP算法全称As-Projective-As-Possible Image Stitching，也是图像拼接算法的一种。

APAP算法流程如下：

1. SIFT得到两幅图像的匹配点对
2. 通过RANSAC剔除外点，得到N对内点
3. 利用DLT和SVD计算全局单应性
4. 将源图划分网格，取网格中心点，计算每个中心点和源图上内点之间的欧式距离和权重
5. 将权重放到DLT算法的A矩阵中，构建成新的W*A矩阵，重新SVD分解，自然就得到了当前网格的局部单应性矩阵
6. 遍历每个网格，利用局部单应性矩阵映射到全景画布上，就得到了APAP变换后的源图
7. 最后就是进行拼接线的加权融合

如果两幅待拼接的图像只有旋转变换或者拍摄同一平面的场景，那么原先的公式$\widetilde{x^{{}^{\prime }}}=H\tilde{x}$可以计算出准确的投影矩阵。如果不能满足这个条件，那计算出来的变换矩阵就会存在误差。矩阵存在误差的结果就是，拼接的结果在两张图像的重叠区域会出现鬼影。为了解决这个问题，论文的中心思想是对每个像素点x使用局部单一映射矩阵。
$$\widetilde{x_{\ast }^{{}^{\prime }}}=H_{\ast }\widetilde{x_{\ast }}$$
每次的局部单应性矩阵都要加上权重w

权重会随着像素x变化，其中，σ 是比例系数，$x_i$是原始图像的匹配的特征点的坐标。从公式可以看出，离像素点x越近的特征点的权重越高。同样地，随着x在图像域连续平移，计算的变换矩阵H也是平缓变化的。

改写成：

太远的点权重不稳定过于接近0，通过γ抵消权重

3.5. 寻找最佳拼接缝（Seam-finding）,最大流最小割算法

参考
最小割集的流量之和=网络中的最大流
每个重叠区域都会包含Patch A和Patch B对应的像素点，这个时候采用最小割算法，将重叠部分以图的形式表示，边的权重计算方式如上表示，最后找到权重最小的割线作为拼接缝。最大流算法的结果=最小割算法的结果。

将匹配点最为接近的点作为拼接缝上的点，最终构造出拼接缝。

3.6. multi-band blending算法

参考1
参考2
在通过图割算法找到拼接缝后，由于图像光照、噪声等因素，会使拼接缝很突兀。而multi-band blending算法通过图像融合消除边缘痕迹。融合目的在于拼缝消除, Multi-Band能够达到比较好的融合效果，但是效率低，采用Laplacian（拉普拉斯）金字塔，通过对相邻两层的高斯金字塔进行差分，将原图分解成不同尺度的子图，对每一个之图进行加权平均，得到每一层的融合结果，最后进行金字塔的反向重建，得到最终融合效果过程。
思想：采用的方法是直接对带拼接的两个图片进行拉普拉斯金字塔分解，后一半对前一半进行融合。
步骤：
（1）首先计算当前待拼接图像和已合成图像的重叠部分。

（2）两幅图像的融合：分别构建图像A、B的高斯金字塔和拉普拉斯金字塔，然后进行加权融合。

（3）对加权后的拉普拉斯金字塔进行重构。

5. 代码实现

from pylab import *
from numpy import *
from PIL import Image
from PCV.geometry import homography, warp
from PCV.localdescriptors import sift
# from scipy.spatial import Delaunay


# 设置数据文件夹的路径
featname = ['D:\study\计算机视觉\code\code3\image' + str(i + 1) + '.sift' for i in range(3)]
imname = ['D:\study\计算机视觉\code\code3\image' + str(i + 1) + '.jpg' for i in range(3)]

# 提取特征并匹配使用sift算法
l = {}
d = {}
for i in range(3):
    sift.process_image(imname[i], featname[i])
    l[i], d[i] = sift.read_features_from_file(featname[i])

matches = {}
for i in range(2):
    matches[i] = sift.match(d[i + 1], d[i])

# 可视化匹配
for i in range(2):
    im1 = array(Image.open(imname[i]))
    im2 = array(Image.open(imname[i + 1]))
    figure()
    sift.plot_matches(im2, im1, l[i + 1], l[i], matches[i], show_below=True)


# 将匹配转换成齐次坐标点的函数
def convert_points(j):
    ndx = matches[j].nonzero()[0]
    fp = homography.make_homog(l[j + 1][ndx, :2].T)
    ndx2 = [int(matches[j][i]) for i in ndx]
    tp = homography.make_homog(l[j][ndx2, :2].T)

    # switch x and y - TODO this should move elsewhere
    fp = vstack([fp[1], fp[0], fp[2]])
    tp = vstack([tp[1], tp[0], tp[2]])
    return fp, tp


# 估计单应性矩阵
model = homography.RansacModel()

fp, tp = convert_points(1)
H_12 = homography.H_from_ransac(fp, tp, model)[0]  # im 1 to 2

fp, tp = convert_points(0)
H_01 = homography.H_from_ransac(fp, tp, model)[0]  # im 0 to 1

# 扭曲图像
delta = 500  # for padding and translation用于填充和平移

im1 = array(Image.open(imname[1]), "uint8")
im2 = array(Image.open(imname[2]), "uint8")
im_12 = warp.panorama(H_12, im1, im2, delta, delta)

im1 = array(Image.open(imname[0]), "f")
im_02 = warp.panorama(dot(H_12, H_01), im1, im_12, delta, delta)

figure()
imshow(array(im_02, "uint8"))
axis('off')
show()

结果：

6. 问题解决

No module named 'matplotlib.delaunay’问题解决见博客
解决方法