线性求解单应矩阵 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}=Hmathbf{x}_{i}$

 

---

 

方法：直接线性变换DLT算法

我们首先讨论由给定 2D 到 2D 的四组点对应$ x_{i}   leftrightarrow   x_{i}^{'}$,确定H的一种简单的方法是线性算法$x^{'}_{i}=Hx_{i}$，这是一个齐次方程,所以$mathbf{x}_{i}^{'} $和$Hmathbf{x}^{'}$不相等,它们有相同的方向,但是在大小上相差一个非零因子.

使用叉乘表示: $mathbf{x}_{2}^{wedge}Hmathbf{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&1end{pmatrix}^T$  ,     $x_2=begin{pmatrix}u_2&v_2&1end{pmatrix}^T$

 

其中 $x_2^{wedge}=begin{pmatrix}0&-1&v_2\1&0&-u_2\-v_1&u_2&0end{pmatrix}$

叉乘表示: $mathbf{x}_{2}^{'wedge}Hmathbf{x}^{'}_{2}=0=begin{pmatrix}v_2mathbf{-h}^{3T} mathbf{x}_{1}- mathbf{h}^{2T} mathbf{x}_{1}  \mathbf{h}^{1T} mathbf{x}_{1}-u_2mathbf{h}^{3T} mathbf{x}_{1}  \u_2mathbf{h}^{2T} mathbf{x}_{1}-v_2mathbf{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}$ 多于四对,方程$Amathbf{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的最后一列
    }