线性求解单应矩阵 Homography
定义：
2D单应：给定图像$\mathbb{P}^{2}$中的特征点集$\mathbf{x}_i$和另一幅图像在$\mathbb{P}^{2}$ 中对应的特征点集$\mathbf{x}_{i}^{'}$， 将$\mathbf{x}_i$映射到$\mathbf{x}^{'}_{i}$的射影变换。在实际情况中，点$\mathbf{x}_{i}$和$\mathbf{x}^{'}_{i}$是两幅图像上的点，每幅图像都视为一张射影平面$\mathbb{P}^{2}$
$\mathbf{x}^{'}_{i}=H\mathbf{x}_{i}$
方法：直接线性变换DLT算法
我们首先讨论由给定 2D 到 2D 的四组点对应$ x_{i} \leftrightarrow x_{i}^{'}$,确定H的一种简单的方法是线性算法$x^{'}_{i}=Hx_{i}$，这是一个齐次方程,所以$\mathbf{x}_{i}^{'} $和$H\mathbf{x}^{'}$不相等,它们有相同的方向,但是在大小上相差一个非零因子.
使用叉乘表示: $\mathbf{x}_{2}^{\wedge}H\mathbf{x}_{1}=0$
将矩阵H的第j行记为$\mathbf{h}^{jT}$得
$$H \mathbf{x}_1=\begin{pmatrix}\mathbf{h}^{1T} \mathbf{x}_{1} \\ \mathbf{h}^{2T} \mathbf{x}_{1} \\ \mathbf{h}^{3T} \mathbf{x}_{1} \\ \end{pmatrix}$$
$\mathbf{h}$ 默认是列排列
记
$x_1=\begin{pmatrix}u_1&v_1&1\end{pmatrix}^T$ , $x_2=\begin{pmatrix}u_2&v_2&1\end{pmatrix}^T$
其中 $x_2^{\wedge}=\begin{pmatrix}0&-1&v_2\\1&0&-u_2\\-v_1&u_2&0\end{pmatrix}$
叉乘表示: $\mathbf{x}_{2}^{'\wedge}H\mathbf{x}^{'}_{2}=0=\begin{pmatrix}v_2\mathbf{-h}^{3T} \mathbf{x}_{1}- \mathbf{h}^{2T} \mathbf{x}_{1} \\\mathbf{h}^{1T} \mathbf{x}_{1}-u_2\mathbf{h}^{3T} \mathbf{x}_{1} \\u_2\mathbf{h}^{2T} \mathbf{x}_{1}-v_2\mathbf{h}^{1T} \mathbf{x}_{1} \\ \end{pmatrix}$
将h 抽出为八行向量得 A$\bf{h}$=0;
$$\begin{pmatrix}\mathbf{0}^{T}&-\mathbf{x}^{T}_{1}&u_{2}\mathbf{x}^{T}_{1} \\\mathbf{x}^{T}_{1} &\mathbf{0}^{T}&-v_{2}\mathbf{x}^{T}_{1}\\v_{2}\mathbf{x}^{T}_{1}&u_{2}\mathbf{x}^{T}_{1}&\mathbf{0}^{T}\\\end{pmatrix} \begin{pmatrix} \mathbf{h}^{1} \\\mathbf{h}^{2} \\\mathbf{h}^{3} \end{pmatrix}=0$$
上式虽然有三个方程但是只有两个是线性无关的,所以每对点可以取出两个方程;
$$\begin{pmatrix}\mathbf{0}^{T}&-\mathbf{x}^{T}_{1}&u_{2}\mathbf{x}^{T}_{1} \\\mathbf{x}^{T}_{1} &\mathbf{0}^{T}&-v_{2}\mathbf{x}^{T}_{1}\end{pmatrix} \begin{pmatrix} \mathbf{h}^{1} \\\mathbf{h}^{2} \\\mathbf{h}^{3} \end{pmatrix}=0$$
A是2x9矩阵,四组点后A为8x9
求解H
H是自由度为8的矩阵. 每对点得到两个关于H元素的两个线性无关的方程,给定四组点得到八个方程即可得到H. 由于A为8x9秩为八,因此A仅有一维零空间,从而存在只相差一个非零因子意义下的解h.但是单应矩阵H只能够确定到相差一个尺度,因此可以通过范数来对H元素进行选择 如||$\mathbf{h}$||=1.
超定解
若给出的点$\mathbf{x}_{i}$和$\mathbf{x}^{'}_{i}$ 多于四对,方程$A\mathbf{h}=0&是超定的.
(1) 如果不存在噪声,那么A的秩为八且有一维零空间,并且存在精确解$\bf{h}$.
(2)如果存在噪声,那么方程除零解外不存在精确解,但是可以试图寻找近似解,该解是$A^{T}A$的最小特征值对应的单位特征向量,即该解是A最小奇异值的单位向量.
代码
/** * @brief 从特征点匹配求homography（normalized DLT） * * @param vP1 归一化后的点, in reference frame * @param vP2 归一化后的点, in current frame * @return 单应矩阵 * @see Multiple View Geometry in Computer Vision - Algorithm 4.2 p109 */ cv::Mat Initializer::ComputeH21(const vector<cv::Point2f> &vP1, const vector<cv::Point2f> &vP2) { const int N = vP1.size(); cv::Mat A(2 * N, 9, CV_32F); // 2N*9 N=8 for (int i = 0; i < N; i++) //16个方程 ======== 14讲上用8个方程也可以 ===前面找到了16个点顺手一起用了 // ==用16个有部分是线性相关的 超过4组点 Ah=0 是超定的 // 存在噪声 超定方程组不存在精确解 可以找到近似解 # h除零外 # h为 使Ah 最小化 // 不存在噪声 超定方程组只有零解 { const float u1 = vP1[i].x; const float v1 = vP1[i].y; const float u2 = vP2[i].x; const float v2 = vP2[i].y; //第一行 A.at<float>(2 * i, 0) = 0.0; A.at<float>(2 * i, 1) = 0.0; A.at<float>(2 * i, 2) = 0.0; A.at<float>(2 * i, 3) = -u1; A.at<float>(2 * i, 4) = -v1; A.at<float>(2 * i, 5) = -1; A.at<float>(2 * i, 6) = v2 * u1; A.at<float>(2 * i, 7) = v2 * v1; A.at<float>(2 * i, 8) = v2; //第二行 A.at<float>(2 * i + 1, 0) = u1; A.at<float>(2 * i + 1, 1) = v1; A.at<float>(2 * i + 1, 2) = 1; A.at<float>(2 * i + 1, 3) = 0.0; A.at<float>(2 * i + 1, 4) = 0.0; A.at<float>(2 * i + 1, 5) = 0.0; A.at<float>(2 * i + 1, 6) = -u2 * u1; A.at<float>(2 * i + 1, 7) = -u2 * v1; A.at<float>(2 * i + 1, 8) = -u2; } cv::Mat u, w, vt; cv::SVDecomp(A, w, u, vt, cv::SVD::MODIFY_A | cv::SVD::FULL_UV); return vt.row(8).reshape(0, 3); // vt的最后一行 即v的最后一列 }
