基于双目视觉的三维重建C++实战

自摄影开始以来，生成图像的过程并没有改变。来自观察场景的光线由相机通过前光圈（镜头）捕获，将光线照射到镜头后面的图像平面上。该过程如下图所示：

在上图中，do从镜头到被观察物体的距离，di是镜头到像平面的距离。f因此，它将成为镜头的焦距。薄透镜方程之间的关系如下所示：

现在让我们来看看现实世界 3 如何投影维对象？ 2 维平面（照片）的过程。我们理解这一点的最好方法是看看相机是如何工作的。

相机可视为将军 3-D 世界映射到 2-D 图像的功能。让我们以最简单的相机模型为例，即针孔相机模型，这是人类历史上更古老的摄影机制。以下是针孔摄像头的工作图：

我们可以从这张图中得出：

图像的大小自然是由物体形成的hi从物体到相机的距离do成反比。此外，位于 (X, Y, Z) 位置的 3-D 投影到场景点 (x,y) 在图像平面上，其中 (x,y) = (fX/Z, fY/Z)。其中 Z 坐标是指点的深度，这在上一张图像中完成。齐次坐标可用于整个相机配置和符号 系列用简单的矩阵来描述。

当相机生成世界的投影图像时，投影几何被用作现实世界中物体几何、旋转和变换的代数。

齐次坐标是射影几何中使用的坐标系统。即使我们能在欧几里得空间中表达现实世界中的对象(或 3-D 空间中的任何点)的位置，但必须执行的任何变化或旋转都必须在齐次坐标空间中执行，然后返回。让我们来看看使用齐次坐标的优点：

• 涉及齐次坐标的公式通常比笛卡尔世界简单。
• 有限坐标可以表示无限远处的点。
• 单个矩阵可以代表相机和世界之间可能发生的所有可能的保护性转换。

在齐次坐标空间中，2-D点用三个向量表示，3-D点用四个向量表示：

在上述方程中，第一个F符号矩阵称为内参矩阵（或通常称为内参矩阵）。这里的内部矩阵现在只包含焦距（f），在本教程之前，我们将研究更多这个矩阵的参数。

具有 r 和 t 符号的第二个矩阵称为外部参数矩阵（或通常称为外部矩阵）。矩阵中的元素表示相机的旋转和平移参数（即相机在现实世界中的位置和方式）。

因此，这些内外矩阵可以为我们提供图像 (x,y) 点和现实世界 (X, Y, Z) 点之间的关系。这是基于给定相机的内外参数 3-D 投影到场景点 2-D 平面上的方法。

现在我们对射影几何和相机模型有了足够的了解，是时候介绍计算机视觉几何中最重要的元素之一了。

现在我们知道该怎么办了 3-D 世界上的点投影到相机的图像平面上。我们将研究显示两个图像在同一场景中的投影关系。当这两个相机被刚性基线分开时，我们使用术语进行三维视觉。考虑到两个针孔相机观察一个给定的场景点共享相同的基线，如下图所示：

从上图中，世界点X的图像位于图像平面上x，现在这个x可以位于 3-D 这条线在空间中的任何位置。这意味着如果我们想在另一个图像中找到相同的点x，我们需要沿着这条线搜索第二幅图像上的投影。

从x绘制的假想线称为x极线。这条核线带来了一个基本的约束，即给定点的匹配在另一个视图中必须位于这条线上。这意味着如果你想从第二张图像中的第一张图像中找到它x ，你必须在第二张图像上找到x的核线。这些极线可以表示两个视图之间的几何形状。这里需要注意的是，所有的核线总是通过一个点。该点对应于一个摄像头中心到另一个摄像头中心的投影，称为极点。

我们可以将基本矩阵F视为将一个视图中的二维图像点映射到另一个图像视图中的核线的矩阵。通过解决一组方程，可以估计图像之间的基本矩阵，这涉及到两个图像之间一定数量的已知匹配点。这种匹配的最小数量是7，最佳数量是8。然后，对于图像中的一个点，基本矩阵给出了应该在另一个视图中找到相应点的方程。

如果一个点(x,y)一个点的对应点是(x’,y)，两个图像平面之间的基本矩阵是F，在齐次坐标中，我们必须有以下方程：

这个方程表达了两个对应点之间的关系，称为对极约束。

当两个相机观察同一个场景时，它们看到相同的物体，但从不同的角度。C 和 Python 中都有像 OpenCV 这样的图书馆为我们提供了一个特征探测器，它们可以在图像中找到一些带有描述符的点。他们认为这些点是图像中唯一的。如果给出同一场景的另一个图像，可以找到这些点。但实际上并不能保证通过比较检测到的特征点的描述符(如 SIFT、ORB 等)两幅图像之间获得的匹配集准确真实。这就是为什么引入是基于RANSAC基本矩阵估计方法(随机采样共识)策略。

RANSAC 背后的想法是从给定的数据点中随机选择一些数据点，并且只使用这些数据点进行估计。在我们的基本矩阵的例子中，所选点的数量应该是估计数学实体所需的最小点。基本矩阵一旦从这八个随机匹配中估计，所有其他匹配都将测试我们讨论的极限约束。基本矩阵的支持集是这些匹配形成计算的。

支持集越大，计算出的矩阵正确的概率就越高。如果随机选择的匹配之一是不正确的匹配，那么计算出的基本矩阵将是不正确的，集预计将非常小。这个过程重复了很多次，最后，最大支持集的矩阵将被保留为最有可能的矩阵。

人类进化成有两只眼睛的物种的原因是我们能感知到深度。当我们以类似的方式在机器中组织相机时，它被称为三维视觉。三维视觉系统通常由两个并排相机组成，观察相同的场景。下图显示了具有理想配置的三维设备的设置，完美对齐。

在上图所示的理想配置下，相机仅通过水平移分离，因此所有核线都是水平的。这意味着相应的点有相同的y坐标，搜索减少到一维线。当相机被如此纯水平移分离时，第二个相机的投影方程将变为：

查看下图，这个等式会更有意义，这是数码相机的一般情况：

其中点(uo, vo)像平面的像素位置是通过镜头主点线穿过的。这里我们得到了一种关系:

#include <opencv2/opencv.hpp> #include "openc2/xfeatures2d.hpp"
using namespace std;
using namespace cv;
int main(){
cv::Mat img1, img2;
img1 = cv::imread("imR.png",cv::IMREAD_GRAYSCALE);
img2 = cv::imread("imL.png",cv::IMREAD_GRAYSCALE);

我们做的第一件事是包含来自 opencv 和 opencv-contrib 的所需库，我要求你在开始本节之前构建它们。在main()函数中，我们初始化了cv:Mat数据类型的两个变量，它是 opencv 库的成员函数，Mat数据类型可以通过动态分配内存来保存任意大小的向量，尤其是图像。然后使用cv::imread()我们将图像导入mat数据类型的img1和img2中。cv::IMREAD_GRAYSCALE参数将图像导入为灰度。

// Define keypoints vector
std::vector<cv::KeyPoint> keypoints1, keypoints2;
// Define feature detector
cv::Ptr<cv::Feature2D> ptrFeature2D = cv::xfeatures2d::SIFT::create(74);
// Keypoint detection
ptrFeature2D->detect(img1,keypoints1);
ptrFeature2D->detect(img2,keypoints2);
// Extract the descriptor
cv::Mat descriptors1;
cv::Mat descriptors2;
ptrFeature2D->compute(img1,keypoints1,descriptors1);
ptrFeature2D->compute(img2,keypoints2,descriptors2);

在这里，我们使用 opencv 的 SIFT 特征检测器来从图像中提取所需的特征点。如果你想了解有关这些特征检测器如何工作的更多信息，请访问此链接。我们上面获得的描述符描述了提取的每个点，这个描述用于在另一个图像中找到它：

// Construction of the matcher
cv::BFMatcher matcher(cv::NORM_L2);
// Match the two image descriptors
std::vector<cv::DMatch> outputMatches;
matcher.match(descriptors1,descriptors2, outputMatches);

BFMatcher获取第一组中一个特征的描述符，并使用一些阈值距离计算与第二组中的所有其他特征匹配，并返回最接近的一个。我们将 BFMatches 返回的所有匹配项存储在vectorcv::DMatch类型的输出匹配变量中。

// Convert keypoints into Point2f
std::vector<cv::Point2f> points1, points2;
for (std::vector<cv::DMatch>::const_iterator it=   outputMatches.begin(); it!= outputMatches.end(); ++it) {
    // Get the position of left keypoints
    points1.push_back(keypoints1[it->queryIdx].pt);
    // Get the position of right keypoints
    points2.push_back(keypoints2[it->trainIdx].pt);
     }

获取的关键点首先需要转换为cv::Point2f类型，以便与cv::findFundamentalMat一起使用，我们将使用该函数使用我们抽象的这些特征点来计算基本矩阵。两个结果向量Points1和Points2包含两个图像中的对应点坐标。

std::vector<uchar> inliers(points1.size(),0);
cv::Mat fundamental= cv::findFundamentalMat(
  points1,points2, // matching points
     inliers,         // match status (inlier or outlier)  
     cv::FM_RANSAC,   // RANSAC method
     1.0,        // distance to epipolar line
     0.98);     // confidence probability
cout<<fundamental; //include this for seeing fundamental matrix

最后，我们调用了 cv::findFundamentalMat。

// Compute homographic rectification
cv::Mat h1, h2;
cv::stereoRectifyUncalibrated(points1, points2, fundamental,
 img1.size(), h1, h2);
// Rectify the images through warping
cv::Mat rectified1;
cv::warpPerspective(img1, rectified1, h1, img1.size());
cv::Mat rectified2;
cv::warpPerspective(img2, rectified2, h2, img1.size());

正如我之前在教程中解释的那样，在实际世界中，获得理想的相机配置而没有任何错误是非常困难的，因此 opencv 提供了一个校正功能，该功能应用单应变换将每个相机的图像平面投影到完美对齐的虚拟平面上. 这种变换是根据一组匹配点和基本矩阵计算得出的。

// Compute disparity
cv::Mat disparity;
cv::Ptr<cv::StereoMatcher> pStereo = cv::StereoSGBM::create(0,               32,5);
pStereo->compute(rectified1, rectified2, disparity);
cv::imwrite("disparity.jpg", disparity);

最后，我们计算了视差图。从下图中，较暗的像素代表离相机较近的物体，较亮的像素代表远离相机的物体。你在输出视差图中看到的白色像素噪声可以使用一些我不会在本教程中介绍的过滤器来去除。

现在我们已经成功地从给定的立体对中获得了深度图。现在让我们尝试使用 opencv 中名为 3D-Viz 的工具将获得的 2-D 图像点重新投影到 3-D 空间，该工具将帮助我们渲染 3-D 点云。

但是这一次，我们不是从给定的图像点估计一个基本矩阵，而是使用一个基本矩阵来投影这些点。

本质矩阵可以看作是基础矩阵，但用于校准的相机。我们也可以将其称为对基础矩阵的专业化，其中矩阵是使用校准的相机计算的，这意味着我们必须首先获取有关我们在世界上的相机的知识。

因此，为了让我们估计本质矩阵，我们首先需要相机的内在矩阵（表示给定相机的光学中心和焦距的矩阵）。让我们看一下下面的等式：

这里，从第一个矩阵中，fx和fy代表相机的焦距，(uo, vo)是主点。这是内在矩阵，我们的目标是估计它。

这种寻找不同相机参数的过程称为相机校准。我们显然可以使用相机制造商提供的规格，但是对于我们将要做的 3-D 重建等任务，这些规格不够准确。因此，我们将执行我们自己的相机校准。

这个想法是向相机显示一组场景点，这些点我们知道它们在现实世界中的实际 3-D 位置，然后观察这些点在获得的图像平面上的投影位置。有了足够数量的 3-D 点和相关的 2-D 图像点，我们就可以从投影方程中抽象出精确的相机参数。

做到这一点的一种方法是从不同的视点拍摄一组世界的 3-D 点及其已知 3-D 位置的多张图像。我们将使用 opencv 的校准方法，其中一种方法将棋盘图像作为输入，并返回所有存在的角。我们可以自由假设板位于 Z=0，X 和 Y 轴与网格很好地对齐。我们将在下面的部分中了解 OpenCV 的这些校准功能是如何工作的。

让我们首先创建三个函数，我们将在 main 函数中使用它们。这三个功能将

• addChessBoardPoints() //返回给定棋盘图像的角点
• calibrate() // 从提取的点返回内在矩阵
• triangulate() //返回重建点的 3-D 坐标

#include "CameraCalibrator.h"
#include <opencv2/opencv.hpp>
#include "opencv2/xfeatures2d.hpp"
using namespace std;
using namespace cv;
std::vector<cv::Mat> rvecs, tvecs;
// Open chessboard images and extract corner points
int CameraCalibrator::addChessboardPoints(
         const std::vector<std::string>& filelist, 
         cv::Size & boardSize) {
// the points on the chessboard
std::vector<cv::Point2f> imageCorners;
std::vector<cv::Point3f> objectCorners;
// 3D Scene Points:
// Initialize the chessboard corners 
// in the chessboard reference frame
// The corners are at 3D location (X,Y,Z)= (i,j,0)
for (int i=0; i<boardSize.height; i++) {
  for (int j=0; j<boardSize.width; j++) {
objectCorners.push_back(cv::Point3f(i, j, 0.0f));
   }
 }
// 2D Image points:
cv::Mat image; // to contain chessboard image
int successes = 0;
// for all viewpoints
for (int i=0; i<filelist.size(); i++) {
// Open the image
        image = cv::imread(filelist[i],0);
// Get the chessboard corners
        bool found = cv::findChessboardCorners(
                        image, boardSize, imageCorners);
// Get subpixel accuracy on the corners
        cv::cornerSubPix(image, imageCorners, 
                  cv::Size(5,5), 
                  cv::Size(-1,-1), 
      cv::TermCriteria(cv::TermCriteria::MAX_ITER +
                          cv::TermCriteria::EPS, 
             30,    // max number of iterations 
             0.1));     // min accuracy
// If we have a good board, add it to our data
      if (imageCorners.size() == boardSize.area()) {
// Add image and scene points from one view
            addPoints(imageCorners, objectCorners);
            successes++;
          }
//Draw the corners
        cv::drawChessboardCorners(image, boardSize, imageCorners, found);
        cv::imshow("Corners on Chessboard", image);
        cv::waitKey(100);
    }
return successes;
}

在上面的代码中，可以观察到我们包含了一个头文件“CameraCalibrator.h”，它将包含该文件的所有函数声明和变量初始化。可以通过访问此链接在我的 Github 上下载本教程中的标题以及所有其他文件。

我们的函数利用了 opencv 的findChessBoardCorners()函数，该函数将图像位置数组（数组必须包含每个棋盘图像的位置）和棋盘尺寸（你应该输入棋盘中水平和垂直角的数量）作为输入参数并返回给我们一个包含角点位置的向量。

double CameraCalibrator::calibrate(cv::Size &imageSize)
{
  // undistorter must be reinitialized
  mustInitUndistort= true;
// start calibration
  return 
     calibrateCamera(objectPoints, // the 3D points
          imagePoints,  // the image points
          imageSize,    // image size
          cameraMatrix, // output camera matrix
          distCoeffs,   // output distortion matrix
          rvecs, tvecs, // Rs, Ts 
          flag);        // set options
}

在这个函数中，我们使用了calibrateCamera()函数，它获取我们上面获得的 3-D 点和图像点，并返回给我们固有矩阵、旋转向量（描述相机相对于场景点的旋转）和平移矩阵（描述相机相对于场景点的位置）。

cv::Vec3d CameraCalibrator::triangulate(const cv::Mat &p1, const cv::Mat &p2, const cv::Vec2d &u1, const cv::Vec2d &u2) {
// system of equations assuming image=[u,v] and X=[x,y,z,1]
  // from u(p3.X)= p1.X and v(p3.X)=p2.X
  cv::Matx43d A(u1(0)*p1.at<double>(2, 0) - p1.at<double>(0, 0),
  u1(0)*p1.at<double>(2, 1) - p1.at<double>(0, 1),
  u1(0)*p1.at<double>(2, 2) - p1.at<double>(0, 2),
  u1(1)*p1.at<double>(2, 0) - p1.at<double>(1, 0),
  u1(1)*p1.at<double>(2, 1) - p1.at<double>(1, 1),
  u1(1)*p1.at<double>(2, 2) - p1.at<double>(1, 2),
  u2(0)*p2.at<double>(2, 0) - p2.at<double>(0, 0),
  u2(0)*p2.at<double>(2, 1) - p2.at<double>(0, 1),
  u2(0)*p2.at<double>(2, 2) - p2.at<double>(0, 2),
  u2(1)*p2.at<double>(2, 0) - p2.at<double>(1, 0),
  u2(1)*p2.at<double>(2, 1) - p2.at<double>(1, 1),
  u2(1)*p2.at<double>(2, 2) - p2.at<double>(1, 2));
cv::Matx41d B(p1.at<double>(0, 3) - u1(0)*p1.at<double>(2,3),
                p1.at<double>(1, 3) - u1(1)*p1.at<double>(2,3),
                p2.at<double>(0, 3) - u2(0)*p2.at<double>(2,3),
                p2.at<double>(1, 3) - u2(1)*p2.at<double>(2,3));
// X contains the 3D coordinate of the reconstructed point
  cv::Vec3d X;
  // solve AX=B
  cv::solve(A, B, X, cv::DECOMP_SVD);
  return X;
}

上述函数采用可以使用先前函数的固有矩阵获得的投影矩阵和归一化图像点，并返回上述点的 3-D 坐标。

下面是从立体对进行 3-D 重建的完整代码。此代码需要至少 25 到 30 张棋盘图像，这些图像来自你拍摄立体对图像的同一台相机。为了首先在你的 PC 中运行此代码，请克隆我的 GitHub 存储库，将双目对替换为自己的，并将棋盘图像位置数组替换为自己的数组，然后构建和编译。我正在将示例棋盘图像上传到我的 GitHub 以供你参考，你必须拍摄大约 30 张这样的图像并在代码中提及。

int main(){
cout<<"compiled"<<endl;
const std::vector<std::string> files = {"boards/1.jpg"......};
  cv::Size board_size(7,7);
CameraCalibrator cal;
  cal.addChessboardPoints(files, board_size);
cv::Mat img = cv::imread("boards/1.jpg");
cv::Size img_size = img.size();
  cal.calibrate(img_size);
  cout<<cameraMatrix<<endl;
cv::Mat image1 = cv::imread("imR.png");
  cv::Mat image2 = cv::imread("imL.png");
// vector of keypoints and descriptors
  std::vector<cv::KeyPoint> keypoints1;
  std::vector<cv::KeyPoint> keypoints2;
  cv::Mat descriptors1, descriptors2;
// Construction of the SIFT feature detector
  cv::Ptr<cv::Feature2D> ptrFeature2D = cv::xfeatures2d::SIFT::create(10000);
// Detection of the SIFT features and associated descriptors
  ptrFeature2D->detectAndCompute(image1, cv::noArray(), keypoints1, descriptors1);
  ptrFeature2D->detectAndCompute(image2, cv::noArray(), keypoints2, descriptors2);
// Match the two image descriptors
  // Construction of the matcher with crosscheck
  cv::BFMatcher matcher(cv::NORM_L2, true);
  std::vector<cv::DMatch> matches;
  matcher.match(descriptors1, descriptors2, matches);
cv::Mat matchImage;
cv::namedWindow("img1");
  cv::drawMatches(image1, keypoints1, image2, keypoints2, matches, matchImage, Scalar::all(-1), Scalar::all(-1), vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
  cv::imwrite("matches.jpg", matchImage);
// Convert keypoints into Point2f
  std::vector<cv::Point2f> points1, points2;
for (std::vector<cv::DMatch>::const_iterator it = matches.begin(); it != matches.end(); ++it) {
    // Get the position of left keypoints
    float x = keypoints1[it->queryIdx].pt.x;
    float y = keypoints1[it->queryIdx].pt.y;
    points1.push_back(cv::Point2f(x, y));
    // Get the position of right keypoints
    x = keypoints2[it->trainIdx].pt.x;
    y = keypoints2[it->trainIdx].pt.y;
    points2.push_back(cv::Point2f(x, y));
  }
// Find the essential between image 1 and image 2
  cv::Mat inliers;
  cv::Mat essential = cv::findEssentialMat(points1, points2, cameraMatrix, cv::RANSAC, 0.9, 1.0, inliers);
cout<<essential<<endl;
// recover relative camera pose from essential matrix
  cv::Mat rotation, translation;
  cv::recoverPose(essential, points1, points2, cameraMatrix, rotation, translation, inliers);
  cout<<rotation<<endl;
  cout<<translation<<endl;
// compose projection matrix from R,T
  cv::Mat projection2(3, 4, CV_64F); // the 3x4 projection matrix
  rotation.copyTo(projection2(cv::Rect(0, 0, 3, 3)));
  translation.copyTo(projection2.colRange(3, 4));
  // compose generic projection matrix
  cv::Mat projection1(3, 4, CV_64F, 0.); // the 3x4 projection matrix
  cv::Mat diag(cv::Mat::eye(3, 3, CV_64F));
  diag.copyTo(projection1(cv::Rect(0, 0, 3, 3)));
  // to contain the inliers
  std::vector<cv::Vec2d> inlierPts1;
  std::vector<cv::Vec2d> inlierPts2;
  // create inliers input point vector for triangulation
  int j(0);
  for (int i = 0; i < inliers.rows; i++) {
    if (inliers.at<uchar>(i)) {
      inlierPts1.push_back(cv::Vec2d(points1[i].x, points1[i].y));
      inlierPts2.push_back(cv::Vec2d(points2[i].x, points2[i].y));
    }
  }
  // undistort and normalize the image points
  std::vector<cv::Vec2d> points1u;
  cv::undistortPoints(inlierPts1, points1u, cameraMatrix, distCoeffs);
  std::vector<cv::Vec2d> points2u;
  cv::undistortPoints(inlierPts2, points2u, cameraMatrix, distCoeffs);
// Triangulation
  std::vector<cv::Vec3d> points3D;
  cal.triangulate(projection1, projection2, points1u, points2u, points3D);
cout<<"3D points :"<<points3D.size()<<endl;
viz::Viz3d window; //creating a Viz window
//Displaying the Coordinate Origin (0,0,0)
  window.showWidget("coordinate", viz::WCoordinateSystem());
window.setBackgroundColor(cv::viz::Color::black());
//Displaying the 3D points in green
  window.showWidget("points", viz::WCloud(points3D, viz::Color::green()));
  window.spin();
}

我知道 代码显示是一团糟，尤其是对于 C++，所以我建议你去我的GitHub并理解上面的代码。

对我来说，给定对的输出如下所示，可以通过调整特征检测器及其类型来改进。

任何有兴趣深入学习这些概念的人，我都会在下面推荐这本书，我认为这本书是计算机视觉几何的圣经。这也是本教程的参考书。

计算机视觉中的多视图几何第 2 版- Richard Hartley 和 Andrew Zisserman。

原文链接：三维重建C++实战 — BimAnt