图像拼接技术就是将数张有重叠部分的图像（可能是不同时间、不同视角或者不同传感器获得的）拼成一幅无缝的全景图或高分辨率图像的技术。图像拼接在医学成像、计算机视觉、卫星数据、军事目标自动识别等领域具有重要意义。图像拼接的输出是两个输入图像的并集。
图像配准（image alignment）和图像融合是图像拼接的两个关键技术。图像配准是图像融合的基础,而且图像配准算法的计算量一般非常大,因此图像拼接技术的发展很大程度上取决于图像配准技术的创新。早期的图像配准技术主要采用点匹配法,这类方法速度慢、精度低,而且常常需要人工选取初始匹配点,无法适应大数据量图像的融合。图像拼接的方法很多,不同的算法步骤会有一定差异,但大致的过程是相同的。

全景图拼接步骤

将SIFT应用到图像拼接上，根据特征点匹配的方式，则利用这些匹配的点来估算单应矩阵使用RANSAC算法，也就是把其中一张通过个关联性和另一张匹配的方法。通过单应矩阵H，可以将原图像中任意像素点坐标转换为新坐标点，转换后的图像即为适合拼接的结果图像。
可以简单分为以下几步:
1.根据给定图像/集，实现特征匹配。
2.通过匹配特征计算图像之间的变换结构。
3.利用图像变换结构，实现图像映射。
4.针对叠加后的图像，采用APAP之类的算法，对齐特征点。(图像配准)
5.通过图割方法，自动选取拼接缝。
6.根据multi-band blending策略实现融合。
?

全景图拼接原理

RANSAC算法

RANSAC[1] (随机抽样一致)是一种迭代算法，该算法从一组包含“外点(outlier)”的观测数据中估计数学模型的参数。“外点”指观测数据中的无效数据，通常为噪声或错误数据，比如图像匹配中的误匹配点和曲线拟合中的离群点。与“外点”相对应的是“内点(inlier)”，即用来估计模型参数的有效数据。因此，RANSAC也是一种“外点”检测算法。此外，RANSAC算法是一种非确定算法，它只能在一定概率下产生可信的结果，当迭代次数增加时，准确的概率也会增加。

RANSAC算法是用来找到正确模型来拟合带有噪声数据的迭代方法。
基本思想：数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摈弃噪声点。

RANSAC的基本思想和算法流程如下：

随机采样K个点，K是求解模型参数的最少点个数
使用K个点估计模型参数
计算剩余点到估计模型的距离，距离小于阈值则为内点，统计内点的数目
重复步骤1~3，重复次数M且保留数目最多的内点
使用所有的内点重新估计模型
?

?图像配准

图像配准是对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。为了能够进行图像对比和更精细的图像分析，图像配准是一步非常重要的操作。图像配准的方法有很多，这里以APAP算法为例：

提取两张图片的sift特征点
对两张图片的特征点进行匹配
匹配后，仍有很多错误点，此时用RANSAC进行特征点对的筛选。筛选后的特征点基本能够一一对应。
使用DLT算法，将剩下的特征点对进行透视变换矩阵的估计。
因为得到的透视变换矩阵是基于全局特征点对进行的，即一个刚性的单应性矩阵完成配准。为提高配准的精度，Apap将图像切割成无数多个小方块，对每个小方块的变换矩阵逐一估计。
?

图割方法

最大流最小割算法原理，
1.最小割问题
一个有向图，并有一个源顶点（source vertex）和目标顶点（target vertex）.边的权值为正,又称之为容量（capacity）。如下图

在这里插入图片描述

一个st-cut(简称割cut)会把有向图的顶点分成两个不相交的集合，其中s在一个集合中，t在另外一个集合中。
这个割的容量（capacity of the cut）就是A到B所有边的容量和。注意这里不包含B到A的。最小割问题就是要找到割容量最小的情况。
2.最大流问题
跟mincut问题类似，maxflow要处理的情况也是一个有向图，并有一个原顶点（source vertex）和目标（target vertex），边的权值为正,又称之为容量（capacity）。
（1）初始化，所有边的flow都初始化为0。
（2）沿着增广路径增加flow。增广路径是一条从s到t的无向路径，但也有些条件，可以经过没有满容量的前向路径（s到t）或者是不为空的反向路径(t->s)。
?

图像融合

图像拼接之后可以发现，在拼接的交界处有明显的衔接痕迹，存在边缘效应，因为光照色泽的原因使得图片交界处的过渡很糟糕，所以需要特定的处理解决这种不自然。那么这时候可以采用blending方法。multi-band blending是目前图像融和方面比较好的方法。
原理：

建立两幅图像的拉普拉斯金字塔
求高斯金字塔（掩模金字塔-为了拼接左右两幅图像）因为其具有尺度不变性
进行拼接blendLapPyrs() ; 在每一层上将左右laplacian图像直接拼起来得结果金字塔resultLapPyr
重建图像: 从最高层结果图

将左右laplacian图像拼成的resultLapPyr金字塔中每一层，从上到下插值放大并和下一层相加，即得blend图像结果（reconstructImgFromLapPyramid）
且我们可以将拉普拉斯金字塔理解为高斯金字塔的逆形式。

APAP算法

在图像拼接融合的过程中，受客观因素的影响，拼接融合后的图像可能会存在“鬼影现象”以及图像间过度不连续等问题。下图就是图像拼接的一种“鬼影现象”。而APAP算法可以解决鬼影现象。

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

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

multi-band bleing算法

在找完拼接缝后，由于图像噪声、光照、曝光度、模型匹配误差等因素，直接进行图像合成会在图像重叠区域的拼接处出现比较明显的边痕迹。

这些边痕迹需要使用图像融合算法来消除。这里介绍一种方法—multi-band bleing
思想：
采用的方法是直接对带拼接的两个图片进行拉普拉斯金字塔分解，后一半对前一半进行融合
步骤：
首先计算当前待拼接图像和已合成图像的重叠部分。并对图像A、B 重叠部分进行高斯金字塔和拉普拉斯金字塔分解

G0为原始图像，G1表示对G0做reduce操作。Reduce操作定义如下：?

对G1进行扩展后与G0相减，可以得到拉普拉斯金字塔的第一层L0。同理，拉普拉斯金字塔的L2、L3等层也可以按照这种方法来计算。
两幅图像的融合过程：分别构建图像A、B的高斯金字塔和拉普拉斯金字塔，然后进行加权融合。
对加权后的拉普拉斯金字塔进行重构

全景图拼接代码实现?

估计出单应性矩阵H后，我们需要将所有图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心平面（否则，需要进行大量变形）。一种方法是创建一个很大的图像，比如将平面中全部填充为0，使其和中心图像平行，然后将所有的图像扭曲到上面。其基本步骤如下：

（1）实现像素间的映射（计算像素和和单应性矩阵H相乘，然后对齐次坐标进行归一化）

（2）判断图像填补位置（查看H中的平移量，左边为正，右边为负）

（3）在单应性矩阵中加入平移量，进行alpha图填充

图片集为五张的拼接

# -*- codeing =utf-8 -*-
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')


# 设置数据文件夹的路径
featname = ['D:\\全景拼接\\' + str(i + 1) + '.sift' for i in range(5)]
imname = ['D:\\全景拼接\\' + 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]  # im 1 to 2

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

tp, fp = convert_points(2)  # NB: reverse order
H_32 = homography.H_from_ransac(fp, tp, model)[0]  # im 3 to 2

tp, fp = convert_points(3)  # NB: reverse order
H_43 = homography.H_from_ransac(fp, tp, model)[0]  # im 4 to 3

# 扭曲图像
delta = 2000  # 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)

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()

?图片集为三张的拼接

# -*- coding: utf-8 -*-
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
"""
This is the panorama example from section 3.3.
"""
# set paths to data folder
featname = ['D:\\全景拼接\\' + str(i + 1) + '.sift' for i in range(3)]
imname = ['D:\\全景拼接\\' + str(i + 1) + '.jpg' for i in range(3)]
# extract features and match
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])
    print(matches)
# visualize the matches (Figure 3-11 in the book)
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)
# function to convert the matches to hom. points
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
# estimate the homographies
model = homography.RansacModel()
fp,tp = convert_points(0)
H_01 = homography.H_from_ransac(fp,tp,model)[0] #im 0 to 1
tp,fp = convert_points(1) #NB: reverse order
H_21 = homography.H_from_ransac(fp,tp,model)[0] #im 2 to 1
# warp the images
delta = 2000 # for padding and translation
im1 = array(Image.open(imname[0]), "uint8")
im2 = array(Image.open(imname[1]), "uint8")
im_01 = warp.panorama(H_01,im1,im2,delta,delta)
im1 = array(Image.open(imname[2]), "f")
im_21 = warp.panorama(H_21,im1,im_01,delta,delta)
figure()
imshow(array(im_21, "uint8"))
axis('off')
savefig("result.png",dpi=300)
show()

代码调试

报错1

Traceback (most recent call last):
  File "F:PJ.py", line 7, in <module>
    from PCV.geometry import homography, warp
  File "F:warp.py", line 1, in <module>
    import matplotlib.delaunay as md 
ModuleNotFoundError: No module named 'matplotlib.delaunay'
[Finished in 4.8s]

解决方法

因为PCV下面的warp.py里面的matplotlib.delaunay不再被使用了，所以把它换成一个相同功能的就可以：
①：把import matplotlib.delaunay as md 换成from scipy.spatial import Delaunay
②：warp.py里面的centers,edges,tri,neighbors = md.delaunay(x,y)换成tri= Delaunay(np.c_[x,y]).simplices就好啦

报错2

runtimewarning：invalid value encountered in divide

解决方法

runtimewarning:在divide中遇到无效值通过网上查阅资料，了解到只要加一句代码:np.seterr(invalid=‘ignore’)? 就可以忽略报错，继续运行得出结果? ?

运行代码拼接结果? ? ? ? ? ?

这是由于照片的排列顺序是要求从右往左的，但是我一开始是从左往右排序，就会导致这样的结果。我们的照片集应该按照从右到左编号，因为匹配是从最右边的照片开始计算。?

数据集1（定点多角度拍摄的室外场景）?
特征匹配结果：?

拼接结果图?

?可以看到，在景深差距小，只有拍摄角度变化下的图片，拼接效果还是很好的，图中操场的线扭曲不大，几乎对齐。不过拼接缝非常明显，这是拍摄照片时的曝光不一致导致的。

数据集2（水平平移拍摄的室外场景）?

特征匹配结果：??

??拼接结果图?

?可以看到拼接效果并不好，这是由于景深不同，拍摄的景物过多，平移拍摄过程中远处建筑群有的被遮挡，有的显现造成的。

数据集3（水平平移拍摄的室内场景）?

特征匹配结果：??

?拼接结果图?

?室内水平平移拍摄的照片的拼接效果非常好，几乎看不到扭曲变形的地方，拼接缝也看不到。

总结

在本次实验过程中得到如下总结：对比固定点位拍摄和视差变化大的场景可以得到结论：固定点位拍照的图像在进行拼接的时候没有发生鬼影现象，但较远处的拼接有略微偏差，而差变化大的场景在进行拼接时，其对于远景的拼接效果极好，细节处的拼接也处理得很到位。两种场景各有优劣。
分别对比室内和室外场景得拼接结果可以得到结论：室外场景的拼接效果较好，其特征点丰富，在拼接时参考性很大；室内场景单一，特征点也较少，拼接的效果也较为不好。
实验中，拼接图像的像素会影响运行速度，所以在拍摄图像上传到电脑后要进行处理，使得像素大小适当，这样对于运行时间会有较大的改变。
代码中有个参数delta。这个参数是针对你拍摄图像时，你相对平移的距离的变量，当你拍摄近景时候，这个参数尽量该小，远景相反。
dpi是图像精细度的变量，可以通过修改来改变图像的分辨率，数值越大表示图像越精细，即分辨率越高。