计算机视觉学习blog（四）

您所在的位置：网站首页 › 什么是全景图和半景图的区别呢 › 计算机视觉学习blog（四）

计算机视觉学习blog（四）

#计算机视觉学习blog（四）| 来源: 网络整理| 查看: 265

创建全景图前言一、全景拼接原理1. 提取图像的特征和匹配2. 将匹配转化成齐次坐标3. 估计单应性矩阵（RANSAC算法）4. 拼接图像二、实验结果及分析1. 全代码2. 实现2.1 第一组2.2 第二组2.3 第三组

前言

在同一位置拍摄的两幅或者多幅图像是单应性相关的。我们经常使用该约束将很多图像缝补起来，拼成一个大的图像来创建图像。接下来，我们将探讨如何创建全景图像。

一、全景拼接原理

全景拼接的基础流程如下：

（1）针对同一场景拍摄系列图像

（2）提取图像的特征和匹配

（3）将匹配转化成齐次坐标点

（4）估计单应性矩阵

（5）拼接图像

1. 提取图像的特征和匹配

SIFT算法实现特征匹配主要有三个流程 1、提取关键点； 2、对关键点附加详细的信息（局部特征），即描述符； 3、通过特征点（附带上特征向量的关键点）的两两比较找出相互匹配的若干对特征点，建立景物间的对应关系。

2. 将匹配转化成齐次坐标 3. 估计单应性矩阵（RANSAC算法）

RANSAC 是“RANdom SAmple Consensus”（随机一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC 基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。即从一组数据集（包含噪声点的数据集）中，能够从中挑选出正确的点，获取能够正确拟合的参数模型。

RANSAC的求解过程如下：

（1）随机选择4对匹配特征对应点对（2）根据DLT计算单应性矩阵H（唯一解）（3）对所有匹配点，计算映射误差（对每个对应点使用该单应性矩阵，然后返回平方距离之和）（4）根据误差阈值，确定正常值inliers（小于阈值的点看做正确点，否则认为是错误的点）重复步骤（1）-（4）（5）针对最大inliers集合，重新计算单应性矩阵H

RANSC其实是DLT（Direct Linear Transformation，直接线性变换）的优化。DLT算法是将给定的4个或更多个对应点方程的系数堆叠到一个矩阵中，使用SVD（奇异值分解）算法求得H的最小二乘解。但是，值得注意的是并非所有的对应点对都是正确的

我们使用 RANSAC 算法来求解单应性矩阵，首先需要将下面模型类添加到homography.py文件中：

class RansacModel(object): """ 用于测试单应性矩阵的类，其中单应性矩阵是由网站http://www.scipy.org/Cookbook/RANSAC 上的ransac.py 计算出来的""" 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) # 归一化齐次坐标 # for i in range(3): # fp_transformed[i] /= fp_transformed[2] fp_transformed = normalize(fp_transformed) # 返回每个点的误差 return sqrt(sum((tp - fp_transformed) ** 2, axis=0))

可以看到，这个类包含 fit() 方法。该方法仅仅接受由 ransac.py 选择的４个对应点对（data 中的前 4 个点对），然后拟合一个单应性矩阵。记住，4 个点对是计算单应性矩阵所需的最少数目。由于 get_error() 方法对每个对应点对使用该单应性矩阵，然后返回相应的平方距离之和，因此 RANSAC 算法能够判定哪些点对是正确的，哪些是错误的。在实际中，我们需要在距离上使用一个阈值来决定哪些单应性矩阵是合理的。为了方便使用，将下面的函数添加到 homography.py 文件中：

def H_from_ransac(fp, tp, model, maxiter=1000, match_theshold=10): """ 使用RANSAC 稳健性估计点对应间的单应性矩阵H（ransac.py 为从 http://www.scipy.org/Cookbook/RANSAC 下载的版本） # 输入：齐次坐标表示的点fp，tp（3×n 的数组）""" 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']

该函数允许提供阈值和最小期望的点对数目。最重要的参数是最大迭代次数：程序退出太早可能得到一个坏解；迭代次数太多会占有太多时间。函数的返回结果为单应性矩阵和对应该单应性矩阵的正确点对。

# 将匹配转换成齐次坐标点的函数 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] # im1 到im2 的单应性矩阵 fp, tp = convert_points(0) H_01 = homography.H_from_ransac(fp, tp, model)[0] # im0 到im1 的单应性矩阵 tp, fp = convert_points(2) # 注意：点是反序的 H_32 = homography.H_from_ransac(fp, tp, model)[0] # im3 到im2 的单应性矩阵 tp, fp = convert_points(3) # 注意：点是反序的 H_43 = homography.H_from_ransac(fp, tp, model)[0] # im4 到im3 的单应性矩阵在每个图像对中，由于匹配是从最右边的图像计算出来的，所以我们将对应的顺序进行了颠倒，使其从左边图像开始扭曲。图像编号顺序应逆序编号。 4. 拼接图像

估计出单应性矩阵H后，我们需要将所有图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心平面（否则，需要进行大量变形）。一种方法是创建一个很大的图像，比如将平面中全部填充为0，使其和中心图像平行，然后将所有的图像扭曲到上面。其基本步骤如下：（1）实现像素间的映射（计算像素和和单应性矩阵H相乘，然后对齐次坐标进行归一化）（2）判断图像填补位置（查看H中的平移量，左边为正，右边为负）（3）在单应性矩阵中加入平移量，进行alpha图填充

由于我们所有的图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充 0，以便为扭曲的图像腾出空间。将下面的代码添加到 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) 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

对于通用的 geometric_transform() 函数，我们需要指定能够描述像素到像素间映射的函数。在这个例子中，transf() 函数就是该指定的函数。该函数通过将像素和 H 相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看 H 中的平移量，我们可以决定应该将该图像填补到左边还是右边。当该图像填补到左边时，由于目标图像中点的坐标也变化了，所以在“左边”情况中，需要在单应性矩阵中加入平移。简单起见，我们同样使用0像素的技巧来寻找alpha图。函数如下：

# 扭曲图像 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) 二、实验结果及分析 1. 全代码 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') """ This is the panorama example from section 3.3. """ # 设置数据文件夹的路径 featname = ['D:/JMU/大三下/计算机视觉/实验/全景/pic/pic2/' + str(i + 1) + '.sift' for i in range(5)] imname = ['D:/JMU/大三下/计算机视觉/实验/全景/pic/pic2/' + 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] # im1 到im2 的单应性矩阵 fp, tp = convert_points(0) H_01 = homography.H_from_ransac(fp, tp, model)[0] # im0 到im1 的单应性矩阵 tp, fp = convert_points(2) # 注意：点是反序的 H_32 = homography.H_from_ransac(fp, tp, model)[0] # im3 到im2 的单应性矩阵 tp, fp = convert_points(3) # 注意：点是反序的 H_43 = homography.H_from_ransac(fp, tp, model)[0] # im4 到im3 的单应性矩阵 # 扭曲图像 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() 2. 实现

以下图片按从右到左拍摄顺序编号

2.1 第一组

测试图：

阴天室外