单应性变换 Homography Estimation

最新推荐文章于 2024-10-15 17:19:27 发布

居愉心

最新推荐文章于 2024-10-15 17:19:27 发布

阅读量2.0k

点赞数 2

CC 4.0 BY-SA版权

分类专栏：结构光三维重建文章标签：计算机视觉

本文链接：https://blog.youkuaiyun.com/weixin_43956164/article/details/119053714

结构光三维重建专栏收录该内容

22 篇文章

订阅专栏

本文介绍如何通过至少4对同名点计算单应性变换矩阵，实现同一平面在不同角度下的图像坐标转换。利用SVD分解求解未知参数。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

应用背景：

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

原理介绍

角度A下图像中一点的坐标 $u_1, v_1)$ , 对应的在角度B下的坐标为 $u_2,v_2)$ , 两个点之间的单应性变换可以表示为: $zc2[u2v21]=[H11H12H13H21H22H23H31H32H33]zc1[u1v11]z_{c2}\begin{bmatrix} u_2\\v_2\\1 \end{bmatrix} = \begin{bmatrix} H_{11}&H_{12}&H_{13}\\ H_{21}&H_{22}&H_{23}\\H_{31}&H_{32}&H_{33}\end{bmatrix}z_{c1}\begin{bmatrix} u_1\\v_1\\1 \end{bmatrix}$ 根据矩阵展开可以得到两个等式，即 $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对以上的点对组合才能解算。
将上式的左边移到右边，可以得到类似 $A x = 0$ 的矩阵形式，利用SVD分解来求解 $x=(H_{11},H_{12},H_{13},H_{21},H_{22},H_{23},H_{31},H_{32},H_{33})$ ，其中 $A=\begin{bmatrix} u_1&v_1&1&0&0&0&-u_1*u_2&-v_1*u_2&-u_2 \\ 0 & 0&0&u_1&v_1&1&-u_1*v_2&-v_1*v_2&-v_2 \\ ...\\...\end{bmatrix}$

代码实现

使用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;