基于OpenCV的全景拼接

背景介绍

在同一位置拍摄的两幅或多幅图像是单应性相关的。我们可以使用该约束将很多图像拼接起来，拼成一幅大的图像来创建全景图像。其步骤总结起来就两个步骤：
1.利用sift算法找出两种图片的相似点，计算变换矩阵（单应性矩阵）。
2.变换一张图片到另一种图片上合适的位置，并重新计算重叠区域的像素值。

基本原理

1.单应性矩阵

定义：在计算机视觉领域，空间同一平面的任意两幅图像被单应矩阵联系着（假设在针孔相机模型中），即一个相机拍摄空间同一平面的两张图像，这两张图像之间的映射关系可以用单应矩阵表示。
在两视几何中，也可以这样理解，两架相机拍同一空间上得到两幅图像A、B，其中图像A到图像B存在一种变换，而且这种变换是一一对应的关系，这个变换矩阵用单应矩阵表示。OpenCV中可以用函数findHomography计算得到单应矩阵H。
要实现两张图片的简单拼接，其实只需找出两张图片中相似的点 (至少四个，因为 homography 矩阵的计算需要至少四个点)，计算一张图片可以变换到另一张图片的变换矩阵 (homography 单应性矩阵)，用这个矩阵把那张图片变换后放到另一张图片相应的位置 ( 就是相当于把两张图片中定好的四个相似的点給重合在一起)。如此，就可以实现简单的全景拼接。当然，因为拼合之后图片会重叠在一起，所以需要重新计算图片重叠部分的像素值，否则结果会很难看。

2.RANSAC算法

RANSAC(Random Sample Consensus)即随机采样一致性，该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC的基本思想在于，找到正确数据点的同时摒弃噪声点。

3.利用RANSAC算法求解单应性矩阵

虽然SIFT是具有很强稳健性的描述子，当这方法仍远非完美，还会存在一些错误的匹配。而单应性矩阵需要选取4对特征点计算，万一选中了不正确的匹配点，那么计算的单应性矩阵肯定是不正确的。因此，为了提高计算结果的鲁棒性，我们下一步就是要把这些不正确的匹配点给剔除掉，获得正确的单应性矩阵。在这里使用了RANSAC算法：随机抽取不同的4对特征匹配坐标（在图1中随机抽取4个特征坐标，以及这4个特征坐标在图2中匹配的4个特征坐标，组成4对特征匹配坐标），利用这4对特征匹配坐标计算出矩阵H1（3x3的一个矩阵，图2经过矩阵变换后，可以把图2映射到图1的坐标空间中，再将图2进行简单的平移即可与图1实现无缝拼接），再将图2中所有特征匹配点经过该透视矩阵H1映射到图1的坐标空间，然后与图1匹配点实际坐标求欧氏距离（就是为了验证计算出来的这个H1矩阵是否满足绝大多数特征匹配点）；之后重复上面内容，再随机抽取不同的四组特征匹配坐标，再计算透视矩阵H2，再求欧式距离，如此重复多次。最后以欧式距离最小的那个透视矩阵（表示这个特征矩阵H满足最多的特征匹配点，它最优秀）作为最终计算结果。

4.图片融合

在用计算出的变换矩阵对其中一张图做变换，然后把变换的图片与另一张图片重叠在一起，并重新计算重叠区域新的像素值。对于计算重叠区域的像素值，其实可以有多种方法去实现一个好的融合效果，这里就用最简单粗暴的但效果也不错的方式。直白来说就是实现一个图像的线性渐变，对于重叠的区域，靠近左边的部分，让左边图像内容显示的多一些，靠近右边的部分，让右边图像的内容显示的多一些。用公式表示就是，假设 alpha 表示像素点横坐标到左右重叠区域边界横坐标的距离，新的像素值就为 newpixel = 左图像素值 × (1 - alpha) + 右图像素值 × alpha 。这样就可以实现一个简单的融合效果。

代码及运行结果

1.两图拼接有融合

该代码实现两种图片的拼接，对拼接处有模糊处理。

import numpy as np
import cv2 as cv
from matplotlib import pyplot as pltif __name__ == '__main__':top, bot, left, right = 100, 100, 0, 500img1 = cv.imread('C:\\Users\DELL\Desktop\PCV\jmu\lib/lib1.jpg')img2 = cv.imread('C:\\Users\DELL\Desktop\PCV\jmu\lib/lib2.jpg')srcImg = cv.copyMakeBorder(img1, top, bot, left, right, cv.BORDER_CONSTANT, value=(0, 0, 0))testImg = cv.copyMakeBorder(img2, top, bot, left, right, cv.BORDER_CONSTANT, value=(0, 0, 0))img1gray = cv.cvtColor(srcImg, cv.COLOR_BGR2GRAY)img2gray = cv.cvtColor(testImg, cv.COLOR_BGR2GRAY)sift = cv.xfeatures2d_SIFT().create()# find the keypoints and descriptors with SIFTkp1, des1 = sift.detectAndCompute(img1gray, None)kp2, des2 = sift.detectAndCompute(img2gray, None)# FLANN parametersFLANN_INDEX_KDTREE = 1index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)search_params = dict(checks=50)flann = cv.FlannBasedMatcher(index_params, search_params)matches = flann.knnMatch(des1, des2, k=2)# Need to draw only good matches, so create a maskmatchesMask = [[0, 0] for i in range(len(matches))]good = []pts1 = []pts2 = []# ratio test as per Lowe's paperfor i, (m, n) in enumerate(matches):if m.distance < 0.7*n.distance:good.append(m)pts2.append(kp2[m.trainIdx].pt)pts1.append(kp1[m.queryIdx].pt)matchesMask[i] = [1, 0]draw_params = dict(matchColor=(0, 255, 0),singlePointColor=(255, 0, 0),matchesMask=matchesMask,flags=0)img3 = cv.drawMatchesKnn(img1gray, kp1, img2gray, kp2, matches, None, **draw_params)# plt.imshow(img3, ), plt.show()rows, cols = srcImg.shape[:2]MIN_MATCH_COUNT = 10if len(good) > MIN_MATCH_COUNT:src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)dst_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)M, mask = cv.findHomography(src_pts, dst_pts, cv.RANSAC, 5.0)warpImg = cv.warpPerspective(testImg, np.array(M), (testImg.shape[1], testImg.shape[0]), flags=cv.WARP_INVERSE_MAP)for col in range(0, cols):if srcImg[:, col].any() and warpImg[:, col].any():left = colbreakfor col in range(cols-1, 0, -1):if srcImg[:, col].any() and warpImg[:, col].any():right = colbreakres = np.zeros([rows, cols, 3], np.uint8)for row in range(0, rows):for col in range(0, cols):if not srcImg[row, col].any():res[row, col] = warpImg[row, col]elif not warpImg[row, col].any():res[row, col] = srcImg[row, col]else:srcImgLen = float(abs(col - left))testImgLen = float(abs(col - right))alpha = srcImgLen / (srcImgLen + testImgLen)res[row, col] = np.clip(srcImg[row, col] * (1-alpha) + warpImg[row, col] * alpha, 0, 255)# opencv is bgr, matplotlib is rgbres = cv.cvtColor(res, cv.COLOR_BGR2RGB)# show the resultplt.figure()plt.imshow(res)plt.show()else:print("Not enough matches are found - {}/{}".format(len(good), MIN_MATCH_COUNT))matchesMask = None

运行结果

可以看出无论是在室内场景还是在室外，都能很好地进行拼接，模糊处理也消除了拼接边界。唯一不足的是无法进行多图拼接，下面介绍多图拼接。

多图拼接

基于RANSAC的实现方法

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 = ['C:\\Users\DELL\Desktop\PCV\jmu\panorama/z0'+str(i+1)+'.sift' for i in range(5)]
imname = ['C:\\Users\DELL\Desktop\PCV\jmu\panorama/z0'+str(i+1)+'.jpg' for i in range(5)]# extract features and match
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])# visualize the matches (Figure 3-11 in the book)
'''
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)
'''# 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 elsewherefp = 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(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    # warp the images
delta = 500 # for padding and translationim1 = 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()

调用warp.py中相关的函数进行单应性变换

def panorama(H,fromim,toim,padding=2400,delta=2400):""" Create horizontal panorama by blending two images using a homography H (preferably estimated using RANSAC).The result is an image with the same height as toim. 'padding' specifies number of fill pixels and 'delta' additional translation. """ # check if images are grayscale or coloris_color = len(fromim.shape) == 3# homography transformation for 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 is to the rightprint('warp - right')# transform fromimif is_color:# pad the destination image with zeros to the righttoim_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:# pad the destination image with zeros to the righttoim_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')# add translation to compensate for padding to the leftH_delta = array([[1,0,0],[0,1,-delta],[0,0,1]])H = dot(H,H_delta)# transform fromimif is_color:# pad the destination image with zeros to the lefttoim_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:# pad the destination image with zeros to the lefttoim_t = hstack((zeros((toim.shape[0],padding)),toim))fromim_t = ndimage.geometric_transform(fromim,transf,(toim.shape[0],toim.shape[1]+padding))# blend and return (put fromim above toim)if is_color:# all non black pixelsalpha = ((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

transf()函数通过将像素和H相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看H中的平移量，判断该图像填补到左边还是右边。当该图像填充到左边时，由于目标图像中心点的坐标也变化了，所以在左边的情况中，需要在单应性矩阵中加入平移。

运行结果

总的全景图，除了主建筑有稍微偏差，其他景观拼接无明显差异。

左边放大部分，拼接效果还行，不过中间有条缝还是能看见：

右边放大部分，拼接效果一般，有2处明显缝隙：

再拼接室内图片，发现效果很差，估计是平衡值没调好，但调了很久也没见多大改善，只好将就。拼接图缝隙明显，角度也不是很对，估计单应性变换那边出了问题，有待解决。

总结

室外场景的拼接效果比室内的要好很多，多图拼接对图片的要求也比较高，差异性大或太小（几乎相同）的拼接效果都很差。而且如果拍摄地点变动过大，还可能报错 ValueError: did not meet fit acceptance criteria，匹配值为0。相比之下，第一个算法两图拼接就对图片要求低一些。