实验三 —— 图像映射与全景拼接

---

前言

图像映射的基本流程

① 针对某个场景拍摄多张/序列图像
② 计算第二张图像与第一张图像之间的变换关
系
③ 将第二张图像叠加到第一张图像的坐标系中
④ 变换后的融合/合成
⑤ 在多图场景中，重复上述过程

---

一、计算变换结构

单应性变换

单应性变换是将一个平面内的点映射到另外一个平面内的二维投影变换。

一般传统方法估计单应性变换矩阵，需要经过以下四个步骤

• 提取每张图SIFT特征
• 提取每个特征点对应的描述子
• 通过匹配特征点描述子，找到两张图中匹配的特征点对（可能存在错误匹配）
• 使用RANSAC算法剔除错误匹配
• 求解方程组，计算Homograph单应性变换矩阵

在这里插入图片描述

参数求解单应性变换

然后应用最小二乘法解最小特征值对应的特征向量

二、全景拼接

1.RANSAC

RANSAC是用来找到正确模型来拟合带有噪声数据的迭代方法，给定一个模型，例如点集之间的单应性矩阵，RANSAC基本思想：数据中包含正确的点和噪声点，合理的模型应该能够描述正确数据点的同时摒弃噪声点。

2.单应性矩阵估计

代码如下（示例）：

# sift特征自动找到匹配对应
import sift
featname = ['Univ' + str(i+1) '.sift' for i in range(5)]
imname = ['Univ' + str(i+1) '.jpg' for i in range(5)] 
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.matches(d[i+1], d[i])

# RANSAC算法来求解单应性矩阵
def __init__(self, debug=False):
	self.debug = debug
	
def fit(self, data):
	""" 计算选取4个单应性矩阵 """
	# 转置，计算单应性矩阵
	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[i] /= dot(H,fp)
	# 归一化齐次坐标
	for i in range(3):
		fp_transformed[i] /= fp_transformed[2]
	# 返回每个点的误差
	return sqrt(sum(tp-fp_transformed)**2, axis = 0))
	
# 估计单应性矩阵
model = homography.RansacModel()
# im0到im1的单应性矩阵
fp, tp = convert_points(0)
H_01 = homography.H_from_ransac(fp,tp,model)[0] 
# im1到im2的单应性矩阵
fp, tp = convert_points(1)
H_12 = homography.H_from_ransac(fp,tp,model)[0]
# im3到im2的单应性矩阵
fp, tp = convert_points(2)
H_32 = homography.H_from_ransac(fp,tp,model)[0]
# im4到im3的单应性矩阵
fp, tp = convert_points(0)
H_43 = homography.H_from_ransac(fp,tp,model)[0]

3. 拼接图像

估计出图像的单应性矩阵（使用RANSAC算法），现在需要拼接图像

def panorama(H, fromin, toim, padding=2400,delta=2400):
	
	# 检查图像是灰度图还是彩色图
	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:
		print('warp - right')
		if is_color:
			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):
				fromin_t[:,:col] = ndimage.geometric_transform(fromin[:,:,col,
							transf,(toim.shape[0], toim.shape[1]+padding))
			else:
				toim_t = hstack((toim,zeros((toim.shape[0],padding))))
				fromim_t = ndimage.geometric_transform(fromin,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)
				if is_color:
				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] = ndiamge.geometric_transform(fromim[:,:,col],
						transf,(toim.shape[0],toim.shape[1]+padding))
				else:
				toim_t = hstack((zeros((toim.shape[0],padding)),toim))
				fromim_t = ndimagegeometric_transform(fromim,transf,(toim.shape[0],toim.shape[1]+padding))
				
				# 协调后返回
				if is_color:
					alpha = ((fromim_t[:,:,0]*fromim_t[:,:,1],*fromim_t[:,:,2])>0
					for col in range(3):
						toim_t[:,:,col] = fromin_t[:,:,1]*alpha + toim_t[:,:,col]*(1-alpha)
				else:
					alpha = (fromim_t>0)
					toim_t = fromim_t*alpha + toim_t*(1-alpha)
				return toim_t