注意：本文中的代码必须使用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]$，那么我们可以得到下面两个等式

$s_1x_1 = K(R_1X + T_1) \\ s_2x_2 = K(R_2X + T_2)$

由于K是可逆矩阵，两式坐乘K的逆，有

$s_1K^{-1}x_1 = R_1X + T_1 \\ s_2K^{-1}x_2 = R_2X + T_2$

设$K^{-1}x_1 = x_1^{'}$，$K^{-1}x_2 = x_2^{'}$，则有

$s_1x_1^{'} = R_1X + T_1 \\ s_2x_2^{'} = R_2X + T_2$

我们一般称$x_1^{'}$和$x_2^{'}$为归一化后的像坐标，它们和图像的大小没有关系，且原点位于图像中心。

由于世界坐标系可以任意选择，我们将世界坐标系选为第一个相机的相机坐标系，这时$R_1 = I,\ T_1 = 0$。上式则变为

$s_1x_1^{'} = X \\ s_2x_2^{'} = R_2X + T_2$

将第一式带入第二式，有

$s_2x_2^{'} = s_1R_2x_1^{'} + T_2$

$x_2^{'}$和$T_2$都是三维向量，它们做外积（叉积）之后得到另外一个三维向量$\widehat{T_2}x_2^{'}$（其中$\widehat{T_2}$为外积的矩阵形式，$\widehat{T_2}x_2^{'}$代表$T_2\times x_2^{'}$），且该向量垂直于$x_2^{'}$和$T_2$，再用该向量对等式两边做内积，有

$0 = s_1(\widehat{T_2}x_2^{'})^TR_2x_1^{'}$

即

$x_2^{'}\widehat{T_2}R_2x_1^{'} = 0$

令$E = \widehat{T_2}R_2$ 有

$x_2^{'}Ex_1^{'} = 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_2x_2 = K(R_2X + T_2)$

这个等式中有两个未知量，分别是$s_2$和$X$。用$x_2$对等式两边做外积，可以消去$s_2$，得

$0 = \widehat{x_2}K(R_2X+T_2)$

整理一下可以得到一个关于空间坐标X的线性方程

$\widehat{x_2}KR_2X = -\widehat{x_2}KT_2$

上面的方程不能直接取逆求解，因此化为其次方程

$\widehat{x_2}K(R_2\ \ T)\left(\begin{matrix}X \\ 1\end{matrix}\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文件并显示三维结构。
代码下载

reference: https://zhuanlan.zhihu.com/p/114047559

（1）传统的三维重建技术



（2）基于深度学习的三维重建算法研究主要有三种：

在传统三维重建算法中引入深度学习方法进行改进；
深度学习重建算法和传统三维重建算法进行融合，优势互补；
模仿动物视觉，直接利用深度学习算法进行三维重建，包括基于体素、基于点云和基于网格。

前言

       三维重建技术通过深度数据获取、预处理、点云配准与融合、生成表面等过程，把真实场景刻画成符合计算机逻辑表达的数学模型。这种模型可以对如文物保护、游戏开发、建筑设计、临床医学等研究起到辅助的作用。

1.1 研究背景及意义

人类通过双眼来探索与发现世界。人类接收外部信息的方式中，有不到三成来自于听觉、触觉、嗅觉等感受器官，而超过七成、最丰富、最复杂的信息则通过视觉[1]进行感知的。计算机视觉便是一种探索给计算机装备眼睛（摄像头）与大脑（算法）的技术，以使计算机能够自主独立的控制行为、解决问题，同时感知、理解、分析外部环境。20世纪60年代，计算机视觉得到了最初的发展。该阶段的研究重心主要体现在如何从二维图像中恢复出如立方体、圆柱体等立体化的三维形状，解释各个物体的空间位置关系。1982年，David Marr[2]从信息处理的角度对数学、神经生理学、计算机图形学等学科的研究成果进行了归纳总结，并在此基础上提出了一系列计算机视觉理论。得益于这个完整明确的理论体系，计算机视觉得到了蓬勃的发展。它的核心思想是从二维图像恢复三维结构。图1-1展示的是经典Marr视觉信息处理过程。 




 

Fig.1-1 Process of Marr visual information 

随着科学技术的日新月异，计算机视觉的应用日益受到各行业的关注和重视，如设备检测与监视、医学图像处理、文物保护[3]、机器人视觉、自动导航、工业产品外观设计与生产等领域。计算机视觉技术为人们带来了机遇，也带来了挑战。三维重建作为计算机视觉技术中最为最为热门的研究方向之一，涉及到包括图像处理、立体视觉、模式识别等多个学科体系。利用计算机建立表达现实客观景物的三维模型，并以此来满足生产和生活的需要。随着工业化进程的不断发展，多种技术的实现均有赖于目标物体三维信息的获取。三维重建现已被广发的应用于生活和科研工作中，特别是在医学治疗、文物保护、游戏开发、工业设计、航天航海等方面，展现出了极强的生命力和影响力。

1.2 三维重建技术简介

三维重建技术的重点在于如何获取目标场景或物体的深度信息。在景物深度信息已知的条件下，只需要经过点云数据[4]的配准及融合，即可实现景物的三维重建。基于三维重建模型的深层次应用研究也可以随即展开。人们按照被动式测量与主动式测量[5]对目标物体深度信息的获取方法进行了分类，下面对这两种方式进行相应的介绍。

1.2.1 被动式三维重建技术

被动式一般利用周围环境如自然光的反射，使用相机获取图像，然后通过特定算法计算得到物体的立体空间信息。主要有以下三种方法：

1.纹理恢复形状法

各种物体表面具有不同的纹理信息，这种信息由纹理元组成，根据纹理元可以确定表面方向，从而恢复出相应的三维表面。这种方法称为纹理恢复形状法[6] (Shape From Texture，SFT)。 

纹理法的基本理论为：作为图像视野中不断重复的视觉基元，纹理元覆盖在各个位置和方向上。当某个布满纹理元的物体被投射在平面上时，其相应的纹理元也会发生弯折与变化。例如透视收缩变形使与图像平面夹角越小的纹理元越长，投影变形会使离图像平面越近的纹理元越大。通过对图像的测量来获取变形，进而根据变形后的纹理元，逆向计算出深度数据。SFT对物体表面纹理信息的要求严苛，需要了解成像投影中纹理元的畸变信息，应用范围较窄，只适合纹理特性确定等某些特殊情形。所有在实际使用中较为少见。

2.阴影恢复形状法

SFS[7] (Shape From Shading，从阴影恢复形状)法也是一种较为常用的方法。考虑到图像的阴影边界包含了图像的轮廓特征信息，因此能够利用不同光照条件下的图像的明暗程度与阴影来计算物体表面的深度信息，并以反射光照模型进行三维重建。需要注意的是，像素点的亮度受到包括光源指标、摄像机参数、目标表面材质等的制约。 

阴影恢复形状法的应用范围比较广泛，可以恢复除镜面外的各种物体的三维模型。缺点体现在过程多为数学计算、重建结果不够精细，另外不能忽视的是，SFS法需要准确的光源参数，包括位置与方向信息。这就导致其无法应用于诸如露天场景等具有复杂光线的情形中。

3.立体视觉法

立体视觉法[8]（Multi-View Stereo，MVS)是另外一种常用的三维重建方法。主要包括直接利用测距器获取程距信息、通过一幅图像推测三维信息和利用不同视点上的两幅或多幅图像恢复三维信息等三种方式。通过模拟人类视觉系统，基于视差原理获取图像对应点之间的位置偏差，恢复出三维信息。S.T.Barnard[9]等人对20世纪70年代到80年代之间出现的三维重建的算法和评价体系做了概述。到了80年代中后期，出现了更多、更深层次的视觉原理，包括立体测量方法和深度传感器等，极大的促进了相关学科的发展。新兴方法可以直接获取景物的三维信息，极大的节省了物力与人力成本。U.R.Dhond[10]等人提出了基于层次处理的三目立体约束方法。二十世纪90年代末，涌现出诸如图像匹配的前沿算法、遮挡处理算法等。M.Z.Brown[11]等人总结了2000年到2010年间的三维视觉发展的总体概况，包括遮挡、配准和效率等的相关分析。 

双目立体视觉重建，在实际应用情况优于其他基于视觉的三维重建方法，也逐渐出现在一部分商业化产品上; 不足的是运算量仍然偏大，而且在基线距离较大的情况下重建效果明显降低 。 

代表文章：AKIMOIO T Automatic creation of 3D facial models 1993 

CHEN C L Visual binocular vison systems to solid model reconstruction2007 

作为计算机视觉的关键技术之一，立体视觉法也其弊端。例如，立体视觉需要假设空间的平面是正平面，而实际情况却与此相差甚远。除此之外，匹配还存在歧义性：对于一幅图像上的某些特征点，另外的图像可能存在若干个与之相似的特征点。那么如何选取最适配的匹配点，显得较为棘手。如图1-2所示，展示了Middlebury[16]数据集中Teddy和Cones场景的基准彩色图像、标准视差以及通过Graph Cuts[17]算法获取的立体匹配视差估计结果。虽然视差结果体现出了景物的三维位置关系，但是某些像素点的视差与标准值仍有细微的差距。除此之外，对于如相机运动参数的确定、大型场景重建需要获取多帧图像等问题，也极大的影响了立体视觉的深层次应用。 




 

 




 

 

参考：立体匹配导论

1.2.2 主动式三维重建技术

主动式是指利用如激光、声波、电磁波等光源或能量源发射至目标物体，通过接收返回的光波来获取物体的深度信息。主动测距有莫尔条纹法、飞行时间法、结构光法和三角测距法等四种方法。

1.莫尔条纹法

莫尔条纹在生活中比较常见，如两层薄薄的丝绸重叠在一起，即可以看到不规则的莫尔(Morie)条纹；微风的吹动窗纱时，条纹亦随之运动。莫尔条纹法[18]起源于18世纪的法国，是一项古老又现代的测量方法。基本原理是将两块等间隔排列的直线簇或曲线簇图案重叠起来，以非常小的角度进行相对运动来形成莫尔条纹。如图1-3所示，在主光栅与指示光栅的交叉重合处，因光线的透射与遮挡而产生不同的明暗带，即莫尔条纹。莫尔条纹随着光栅的左右平移而发生垂直位移，此时产生的条纹相位信息体现了待测物体表面的深度信息，再通过逆向的解调函数，实现深度信息的恢复。这种方法具有精度高、实时性强的优点，但是其对光照较为敏感，抗干扰能力弱。 




 

提出：WIKTIN recovering surface shape and orientation from texture (1987)（被引用454 次）。 

发展：Warren 2010 对 wiktin 方法进行改进使用了透视投影； 

Liboy 2006 给出了在纹理单元结构发生改变的情况下的重建方法。 

优点：精度高，对光照和噪声不敏感。 

缺点：只应用于具有规则纹理的物体。

2.飞行时间法

飞行时间法[19] (Time of Flight，ToF)指的是在光速及声速一定的前提下，通过测量发射信号与接收信号的飞行时间间隔来获得距离的方法。这种信号可以是超声波，也可以是红外线等。飞行时间法相较于立体视觉法而言，具有不受基线长度限制、与纹理无关、成像速度快等特点。但是其也有一定的缺点。首先，ToF相机的分辨率非常低。例如图1-4所示，当今分辨率最高的PMD Camcube 2.0 相机，也仅为204×204像素；其次，ToF相机容易受到环境因素的影响，如混合像素、外界光源等，导致景物深度不准确；最后，系统误差与随机误差对测量结果的影响很大，需要进行后期数据处理，主要体现在场景像素点的位置重合上。值得注意的是，ToF相机的售价达到了数万美元，受众较窄。 




 

 

Fig.1-4 SR4000 ToF camera

3.结构光法

结构光法[20](Structured Light)通过向表面光滑无特征的物体发射具有特征点的光线，依据光源中的立体信息辅助提取物体的深度信息。具体的过程包括两个步骤，首先利用激光投影仪向目标物体投射可编码的光束，生成特征点；然后根据投射模式与投射光的几何图案，通过三角测量原理计算摄像机光心与特征点之间的距离，由此便可获取生成特征点的深度信息，实现模型重建。这种可编码的光束就是结构光，包括各种特定样式的点、线、面等图案。结构光法解决了物体表面平坦、纹理单一、灰度变化缓慢等问题。因为实现简单且精度较高，所以结构光法的应用非常广泛，目前已有多家公司生产了以结构光技术为基础的硬件设备，如PrimeSense公司的Prime Sensor、微软公司的Kinect和华硕公司的Xtion PRO LIVE等产品[21]。图1-5展示了利用结构光技术采集文物三维信息的场景。 

提出：Woodham 对 SFS 进行改进（1980 年）：photometric method for determining surface orientation from multiple images (该文章被引用了 891 次) 

发展：Noakes ：非线性与噪声减除 2003 年； 

Horocitz ：梯度场合控制点 2004 年； 

Tang ： 可信度传递与马尔科夫随机场 2005 年； 

Basri ： 光源条件未知情况下的三维重建 2007 年； 

Sun ：非朗伯特 2007 年； 

Hernandez ： 彩色光线进行重建方法 2007 年； 

Shi ： 自标定的光度立体视觉法 2010 年。 




 

4.三角测距法

三角测距法[22]是一种非接触式的测距方法，以三角测量原理为基础。红外设备以一定的角度向物体投射红外线，光遇到物体后发生反射并被CCD（Charge-coupled Device，电荷耦合元件）图像传感器所检测。随着目标物体的移动，此时获取的反射光线也会产生相应的偏移值。根据发射角度、偏移距离、中心矩值和位置关系，便能计算出发射器到物体之间的距离。三角测距法在军工测量、地形勘探等领域中应用广泛。

参考文献

[1] Szeliski R. Computer vision: algorithms and applications[M]. Berlin: Springer, 2010. 

[2] D. Marr, et al. A Computational Theory of Human Stereo Vision. Proc.R.Soc.Lond. 1979, B.204:301-328. 

[3] Levoy, M. Pulli, et al. The Digital Michelangelo Project:3D Scanning of Large Statues. Proc.SIGGRAPH,2000. 

[4] Anand A, Koppula H S, Joachims T, et al. Contextually guided semantic labeling and search for three-dimensional point clouds[J]. The International Journal of Robotics Research, 2013, 32(1):19-34. 

[5] Mada S K, Smith M L, Smith L N, et al. Overview of passive and active vision techniques for hand-held 3D data acquisition [C]//Opto Ireland. International Society for Optics and Photonics, 2003: 16-27. 

[6] D. A. Forsyth, J. Ponce, Computer Vision: A Modern Approach. Prentice Hall 2001 

[7] Horn B. K. P. Shape from shading: a method for obtaining the shape of a smooth opaque object from one view. PhD thesis, Department of Electrical Engineering, MIT, Cambridge. 1970. 

[8] Ikeuchi K. Determining surface orientations of specular surfaces by using the photometric stereo method [J]. Pattern Analysis and Machine Intelligence, IEEE Transactions, 1981, (6): 661-669. 

[9] S. T. Barnard, M. A. Fisehler. Computational Stereo[J].ACM Computing Surveys. 1982, Vol.14:553-572. 

[10] U. R. Dhond, J. K. Aggarval. Struct from Stereo—A Review [J]. IEEE Trans. Systems, Man, and Cybemeties.1989, Vol.19: 1489-1510.

基于视觉的三维重建，指的是通过摄像机获取场景物体的数据图像，并对此图像进行分析处理，再结合计算机视觉知识推导出现实环境中物体的三维信息。

1. 相关概念

（1）彩色图像与深度图像

彩色图像也叫作RGB图像，R、G、B三个分量对应于红、绿、蓝三个通道的颜色，它们的叠加组成了图像像素的不同灰度级。RGB颜色空间是构成多彩现实世界的基础。深度图像又被称为距离图像，与灰度图像中像素点存储亮度值不同，其像素点存储的是该点到相机的距离，即深度值。图2-1表示深度图像与灰度图像之间的关系。




 


Fig.2-1 The depth image and gray image

深度值指的目标物体与测量器材之间的距离。由于深度值的大小只与距离有关，而与环境、光线、方向等因素无关，所以深度图像能够真实准确的体现景物的几何深度信息。通过建立物体的空间模型，能够为深层次的计算机视觉应用提供更坚实的基础。




 


Fig.2-2 Color image and depth image of the characters

（2）PCL

PCL（Point Cloud Library，点云库）是由斯坦福大学的Dr.Radu等学者基于ROS(Robot Operating System，机器人操作系统)下开发与维护的开源项目，最初被用来辅助机器人传感、认知和驱动等领域的开发。2011年PCL正式向公众开放。随着对三维点云算法的加入与扩充，PCL逐步发展为免费、开源、大规模、跨平台的C++编程库。

PCL框架包括很多先进的算法和典型的数据结构，如滤波、分割、配准、识别、追踪、可视化、模型拟合、表面重建等诸多功能。能够在各种操作系统和大部分嵌入式系统上运行，具有较强的软件可移植性。鉴于PCL的应用范围非常广，专家学者们对点云库的更新维护也非常及时。PCL的发展时至今日，已经来到了1.7.0版本。相较于早期的版本，加入了更多新鲜、实用、有趣的功能，为点云数据的利用提供了模块化、标准化的解决方案。再通过诸如图形处理器、共享存储并行编程、统一计算设备架构等领先的高性能技术，提升PCL相关进程的速率，实现实时性的应用开发。

在算法方面，PCL是一套包括数据滤波、点云配准、表面生成、图像分割和定位搜索等一系列处理点云数据的算法。基于不同类型区分每一套算法，以此把整合所有三维重建流水线功能，保证每套算法的紧凑性、可重用性与可执行性。例如PCL中实现管道运算的接口流程：

①创建处理对象，例如滤波、特征估计、图像分割等； 

②通过setInputCloud输入初始点云数据，进入处理模块； 

③设置算法相关参数； 

④调用不同功能的函数实现运算，并输出结果。

为了实现模块化的应用与开发，PCL被细分成多组独立的代码集合。因此便可方便快捷的应用于嵌入式系统中，实现可移植的单独编译。如下列举了部分常用的算法模块：

libpcl I/O：完成数据的输入、输出过程，如点云数据的读写； 

libpcl filters：完成数据采样、特征提取、参数拟合等过程； 

libpcl register：完成深度图像的配准过程，例如迭代最近点算法； 

libpcl surface：完成三维模型的表面生成过程，包括三角网格化、表面平滑等。

此类常用的算法模块均具有回归测试功能，以确保使用过程中没有引进错误。测试一般由专门的机构负责编写用例库。检测到回归错误时，会立即将消息反馈给相应的作者。因此能提升PCL和整个系统的安全稳定性。

（3）点云数据

如图2-3所示，展示了典型的点云数据（Point Cloud Data，PCD）模型。




 

点云数据通常出现在逆向工程中，是由测距设备获取的物体表面的信息集合。其扫描资料以点的形式进行记录，这些点既可以是三维坐标，也可以是颜色或者光照强度等信息。通常所使用的点云数据一般包括点坐标精度、空间分辨率和表面法向量等内容。点云一般以PCD格式进行保存，这种格式的点云数据可操作性较强，同时能够提高点云配准融合的速度。本文研究的点云数据为非结构化的散乱点云，属于三维重建特有的点云特点。

（4）坐标系 

在三维空间中，所有的点必须以坐标的形式来表示，并且可以在不同的坐标系之间进行转换。首先介绍基本坐标系的概念、计算及相互关系。

①图像坐标系

图像坐标系分为像素和物理两个坐标系种类。数字图像的信息以矩阵形式存储，即一副像素的图像数据存储在维矩阵中。图像像素坐标系以为原点、以像素为基本单位，U、V分别为水平、垂直方向轴。图像物理坐标系以摄像机光轴与图像平面的交点作为原点、以米或毫米为基本单位，其X、Y轴分别与U、V轴平行。图2-4展示的是两种坐标系之间的位置关系：




 

 

Fig.2-4 Image pixel coordinate system and physical coordinate system

令U-V坐标系下的坐标点(u0，v0)，与代表像素点在X轴与Y轴上的物理尺寸。那么图像中的所有像素点在U-V坐标系与在X-Y坐标系下的坐标间有着如式（2-1）表示的关系： 



其中指的是图像坐标系的坐标轴倾斜相交而形成的倾斜因子（Skew Factor）。

②摄像机坐标系 

摄像机坐标系由摄像机的光心及三条、、轴所构成。它的、轴对应平行于图像物理坐标系中的、轴，轴为摄像机的光轴，并与由原点、、轴所组成的平面垂直。如图2-5所示：




 

令摄像机的焦距是f，则图像物理坐标系中的点与摄像机坐标系中的点的关系为：

③世界坐标系

考虑到摄像机位置具有不确定性，因此有必要采用世界坐标系来统一摄像机和物体的坐标关系。世界坐标系由原点及、、三条轴组成。世界坐标与摄像机坐标间有着（2-3）所表达的转换关系：



 





( 23 ) 

其中，是旋转矩阵，代表摄像机在世界坐标系下的指向；是平移向量，代表了摄像机的位置。

2.三维重建流程

本文使用Kinect采集景物的点云数据，经过深度图像增强、点云计算与配准、数据融合、表面生成等步骤，完成对景物的三维重建。




 


Fig.2-6 Flow chart of 3D reconstruction based on depth sensor

图2-6显示的流程表明，对获取到的每一帧深度图像均进行前六步操作，直到处理完若干帧。最后完成纹理映射。下面对每个步骤作详细的说明。

2.1 深度图像的获取

景物的深度图像由Kinect在Windows平台下拍摄获取，同时可以获取其对应的彩色图像。为了获取足够多的图像，需要变换不同的角度来拍摄同一景物，以保证包含景物的全部信息。具体方案既可以是固定Kinect传感器来拍摄旋转平台上的物体；也可以是旋转Kinect传感器来拍摄固定的物体。价格低廉、操作简单的深度传感器设备能够获取实时的景物深度图像，极大的方便了人们的应用。

2.2 预处理

受到设备分辨率等限制，它的深度信息也存在着许多缺点。为了更好的促进后续基于深度图像的应用，必须对深度图像进行去噪和修复等图像增强过程。作为本文的重点问题，具体的处理方法将在第四章进行详细的解释说明。

2.3 点云计算

经过预处理后的深度图像具有二维信息，像素点的值是深度信息，表示物体表面到Kinect传感器之间的直线距离，以毫米为单位。以摄像机成像原理为基础，可以计算出世界坐标系与图像像素坐标系之间具有下式的转换关系：



 



则k值只与有关，而等参数只与摄像机的内部构造有关，所以称为像机的内参数矩阵。以摄像机作为世界坐标系，即，则深度值即为世界坐标系中的值，与之对应的图像坐标就是图像平面的点。

2.4 点云配准

对于多帧通过不同角度拍摄的景物图像，各帧之间包含一定的公共部分。为了利用深度图像进行三维重建，需要对图像进行分析，求解各帧之间的变换参数。深度图像的配准是以场景的公共部分为基准，把不同时间、角度、照度获取的多帧图像叠加匹配到统一的坐标系中。计算出相应的平移向量与旋转矩阵，同时消除冗余信息。点云配准除了会制约三维重建的速度，也会影响到最终模型的精细程度和全局效果。因此必须提升点云配准算法的性能。

三维深度信息的配准按不同的图像输入条件与重建输出需求被分为：粗糙配准、精细配准和全局配准等三类方法。

（1）粗糙配准（Coarse Registration）

粗糙配准研究的是多帧从不同角度采集的深度图像。首先提取两帧图像之间的特征点，这种特征点可以是直线、拐点、曲线曲率等显式特征，也可以是自定义的符号、旋转图形、轴心等类型的特征。

随后根据特征方程实现初步的配准。粗糙配准后的点云和目标点云将处于同一尺度(像素采样间隔)与参考坐标系内，通过自动记录坐标，得到粗匹配初始值。

（2）精细配准（Fine Registration）

精细配准是一种更深层次的配准方法。经过前一步粗配准，得到了变换估计值。将此值作为初始值，在经过不断收敛与迭代的精细配准后，达到更加精准的效果。以经典的由Besl和Mckay[49]提出的ICP（Iterative Closest Point，迭代最近点）算法为例，该算法首先计算初始点云上所有点与目标点云的距离，保证这些点和目标点云的最近点相互对应，同时构造残差平方和的目标函数。

基于最小二乘法对误差函数进行最小化处理，经过反复迭代，直到均方误差小于设定的阈值。ICP算法能够获得精正确无误的配准结果，对自由形态曲面配准问题具有重要意义。另外还有如SAA（Simulate Anneal Arithmetic，模拟退火）算法、GA（Genetic Algorithm，遗传）算法等也有各自的特点与使用范畴。

（3）全局配准（Global Registration）

全局配准是使用整幅图像直接计算转换矩阵。通过对两帧精细配准结果，按照一定的顺序或一次性的进行多帧图像的配准。这两种配准方式分别称为序列配准（Sequential Registration）和同步配准（Simultaneous Registration）。

配准过程中，匹配误差被均匀的分散到各个视角的多帧图像中，达到削减多次迭代引起的累积误差的效果。值得注意的是，虽然全局配准可以减小误差，但是其消耗了较大的内存存储空间，大幅度提升了算法的时间复杂度。

2.5 数据融合

经过配准后的深度信息仍为空间中散乱无序的点云数据，仅能展现景物的部分信息。因此必须对点云数据进行融合处理，以获得更加精细的重建模型。以Kinect传感器的初始位置为原点构造体积网格，网格把点云空间分割成极多的细小立方体，这种立方体叫做体素(Voxel)。通过为所有体素赋予SDF（Signed Distance Field，有效距离场）值，来隐式的模拟表面。

SDF值等于此体素到重建表面的最小距离值。当SDF值大于零，表示该体素在表面前；当SDF小于零时，表示该体素在表面后；当SDF值越接近于零，表示该体素越贴近于场景的真实表面。KinectFusion技术虽然对场景的重建具有高效实时的性能，但是其可重建的空间范围却较小，主要体现在消耗了极大的空间用来存取数目繁多的体素。

为了解决体素占用大量空间的问题，Curless[50]等人提出了TSDF (Truncated Signed Distance Field，截断符号距离场)算法，该方法只存储距真实表面较近的数层体素，而非所有体素。因此能够大幅降低KinectFusion的内存消耗，减少模型冗余点。




 

TSDF算法采用栅格立方体代表三维空间，每个栅格中存放的是其到物体表面的距离。TSDF值的正负分别代表被遮挡面与可见面，而表面上的点则经过零点，如图2-7中左侧展示的是栅格立方体中的某个模型。若有另外的模型进入立方体，则按照下式(2-9)与(2-10)实现融合处理。

其中，指的是此时点云到栅格的距离，是栅格的初始距离，是用来对同一个栅格距离值进行融合的权重。如图2-7中右侧所示，两个权重之和为新的权重。对于KinectFusion算法而言，当前点云的权重值设置为1。

鉴于TSDF算法采用了最小二乘法进行了优化，点云融合时又利用了权重值，所有该算法对点云数据有着明显的降噪功能。

2.6 表面生成

表面生成的目的是为了构造物体的可视等值面，常用体素级方法直接处理原始灰度体数据。Lorensen[51]提出了经典体素级重建算法：MC（Marching Cube，移动立方体）法。移动立方体法首先将数据场中八个位置相邻的数据分别存放在一个四面体体元的八个顶点处。对于一个边界体素上一条棱边的两个端点而言，当其值一个大于给定的常数T，另一个小于T时，则这条棱边上一定有等值面的一个顶点。

然后计算该体元中十二条棱和等值面的交点，并构造体元中的三角面片，所有的三角面片把体元分成了等值面内与等值面外两块区域。最后连接此数据场中的所有体元的三角面片，构成等值面。合并所有立方体的等值面便可生成完整的三维表面。

3 性能优化

Kinect等深度传感器的出现，不仅给娱乐应用带来了变革，同样对科学研究提供了新的方向。尤其是在三维重建领域。然而由于三维重建过程涉及到大量密集的点云数据处理，计算量巨大，所以对系统进行相应的性能优化显得非常的重要。本文采用基于GPU（Graphic Processing Unit，图形处理器）并行运算功能，以提高整体的运行效率。

NVIDIA公司于1999年提出了GPU概念。在这十几年间，依靠硬件行业的改革创新，芯片上晶体管数量持续增多，GPU性能以半年翻一番的速度成倍提升。GPU的浮点运算能力远超CPU上百倍，却具有非常低的能耗，极具性价比。因GPU不仅广泛应用于图形图像处理中，也在如视频处理、石油勘探、生物化学、卫星遥感数据分析、气象预报、数据挖掘等方面崭露头角。

作为GPU的提出者，NVIDIA公司一直致力于GPU性能提升的研究工作，并在2007年推出了CUDA架构。CUDA（Compute Unified Device Architecture，统一计算设备架构）是一种并行计算程序架构。在CUDA的支持下，使用者可以编写程序以利用NVIDIA系列GPU完成大规模并行计算。GPU在CUDA中被用作通用计算设备，而不只是处理图像。在CUDA中，将计算机CPU称为主机（Host），GPU称为设备（Device）。

主机端和设备端都有程序运行，主机端主要完成程序的流程与串行计算模块，而设备端则专门处理并行计算。其中，设备端的并行计算过程被记录在Kernel内核函数中，主机端可以从Kernel函数入口执行并行计算的调用功能。在此过程中，虽然Kernel函数执行同一代码，但却处理着不同的数据内容。

Kernel函数采用扩展的C语言来编程，称为CUDAC语言。需要注意的是，并不是所有的运算都可以采用CUDA并行计算。只有独立性的计算，如矩阵的加减，因为只涉及到对应下标的元素的加减，不同下标元素毫无关联，所以适用于并行计算；而对于如阶乘的计算则必须对所有数累积相乘，故无法采用并行计算。

CUDA具有线程（Thread）、程序块（Block）、网格（Grid）三级架构，计算过程一般由单一的网格完成，网格被平均分成多个程序块，每个程序块又由多个线程组成，最终由单个线程完成每个基本运算，如图2-8所示。 




 

为了更深入的理解CUDA模型的计算过程，这里以前一章中提到的公式（2-11）为例，计算某点的深度值与三维坐标之间的转换： 



 

上式中的表示深度值，内参数矩阵是已知量，是该点的坐标。可以发现这个点的转换过程与其他点转换过程是相互独立的，所以整幅图像中各点的坐标转换能够并行执行。这种并行计算可以大幅提升整体计算的速率。例如，利用一个网格来计算一幅像素的深度图像到三维坐标的转换，只需要将此网格均分成块，每块包括个线程，每个线程分别操作一个像素点，便可以便捷的完成所有的坐标转换运算。

通过GPU的并行计算，三维重建性能得到了大幅的提升，实现了实时的输入输出。对于Kinect在实际生产生活中的应用奠定了基础。

小结

首先介绍了与三维重建相关的基本概念，包括深度图像、点云数据、四种坐标系及其之间的转换关系等。