Python计算机视觉——全景图像拼接

---

一、基本原理

  将两幅或多幅具有重叠区域的图像，通过特征匹配将具有相同的特征点（SIFT特征点）的图像拼接在一起，将来自多个不同视角拍摄的图像变换到同一视角下，拼接成一张宽视野图像

图像拼接基础流程

1. 针对某个场景拍摄多张/序列图像
2. 计算第二张图像与第一张图像之间的变换关系
   • 提取特征点，生成描述符
   • 特征匹配
3. 将第二张图像叠加到第一张图像的坐标系中
   • 图像映射、图像扭曲
4. 变换后的融合/合成
5. 在多图场景中，重复上述过程

1.1 图像拼接的几何原理

全景融合的 3D 几何解释

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

1.2 2D 图像变换

给定变换模型 x’= h(x) ，以及输入图像 f(x), 如何根据f(x)计算变换后的图像 g(x’) = f(h(x))?

一种方法是前向映射，一种方法是逆向映射。

2D 图像变换：前向映射

• 对于 $f(x)$ 中的每个像素 x，根据变换模型计算相应的映射坐标$x^{\prime }=h(x)$，并将x的像素值赋给 $g(x^{\prime })$

2D 图像变换：逆向映射

• 对于 $g(x^{\prime })$中的每个像素 $x^{\prime }$，根据变换模型计算相应的映射坐标$x=h^{-1}(x^{\prime })$，并将x的像素值赋给$g(x^{\prime })$

前向映射	逆向映射

如果像素落在两个像素之间会发生什么?

前向映射：近邻插值	逆向映射：线性插值/双线性插值

主要插值方法:

• 最近邻插值
• 线性插值
• 双线性插值
• 三线性插值

图像映射流程

1. 针对两张/多张图像提取特征
2. 特征匹配
3. 根据图像变换特点，选取合适的变换结构
4. 根据DLT等方法计算变换结构
5. 采用正向/逆向映射，利用插值方式实现图像映射变换

1.3 直线拟合

  给定若干二维空间中的点，求直线 $y=ax+b$ ，使得该直线对空间点的拟合误差最小。

如何选择该直线呢？

直线确定步骤：

• 随机选择两个点
• 根据该点构造直线
• 给定阈值，计算 inliers 数量

  我们在空间中随机寻找两个点，并以这两点构造出一条直线，然后计算在这条线上的对应 inliers 的数量，当 inliers 值达到我们需要的阈值，我们就得到了一条拟合直线。

  用这种方法我们可以拟合出一条直线，那么圆可以拟合吗？复杂的曲线方程可以拟合吗？

  对于圆，我们同样可以用该方法，两点确定一条直线，三点可以确定一个圆，我们随机确定三个点，然后确定一个圆的方程，然后计算在圆上的 inliers 值，之后的过程就与直线的一致了。

  对于复杂的曲线方程，我们需要先给定一个方程式（例如$y=ax+b{x}^2+c{x}^3...+n$），这个方程式的参数是我们要求的，求出一个解，使得在 inliers 值尽可能多，或者是求出的曲线与点之间是最逼近。

1.4 RANSAC算法

  RANSAC 是“RANdom SAmple Consensus”（随机一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC 基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。 即从一组数据集（包含噪声点的数据集）中，能够从中挑选出正确的点，获取能够正确拟合的参数模型。

使用 RANSAC 算法用一条直线来拟合包含噪声数据点集

ransac.py文件可以在GitHub上下载： ransac.py下载路径

  图像匹配中会出现匹配异常的情况，原因可能是噪声点的影响，使得对应点匹配出现问题；也可能是图像本身因为拍摄角度、曝光、缩放等因素导致的问题。

  实际上，SIFT 是具有很强稳健性的描述子，能够比其他描述子，例如图像块相关的 Harris 角点，产生更少的错误的匹配。但是该方法仍然远非完美。

1.4.1 RANSAC 求解单应矩阵

RANSAC loop:

1. 随机选择四对匹配特征
2. 根据直接线性变换解法DLT计算单应矩阵 H (唯一解)
3. 对所有匹配点，计算映射误差 $\epsilon =||pi’,Hpi||$
4. 根据误差阈值，确定 inliers（例如3-5像素）
5. 针对最大inliers集合，重新计算单应矩阵 H

  我们使用 RANSAC 算法来求解单应性矩阵，首先需要将下面模型类添加到 homography.py 文件中：

class RansacModel(object):
    """ 用于测试单应性矩阵的类，其中单应性矩阵是由网站http://www.scipy.org/Cookbook/RANSAC 上的ransac.py 计算出来的"""

    def __init__(self, debug=False):
        self.debug = debug

    def fit(self, data):
        """ 计算选取的4 个对应的单应性矩阵 """

        # 将其转置，来调用H_from_points() 计算单应性矩阵
        data = data.T

        # 映射的起始点
        fp = data[:3, :4]
        # 映射的目标点
        tp = data[3:, :4]

        # 计算单应性矩阵，然后返回
        return H_from_points(fp, tp)

    def get_error(self, data, H):
        """ 对所有的对应计算单应性矩阵，然后对每个变换后的点，返回相应的误差"""

        data = data.T

        # 映射的起始点
        fp = data[:3]
        # 映射的目标点
        tp = data[3:]

        # 变换fp
        fp_transformed = dot(H, fp)

        # 归一化齐次坐标
        # for i in range(3):
        #     fp_transformed[i] /= fp_transformed[2]
        fp_transformed = normalize(fp_transformed)

        # 返回每个点的误差
        return sqrt(sum((tp - fp_transformed) ** 2, axis=0))

  可以看到，这个类包含 fit() 方法。该方法仅仅接受由 ransac.py 选择的４个对应点对（data 中的前 4 个点对），然后拟合一个单应性矩阵。记住，4 个点对是计算单应性矩阵所需的最少数目。由于 get_error() 方法对每个对应点对使用该单应性矩阵，然后返回相应的平方距离之和，因此 RANSAC 算法能够判定哪些点对是正确的，哪些是错误的。在实际中，我们需要在距离上使用一个阈值来决定哪些单应性矩阵是合理的。为了方便使用，将下面的函数添加到 homography.py 文件中：

def H_from_ransac(fp, tp, model, maxiter=1000, match_theshold=10):
    """ 使用RANSAC 稳健性估计点对应间的单应性矩阵H（ransac.py 为从 http://www.scipy.org/Cookbook/RANSAC 下载的版本）

        # 输入：齐次坐标表示的点fp，tp（3×n 的数组）"""

    from PCV.tools import ransac

    # 对应点组
    data = vstack((fp, tp))

    # 计算H，并返回
    H, ransac_data = ransac.ransac(data.T, model, 4, maxiter, match_theshold, 10, return_all=True)
    return H, ransac_data['inliers']

  该函数同样允许提供阈值和最小期望的点对数目。最重要的参数是最大迭代次数：程序退出太早可能得到一个坏解；迭代次数太多会占用太多时间。函数的返回结果为单应性矩阵和对应该单应性矩阵的正确点对。

1.5 APAP算法

  在图像拼接融合的过程中，拼接融合后的图像可能会存在“鬼影现象”，即图像叠加后出现重影的现象。如下图所示：

  有问题，也有相对应的解决方法，2013年，Julio Zaragoza等人发表了一种新的图像配准算法APAP（As-Projective-As-Possible Image Stitching with Moving DLT），解决鬼影现象可以采用APAP算法。在拼接时，将图像分成很多个规则的小网格，分别计算局部的单应性变换矩阵，使整体的单应性变换误差降低。

  可以看出，使用 APAP算法后的拼接图像，“鬼影”的现象明显消失了。

APAP算法流程：

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

  Apap虽然能够较好地完成配准，但非常依赖于特征点对。若图像高频信息较少，特征点对过少，配准将完全失效，并且对大尺度的图像进行配准，其效果也不是很好，一切都决定于特征点对的数量。

1.6 寻找最佳拼接缝（Seam Finding）

  在两张图像的重叠区域选择一条最佳拼接缝，拼接缝左边的图像为第一张图像的内容，拼接缝右边的图像为第二张图像的内容，并且拼接缝上两幅图像的差异应该尽可能小，这样也可以有效的解决“鬼影”现象，而且拼接的效果也会更好。

  那么，如何找到最佳的拼接缝？这就需要使用到 Seam Finding（寻找最佳拼接缝）算法。算法主要原理运用了最大流最小割的思想。找到两张图像中差异最小的一条线，进行切割拼接，使得拼接缝两边的图像差异最小。

1.6.1 最大流问题

  最大流问题(maximum flow problem)，一种组合最优化问题。就是要讨论如何充分利用装置的能力，使得运输的流量最大，以取得最好的效果。

最大流问题归纳如下：

1. 所有流经网络（有向且流通）的流都起源于同一点，称为发点（源），终止与另一个点，称为收点（汇）。
2. 剩余的其他点称为中间点。
3. 流的方向由箭头表明，弧的容量就是允许的最大流量。在发点，所有的流都从这一点发出。在收点，所有的流都指向这一点。
4. 问题的目标是使从发点到收点的总流量达到最大。该值可以由两种等价的方式来衡量，即发点的流出量或收点的流入量。

1.6.2 最大流问题的增广链算法

  如何解决最大流问题，我们一般使用增广链算法，通过不断寻找网络中的增广链，直到增广链不存在了，从而确定了最大流。

最大流问题的增广链算法

1. 通过寻找剩余网络中从出发点到收点的正向链中每个弧上都有非零剩余容量的链而找到增广链（如果不存在增广链，该网络已达到最大流）。
2. 找出增广链中弧的最小剩余容量 $c^{\ast }$，这就是该增广链的剩余容量，该增广链中增加流量为 $c^{\ast }$的流。
3. 在增广链的每个正向弧的剩余容量中减去 $c^{\ast }$，而在每个反向弧的剩余容量中加上 $c^{\ast }$。返回步骤1。

  在进行步骤1时，往往有很多增广链可供选择。在用这种方法解决大规模问题时，增广链的正确选择对解决问题的效率是至关重要的。

1.6.3 最大流最小割定理

  对于任意一个只包含一个发点和收点的网络，从发点到收点的最大可行流等于网络所有割的最小割值。

二、图像拼接具体实现

  估计出图像间的单应性矩阵（使用 RANSAC 算法），现在我们需要将所有的图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心图像平面（否则，需要进行大量变形）。一种方法是创建一个很大的图像，比如图像中全部填充 0，使其和中心图像平行，然后将所有的图像扭曲到上面。由于我们所有的图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充 0，以便为扭曲的图像腾出空间。将下面的代码添加到 warp.py 文件中：

def panorama(H, fromim, toim, padding=2400, delta=2400):
    """ 使用单应性矩阵H（使用RANSAC 健壮性估计得出），协调两幅图像，创建水平全景图像。结果
    为一幅和toim 具有相同高度的图像。padding 指定填充像素的数目，delta 指定额外的平移量"""

    # 检查图像是灰度图像，还是彩色图像
    is_color = len(fromim.shape) == 3

    # 用于geometric_transform() 的单应性变换
    def transf(p):
        p2 = dot(H, [p[0], p[1], 1])
        return (p2[0] / p2[2], p2[1] / p2[2])

    if H[1, 2] < 0:  # fromim 在右边
        print('warp - right')
        # 变换fromim
        if is_color:
            # 在目标图像的右边填充0
            toim_t = hstack((toim, zeros((toim.shape[0], padding, 3))))
            fromim_t = zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
            for col in range(3):
                fromim_t[:, :, col] = ndimage.geometric_transform(fromim[:, :, col],
                                                                  transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # 在目标图像的右边填充0
            toim_t = hstack((toim, zeros((toim.shape[0], padding))))
            fromim_t = ndimage.geometric_transform(fromim, transf,
                                                   (toim.shape[0], toim.shape[1] + padding))
    else:
        print('warp - left')
        # 为了补偿填充效果，在左边加入平移量
        H_delta = array([[1, 0, 0], [0, 1, -delta], [0, 0, 1]])
        H = dot(H, H_delta)
        # fromim 变换
        if is_color:
            # 在目标图像的左边填充0
            toim_t = hstack((zeros((toim.shape[0], padding, 3)), toim))
            fromim_t = zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
            for col in range(3):
                fromim_t[:, :, col] = ndimage.geometric_transform(fromim[:, :, col],
                                                                  transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # 在目标图像的左边填充0
            toim_t = hstack((zeros((toim.shape[0], padding)), toim))
            fromim_t = ndimage.geometric_transform(fromim,
                                                   transf, (toim.shape[0], toim.shape[1] + padding))

    # 协调后返回（将fromim 放置在toim 上）
    if is_color:
        # 所有非黑色像素
        alpha = ((fromim_t[:, :, 0] * fromim_t[:, :, 1] * fromim_t[:, :, 2]) > 0)
        for col in range(3):
            toim_t[:, :, col] = fromim_t[:, :, col] * alpha + toim_t[:, :, col] * (1 - alpha)
    else:
        alpha = (fromim_t > 0)
        toim_t = fromim_t * alpha + toim_t * (1 - alpha)

    return toim_t

  对于通用的 geometric_transform() 函数，我们需要指定能够描述像素到像素间映射的函数。在这个例子中，transf() 函数就是该指定的函数。该函数通过将像素和 H 相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看 H 中的平移量， 我们可以决定应该将该图像填补到左边还是右边。当该图像填补到左边时，由于目标图像中点的坐标也变化了，所以在“左边”情况中，需要在单应性矩阵中加入平移。

2.1 运行代码

from pylab import *
from numpy import *
from PIL import Image

# If you have PCV installed, these imports should work
from PCV.geometry import homography, warp
from PCV.localdescriptors import sift

np.seterr(invalid='ignore')

"""
This is the panorama example from section 3.3.
"""

# 设置数据文件夹的路径
featname = ['image/jmu_lu' + str(i + 1) + '.sift' for i in range(5)]
imname = ['image/jmu_lu' + str(i + 1) + '.jpg' for i in range(5)]

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

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

# 可视化匹配
for i in range(4):
    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]  # im1 到im2 的单应性矩阵

fp, tp = convert_points(0)
H_01 = homography.H_from_ransac(fp, tp, model)[0]  # im0 到im1 的单应性矩阵

tp, fp = convert_points(2)  # 注意：点是反序的
H_32 = homography.H_from_ransac(fp, tp, model)[0]  # im3 到im2 的单应性矩阵

tp, fp = convert_points(3)  # 注意：点是反序的
H_43 = homography.H_from_ransac(fp, tp, model)[0]  # im4 到im3 的单应性矩阵

# 扭曲图像
delta = 2000  # 用于填充和平移

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)

im1 = array(Image.open(imname[3]), "f")
im_32 = warp.panorama(H_32, im1, im_02, delta, delta)

im1 = array(Image.open(imname[4]), "f")
im_42 = warp.panorama(dot(H_32, H_43), im1, im_32, delta, 2 * delta)

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

2.2 运行结果

2.2.1 第一组

图片集：

运行结果：

  从这组图片的拼接结果看，还是可以明显看出拼接缝，因为第三张图片的曝光程度和其它几张图片的曝光程度明显不同，使得拼接后的全景图可以明显的看出拼接的痕迹。

2.2.2 第二组

图片集：

运行结果：

  从这组图片的拼接情况来看，也是能够很明显的感觉到曝光参数不同导致的拼接差异明显。

2.2.3 第三组

图片集：

运行结果：

  这一组图片比起第一组，整体的视野更加宽广，但还是可以明显看出拼接缝，在拼接缝两边的图像并没有很好的匹配上，这可能是因为拍摄的角度出现偏差，使得拼接后的景物对不上。

2.2.4 第四组

图片集：

运行结果：

  这组图片从匹配结果来看，拼接的结果很差，拼接后的图像景物也没有对上。首先因为拍摄的原因，使得图片集有很多树木挡住了后面的建筑物，在匹配过程中，算法检测到的匹配点也比较少，使得最后的拼接效果很差。这也说明了 Apap算法对特征点对的依赖，若图像高频信息较少，特征点对过少，配准将完全失效，并且对大尺度的图像进行配准，其效果也不是很好，一切都决定于特征点对的数量。

2.2.5 第五组

图片集：

运行结果：

  这一组图片的拼接的整体效果看起来还可以，但也是能够很明显的看出拼接缝。

2.3 结果分析

  从上面几组图片的情况看，除了第四组的图片因为被树木遮挡了许多特征点对，导致匹配拼接的效果很差，其他几组的拼接效果都还可以，但是也都能够明显的看出拼接缝。

  首先一个原因是因为图片的曝光程度不同，第一组、第二组的情况可以展现出，因为曝光参数的不同，导致拼接后的全景图可以明显的看出明暗差异；也可能是因为图像景深的复杂程度导致的，第四组图片就是因为树木遮挡的后面建筑物的特征，使得算法在匹配时得到的特征点对不是很多，导致匹配拼接的结果很差。而另外几组相对的景深比较单一，没有太多复杂的背景，最后得到的拼接全景图效果就比较好。

  总体来说，用 Apap算法对景深落差小的图像进行拼接，能够有比较好的拼接效果；但是对于景深比较复杂的图像，效果相对会更差；但是都会有比较明显的拼接缝，拼接缝两边的景物也可能出现对不上的问题。同时，Apap算法对特征点对非常的依赖，若图像高频信息较少，特征点对过少，配准将完全失效，并且对大尺度的图像进行配准，其效果也不是很好。所以特征点对的数量是决定 Apap算法效果的一个重要因素。

图片集的编号顺序要从右到左顺序编号，同时图片集的大小要保持一致。