用代码解释单应性的基本概念

目录

• 引言
  • 基本理论
    • 什么是单应矩阵？
    • 单应变换如何有用？
  • 演示代码
    • 演示 1：来自共面点的姿态估计
    • 演示 2：透视校正
    • 演示 3：来自相机位移的单应性
      • 练习
    • 演示 4：分解单应矩阵
    • 演示 5：来自旋转相机的基本全景拼接
  • 其他参考资料

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

兼容性	OpenCV >= 3.0

引言

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

• 计算机视觉中的多视图几何，Richard Hartley 和 Andrew Zisserman，[116]（一些示例章节可在此处获得 这里，CVPR 教程可在此处获得 这里）
• 3D 视觉入门：从图像到几何模型，Yi Ma、Stefano Soatto、Jana Kosecka 和 S. Shankar Sastry，[177]（计算机视觉书籍讲义可在此处获得 这里）
• 计算机视觉：算法与应用，Richard Szeliski，[260]（电子版可在此处获得 这里）
• 基于视觉的控制中单应性分解的更深入理解，Ezio Malis，Manuel Vargas，[180]（开放访问 这里）
• 增强现实的姿态估计：一个动手调查，Eric Marchand，Hideaki Uchiyama，Fabien Spindler，[182]（开放访问 这里）

本教程代码可在此处找到 C++、Python、Java。本教程中使用的图像可在此处找到 这里 (left*.jpg)。

基本理论

什么是单应矩阵？

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

$$s\left[\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ 1\end{array}\right]=\mathbf{H}\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$

单应矩阵是一个 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{array}{rl}\mathbf{X}& =(X,Y,0,1)\\ \mathbf{x}& =\mathbf{P}\mathbf{X}\\ & =\mathbf{K}\left[{\mathbf{r}}_{\mathbf{1}}{\mathbf{r}}_{\mathbf{2}}{\mathbf{r}}_{\mathbf{3}}\mathbf{t}\right]\left(\begin{array}{c}X\\ Y\\ 0\\ 1\end{array}\right)\\ & =\mathbf{K}\left[{\mathbf{r}}_{\mathbf{1}}{\mathbf{r}}_{\mathbf{2}}\mathbf{t}\right]\left(\begin{array}{c}X\\ Y\\ 1\end{array}\right)\\ & =\mathbf{H}\left(\begin{array}{c}X\\ Y\\ 1\end{array}\right)\end{array}$$

$$\begin{array}{rl}\mathbf{H}& =\lambda \mathbf{K}\left[{\mathbf{r}}_{\mathbf{1}}{\mathbf{r}}_{\mathbf{2}}\mathbf{t}\right]\\ {\mathbf{K}}^{-1}\mathbf{H}& =\lambda \left[{\mathbf{r}}_{\mathbf{1}}{\mathbf{r}}_{\mathbf{2}}\mathbf{t}\right]\\ \mathbf{P}& =\mathbf{K}\left[{\mathbf{r}}_{\mathbf{1}}{\mathbf{r}}_{\mathbf{2}}({\mathbf{r}}_{\mathbf{1}}\times {\mathbf{r}}_{\mathbf{2}})\mathbf{t}\right]\end{array}$$

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

要获得正确的旋转矩阵（具有旋转矩阵的属性），一种解决方案是应用极分解或旋转矩阵的正交化（参见6或7或8或9 获取更多信息）

cout << "R (极分解前):\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 (极分解后):\n" << R << "\ndet(R): " << determinant(R) << endl;

为了检查结果，将使用估计的相机姿态将物体坐标系投影到图像中。

演示 2：透视校正

C++PythonJava

在这个例子中，源图像将被转换为所需的透视图，方法是计算将源点映射到所需点的单应性。下图显示了源图像（左）和我们想要将其转换为所需棋盘视图的棋盘视图（右）。

源视图和目标视图

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

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++

    cv::Mat img_draw_matches;
    cv::hconcat(img1, img2, img_draw_matches);
    for (size_t i = 0; i < corners1.size(); i++)
    {
        cv::Mat pt1 = (cv::Mat_<double>(3,1) << corners1[i].x, corners1[i].y, 1);
        cv::Mat pt2 = H * pt1;
pt2 /= pt2.at<double>(2);
 
        cv::Point end( (int) (img1.cols + pt2.at<double>(0)), (int) pt2.at<double>(1) );
        cv::line(img_draw_matches, corners1[i], end, randomColor(rng), 2);
    }
 
    cv::imshow("Draw matches", img_draw_matches);
    cv::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：来自相机位移的单应性

单应性描述了两个平面之间的变换关系，可以从中获取相应的相机位移，从而实现从第一个平面视图到第二个平面视图的转换（更多信息请参见[180]）。在详细介绍如何根据相机位移计算单应性之前，需要先回顾一下相机姿态和齐次变换。

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

$$\begin{array}{rl}s\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]& =\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{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{array}\right]\left[\begin{array}{c}{X}_o\\ Y_o\\ Z_o\\ 1\end{array}\right]\\ & =\mathbf{K}{}^c{\mathbf{M}}_o\left[\begin{array}{c}{X}_o\\ Y_o\\ Z_o\\ 1\end{array}\right]\end{array}$$

$\mathbf{K}$ 是内参矩阵， ${}^c{\mathbf{M}}_o$ 是相机姿态。cv::solvePnP 的输出正是如此：rvec是罗德里格斯旋转向量，tvec是平移向量。

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

$$\begin{array}{rl}\left[\begin{array}{c}{X}_c\\ Y_c\\ Z_c\\ 1\end{array}\right]& ={}^c{\mathbf{M}}_o\left[\begin{array}{c}{X}_o\\ Y_o\\ Z_o\\ 1\end{array}\right]\\ & =\left[\begin{array}{cc}{}^c{\mathbf{R}}_o& {}^c{\mathbf{t}}_o\\ 0_{1\times 3}& 1\end{array}\right]\left[\begin{array}{c}{X}_o\\ Y_o\\ Z_o\\ 1\end{array}\right]\\ & =\left[\begin{array}{cccc}{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{array}\right]\left[\begin{array}{c}{X}_o\\ Y_o\\ Z_o\\ 1\end{array}\right]\end{array}$$

使用矩阵乘法可以轻松地将以一个坐标系表示的点转换为另一个坐标系。

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

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

$${}^{c_2}{\mathbf{M}}_{c_1}={}^{c_2}{\mathbf{M}}_o\cdot {}^o{\mathbf{M}}_{c_1}={}^{c_2}{\mathbf{M}}_o\cdot {\left({}^{c_1}{\mathbf{M}}_o\right)}^{-1}=\left[\begin{array}{cc}{}^{c_2}{\mathbf{R}}_o& {}^{c_2}{\mathbf{t}}_o\\ 0_{3\times 1}& 1\end{array}\right]\cdot \left[\begin{array}{cc}{}^{c_1}{\mathbf{R}}_o^T& -{}^{c_1}{\mathbf{R}}_o^T\cdot {}^{c_1}{\mathbf{t}}_o\\ 0_{1\times 3}& 1\end{array}\right]$$

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

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;
 
    cv::Mat rvec1, tvec1;
    cv::solvePnP(objectPoints, corners1, cameraMatrix, distCoeffs, rvec1, tvec1);
    cv::Mat rvec2, tvec2;
    cv::solvePnP(objectPoints, corners2, cameraMatrix, distCoeffs, rvec2, tvec2);

可以使用上面的公式根据相机姿态计算相机位移。

void computeC2MC1(const cv::Mat &R1, const cv::Mat &tvec1, const cv::Mat &R2, const cv::Mat &tvec2,
                  cv::Mat &R_1to2, cv::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 derivative work: Appoose (Homography-transl.svg) [CC BY 3.0 (http://creativecommons.org/licenses/by/3.0)], via Wikimedia Commons

在该图中，n 是平面的法向量，d 是相机坐标系与平面沿平面法向量的距离。从相机位移计算单应性的公式为：

$${}^2{\mathbf{H}}_1={}^2{\mathbf{R}}_1-\frac{{}^2{\mathbf{t}}_1\cdot {}^1{\mathbf{n}}^{\mathrm{\top }}}{{}^{1}d}$$

其中 ${}^2{\mathbf{H}}_1$ 是将第一个相机坐标系中的点映射到第二个相机坐标系中对应点的单应性矩阵，${}^2{\mathbf{R}}_1={}^{c_2}{\mathbf{R}}_o\cdot {}^{c_1}{\mathbf{R}}_o^{\mathrm{\top }}$ 是表示两个相机坐标系之间旋转的旋转矩阵，${}^2{\mathbf{t}}_1={}^{c_2}{\mathbf{R}}_o\cdot (-{}^{c_1}{\mathbf{R}}_o^{\mathrm{\top }}\cdot {}^{c_1}{\mathbf{t}}_o)+{}^{c_2}{\mathbf{t}}_o$ 是两个相机坐标系之间的平移向量。

这里法向量n是在相机坐标系1中表达的平面法向量，可以通过两个向量的叉积计算（使用位于平面上的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);

可以使用内参矩阵$\mathbf{\text{K}}$（参见[180]），从欧几里得单应性矩阵$\mathbf{\text{H}}$计算投影单应性矩阵$\mathbf{\text{G}}$，这里假设两个平面视图之间使用相同的相机。

$$\mathbf{\text{G}}=\gamma \mathbf{\text{K}}\mathbf{\text{H}}{\mathbf{\text{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={}^2{\mathbf{R}}_1+\frac{{}^2{\mathbf{t}}_1\cdot {}^1{\mathbf{n}}^{\mathrm{\top }}}{{}^{1}d}$$

    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 << "从相机位移计算出的单应性矩阵分解结果:" << 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 << "解 " << i << ":" << endl;
cout << "单应性矩阵分解得到的rvec: " << rvec_decomp.t() << endl;
cout << "相机位移得到的rvec: " << rvec_1to2.t() << endl;
cout << "单应性矩阵分解得到的tvec: " << ts_decomp[i].t() << " 并按d缩放: " << factor_d1 * ts_decomp[i].t() << endl;
cout << "相机位移得到的tvec: " << t_1to2.t() << endl;
cout << "单应性矩阵分解得到的平面法向量: " << normals_decomp[i].t() << endl;
cout << "相机姿态1处的平面法向量: " << 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={}^{c_2}{\mathbf{M}}_{c_1}\cdot {}^{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：来自旋转相机的基本全景拼接

C++PythonJava

注意

此示例旨在说明基于相机的纯旋转运动的图像拼接概念，不应将其用于拼接全景图像。拼接模块 提供了一个完整的图像拼接流程。

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

然后可以使用旋转变换和相机内参计算单应性，例如（参见 10）

$$s\left[\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ 1\end{array}\right]=\mathbf{K}\mathbf{R}{\mathbf{K}}^{\mathbf{-}\mathbf{1}}\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$

为了说明这一点，我们使用了 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++

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

Python

img_stitch = cv::warpPerspective(img2, H, cv::Size(img2.cols*2, img2.rows));
img_stitch(cv::Range(0, img1.rows), cv::Range(0, img1.cols)) = img1;

Java

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

生成的图像为

更多参考文献

• 1. 第16讲：平面单应性，Robert Collins
• 2. 二维投影变换（单应性），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相机矩阵