定義：
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的最后一列 }