计算机视觉 - python实现求解homography矩阵

计算机视觉 - python实现求解homography矩阵

要求

给定$P_i\equiv [x_i,y_i,1]$ 和$P_i^{{}^{\prime }}\equiv [x_i^{{}^{\prime }},y_i^{{}^{\prime }},1]$，需要计算一个齐次变换矩阵$H\in R^3$满足$P_i^{{}^{\prime }}\equiv H{P}_i$。

而大部分的数据集不能准确地计算出$H$，所以可以计算出最优的$H$

计算方法如下：$\underset{||H|{|}_2^2=1}{\mathrm{argmin}}||Ah||,h\in R^9,A\in R^{2N\times 9}$

使$Ah$的模最小的$h$的参数则为所求的$H$的九个参数

分析与求解

首先，我们来看看这里的$A$是什么，
$$A=\left[\begin{array}{ccccccccc}-x_1& -y_1& -1& 0& 0& 0& x_1{x}_1^{{}^{\prime }}& y_1{x}_1^{{}^{\prime }}& x_1^{{}^{\prime }}\\ 0& 0& 0& -x_1& -y_1& -1& x_1{y}_1^{{}^{\prime }}& y_1{y}_1^{{}^{\prime }}& y_1^{{}^{\prime }}\\ -x_2& -y_2& -1& 0& 0& 0& x_2{x}_2^{{}^{\prime }}& y_2{x}_2^{{}^{\prime }}& x_2^{{}^{\prime }}\\ 0& 0& 0& -x_2& -y_2& -1& x_2{y}_2^{{}^{\prime }}& y_2{y}_2^{{}^{\prime }}& y_2^{{}^{\prime }}\\ -x_3& -y_3& -1& 0& 0& 0& x_3{x}_3^{{}^{\prime }}& y_3{x}_3^{{}^{\prime }}& x_3^{{}^{\prime }}\\ 0& 0& 0& -x_3& -y_3& -1& x_3{y}_3^{{}^{\prime }}& y_3{y}_3^{{}^{\prime }}& y_3^{{}^{\prime }}\\ -x_4& -y_4& -1& 0& 0& 0& x_4{x}_4^{{}^{\prime }}& y_4{x}_4^{{}^{\prime }}& x_4^{{}^{\prime }}\\ 0& 0& 0& -x_4& -y_4& -1& x_4{y}_4^{{}^{\prime }}& y_4{y}_4^{{}^{\prime }}& y_4^{{}^{\prime }}\\ ...& & & & & & & & \end{array}\right]$$
矩阵$A$就是将数据点按照上述规则进行构造的矩阵

我们再来看看h
$$h={\left[\begin{array}{cccccccccc}{H}_{00}& H_{01}& H_{02}& H_{10}& H_{11}& H_{12}& H_{20}& H_{21}& H_{22}& \end{array}\right]}^T$$
这样，问题就转化为了齐次方程线性最小二乘的求解问题，我们期望找到$Ah=0$中$h$不等于零的解。由于该方程的特殊形式我们发现：对于$h$不等于零的解，我们乘上任意的尺度因子$k$，使解变为$kh$，都不会改变最终结果，因此我们可以将问题转化为求解$h$，使得$||Ah||$值最小并且$||h||=1$。

我们可以采取奇异值分解的方法，对A进行SVD分解有：
$$||Ah||=||UD{V}^{T}h||=||D{V}^{T}h||$$
令$||h||=||V^{T}h||$，可将问题进一步转化为，最小化$||D{V}^{T}h||$，且$||V^{T}h||=1$。

又令$y=V^{T}h$，有：$||Dy||$ 且 $||y||=1$

对于当前问题，A列满秩，则D是对角阵，V是方阵

在奇异值分解中D的对角线元素是递减排列的，那么只需取$=(0,0,0,...,0,1)$，则$||Ah||=||Dy||={\sigma }_n$，${\sigma }_n$是A最小的奇异值

此时：$h=Vy$，

即$h$的值为矩阵$V$的最后一列，即为矩阵$V^T$的最后一行。

这样，我们就可以解决问题了：先构造出A，然后调用scipy库中的svd函数进行奇异值分解，然后利用h构造出H即可，代码如下：

注：

1. 最好不要用numpy库里的svd函数，当数据量较大时，会显示init_dgesdd failed init报错
2. scipy库里的svd函数的第三个返回值为$V^T$，而不是$V$，因此是取它的最后一行，而不是最后一列
3. 可视化时，注意是H左乘齐次坐标$[x_i,y_i,1]$

具体代码如下：

import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as linalg


def fit_homography(XY):
	'''
	Given a set of N correspondences XY of the form [x,y,x',y'],
	fit a homography from [x,y,1] to [x',y',1].
	
	Input - XY: an array with size(N,4), each row contains two
			points in the form [x_i, y_i, x'_i, y'_i] (1,4)
	Output -H: a (3,3) homography matrix that (if the correspondences can be
			described by a homography) satisfies [x',y',1]^T === H [x,y,1]^T

	'''

	A = np.zeros(shape=(XY.shape[0] * 2,9),dtype=np.float32)
	for i in range(XY.shape[0]):
		A[2 * i,0] = -1.0 * XY[i,0]
		A[2 * i,1] = -1.0 * XY[i,1]
		A[2 * i,2] = -1
		A[2 * i,6] = XY[i,0] * XY[i,2]
		A[2 * i,7] = XY[i,1] * XY[i,2]
		A[2 * i,8] = XY[i,2]
		A[2 * i + 1,3] = -1.0 * XY[i,0]
		A[2 * i + 1,4] = -1.0 * XY[i,1]
		A[2 * i + 1,5] = -1
		A[2 * i + 1,6] = XY[i,0] * XY[i,3]
		A[2 * i + 1,7] = XY[i,1] * XY[i,3]
		A[2 * i + 1,8] = XY[i,3]
	_, _, vt = linalg.svd(A)
	H = vt[-1].reshape(3,3) # 注意是-1行，不是-1列
	H = H / H[2,2] # 第三行第三列为1
	return H 

if __name__ == "__main__":
    file = 'file.npy'
	data = np.load(file) # 数据每一行为[x,y,x',y']

	plt.clf()
	plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
	plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
	xy = data[:,0:2]
	xy_pie = data[:,2:]
    # 添加一列1,转化为齐次坐标
    A = np.c_[xy,np.ones(len(xy))] 
	H = fit_homography(data)
	print(H)
	xy_pred = np.dot(H, A.T).T
    # 转化为二维坐标
	xy_pred[:,0] = xy_pred[:,0] / xy_pred[:,2]
	xy_pred[:,1] = xy_pred[:,1] / xy_pred[:,2]
    plt.scatter(xy[:,0], xy[:,1], 1, color='orange', label='x,y')
	plt.scatter(xy_pie[:,0], xy_pie[:,1], 1, color='blue', label='x\',y\'')
	plt.scatter(xy_pred[:,0], xy_pred[:,1], 1, color='red', label='x_pred,y_pred')
	plt.legend()
	plt.show()