使用RANSAC算法实现全景拼接
作者：互联网
一、全景拼接的原理
1、RANSAC算法介绍
RANSAC(Random Sample Consensus)即随机采样一致性，该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC的作用就在于，找到正确数据点的同时摒弃噪声点。
2、使用RANSAC算法来求解单应性矩阵
在进行图像拼接时，我们首先要解决的是找到图像之间的匹配的对应点。通常我们采用SIFT算法来实现特征点的自动匹配，SIFT算法的具体内容参照我的另一篇博客：https://blog.csdn.net/sinat_37751993/article/details/88578410 SIFT是具有很强稳健性的描述子，比起图像块相关的Harris角点，它能产生更少的错误的匹配，但仍然还是存在错误的对应点。所以需要用RANSAC算法，对SIFT算法产生的128维特征描述符进行剔除误匹配点。
由直线的知识点可知，两点可以确定一条直线，所以可以随机的在数据点集中选择两点，从而确定一条直线。然后通过设置给定的阈值，计算在直线两旁的符合阈值范围的点，统计点的个数inliers。inliers最多的点集所在的直线，就是我们要选取的最佳直线。
RANSAC算法就是在一原理的基础上，进行的改进，从而根据阈值，剔除错误的匹配点。首先，从已求得的匹配点对中抽取几对匹配点，计算变换矩阵。然后对所有匹配点，计算映射误差。接着根据误差阈值，确定inliers。最后针对最大inliers集合，重新计算单应矩阵 H。
3、拼接图像
在估计出图像的单应性矩阵后，我们需要将所有的图像扭曲到一个公共的图像平面上。通常我们以中心图像平面作为公共平面，然后将中心图像左边或者右边的区域填充为0，以便为扭曲的图像腾出空间。
二、源代码
1、使用SIFT特征检测，匹配对应点
主函数：
# -*- coding: cp936 -*-
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 = ['n'+str(i+1)+'.sift' for i in range(5)]
imname = ['n'+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)
调用sift.py中的相应函数
def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):
""" process an image and save the results in a file"""
if imagename[-3:] != 'pgm':
#create a pgm file
im = Image.open(imagename).convert('L')
im.save('tmp.pgm')
imagename = 'tmp.pgm'
cmmd = str("F:\win64vlfeat\sift.exe "+imagename+" --output="+resultname+
" "+params)
os.system(cmmd)
print 'processed', imagename, 'to', resultname
def read_features_from_file(filename):
""" read feature properties and return in matrix form"""
f = loadtxt(filename)
return f[:,:4],f[:,4:] # feature locations, descriptors
def match(desc1,desc2):
""" for each descriptor in the first image,
select its match in the second image.
input: desc1 (descriptors for the first image),
desc2 (same for second image). """
desc1 = array([d/linalg.norm(d) for d in desc1])
desc2 = array([d/linalg.norm(d) for d in desc2])
dist_ratio = 0.6
desc1_size = desc1.shape
matchscores = zeros((desc1_size[0],1))
desc2t = desc2.T #precompute matrix transpose
for i in range(desc1_size[0]):
dotprods = dot(desc1[i,:],desc2t) #vector of dot products
dotprods = 0.9999*dotprods
#inverse cosine and sort, return index for features in second image
indx = argsort(arccos(dotprods))
#check if nearest neighbor has angle less than dist_ratio times 2nd
if arccos(dotprods)[indx[0]] < dist_ratio * arccos(dotprods)[indx[1]]:
matchscores[i] = int(indx[0])
return matchscores
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
""" show a figure with lines joining the accepted matches
input: im1,im2 (images as arrays), locs1,locs2 (location of features),
matchscores (as output from 'match'), show_below (if images should be shown below). """
im3 = appendimages(im1,im2)
if show_below:
im3 = vstack((im3,im3))
# show image
imshow(im3)
# draw lines for matches
cols1 = im1.shape[1]
for i in range(len(matchscores)):
if matchscores[i] > 0:
plot([locs1[i,0], locs2[matchscores[i,0],0]+cols1], [locs1[i,1], locs2[matchscores[i,0],1]], 'c')
axis('off')
2、使用RANSAC算法来求解单应性矩阵
主函数：
# 将匹配转换成齐次坐标点的函数
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) #注意: 点是反序的
H_32 = homography.H_from_ransac(fp,tp,model)[0] #im 3 to 2
tp,fp = convert_points(3) #注意: 点是反序的
H_43 = homography.H_from_ransac(fp,tp,model)[0] #im 4 to 3
调用homegraphy.py中相应的函数
from numpy import *
from scipy import ndimage
class RansacModel(object):
""" Class for testing homography fit with ransac.py from
http://www.scipy.org/Cookbook/RANSAC"""
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)
# 归一化齐次坐标
fp_transformed = normalize(fp_transformed)
# 返回每个点的误差
return sqrt( sum((tp-fp_transformed)**2,axis=0) )
fit()方法只接受由ransac.py选择的4个对应点对，然后拟合一个单应性矩阵。其中4个点对是计算单应性矩阵所需的最少数目。
get_error()方法对每个对应点对使用该单应性矩阵，然后返回相应的平方距离之和，通过计算其误差，从而判断哪些点的正确的，哪些是错误的。
def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):
""" 使用RANSAC稳健性估计点对应间的单应性矩阵H (ransac.py from
http://www.scipy.org/Cookbook/RANSAC).
input: fp,tp (3*n arrays) points in hom. coordinates. """
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']
在实际中，我们需要在距离上使用一个阈值来决定哪些单应性矩阵是合理的。该函数提供了阈值和最小期望值的点对数目，还有就是最大迭代次数。迭代次数在求解的过程中起着非常重要的作用，若迭代次数太少，程序太早退出，可能还没找到正确的解，若迭代次数太多，使得运行太久，效率降低。该函数返回的结果为单应性矩阵和对应该单应性矩阵的正确点对。
3、扭曲图像
主函数：
# 扭曲图像
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()
调用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))
# 协调后返回（将fromin放置在tomi上）
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
transf()函数通过将像素和H相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看H中的平移量，判断该图像填补到左边还是右边。当该图像填充到左边时，由于目标图像中心点的坐标也变化了，所以在左边的情况中，需要在单应性矩阵中加入平移。
三、针对不同场景做全景拼接（测试图片的拍摄均来自集美大学）
1、室内场景
2、室外景深落差较小的场景
3、室外景深落差较大的场景
四、实验结果分析
通过三种不同场景的测试，可以看出：
室内场景虽然大致上都拼接出来了，但总体的拼接图像在视觉上有点扭曲，且右边的细节处理的不够好，扭曲的比较严重。
室外景深落差较小的场景，拼接效果比较理想，但由于图像曝光度的不同，导致在图像的边界上存在边缘效应，颜色差异明显，这也是该算法需要改进的地方。
室外景深落差较大的场景，远处的嘉庚图书馆拼接的比较完整，但美中不足的就是，近处草坪上的石阶，发生了扭曲。
这边补充一下拍摄测试图像的注意点：为了拼接出效果比较好的图像，在保证有相同匹配点的情况下，拍摄图像的间隔尽可能大一些。且一定要站在同一点，水平移动手机进行拍摄，就像拍摄全景图那样。若人拍摄的位置发生移动的话，算法可能就会因为找不到正确的点对而报错。如下所示：
且最好不要拍摄那种对称的建筑物，且两边的的特征点长的几乎一样的。这样会使算法的匹配出现失误。
最后在给图像编号进行测试时，切记，一定要从右往左进行编号，因为我们的算法的匹配是从最右边的图像计算出来的，代码中有一步骤是将对应的顺序进行颠倒，使其从左边图像开始进行扭曲。
标签：RANSAC,fp,单应性,矩阵,tp,算法,拼接,图像,data 来源： https://blog.csdn.net/sinat_37751993/article/details/88919285