图像配准（一）——仿射（实战学习记录）

参考内容：ImageRegistration

				[RANSAC算法——看完保证你理解](https://blog.csdn.net/zhoucoolqi/article/details/105497572)

RANSAC：随机抽样一致算法
采用迭代的方式从一组包含离群的被观测数据中估算出数学模型的参数
相较于最小二乘算法：其融合了剔除不合格数据的思想；
通用的RANSAC算法的工作流程如下：
输入：
data – 一组观测数据组.
model – 拟合模型（例如线性、二次曲线等等）.
n – 用于拟合的最小数据组数.
k – 算法规定的最大遍历次数.
t – 数据和模型匹配程度的阈值，在t范围内即inliers，在范围外即outliers.
d – 表示模型合适的最小数据组数.
返回：
bestFit – 一组最匹配的模型参数，即model的参数

待配准图像如下：

目标图像：

代码详解学习：

# 读取图像
source_path = 'mona_source.png'
target_path = 'mona_target.jpg'
img_source = cv2.imread(source_path, 1)
img_target = cv2.imread(target_path, 1)

# 提取SIFT关键点与描述符
def extract_SIFT(img):
	img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
	sift = cv2.SIFT_create()
	kp, desc = sift.detectAndCompute(img_gray, None)
	kp = np.array([p.pt for p in kp]).T
	return kp, desc

# 提取两张图像中的关键点与SIFT描述符
kp_s, desc_s = extract_SIFT(img_source)
kp_t, desc_t = extract_SIFT(img_target)

# 匹配源图像和目标图像的SIFT描述符
def match_SIFT(desc_s, desc_t):
	# 匹配描述符并获得两个最佳匹配
	bf = cv2.BFMatcher()
	matches = bf.knnMatch(desc_s, desc_t, k = 2)
	fit_pos = np.array([], dtype=np.int32).reshape((0, 2))

	matches_num = len(matches)
	for i in range(matches_num):
		# 如果比例id小于0.8，则获取良好匹配
		if matches[i][0].distance <= 0.8 * matches[i][1].distance:
			temp = np.array([matches[i][0].queryIdx, matches[i][0].trainIdx])
			fit_pos = np.vstack((fit_pos, temp))
	return fit_pos

# 获得最佳匹配
fit_pos = match_SIFT(desc_s, desc_t)

# 由对应点计算仿射变换矩阵
def affine_matrix(kp_s, kp_t, fit_pos):
	# 从所有关键点中提取相应的点
	kp_s = kp_s[:, fit_pos[:, 0]]
	kp_t = kp_t[:, fit_pos[:, 1]]
	# 应用RANSAC查找大多数内层
	_,_, inliers = ransac_fit(kp_s, kp_t)
	
	# 从所有关键点提取所有内点
	kp_s = kp_s[:, inliers[0]]
	kp_t = kp_t[:, inliers[0]]

	# 使用所有的内点来估计变换矩阵
	A, t = estimate_affine(kp_s, kp_t)
	M = np.hstack((A, t))
	return M

在affine_matrix函数中涉及到了ransac_fit与estimate_affine函数，进行解释为：

# 
def ransac_fit(pts_s, pts_t):
	# 初始化内部器的数量
	inliers_num = 0
	# 初始化仿射变换A和t，和一个存储内部容器索引的向量
	A = None
	t = None
	inliers = None
	# 内点数量设置为2000
	for i in range(2000):
		idx = np.random.randint(0, pts_s.shape[1],(3, 1))
		# 通过这些点来估计仿射变换
		A_tmp, t_tmp = estimate_affine(pts_s[:, idx], pts_t[:, idx])

		# 应用估计的变换计算残差
		residual = residual_lengths(A_tmp, t_tmp, pts_s, pts_t)

		if not(residual is None):
			inliers_tmp = np.where(resudual < 1)
			# 获取内层数
			inliers_num_tmp = len(inliers_tmp[0])
			# 设置仿射变换和索引
			# 在一个迭代中有最多的内部器
			if inliers_num_tmp > inliers_num:
				inliers_num = inliers_num_tmp
				inliers = inliers_tmp
				A = A_tmp
				t = t_tmp
		else:
			pass
	return A, t, inliers

# 通过给定的点估计仿射变换
def estimate_affine(pts_s, pts_t):
	# 获取对应点的个数
	pts_num = pts_s.shape[1]
	# 初始化矩阵M，由于仿射变换，M有6列
	M = np.zeros((2 * pts_num, 6))
	for i in range(pts_num):
		temp = [[pts_s[0, i], pts_s[1, i], 0, 0, 1, 0],
				[0, 0, pts_s[0, i], pts_s[1, i], 0, 1]]
		M[2 * i:2*i+2, :] = np.array(temp, dtype = object)

	# 形成矩阵b，b包含所有已知的目标点
	b = pts_t.T.reshape((2 * pts_num, 1))
	try:
		theta = np.linalg.lstsq(M, b, recond=None)[0]
	
		# 形成仿射变换
		A = thera[:4].reshape((2, 2))
		t= theta[4:]
	except np.linalg.linalg.LinAlgError:
		A = None
		t = None
	return A, t

继续进行下面的步骤

# 计算仿射变换矩阵M ：本例中维度 2*3
M = affine_matrix(kp_s, kp_t, fit_pos)

# 用仿射将源图像变形为目标图像
def warp_image(source, target, M):
	rows, cols, _ = target.shape
	# 扭曲源图像
	warp = cv2.warpAffine(source, M, (cols, rows))
	# 将扭曲图像与要现实的目标图像合并
	merge = np.uint8(target * 0.5 + warp * 0.5)
	# 展示图像
	cv2.inshow('img', merge)
	cv2.waitKey(0)
	cv2.destroyAllWindows()
	return

# 实现变形并显示
warp_image(img_source, img_target, M)