单应性变换 Homography Estimation

应用背景：

已知同一个平面标定板在两个不同角度下的图像，根据图像上对应点的二维（u, v）坐标求解其变换矩阵。单应性变换是两个平面之间的几何变换，利用至少4对以上的同名点计算变换矩阵。

原理介绍

角度A下图像中一点的坐标$(u_1,v_1)$ , 对应的在角度B下的坐标为$(u_2,v_2)$, 两个点之间的单应性变换可以表示为:$$z_{c2}\left[\begin{array}{c}{u}_2\\ v_2\\ 1\end{array}\right]=\left[\begin{array}{ccc}{H}_{11}& H_{12}& H_{13}\\ H_{21}& H_{22}& H_{23}\\ H_{31}& H_{32}& H_{33}\end{array}\right]z_{c1}\left[\begin{array}{c}{u}_1\\ v_1\\ 1\end{array}\right]$$根据矩阵展开可以得到两个等式，即$$u_2(H_{31}{u}_1+H_{32}{v}_1+H_{33})=H_{11}{u}_1+H_{12}{v}_1+H_{13}$$$$v_2(H_{31}{u}_1+H_{32}{v}_1+H_{33})=H_{21}{u}_1+H_{22}{v}_1+H_{23}$$由此可见，1组对应点能提供2个等式。单应性矩阵虽然是3x3共9个参数，但是变量是8个，上式的左右两边可以同时除以一个变量，去除一个相关变量。对于8个变量需要提供4对以上的点对组合才能解算。
将上式的左边移到右边，可以得到类似$Ax=0$的矩阵形式，利用SVD分解来求解$x=(H_{11},H_{12},H_{13},H_{21},H_{22},H_{23},H_{31},H_{32},H_{33})$，其中$$A=\left[\begin{array}{ccccccccc}{u}_1& v_1& 1& 0& 0& 0& -u_1\ast u_2& -v_1\ast u_2& -u_2\\ 0& 0& 0& u_1& v_1& 1& -u_1\ast v_2& -v_1\ast v_2& -v_2\\ ...& & & & & & & & \\ ...& & & & & & & & \end{array}\right]$$

代码实现

使用EIGEN库的主要计算代码如下，

// 计算单应性矩阵
int nPairSize = vCodedPairPt.size();
if (nPairSize < 4)
{
	cout << "匹配对数少于4对" << endl;
	return 0;
}
MatrixXd MatrixA(2 * nPairSize, 9);
double u1, v1, u2, v2;
int index = 0;
for (int i = 0; i < nPairSize; i++)
{
	u1 = vCodedPairPt[i].first->pt.dX;
	v1 = vCodedPairPt[i].first->pt.dY;
	u2 = vCodedPairPt[i].second->pt.dX;
	v2 = vCodedPairPt[i].second->pt.dY;
	index = i * 2;
	MatrixA(index, 0) = u1,     	 MatrixA(index, 1) = v1, 		  MatrixA(index, 2) = 1;
	MatrixA(index, 3) = 0,  		 MatrixA(index, 4) = 0, 		  MatrixA(index, 5) = 0;
	MatrixA(index, 6) = -u1*u2,      MatrixA(index, 7) = -v1*u2,      MatrixA(index, 8) = -u2;
	MatrixA(index + 1, 0) = 0,       MatrixA(index + 1, 1) = 0,       MatrixA(index + 1, 2) = 0;
	MatrixA(index + 1, 3) = u1,      MatrixA(index + 1, 4) = v1,      MatrixA(index + 1, 5) = 1;
	MatrixA(index + 1, 6) = -u1*v2,  MatrixA(index + 1, 7) = -v1*v2,  MatrixA(index + 1, 8) = -v2;
}
//cout << "MatrixA: " << MatrixA << endl;
JacobiSVD<MatrixXd> svd(MatrixA, ComputeThinU | ComputeThinV);
MatrixXd MatrixU;
MatrixU = svd.matrixV();
int nCols = MatrixU.cols();
VectorXd vSolution = MatrixU.col(nCols - 1);
Matrix3d MatrixH = Map<Matrix3d>(vSolution.data()).transpose();
cout << "MatrixH: " << MatrixH << endl;