基于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 plt

if __name__ == '__main__':
    top, bot, left, right = 100, 100, 0, 500
    img1 = 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 SIFT
    kp1, des1 = sift.detectAndCompute(img1gray, None)
    kp2, des2 = sift.detectAndCompute(img2gray, None)
    # FLANN parameters
    FLANN_INDEX_KDTREE = 1
    index_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 mask
    matchesMask = [[0, 0] for i in range(len(matches))]

    good = []
    pts1 = []
    pts2 = []
    # ratio test as per Lowe's paper
    for 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 = 10
    if 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 = col
                break
        for col in range(cols-1, 0, -1):
            if srcImg[:, col].any() and warpImg[:, col].any():
                right = col
                break

        res = 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 rgb
        res = cv.cvtColor(res, cv.COLOR_BGR2RGB)
        # show the result
        plt.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 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(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 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()

调用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 color
    is_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 right
        print('warp - right')
        # transform fromim
        if is_color:
            # pad the destination image with zeros to the right
            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:
            # pad the destination image with zeros to the right
            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')
        # add translation to compensate for padding to the left
        H_delta = array([[1,0,0],[0,1,-delta],[0,0,1]])
        H = dot(H,H_delta)
        # transform fromim
        if is_color:
            # pad the destination image with zeros to the left
            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:
            # pad the destination image with zeros to the left
            toim_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 pixels
        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

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

运行结果

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

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

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

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

总结

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