[计算机视觉多视图几何] -- Homography

Homography求解与基于eigen库的c++实现

---

---

前言

在使用标志物识别与跟踪技术中，利用homography来描述两个视图中标志物平面之间的映射关系。

一、理论推导

当两个视图中拥有相同的一个平面时，可以根据平面来估计视图之间的映射关系，并称之为单应性矩阵。

在世界中一个点$P[X,Y,Z]$，在两张图像中的位置分别为$p_1,p_2$，则有下列关系：
$${\mathbf{n}}^T\mathbf{P}+d=0$$

$$-\frac{{\mathbf{n}}^T\mathbf{P}}{d}=1$$

同一点在两张视图中的映射关系：
$$\mathbf{p}1=K(\mathbf{P})$$

$${\mathbf{p}}_2=K(\mathbf{R}\mathbf{P}+\mathbf{t})$$

$${\mathbf{p}}_2=K(\mathbf{R}\mathbf{P}-\mathbf{t}\frac{{\mathbf{n}}^T\mathbf{P}}{d})$$

$${\mathbf{p}}_2=\mathbf{K}(\mathbf{R}-\mathbf{t}\frac{{\mathbf{n}}^T}{d})\mathbf{P}$$

$${\mathbf{p}}_2\sim \mathbf{K}(\mathbf{R}-\mathbf{t}\frac{{\mathbf{n}}^T}{d}){\mathbf{K}}^{-1}{\mathbf{p}}_1$$

所以，化简得到点的映射关系：
$${\mathbf{p}}_2\sim \mathbf{H}{\mathbf{p}}_1$$

$$\left(\begin{array}{c}{u}_2\\ v2\\ 1\end{array}\right)\sim \left(\begin{array}{ccc}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ h_7& h_8& h_9\end{array}\right)\left(\begin{array}{c}{u}_1\\ v_1\\ 1\end{array}\right)$$

$$u_2\sim h_1{u}_1+h_2{v}_1+h_3$$

$$v_2\sim h_4{u}_1+h_5{v}_1+h_6$$

$$1\sim h_7{u}_1+h_8{v}_1+h_9$$

由于缩放因子的存在，添加一个约束，令$h_9=1$，将矩阵展开得：
$$h_1{u}_1+h_2{v}_1+h_3-h_7{u}_1{u}_2-h_8{v}_1{u}_2=u_2$$

$$h_4{u}_1+h_5{v}_1+h_6-h_7{u}_1{v}_2-h_8{v}_1{v}_2=v_2$$

$$\left(\begin{array}{cccccccc}{u}_1& v_1& 1& 0& 0& 0& -u_1{u}_2& -v_1{u}_2\\ 0& 0& 0& u_1& v_1& 1& -u_1{v}_2& -v_1{v}_2\end{array}\right)\left(\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\end{array}\right)=\left(\begin{array}{c}{u}_2\\ v_2\end{array}\right)$$

所以，一对点可以引入两个约束，一共8个未知参数，故需要4对点即可求得映射关系：Ax=B

二、c++实现

1.CMakeLists.txt

cmake_minimum_required(VERSION 3.22.0)
project(homography)

# Eigen
include_directories("/usr/include/eigen3")
add_executable(homography homography.cpp)

2.homography.cpp

void calc_homography(int pts1[4][2], int pts2[4][2], int pts_in[][2], int pts_out[][2])
{
    Eigen::Matrix<float,8,8> A;
    Eigen::Matrix<float,8,1> B;
    for(int i = 0; i<4;++i) {
        A(i*2, 0) = pts1[i][0];
        A(i*2, 1) = pts1[i][1];
        A(i*2, 2) = 1.0;
        A(i*2, 3) = 0.0;
        A(i*2, 4) = 0.0;
        A(i*2, 5) = 0.0;
        A(i*2, 6) = -pts1[i][0] * pts2[i][0];
        A(i*2, 7) = -pts1[i][1] * pts2[i][0];

        A(i*2 + 1, 0) = 0.0;
        A(i*2 + 1, 1) = 0.0;
        A(i*2 + 1, 2) = 0.0;
        A(i*2 + 1, 3) = pts1[i][0];
        A(i*2 + 1, 4) = pts1[i][1];
        A(i*2 + 1, 5) = 1.0;
        A(i*2 + 1, 6) = -pts1[i][0] * pts2[i][1];
        A(i*2 + 1, 7) = -pts1[i][1] * pts2[i][1];
        
        B(i*2, 0) = pts2[i][0];
        B(i*2 + 1, 0) = pts2[i][1];
    }

    std::cout << "求解 Ax=B " << std::endl;
    std::cout << A << "\n" << std::endl;
    std::cout << B.transpose() << "\n" << std::endl;

    Eigen::MatrixXf H = A.colPivHouseholderQr().solve(B);

    Eigen::Matrix<float, 3, 3> X;
    X(0,0) = H(0,0);
    X(0,1) = H(1,0);
    X(0,2) = H(2,0);
    X(1,0) = H(3,0);
    X(1,1) = H(4,0);
    X(1,2) = H(5,0);
    X(2,0) = H(6,0);
    X(2,1) = H(7,0);
    X(2,2) = 1.0;
    std::cout <<"单应性矩阵为" <<std::endl;
    std::cout << X << std::endl;

    Eigen::MatrixXf X_ = X.inverse();
    for(int i=0; i<4; i++)
    {
        Eigen::Matrix<float,3,1> in;
        in << static_cast<float>(pts_in[i][0]), static_cast<float>(pts_in[i][1]), 1;

        Eigen::Matrix<float,3,1> out;
        out = X_ * in;

        pts_out[i][0] = static_cast<int>(out[0]/out[2]);
        pts_out[i][1] = static_cast<int>(out[1]/out[2]);
    }
    std::cout <<"\n" << "N个相机1图像中坐标的点: " << std::endl;
    std::cout << pts_out[0][0] << "   "<< pts_out[0][1] << std::endl;
    std::cout << pts_out[1][0] << "   "<< pts_out[1][1] << std::endl;
    std::cout << pts_out[2][0] << "   "<< pts_out[2][1] << std::endl;
    std::cout << pts_out[3][0] << "   "<< pts_out[3][1] << std::endl;
}

---

总结

注意由于缩放因子引入的约束。