单应性基本概念及代码解释

目录

• 简介
  • 基本理论
    • 什么是单应性矩阵？
    • 单应性变换有何用途？
  • 演示代码
    • 演示 1：共面点的姿态估计
    • 演示 2：透视校正
    • 演示 3：相机位移产生的单应性
      • 练习
    • 演示 4：分解单应性矩阵
    • 演示 5：旋转相机拍摄的全景图像基本拼接
  • 附加参考文献

上一个教程： AKAZE 和 ORB 平面跟踪

兼容性	OpenCV >= 3.0

引言

本教程将通过一些代码演示单应性的基本概念。有关理论的详细解释，请参考计算机视觉课程或计算机视觉书籍，例如：

• 《计算机视觉中的多视图几何》，Richard Hartley 和 Andrew Zisserman，[120]（部分章节可在此处获取：这里，CVPR 教程可在此处获取：这里）
• 《三维视觉邀请：从图像到几何模型》，Yi Ma、Stefano Soatto、Jana Kosecka 和 S. Shankar Sastry，[183]（一份计算机视觉书籍讲义可在此处获取：这里）
• 《计算机视觉：算法与应用》，Richard Szeliski，[266]（电子版可在此处获取：这里）
• 《更深入理解基于视觉控制的单应性分解》，Ezio Malis、Manuel Vargas，[186]（开放访问：这里）
• 《增强现实中的姿态估计：实践调查》，Eric Marchand、Hideaki Uchiyama、Fabien Spindler，[188]（开放访问：这里）

本教程的代码可在以下链接找到：C++，Python，Java。本教程中使用的图像可在此处找到（left*.jpg）。

基本理论

什么是单应性矩阵？

简而言之，平面单应性关联了两个平面之间的变换（直至一个比例因子）

\[ s \begin{bmatrix} x^{'} \\ y^{'} \\ 1 \end{bmatrix} = \mathbf{H} \begin{bmatrix} x \\ y \\ 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} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} \]

单应性矩阵是一个 `3x3` 矩阵，但具有 8 个自由度 (DoF)，因为它是在一个尺度上进行估计的。它通常被归一化（另见1），使得 \( h_{33} = 1 \) 或 \( h_{11}^2 + h_{12}^2 + h_{13}^2 + h_{21}^2 + h_{22}^2 + h_{23}^2 + h_{31}^2 + h_{32}^2 + h_{33}^2 = 1 \)。

以下示例展示了不同类型的变换，但都关联了两个平面之间的变换。

• 一个平面和一个图像平面（图像取自2）

• 由两个摄像机位置观察到的一个平面（图像取自3和2）

• 围绕其投影轴旋转的相机，等同于认为点位于无穷远平面上（图像取自2）

单应性变换有何用途？

• 例如，通过使用标记从共面点估计相机姿态以实现增强现实（参见之前的第一个示例）

• 透视去除/校正（参见之前的第二个示例）

• 全景图像拼接（参见之前的第二个和第三个示例）

演示代码

演示 1：共面点的姿态估计

注意

请注意，从单应性估计相机姿态的代码只是一个示例，如果您想估计平面或任意对象的相机姿态，您应该改用 cv::solvePnP。

单应性可以使用例如直接线性变换（DLT）算法进行估计（更多信息请参见1）。由于对象是平面的，因此在对象坐标系中表示的点与在归一化相机坐标系中表示的图像平面中的投影点之间的变换是单应性。只有当对象是平面时，才能从单应性中检索相机姿态，前提是相机内参已知（参见2或4）。这可以通过使用棋盘格对象和 `findChessboardCorners()` 来获取图像中的角点位置轻松测试。

第一步是检测棋盘格角点，需要棋盘格大小（`patternSize`），此处为 `9x6`

vector<Point2f> corners;
    bool found = findChessboardCorners(img, patternSize, corners);

已知棋盘格方块的大小，可以轻松计算出在对象坐标系中表示的对象点

        for( int i = 0; i < boardSize.height; i++ )
            for( int j = 0; j < boardSize.width; j++ )
corners.push_back(Point3f(float(j*squareSize),
                                          float(i*squareSize), 0));

对于单应性估计部分，必须移除坐标 `Z=0`

vector<Point3f> objectPoints;
calcChessboardCorners(patternSize, squareSize, objectPoints);
vector<Point2f> objectPointsPlanar;
    for (size_t i = 0; i < objectPoints.size(); i++)
    {
objectPointsPlanar.push_back(Point2f(objectPoints[i].x, objectPoints[i].y));
    }

通过角点并使用相机内参和畸变系数应用逆透视变换，可以计算出在归一化相机中表示的图像点

    FileStorage fs( samples::findFile( intrinsicsPath ), FileStorage::READ);
    Mat cameraMatrix, distCoeffs;
fs["camera_matrix"] >> cameraMatrix;
fs["distortion_coefficients"] >> distCoeffs;

vector<Point2f> imagePoints;
    undistortPoints(corners, imagePoints, cameraMatrix, distCoeffs);

然后可以使用以下方法估计单应性

    Mat H = findHomography(objectPointsPlanar, imagePoints);
cout << "H:\n" << H << endl;

从单应性矩阵中检索姿态的一个快速解决方案是（参见5）

    // 归一化以确保 ||c1|| = 1
    double norm = sqrt(H.at<double>(0,0)*H.at<double>(0,0) +
H.at<double>(1,0)*H.at<double>(1,0) +
H.at<double>(2,0)*H.at<double>(2,0));
 
H /= norm;
    Mat c1 = H.col(0);
    Mat c2 = H.col(1);
    Mat c3 = c1.cross(c2);
 
    Mat tvec = H.col(2);
    Mat R(3, 3, CV_64F);
 
    for (int i = 0; i < 3; i++)
    {
R.at<double>(i,0) = c1.at<double>(i,0);
R.at<double>(i,1) = c2.at<double>(i,0);
R.at<double>(i,2) = c3.at<double>(i,0);
    }

\[ \begin{align*} \mathbf{X} &= \left( X, Y, 0, 1 \right ) \\ \mathbf{x} &= \mathbf{P}\mathbf{X} \\ &= \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{r_3} \hspace{0.5em} \mathbf{t} \right ] \begin{pmatrix} X \\ Y \\ 0 \\ 1 \end{pmatrix} \\ &= \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{t} \right ] \begin{pmatrix} X \\ Y \\ 1 \end{pmatrix} \\ &= \mathbf{H} \begin{pmatrix} X \\ Y \\ 1 \end{pmatrix} \end{align*} \]

\[ \begin{align*} \mathbf{H} &= \lambda \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{t} \right ] \\ \mathbf{K}^{-1} \mathbf{H} &= \lambda \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{t} \right ] \\ \mathbf{P} &= \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \left( \mathbf{r_1} \times \mathbf{r_2} \right ) \hspace{0.5em} \mathbf{t} \right ] \end{align*} \]

这是一个快速解决方案（另见2），因为它不能确保生成的旋转矩阵是正交的，并且尺度是通过将第一列归一化为1来粗略估计的。

为了获得一个正确的旋转矩阵（具有旋转矩阵的特性），可以应用极分解，或者对旋转矩阵进行正交化（参见6或7或8或9以获取一些信息）

cout << "R (before polar decomposition):\n" << R << "\ndet(R): " << determinant(R) << endl;
    Mat_<double> W, U, Vt;
    SVDecomp(R, W, U, Vt);
R = U*Vt;
    double det = determinant(R);
    if (det < 0)
    {
Vt.at<double>(2,0) *= -1;
Vt.at<double>(2,1) *= -1;
Vt.at<double>(2,2) *= -1;
 
R = U*Vt;
    }
cout << "R (after polar decomposition):\n" << R << "\ndet(R): " << determinant(R) << endl;

为了检查结果，将以估计的相机姿态投射到图像中的对象坐标系显示出来

演示 2：透视校正

在此示例中，通过计算将源点映射到所需点的单应性，将源图像转换为所需的透视视图。下图显示了源图像（左）和我们想要转换为所需棋盘格视图的棋盘格视图（右）。

源视图和目标视图

第一步是在源图像和目标图像中检测棋盘格角点

C++

vector<Point2f> corners1, corners2;
    bool found1 = findChessboardCorners(img1, patternSize, corners1);
    bool found2 = findChessboardCorners(img2, patternSize, corners2);

Python

ret1, corners1 = cv.findChessboardCorners(img1, patternSize)
ret2, corners2 = cv.findChessboardCorners(img2, patternSize)

Java

MatOfPoint2f corners1 = new MatOfPoint2f(), corners2 = new MatOfPoint2f();
        boolean found1 = Calib3d.findChessboardCorners(img1, new Size(9, 6), corners1 );
        boolean found2 = Calib3d.findChessboardCorners(img2, new Size(9, 6), corners2 );

单应性可以很容易地估计出来

C++

    Mat H = findHomography(corners1, corners2);
cout << "H:\n" << H << endl;

Python

H, _ = cv.findHomography(corners1, corners2)
print(H)

Java

Mat H = new Mat();
H = Calib3d.findHomography(corners1, corners2);
System.out.println(H.dump());

为了将源棋盘格视图扭曲到所需的棋盘格视图，我们使用 cv::warpPerspective

C++

    Mat img1_warp;
    warpPerspective(img1, img1_warp, H, img1.size());

Python

img1_warp = cv.warpPerspective(img1, H, (img1.shape[1], img1.shape[0]))

Java

Mat img1_warp = new Mat();
Imgproc.warpPerspective(img1, img1_warp, H, img1.size());

结果图像是

计算由单应性变换的源角点坐标

C++

    Mat img_draw_matches;
    hconcat(img1, img2, img_draw_matches);
    for (size_t i = 0; i < corners1.size(); i++)
    {
        Mat pt1 = (Mat_<double>(3,1) << corners1[i].x, corners1[i].y, 1);
        Mat pt2 = H * pt1;
pt2 /= pt2.at<double>(2);
 
        Point end( (int) (img1.cols + pt2.at<double>(0)), (int) pt2.at<double>(1) );
        line(img_draw_matches, corners1[i], end, randomColor(rng), 2);
    }
 
    imshow("Draw matches", img_draw_matches);
    waitKey();

Python

img_draw_matches = cv.hconcat([img1, img2])
    for i in range(len(corners1))
pt1 = np.array([corners1[i][0], corners1[i][1], 1])
pt1 = pt1.reshape(3, 1)
pt2 = np.dot(H, pt1)
pt2 = pt2/pt2[2]
end = (int(img1.shape[1] + pt2[0]), int(pt2[1]))
        cv.line(img_draw_matches, tuple([int(j) for j in corners1[i]]), end, randomColor(), 2)
 
    cv.imshow("Draw matches", img_draw_matches)
    cv.waitKey(0)

Java

Mat img_draw_matches = new Mat();
list2.add(img1);
list2.add(img2);
Core.hconcat(list2, img_draw_matches);
        Point []corners1Arr = corners1.toArray();
 
        for (int i = 0 ; i < corners1Arr.length; i++) {
Mat pt1 = new Mat(3, 1, CvType.CV_64FC1), pt2 = new Mat();
pt1.put(0, 0, corners1Arr[i].x, corners1Arr[i].y, 1 );
 
Core.gemm(H, pt1, 1, new Mat(), 0, pt2);
            double[] data = pt2.get(2, 0);
Core.divide(pt2, new Scalar(data[0]), pt2);
 
            double[] data1 =pt2.get(0, 0);
            double[] data2 = pt2.get(1, 0);
            Point end = new Point((int)(img1.cols()+ data1[0]), (int)data2[0]);
Imgproc.line(img_draw_matches, corners1Arr[i], end, RandomColor(), 2);
        }
 
HighGui.imshow("Draw matches", img_draw_matches);
HighGui.waitKey(0);

为了检查计算的正确性，显示匹配线

演示 3：相机位移产生的单应性

单应性关联了两个平面之间的变换，并且可以检索相应的相机位移，该位移允许从第一个平面视图转换到第二个平面视图（更多信息请参见[186]）。在深入了解如何从相机位移计算单应性的细节之前，先回顾一下相机姿态和齐次变换。

函数 cv::solvePnP 允许根据 3D 对象点（在对象坐标系中表示的点）和投影的 2D 图像点（在图像中观察到的对象点）的对应关系来计算相机姿态。需要内参和畸变系数（参见相机标定过程）。

\[ \begin{align*} s \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} &= \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} r_{11} & r_{12} & r_{13} & t_x \\ r_{21} & r_{22} & r_{23} & t_y \\ r_{31} & r_{32} & r_{33} & t_z \end{bmatrix} \begin{bmatrix} X_o \\ Y_o \\ Z_o \\ 1 \end{bmatrix} \\ &= \mathbf{K} \hspace{0.2em} ^{c}\mathbf{M}_o \begin{bmatrix} X_o \\ Y_o \\ Z_o \\ 1 \end{bmatrix} \end{align*} \]

\( \mathbf{K} \) 是内参矩阵，而 \( ^{c}\mathbf{M}_o \) 是相机姿态。 cv::solvePnP 的输出正是这个：`rvec` 是 Rodrigues 旋转向量，`tvec` 是平移向量。

\( ^{c}\mathbf{M}_o \) 可以用齐次形式表示，并允许将以对象坐标系表示的点转换为相机坐标系

\[ \begin{align*} \begin{bmatrix} X_c \\ Y_c \\ Z_c \\ 1 \end{bmatrix} &= \hspace{0.2em} ^{c}\mathbf{M}_o \begin{bmatrix} X_o \\ Y_o \\ Z_o \\ 1 \end{bmatrix} \\ &= \begin{bmatrix} ^{c}\mathbf{R}_o & ^{c}\mathbf{t}_o \\ 0_{1\times3} & 1 \end{bmatrix} \begin{bmatrix} X_o \\ Y_o \\ Z_o \\ 1 \end{bmatrix} \\ &= \begin{bmatrix} r_{11} & r_{12} & r_{13} & t_x \\ r_{21} & r_{22} & r_{23} & t_y \\ r_{31} & r_{32} & r_{33} & t_z \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} X_o \\ Y_o \\ Z_o \\ 1 \end{bmatrix} \end{align*} \]

将一个坐标系中的点变换到另一个坐标系可以通过矩阵乘法轻松完成

• \( ^{c_1}\mathbf{M}_o \) 是相机 1 的相机姿态
• \( ^{c_2}\mathbf{M}_o \) 是相机 2 的相机姿态

将以相机 1 坐标系表示的 3D 点变换到相机 2 坐标系

\[ ^{c_2}\mathbf{M}_{c_1} = \hspace{0.2em} ^{c_2}\mathbf{M}_{o} \cdot \hspace{0.1em} ^{o}\mathbf{M}_{c_1} = \hspace{0.2em} ^{c_2}\mathbf{M}_{o} \cdot \hspace{0.1em} \left( ^{c_1}\mathbf{M}_{o} \right )^{-1} = \begin{bmatrix} ^{c_2}\mathbf{R}_{o} & ^{c_2}\mathbf{t}_{o} \\ 0_{3 \times 1} & 1 \end{bmatrix} \cdot \begin{bmatrix} ^{c_1}\mathbf{R}_{o}^T & - \hspace{0.2em} ^{c_1}\mathbf{R}_{o}^T \cdot \hspace{0.2em} ^{c_1}\mathbf{t}_{o} \\ 0_{1 \times 3} & 1 \end{bmatrix} \]

在这个例子中，我们将计算相对于棋盘格对象在两个相机姿态之间的相机位移。第一步是计算两个图像的相机姿态。

vector<Point2f> corners1, corners2;
    bool found1 = findChessboardCorners(img1, patternSize, corners1);
    bool found2 = findChessboardCorners(img2, patternSize, corners2);
 
    if (!found1 || !found2)
    {
cout << "Error, cannot find the chessboard corners in both images." << endl;
        return;
    }
 
vector<Point3f> objectPoints;
calcChessboardCorners(patternSize, squareSize, objectPoints);
 
    FileStorage fs( samples::findFile( intrinsicsPath ), FileStorage::READ);
    Mat cameraMatrix, distCoeffs;
fs["camera_matrix"] >> cameraMatrix;
fs["distortion_coefficients"] >> distCoeffs;
 
    Mat rvec1, tvec1;
    solvePnP(objectPoints, corners1, cameraMatrix, distCoeffs, rvec1, tvec1);
    Mat rvec2, tvec2;
    solvePnP(objectPoints, corners2, cameraMatrix, distCoeffs, rvec2, tvec2);

相机位移可以使用上述公式从相机姿态计算得出

void computeC2MC1(const Mat &R1, const Mat &tvec1, const Mat &R2, const Mat &tvec2,
                  Mat &R_1to2, Mat &tvec_1to2)
{
    //c2Mc1 = c2Mo * oMc1 = c2Mo * c1Mo.inv()
R_1to2 = R2 * R1.t();
tvec_1to2 = R2 * (-R1.t()*tvec1) + tvec2;
}

从相机位移计算出的与特定平面相关的单应性是

作者：Homography-transl.svg：Per Rosengren 衍生作品：Appoose (Homography-transl.svg) [CC BY 3.0 (http://creativecommons.org/licenses/by/3.0)]，通过维基共享资源

在此图中，`n` 是平面的法向量，`d` 是相机坐标系与平面沿平面法线之间的距离。计算从相机位移得到单应性的公式是

\[ ^{2}\mathbf{H}_{1} = \hspace{0.2em} ^{2}\mathbf{R}_{1} - \hspace{0.1em} \frac{^{2}\mathbf{t}_{1} \cdot \hspace{0.1em} ^{1}\mathbf{n}^\top}{^1d} \]

其中 \( ^{2}\mathbf{H}_{1} \) 是将第一个相机坐标系中的点映射到第二个相机坐标系中相应点的单应性矩阵，\( ^{2}\mathbf{R}_{1} = \hspace{0.2em} ^{c_2}\mathbf{R}_{o} \cdot \hspace{0.1em} ^{c_1}\mathbf{R}_{o}^{\top} \) 是表示两个相机坐标系之间旋转的旋转矩阵，\( ^{2}\mathbf{t}_{1} = \hspace{0.2em} ^{c_2}\mathbf{R}_{o} \cdot \left( - \hspace{0.1em} ^{c_1}\mathbf{R}_{o}^{\top} \cdot \hspace{0.1em} ^{c_1}\mathbf{t}_{o} \right ) + \hspace{0.1em} ^{c_2}\mathbf{t}_{o} \) 是两个相机坐标系之间的平移向量。

这里，法向量 `n` 是在相机坐标系 1 中表示的平面法线，可以计算为 2 个向量的叉积（使用位于平面上的 3 个非共线点）或者在我们的例子中直接使用

    Mat normal = (Mat_<double>(3,1) << 0, 0, 1);
    Mat normal1 = R1*normal;

距离 `d` 可以计算为平面法线与平面上一点的点积，或者通过计算平面方程并使用 D 系数

    Mat origin(3, 1, CV_64F, Scalar(0));
    Mat origin1 = R1*origin + tvec1;
    double d_inv1 = 1.0 / normal1.dot(origin1);

射影单应性矩阵 \( \textbf{G} \) 可以使用内参矩阵 \( \textbf{K} \) 从欧几里得单应性 \( \textbf{H} \) 计算得出（参见[186]），这里假设两个平面视图之间使用相同的相机

\[ \textbf{G} = \gamma \textbf{K} \textbf{H} \textbf{K}^{-1} \]

Mat computeHomography(const Mat &R_1to2, const Mat &tvec_1to2, const double d_inv, const Mat &normal)
{
    Mat homography = R_1to2 + d_inv * tvec_1to2*normal.t();
    return homography;
}

在我们的例子中，棋盘的 Z 轴指向物体内部，而在单应性图中则指向外部。这只是符号问题。

\[ ^{2}\mathbf{H}_{1} = \hspace{0.2em} ^{2}\mathbf{R}_{1} + \hspace{0.1em} \frac{^{2}\mathbf{t}_{1} \cdot \hspace{0.1em} ^{1}\mathbf{n}^\top}{^1d} \]

    Mat homography_euclidean = computeHomography(R_1to2, t_1to2, d_inv1, normal1);
    Mat homography = cameraMatrix * homography_euclidean * cameraMatrix.inv();
 
homography /= homography.at<double>(2,2);
homography_euclidean /= homography_euclidean.at<double>(2,2);

现在我们将比较从相机位移计算的射影单应性与使用 cv::findHomography 估计的单应性

findHomography H
[0.32903393332201, -1.244138808862929, 536.4769088231476;
 0.6969763913334046, -0.08935909072571542, -80.34068504082403;
 0.00040511729592961, -0.001079740100565013, 0.9999999999999999]
 
相机位移产生的单应性
[0.4160569997384721, -1.306889006892538, 553.7055461075881;
 0.7917584252773352, -0.06341244158456338, -108.2770029401219;
 0.0005926357240956578, -0.001020651672127799, 1]

单应性矩阵相似。如果我们比较使用两个单应性矩阵扭曲的图像 1

左图：使用估计的单应性扭曲的图像。右图：使用从相机位移计算的单应性扭曲的图像。

从视觉上看，很难区分从相机位移计算出的单应性与使用 cv::findHomography 函数估计出的单应性所产生的图像结果之间的差异。

练习

此演示向您展示如何从两个相机姿态计算单应性变换。尝试执行相同的操作，但这次计算 N 个中间单应性。而不是计算一个单应性来直接将源图像扭曲到所需的相机视角，而是执行 N 个扭曲操作以查看不同的变换。

您应该得到类似以下的结果

前三张图像显示了源图像在三个不同插值相机视角下的扭曲。第四张图像显示了最终相机视角下的扭曲源图像与所需图像之间的“误差图像”。

演示 4：分解单应性矩阵

OpenCV 3 包含函数 cv::decomposeHomographyMat，它允许将单应性矩阵分解为一组旋转、平移和平面法线。首先，我们将分解从相机位移计算的单应性矩阵。

    Mat homography_euclidean = computeHomography(R_1to2, t_1to2, d_inv1, normal1);
    Mat homography = cameraMatrix * homography_euclidean * cameraMatrix.inv();
 
homography /= homography.at<double>(2,2);
homography_euclidean /= homography_euclidean.at<double>(2,2);

cv::decomposeHomographyMat 的结果是

vector<Mat> Rs_decomp, ts_decomp, normals_decomp;
    int solutions = decomposeHomographyMat(homography, cameraMatrix, Rs_decomp, ts_decomp, normals_decomp);
cout << "Decompose homography matrix computed from the camera displacement:" << endl << endl;
    for (int i = 0; i < solutions; i++)
    {
      double factor_d1 = 1.0 / d_inv1;
      Mat rvec_decomp;
      Rodrigues(Rs_decomp[i], rvec_decomp);
cout << "Solution " << i << ":" << endl;
cout << "rvec from homography decomposition: " << rvec_decomp.t() << endl;
cout << "rvec from camera displacement: " << rvec_1to2.t() << endl;
cout << "tvec from homography decomposition: " << ts_decomp[i].t() << " and scaled by d: " << factor_d1 * ts_decomp[i].t() << endl;
cout << "tvec from camera displacement: " << t_1to2.t() << endl;
cout << "plane normal from homography decomposition: " << normals_decomp[i].t() << endl;
cout << "plane normal at camera 1 pose: " << normal1.t() << endl << endl;
    }

解 0
单应性分解中的 rvec: [-0.0919829920641369, -0.5372581036567992, 1.310868863540717]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [-0.7747961019053186, -0.02751124463434032, -0.6791980037590677]，按 d 缩放后: [-0.1578091561210742, -0.005603443652993778, -0.1383378976078466]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [-0.1973513139420648, 0.6283451996579074, -0.7524857267431757]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
 
解 1
单应性分解中的 rvec: [-0.0919829920641369, -0.5372581036567992, 1.310868863540717]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [0.7747961019053186, 0.02751124463434032, 0.6791980037590677]，按 d 缩放后: [0.1578091561210742, 0.005603443652993778, 0.1383378976078466]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [0.1973513139420648, -0.6283451996579074, 0.7524857267431757]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
 
解 2
单应性分解中的 rvec: [0.1053487907109967, -0.1561929144786397, 1.401356552358475]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [-0.4666552552894618, 0.1050032934770042, -0.913007654671646]，按 d 缩放后: [-0.0950475510338766, 0.02138689274867372, -0.1859598508065552]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [-0.3131715472900788, 0.8421206145721947, -0.4390403768225507]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
 
解 3
单应性分解中的 rvec: [0.1053487907109967, -0.1561929144786397, 1.401356552358475]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [0.4666552552894618, -0.1050032934770042, 0.913007654671646]，按 d 缩放后: [0.0950475510338766, -0.02138689274867372, 0.1859598508065552]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [0.3131715472900788, -0.8421206145721947, 0.4390403768225507]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]

单应性矩阵分解的结果只能恢复到与距离 `d` 对应的尺度因子，因为法线是单位长度的。如您所见，有一个解与计算出的相机位移几乎完全匹配。如文档所述：

如果点对应可用，通过应用正深度约束（所有点必须在相机前方），至少可以使其中两个解无效。

由于分解的结果是相机位移，如果我们有初始相机姿态 \( ^{c_1}\mathbf{M}_{o} \)，我们可以计算当前相机姿态 \( ^{c_2}\mathbf{M}_{o} = \hspace{0.2em} ^{c_2}\mathbf{M}_{c_1} \cdot \hspace{0.1em} ^{c_1}\mathbf{M}_{o} \) 并测试属于平面的 3D 对象点是否投影在相机前方。另一种解决方案可能是，如果我们知道在相机 1 姿态处表示的平面法线，则保留与最近法线匹配的解决方案。

同样的操作，但是使用 cv::findHomography 估计的单应性矩阵

解 0
单应性分解中的 rvec: [0.1552207729599141, -0.152132696119647, 1.323678695078694]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [-0.4482361704818117, 0.02485247635491922, -1.034409687207331]，按 d 缩放后: [-0.09129598307571339, 0.005061910238634657, -0.2106868109173855]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [-0.1384902722707529, 0.9063331452766947, -0.3992250922214516]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
 
解 1
单应性分解中的 rvec: [0.1552207729599141, -0.152132696119647, 1.323678695078694]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [0.4482361704818117, -0.02485247635491922, 1.034409687207331]，按 d 缩放后: [0.09129598307571339, -0.005061910238634657, 0.2106868109173855]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [0.1384902722707529, -0.9063331452766947, 0.3992250922214516]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
 
解 2
单应性分解中的 rvec: [-0.2886605671759886, -0.521049903923871, 1.381242030882511]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [-0.8705961357284295, 0.1353018038908477, -0.7037702049789747]，按 d 缩放后: [-0.177321544550518, 0.02755804196893467, -0.1433427218822783]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [-0.2284582117722427, 0.6009247303964522, -0.7659610393954643]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
 
解 3
单应性分解中的 rvec: [-0.2886605671759886, -0.521049903923871, 1.381242030882511]
相机位移中的 rvec: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
单应性分解中的 tvec: [0.8705961357284295, -0.1353018038908477, 0.7037702049789747]，按 d 缩放后: [0.177321544550518, -0.02755804196893467, 0.1433427218822783]
相机位移中的 tvec: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
单应性分解中的平面法线: [0.2284582117722427, -0.6009247303964522, 0.7659610393954643]
相机 1 姿态处的平面法线: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]

同样，也有一个与计算出的相机位移匹配的解决方案。

演示 5：旋转相机拍摄的全景图像基本拼接

注意

本例旨在说明基于相机纯旋转运动（非平移）的图像拼接概念，不应用于拼接全景图像。stitching 模块提供了一个完整的图像拼接流水线。

单应性变换仅适用于平面结构。但在相机旋转（围绕相机投影轴纯旋转，无平移）的情况下，可以考虑任意世界（参见上文）。

单应性可以根据旋转变换和相机内参计算，公式如下（例如参见10）：

\[ s \begin{bmatrix} x^{'} \\ y^{'} \\ 1 \end{bmatrix} = \bf{K} \hspace{0.1em} \bf{R} \hspace{0.1em} \bf{K}^{-1} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} \]

为了说明这一点，我们使用了 Blender，一款免费的开源 3D 计算机图形软件，生成了两个相机视图，它们之间只有旋转变换。有关如何使用 Blender 检索相机内参和相对于世界的 `3x4` 外参矩阵的更多信息，请参见11（还需要额外的变换才能获得相机和物体坐标系之间的变换）。

下图显示了 Suzanne 模型的两个生成视图，仅进行了旋转变换

通过已知的相机姿态和内参，可以计算出两个视图之间的相对旋转

C++

    Mat R1 = c1Mo(Range(0,3), Range(0,3));
    Mat R2 = c2Mo(Range(0,3), Range(0,3));

Python

R1 = c1Mo[0:3, 0:3]
R2 = c2Mo[0:3, 0:3]

Java

Range rowRange = new Range(0,3);
Range colRange = new Range(0,3);

C++

    //c1Mo * oMc2
    Mat R_2to1 = R1*R2.t();

Python

R2 = R2.transpose()
R_2to1 = np.dot(R1,R2)

Java

        //c1Mo * oMc2
Mat R1 = new Mat(c1Mo, rowRange, colRange);
Mat R2 = new Mat(c2Mo, rowRange, colRange);
Mat R_2to1 = new Mat();
Core.gemm(R1, R2.t(), 1, new Mat(), 0, R_2to1 );

这里，第二张图像将相对于第一张图像进行拼接。单应性可以使用上述公式计算。

C++

    Mat H = cameraMatrix * R_2to1 * cameraMatrix.inv();
H /= H.at<double>(2,2);
cout << "H:\n" << H << endl;

Python

H = cameraMatrix.dot(R_2to1).dot(np.linalg.inv(cameraMatrix))
H = H / H[2][2]

Java

Mat tmp = new Mat(), H = new Mat();
Core.gemm(cameraMatrix, R_2to1, 1, new Mat(), 0, tmp);
Core.gemm(tmp, cameraMatrix.inv(), 1, new Mat(), 0, H);
        Scalar s = new Scalar(H.get(2, 2)[0]);
Core.divide(H, s, H);
System.out.println(H.dump());

拼接操作非常简单

C++

    Mat img_stitch;
    warpPerspective(img2, img_stitch, H, Size(img2.cols*2, img2.rows));
    Mat half = img_stitch(Rect(0, 0, img1.cols, img1.rows));
img1.copyTo(half);

Python

img_stitch = cv.warpPerspective(img2, H, (img2.shape[1]*2, img2.shape[0]))
img_stitch[0:img1.shape[0], 0:img1.shape[1]] = img1

Java

Mat img_stitch = new Mat();
Imgproc.warpPerspective(img2, img_stitch, H, new Size(img2.cols()*2, img2.rows()) );
Mat half = new Mat();
half = new Mat(img_stitch, new Rect(0, 0, img1.cols(), img1.rows()));
img1.copyTo(half);

结果图像是

附加参考文献

• 1. 讲座 16：平面单应性，Robert Collins
• 2. 2D 射影变换（单应性），Christiano Gava，Gabriele Bleser
• 3. 计算机视觉：算法与应用，Richard Szeliski
• 4. 可视化跟踪和平面标记的逐步相机姿态估计
• 5. 从单应性估计姿态
• 6. 极分解（在连续介质力学中）
• 7. 第 3 章 - 3.1.2 从矩阵到旋转 - 定理 3.1（从矩阵 K 最小二乘估计旋转）
• 8. 与奇异值分解的个人访谈，Matan Gavish
• 9. Kabsch 算法，最优旋转矩阵的计算
• 10. 单应性，Gerhard Roth 博士
• 11. 从 Blender 相机获取 3x4 相机矩阵