OpenCV实现SfM（二）：双目三维重建[通俗易懂]

大家好，又见面了，我是你们的朋友全栈君。

注意：本文中的代码必须使用OpenCV3.0或以上版本进行编译，因为很多函数是3.0以后才加入的。
目录：

在三维重建前，我们先研究一下同一点在两个相机中的像的关系。假设在世界坐标系中有一点$p$，坐标为$X$，它在1相机中的像为$x_1$，在2相机中的像为$x_2$（注意$x_1$和$x_2$为齐次坐标，最后一个元素是1），如下图。

设$X$到两个相机像面的垂直距离分别为$s_1$和$s_2$，且这两个相机具有相同的内参矩阵$K$，与世界坐标系之间的变换关系分别为$[R_1{T}_1]$和$[R_2{T}_2]$，那么我们可以得到下面两个等式
$$\begin{gathered} s_1{x}_1=K(R_{1}X+T_1) \\ {s}_2{x}_2=K(R_{2}X+T_2) \end{gathered}$$
由于K是可逆矩阵，两式坐乘K的逆，有
$$\begin{gathered} s_1{K}^{-1}{x}_1=R_{1}X+T_1 \\ {s}_2{K}^{-1}{x}_2=R_{2}X+T_2 \end{gathered}$$
设$K^{-1}{x}_1=x_1^{{}^{\prime }}$，$K^{-1}{x}_2=x_2^{{}^{\prime }}$，则有
$$\begin{gathered} s_1{x}_1^{{}^{\prime }}=R_{1}X+T_1 \\ {s}_2{x}_2^{{}^{\prime }}=R_{2}X+T_2 \end{gathered}$$
我们一般称$x_1^{{}^{\prime }}$和$x_2^{{}^{\prime }}$为归一化后的像坐标，它们和图像的大小没有关系，且原点位于图像中心。
由于世界坐标系可以任意选择，我们将世界坐标系选为第一个相机的相机坐标系，这时$R_1=I,T_1=0$。上式则变为
$$\begin{gathered} s_1{x}_1^{{}^{\prime }}=X \\ {s}_2{x}_2^{{}^{\prime }}=R_{2}X+T_2 \end{gathered}$$
将第一式带入第二式，有
$$s_2{x}_2^{{}^{\prime }}=s_1{R}_2{x}_1^{{}^{\prime }}+T_2$$
$x_2^{{}^{\prime }}$和$T_2$都是三维向量，它们做外积（叉积）之后得到另外一个三维向量$\widehat{T_2}{x}_2^{{}^{\prime }}$（其中$\widehat{T_2}$为外积的矩阵形式，$\widehat{T_2}{x}_2^{{}^{\prime }}$代表$T_2\times x_2^{{}^{\prime }}$），且该向量垂直于$x_2^{{}^{\prime }}$和$T_2$，再用该向量对等式两边做内积，有
$$0=s_1(\widehat{T_2}{x}_2^{{}^{\prime }}{)}^T{R}_2{x}_1^{{}^{\prime }}$$
即
$$x_2^{{}^{\prime }}\widehat{T_2}{R}_2{x}_1^{{}^{\prime }}=0$$
令$E=\widehat{T_2}{R}_2$ 有
$$x_2^{{}^{\prime }}E{x}_1^{{}^{\prime }}=0$$
可以看出，上式是同一点在两个相机中的像所满足的关系，它和点的空间坐标、点到相机的距离均没有关系，我们称之为极线约束，而矩阵$E$则称为关于这两个相机的本征矩阵。如果我们知道两幅图像中的多个对应点（至少5对），则可以通过上式解出矩阵$E$，又由于$E$是由$T_2$和$R_2$构成的，可以从E中分解出$T_2$和$R_2$。
如何从$E$中分解出两个相机的相对变换关系（即$T_2$和$R_2$），背后的数学原理比较复杂，好在OpenCV为我们提供了这样的方法，在此就不谈原理了。

#特征点提取与匹配
从上面的分析可知，要求取两个相机的相对关系，需要两幅图像中的对应点，这就变成的特征点的提取和匹配问题。对于图像差别较大的情况，推荐使用SIFT特征，因为SIFT对旋转、尺度、透视都有较好的鲁棒性。如果差别不大，可以考虑其他更快速的特征，比如SURF、ORB等。
本文中使用SIFT特征，由于OpenCV3.0将SIFT包含在了扩展部分中，所以官网上下载的版本是没有SIFT的，为此需要到这里下载扩展包，并按照里面的说明重新编译OpenCV（哎~真麻烦，-_-!）。如果你使用其他特征，就不必为此辛劳了。
下面的代码负责提取图像特征，并进行匹配。

void extract_features(
vector<string>& image_names,
vector<vector<KeyPoint>>& key_points_for_all,
vector<Mat>& descriptor_for_all,
vector<vector<Vec3b>>& colors_for_all
)
{ 

key_points_for_all.clear();
descriptor_for_all.clear();
Mat image;
//读取图像，获取图像特征点，并保存
Ptr<Feature2D> sift = xfeatures2d::SIFT::create(0, 3, 0.04, 10);
for (auto it = image_names.begin(); it != image_names.end(); ++it)
{ 

image = imread(*it);
if (image.empty()) continue;
vector<KeyPoint> key_points;
Mat descriptor;
//偶尔出现内存分配失败的错误
sift->detectAndCompute(image, noArray(), key_points, descriptor);
//特征点过少，则排除该图像
if (key_points.size() <= 10) continue;
key_points_for_all.push_back(key_points);
descriptor_for_all.push_back(descriptor);
vector<Vec3b> colors(key_points.size());
for (int i = 0; i < key_points.size(); ++i)
{ 

Point2f& p = key_points[i].pt;
colors[i] = image.at<Vec3b>(p.y, p.x);
}
colors_for_all.push_back(colors);
}
}
void match_features(Mat& query, Mat& train, vector<DMatch>& matches)
{ 

vector<vector<DMatch>> knn_matches;
BFMatcher matcher(NORM_L2);
matcher.knnMatch(query, train, knn_matches, 2);
//获取满足Ratio Test的最小匹配的距离
float min_dist = FLT_MAX;
for (int r = 0; r < knn_matches.size(); ++r)
{ 

//Ratio Test
if (knn_matches[r][0].distance > 0.6*knn_matches[r][1].distance)
continue;
float dist = knn_matches[r][0].distance;
if (dist < min_dist) min_dist = dist;
}
matches.clear();
for (size_t r = 0; r < knn_matches.size(); ++r)
{ 

//排除不满足Ratio Test的点和匹配距离过大的点
if (
knn_matches[r][0].distance > 0.6*knn_matches[r][1].distance ||
knn_matches[r][0].distance > 5 * max(min_dist, 10.0f)
)
continue;
//保存匹配点
matches.push_back(knn_matches[r][0]);
}
}

需要重点说明的是，匹配结果往往有很多误匹配，为了排除这些错误，这里使用了Ratio Test方法，即使用KNN算法寻找与该特征最匹配的2个特征，若第一个特征的匹配距离与第二个特征的匹配距离之比小于某一阈值，就接受该匹配，否则视为误匹配。当然，也可以使用Cross Test（交叉验证）方法来排除错误。

得到匹配点后，就可以使用OpenCV3.0中新加入的函数findEssentialMat()来求取本征矩阵了。得到本征矩阵后，再使用另一个函数对本征矩阵进行分解，并返回两相机之间的相对变换R和T。注意这里的T是在第二个相机的坐标系下表示的，也就是说，其方向从第二个相机指向第一个相机（即世界坐标系所在的相机），且它的长度等于1。

bool find_transform(Mat& K, vector<Point2f>& p1, vector<Point2f>& p2, Mat& R, Mat& T, Mat& mask)
{ 

//根据内参矩阵获取相机的焦距和光心坐标（主点坐标）
double focal_length = 0.5*(K.at<double>(0) + K.at<double>(4));
Point2d principle_point(K.at<double>(2), K.at<double>(5));
//根据匹配点求取本征矩阵，使用RANSAC，进一步排除失配点
Mat E = findEssentialMat(p1, p2, focal_length, principle_point, RANSAC, 0.999, 1.0, mask);
if (E.empty()) return false;
double feasible_count = countNonZero(mask);
cout << (int)feasible_count << " -in- " << p1.size() << endl;
//对于RANSAC而言，outlier数量大于50%时，结果是不可靠的
if (feasible_count <= 15 || (feasible_count / p1.size()) < 0.6)
return false;
//分解本征矩阵，获取相对变换
int pass_count = recoverPose(E, p1, p2, R, T, focal_length, principle_point, mask);
//同时位于两个相机前方的点的数量要足够大
if (((double)pass_count) / feasible_count < 0.7)
return false;
return true;
}

#三维重建
现在已经知道了两个相机之间的变换矩阵，还有每一对匹配点的坐标。三维重建就是通过这些已知信息还原匹配点在空间当中的坐标。在前面的推导中，我们有
$$s_2{x}_2=K(R_{2}X+T_2)$$
这个等式中有两个未知量，分别是$s_2$和$X$。用$x_2$对等式两边做外积，可以消去$s_2$，得
$$0=\widehat{x_2}K(R_{2}X+T_2)$$
整理一下可以得到一个关于空间坐标X的线性方程
$$\widehat{x_2}K{R}_{2}X=-\widehat{x_2}K{T}_2$$
上面的方程不能直接取逆求解，因此化为其次方程
$$\widehat{x_2}K(R_{2}T)\left(\begin{array}{c}X\\ 1\end{array}\right)=0$$
用SVD求X左边矩阵的零空间，再将最后一个元素归一化到1，即可求得X。其几何意义相当于分别从两个相机的光心作过$x_1$和$x_2$的延长线，延长线的焦点即为方程的解，如文章最上方的图所示。由于这种方法和三角测距类似，因此这种重建方式也被称为三角化（triangulate）。OpenCV提供了该方法，可以直接使用。

void reconstruct(Mat& K, Mat& R, Mat& T, vector<Point2f>& p1, vector<Point2f>& p2, Mat& structure)
{ 

//两个相机的投影矩阵[R T]，triangulatePoints只支持float型
Mat proj1(3, 4, CV_32FC1);
Mat proj2(3, 4, CV_32FC1);
proj1(Range(0, 3), Range(0, 3)) = Mat::eye(3, 3, CV_32FC1);
proj1.col(3) = Mat::zeros(3, 1, CV_32FC1);
R.convertTo(proj2(Range(0, 3), Range(0, 3)), CV_32FC1);
T.convertTo(proj2.col(3), CV_32FC1);
Mat fK;
K.convertTo(fK, CV_32FC1);
proj1 = fK*proj1;
proj2 = fK*proj2;
//三角化重建
triangulatePoints(proj1, proj2, p1, p2, structure);
}

#测试
我用了下面两幅图像进行测试

得到了着色后的稀疏点云，是否能看出一点轮廓呢？！

图片中的两个彩色坐标系分别代表两个相机的位置。
在接下来的文章中，会将相机的个数推广到任意多个，成为一个真正的SfM系统。

关于源代码的使用
代码是用VS2013写的，OpenCV版本为3.0且包含扩展部分，如果不使用SIFT特征，可以修改源代码，然后使用官方未包含扩展部分的库。软件运行后会将三维重建的结果写入Viewer目录下的structure.yml文件中，在Viewer目录下有一个SfMViewer程序，直接运行即可读取yml文件并显示三维结构。

代码下载