基于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。相比之下,第一个算法两图拼接就对图片要求低一些。