本发明涉及三维重建领域，具体涉及一种快速鲁棒的rgb-d室内三维场景重建方法。
背景技术：
室内场景三维重建是计算机视觉中具有挑战性的研究课题之一，涉及计算机视觉、计算机图形学、模式识别、最优化等多个领域的理论与技术。实现三维重建有多种途径，传统方法是采用激光、雷达等测距传感器或结构光技术来获取场景或物体表面的结构信息进行三维重建，但这些仪器大多价格昂贵并且不易携带，所以应用场合有限。随着计算机视觉技术的发展，研究者们开始研究使用纯视觉的方法进行三维重建，涌现出来了大量有益的研究工作。
消费级rgb-d相机microsoftkinect推出后，人们可以直接利用深度数据比较便捷地进行室内场景三维重建。newcombe等人提出的kinectfusion算法利用kinect来获取图像中各点的深度信息，通过迭代近似最近邻点(iterativeclosestpoint,icp)算法将三维点在当前帧相机坐标系下的坐标与在全局模型中的坐标进行对齐来估计当前帧相机的姿态，再通过曲面隐函数(truncatedsigneddistancefunction,tsdf)迭代进行体数据融合，得到稠密的三维模型。虽然kinect获取深度不受光照条件和纹理丰富程度的影响，但其深度数据范围只有0.5-4m，而且网格模型的位置和大小是固定的，所以只适用于局部、静态的室内场景。
基于消费级rgb-d相机进行室内场景三维重建，一般存在以下几个问题：(1)消费级rgb-d相机获取的深度图像分辨率小、噪声大使得物体表面细节难以保持，而且深度值范围有限不能直接用于完整场景三维重建；(2)室内弱纹理区域几何结构简单、可以提取的点特征比较少，基于icp和特征法的相机定位比较困难；(3)利用装配rgb-d相机的机器人采集室内场景，会存在冗余重复的数据，影响建模速度和精度。
为了进行完整的室内场景三维重建，whelan等人提出了kintinuous算法，是对kinectfusion的进一步扩展。该算法使用shiftingtsdfvolume循环利用显存的方式解决大场景重建时网格模型显存消耗的问题，并通过dbow寻找匹配的关键帧进行闭环检测，最后对位姿图和模型做优化，从而得到大场景三维模型。choi和zeng等人提出了局部-全局配准的思想，先将rgb-d数据流做分段，对每段单独做视觉里程计估计，从两两段间的点云数据中提取描述子寻找匹配进行闭环检测，最后利用优化后的相机姿态进行体数据融合。elasticfusion使用surfel模型来表示三维模型，该算法通过优化deformationgraph的形式，提高重建和位姿估计的精度。bunderfusion基于稀疏特征、密集几何特征和光度匹配的对应关系，来估计相机姿态，可以实时得到高质量的场景模型。但这些算法都需要依赖于高配置的计算资源才能达到实时或耗时比较久。
近年来基于特征法的rgb-d相机定位逐渐引起了关注。orb-slam2是一套基于orb特征进行定位，高精度、高效的稀疏slam(simultaneouslocalizationandmapping)系统，可以实时运行在普通cpus上，但在室内弱纹理区域很容易定位失败。由于室内场景中存在大量的线特征，不少系统基于直线特征进行相机定位以提高在弱纹理区域的相机定位鲁棒性，但是直线特征的提取和匹配都比较耗时。edgevo直接采用边缘匹配，可以非常高效地进行相机跟踪，但定位精度很差。edgeslam在该算法的基础上对初始化做了优化，并加入了闭环控制，提高了定位精度，不过也降低了效率。
基于tsdf的体数据融合方法可以方便地将多幅rgb-d图像融合为完整的三维模型。steinbrucker等人提出了采用八叉树数据结构存储tsdf值的算法，可以高效、实时的进行数据融合。nissner等人采用哈希数据结构来压缩空间，可以实现体数据的实时读取与更新。不过这些算法没有考虑rgb-d相机的噪声特点和数据冗余的情况，并不能完全保证模型表面的几何细节。
技术实现要素：
为了解决现有技术中的上述问题，即为了解决室内三维场景重建效率不能满足需求的问题，本发明提出了一种快速鲁棒的rgb-d室内三维场景重建方法，该方法包括：
采用rgb-d相机实时扫描室内场景；
基于点线融合相机跟踪算法进行相机实时跟踪；
通过检测相机状态对相机轨迹进行约减；
利用约减后的相机轨迹信息对rgb-d图像进行多尺度体数据融合，生成完整的场景三维模型。
在本发明的优选实施方式中，“基于点线融合相机跟踪算法进行相机实时跟踪”包括三个线程：跟踪线程、局部构图线程、闭环检测线程；
所述跟踪线程，基于rgb-d相机采集到的彩色图像与深度图像序列，进行初始姿态估计、跟踪局部地图、关键帧的判断；其中，所述初始状态估计通过松耦合的方式联合特征点跟踪与边缘跟踪，通过最小化能量函数估计相机姿态；
所述局部构图线程，用于对每个新的关键帧完成局部地图构建；
所述闭环检测线程，用于对局部地图构建过程中最后处理的关键帧进行闭环检测，并在确定闭环后进行闭环修正。
在本发明的优选实施方式中，所述初始状态估计中，通过松耦合的方式联合特征点跟踪与边缘跟踪估计相机姿态所采用的最小化能量函数为：
其中，ep是特征点的重投影误差，ee是边缘的几何距离误差，χc表示相邻帧的匹配集合，ρp和ρe是huber函数，λ∈[0,1]是权重因子，b为相邻帧图像间的变换关系。
在本发明的优选实施方式中，特征点的重投影误差ep为：
ep＝x'-kbx
其中，x表示orb特征点的三维坐标，x'表示匹配点的位置，k为相机内参矩阵，b为相邻帧图像间的变换关系；
边缘的几何距离误差ee为：
ee＝d(l',τ(k,b,l))
其中，l表示边缘上像素的坐标，l'表示匹配点的位置，d表示几何距离误差，τ(k,b,l)为边缘上像素的坐标l在相邻帧图像中对应的三维点坐标；
权重因子λ基于场景纹理的丰富程度进行取值：
其中，n表示每帧图像上提取的orb特征数量，nmax和nmin是最大、最小阈值。
在本发明的优选实施方式中，所述初始状态估计中，进行边缘匹配之前设置有像素点深度值修正步骤，包括：
将像素点深度值与预设深度区间进行对比，如果该数值不在所述预设深度区间，则采用ekf进行深度估计，得到新的深度值。
在本发明的优选实施方式中，“通过检测相机状态对相机轨迹进行约减”，其方法为：
计算相邻帧之间相机运动的线速度；
如果线速度低于设定第一阈值则删除该帧位姿，否则计算该帧图像与闭环帧图像的相似度；
如果相似度低于设定第二阈值，则删除该帧位姿；
输出约减后的图像帧位姿。
在本发明的优选实施方式中，“计算该帧图像与闭环帧图像的相似度”，其方法为：
步骤301，根据像素点与空间三维点的投影关系、第i帧深度图像的深度值，计算深度图像上每个像素所对应的空间三维坐标，计算公式如下：
p＝π-1(x,z(x))
其中，xp是图像上的任一像素，对应的深度值和空间三维点分别为z(x)和p，π表示2d-3d投影变换关系；
步骤302，将步骤301计算得到空间三维点p从第i帧相机坐标系做旋转平移变换到世界坐标系下，得到新的三维点坐标q，计算公式为：
q＝bip
其中，bi为相机坐标系与世界坐标系的变换关系，p和q的坐标分别表示为p＝(px,py,pz)，q＝(qx,qy,qz)。
步骤303，将步骤302计算得到三维点q反投影到二维图像平面，得到新的深度图，计算公式为：
其中，xq是q对应的新得到深度图上的像素，(u0,v0)为相机主点，(fu,fv)为相机焦距；
步骤304，分别计算闭环帧深度图有效像素个数n0、第i帧新生成深度图上的有效像素个数ni，按照如下公式计算相似度ρ：
在本发明的优选实施方式中，“利用约减后的相机轨迹信息对rgb-d图像进行多尺度体数据融合”，其方法为：利用tsdf网格模型融合各帧的图像，使用分辨率为m的体素网格来表示三维空间；其中，m为空间被分割的块数。
在本发明的优选实施方式中，每个网格o存储两个值：截断符号距离函数fi(o)及其权重wi(o)；
截断符号距离函数fi(o)的公式为：
其中，fi(o)表示该网格到物体模型表面的距离，正、负分别表示该网格是在表面被遮挡一侧、可见一侧，过零点就是体素表面上的点；[]z表示z轴的取值；φ是截断阈值，k为相机内参矩阵，zi(x)为第i个网格中像素x在深度图对应的深度值；
权重wi(o)的公式为：
其中，δ的取值为体素分辨率m的十分之一。
在本发明的优选实施方式中，在多尺度体数据融合过程中，采用多分辨率八叉树数据结构来存储tsdf值。
在本发明的优选实施方式中，八叉树数据结构中数据块的尺度sl为：
其中，pi(xi,yi,zi)表示空间三维点坐标，该点到光轴的平方距离为
在本发明的优选实施方式中，对于帧数为n的图像序列，进行多尺度体数据融合时，按照下式进行迭代：
f(o)和w(o)表示融合后体素o所对应的tsdf值及其权重。
本发明的有益效果：
本发明中，点线融合的相机跟踪算法可以提高系统的跟踪鲁棒性，轨迹约减算法能有效降低数据融合的速度和模型精度，多尺度体数据融合算法可以高效、准确地得到三维模型；基于cpus的室内场景三维重建系统，能够高效、准确地得到完整室内场景模型，系统具有良好的鲁棒性和扩展性。
附图说明
图1是室内场景三维重建系统示意图；
图2是对彩色图进行orb特征提取和边缘提取的效果比较图；
图3是轨迹约减算法流程示意图；
图4是轨迹约减前后体数据融合的时间比较示意图；
图5是室内场景数据轨迹约减前后的模型细节比较图示意图；
图6是空间三维点与多分辨率八叉树数据结构的关系对应图；
图7是在tumrgb-ddataset数据集上使用本方法得到的三维重建效果图；
图8是对microsoftkinectforwindows采集的场景三维重建的效果比较图。
具体实施方式
下面参照附图来描述本发明的优选实施方式。本领域技术人员应当理解的是，这些实施方式仅仅用于解释本发明的技术原理，并非旨在限制本发明的保护范围。
本发明实施例的快速鲁棒的rgb-d室内三维场景重建方法，如图1所示，包括以下步骤：
采用rgb-d相机实时扫描室内场景；
基于点线融合相机跟踪算法进行相机实时跟踪；
通过检测相机状态对相机轨迹进行约减；
利用约减后的相机轨迹信息对rgb-d图像进行多尺度体数据融合，生成完整的场景三维模型。
根据上述方法，可以将本发明分为相机跟踪阶段、轨迹约减阶段、多尺度体数据融合阶段三个阶段进行详细描述。
1、相机跟踪阶段
相机跟踪阶段分为三个线程：跟踪线程、局部构图线程、闭环检测线程。
(1)跟踪线程，基于rgb-d相机采集到的彩色图像与深度图像序列，进行初始姿态估计、跟踪局部地图、关键帧的判断；其中，所述初始状态估计通过松耦合的方式联合特征点跟踪与边缘跟踪，通过最小化能量函数估计相机姿态。
该线程的输入为rgb-d相机采集到的彩色图像与深度图像序列，输出为初始姿态估计和跟踪局部地图，并根据orb-slam方法的间隔帧数、跟踪点数和视图变化要求，对当前帧进行关键帧的判断。其中，初始姿态估计包括边缘跟踪和特征点跟踪，通过最小化能量函数得到相机姿态。局部地图包含与当前帧有共同地图点的关键帧集合k1，以及在互见图(covisibilitygraph)中与k1相邻的关键帧集合k2。互见图的每个节点是一个关键帧，如果两个关键帧之间有共同看到的地图点，那么两者之间就会生成一条边，共同看到的地图点的个数就是该条边的权值。跟踪局部地图就是寻找与这些关键帧中地图点的匹配，通过motion-onlyba进行位姿优化。跟踪局部地图、关键帧判断的方法在较多技术文献、论文中均有描述，此处不再赘述。
(2)局部构图线程，用于对每个新的关键帧完成局部地图构建。在跟踪线程插入关键帧后，优化局部地图中的点和位姿。在关键帧插入后，会结合局部地图内的另一帧创建新的地图点。局部构图线程的具体方法在较多技术文献、论文中均有描述，此处不再赘述。
(3)闭环检测线程：通过dbow2字典树进行闭环检测，当检测到闭环时，计算闭环帧与当前帧的se(3)变换，并通过位姿图的优化纠正累积误差和地图点的位置。闭环检测线程的具体方法在较多技术文献、论文中均有描述，此处不再赘述。
本发明实施例中为了提高初始化和边缘匹配的效率，我们使用深度图中的深度值进行辅助。由于基于结构光原理的消费级rgb-d相机在结构边缘处采集不到深度，所以可以在进行边缘匹配之前设置有像素点深度值修正步骤，将像素点深度值与预设深度区间进行对比，如果该数值在预设深度区间之内就直接用来做匹配；否则在采用ekf(extendedkalmanfilter)估计新的像素点深度值，再用来做匹配。在该重建系统中预设深度区间为(0.5m，4m)。
2、轨迹约减阶段
轨迹约减阶段输入为相机跟踪阶段得到的所有图像帧的相机姿态，根据约减算法去除冗余帧的位姿，以提高建模精度和建模速度。约减算法流程如图2所示，包括速度检测和相似度检测两部分，分别处理运动缓慢和轨迹重复造成的数据冗余。具体包括以下步骤：
计算相邻帧之间相机运动的线速度；
如果线速度低于设定第一阈值则删除该帧位姿，否则计算该帧图像与闭环帧图像的相似度；
如果相似度低于设定第二阈值，则删除该帧位姿；
输出约减后的图像帧位姿。
3、多尺度体数据融合阶段
多尺度体数据融合阶段利用volumetricmethod框架进行tsdf数据加权融合，tsdf数值的存储采用多分辨率八叉树数据结构，来进一步提高体数据融合效率。并根据rgb-d相机的噪声特点来定义tsdf的参数以实现物体表面几何细节的保持。
本发明实施例的三个阶段分别包含三个核心算法：点线融合的相机跟踪算法、轨迹约减算法、多尺度体数据融合算法，下面将具体展开进行说明。
1、点线融合的相机跟踪算法
点线融合的相机跟踪算法包括：特征点跟踪、边缘跟踪和点线融合跟踪。特征点跟踪我们统一采用orb特征，通过最小化重投影误差来估计相机姿态。边缘跟踪的特点是寻找相邻帧边缘图的匹配函数，通过最小化边缘上匹配像素的几何距离误差来估计相机姿态。点线融合的相机跟踪采用松耦合的方式联合特征点跟踪与边缘跟踪，得到相机姿态。
边缘匹配的目标就是得到相邻帧图像间的变换关系b，我们用六维向量ξ参数化表示该变换，变换函数的推导步骤如下：
(1)通过反投影变换关系，计算第一帧边缘图上像素x＝(u,v)t对应的三维点p坐标如公式(1)所示：
其中，z(x)为像素x在深度图对应的深度值，π表示2d-3d投影变换关系，(u0,v0)为相机主点，(fu,fv)为相机焦距，k为相机内参矩阵：
(2)通过变换得到第二帧图像上对应三维点坐标为b(g(ξ),p，)其中g(ξ)表示se(3)变换。相机观测该点得到的三维点坐标为q＝(qx,qy,qz)，所以我们得到变换关系如公式(2)所示：
联立公式(1)和(2)，最终得到变换函数如公式(3)所示：
τ(k,ξ,x)＝π(b(g(ξ),p))＝π(b(g(ξ),π-1(k,x,z(x))))(3)
通过最小化第一帧图像上像素点和第二帧图像上匹配到的像素点之间的几何距离误差，来求解两帧之间的变换关系，几何距离误差定义如公式(4)所示：
其中，
(3)通过松耦合的方式联合特征点跟踪与边缘跟踪估计相机姿态所采用的最小化能量函数如式(5)所示：
其中，ep是特征点的重投影误差，ee是边缘的几何距离误差，χc表示相邻帧的匹配集合，ρp和ρe是huber函数，λ∈[0,1]是权重因子。
特征点的重投影误差如公式(6)所示：
ep＝x'-kbx(6)
其中，x表示orb特征点的三维坐标，x'表示匹配点的位置。
边缘的几何距离误差如公式(7)所示：
ee＝d(l',τ(k,b,l))(7)
其中，l表示边缘上像素的坐标；l'表示匹配点的位置；d()表示几何距离误差,采用公式(4)计算；τ(k,b,l)为边缘上像素的坐标l在相邻帧图像中对应的三维点坐标，可以通过公式(3)进行计算。
权重因子λ基于场景纹理的丰富程度进行取值，如公式(8)所示：
其中，n表示每帧图像上提取的orb特征数量，nmax和nmin是最大、最小阈值。
为了提高初始化和边缘匹配效率，我们使用深度图中的深度值进行辅助，具体如上述的像素点深度值修正步骤。
2、轨迹约减算法：
通过相机跟踪阶段，可以得到每帧图像对应的相机姿态以及相机坐标系与世界坐标系的变换矩阵bi＝[ri|ti]，其中ri表示旋转关系，ti表示平移关系。对于输入为n帧的图像序列，轨迹约减算法流程为：首先检测每帧图像相对于上帧图像的运动状态，计算相邻帧之间相机运动的线速度；如果线速度低于设定第一阈值则删除该帧位姿，否则计算该帧图像与闭环帧图像的相似度；如果相似度低于设定第二阈值，则删除该帧位姿；最后输出约减后的图像帧位姿。
结合图3对轨迹约减算法流程进行详细说明如下：
步骤201，获取n个图像帧的位姿，定义变量i，进行初始化赋值i＝0；
步骤202，在i值小于n时，计算相邻帧之间相机运动的线速度，否则跳转步骤205；
步骤203，判断线速度是否小于设定第一阈值，若小于则删除该帧位姿并跳转步骤205，否则计算该帧图像与闭环帧图像的相似度；
步骤204，判断相似度是否小于设定第二阈值，若小于则删除该帧位姿并跳转步骤205，否则变量i＝i+1并返回步骤202继续下一个帧图像的判断；
步骤205，最后输出约减后的图像帧位姿。
相邻帧之间相机运动的线速度计算如公式(9)所示：
vi＝ti-ti-1(9)
相似度计算包括四个步骤：
步骤301，根据像素点与空间三维点的投影关系、第i帧深度图像的深度值，计算深度图像上每个像素所对应的空间三维坐标，如公式(10)所示：
p＝π-1(x,z(x))(10)
其中，x是图像上的任一像素，对应的深度值和空间三维点分别为z(x)和p；该公式中π变换隐含了相机参数矩阵k，相比于公式(1)无需在公式中显示k矩阵。
步骤302，将步骤301计算得到空间三维点p从第i帧相机坐标系做旋转平移变换到世界坐标系下，得到新的三维点坐标q，如公式(11)所示：
q＝bip(11)
其中，p和q的坐标分别为p＝(px,py,pz)，q＝(qx,qy,qz)；
步骤303，将步骤302计算得到三维点q反投影到二维图像平面，得到新的深度图，如公式(12)所示：
其中，xq是q对应的新得到深度图上的像素。
步骤304，分别计算闭环帧深度图有效像素个数n0、第i帧新生成深度图上的有效像素个数ni，按照公式(13)计算相似度ρ：
图4是轨迹约减前后体数据融合的时间比较，横轴是图像序列的帧数，纵轴是数据融合用到的时间(秒s)。
图5是室内场景数据轨迹约减前后的模型细节比较图示意图，图中方框标出的区域是存在闭环和缓慢运动的区域。
3、多尺度体数据融合算法：
根据轨迹约减后的相机轨迹信息，利用tsdf网格模型融合各帧的图像，使用分辨率为m的体素网格来表示三维空间。其中，m为空间被分割的块数，每个网格o存储两个值：截断符号距离函数fi(o)及其权重wi(o)。
截断符号距离函数fi(o)的定义如公式(14)所示：
fi(o)＝max{min{φ,|[k-1zi(x)[xt,1]t]z-[o]z|},-φ}(14)
其中，fi(o)表示该网格到物体模型表面的距离，正、负分别表示该网格是在表面被遮挡一侧、可见一侧，过零点就是体素表面上的点；[]z表示z轴的取值，zi(x)为第i个网格中像素x在深度图对应的深度值；φ是预设的截断阈值，其取值与相机参数有关。
权重wi(o)的计算如公式(15)所示：
其中，δ的取值在本实施例中可以为体素分辨率m的十分之一。
为了高效地进行体数据融合，结合消费级rgb-d相机的特点，采用多分辨率八叉树数据结构来存储tsdf值。空间三维点位置与八叉树数据结构的对应关系如图6所示，pi(xi,yi,zi)表示空间三维点坐标，该点到光轴的平方距离为di2＝xi2+yi2，八叉树结构中数据块的尺度sl如公式(16)所示：
图6中示出了四个空间三维点坐标pr、py、pg、pb，对应的到光轴的具体为dr、dy、dg、db，可以看出深度值zb＝zg＞zy＞zr，距离db＞dg＞dy＝dr，得到的尺度sb＞sg＞sy＞sr，同一坐标点的距离和尺度在图中采用相同线型进行表示，重叠部份仅显示尺度较小的线型，如图6右图中不同线型构成的四个立方体。
对于帧数为n的图像序列，进行数据融合时，按照公式(17)进行迭代：
f(o)和w(o)表示融合后体素o所对应的tsdf值及其权重。
图7是在tumrgb-ddataset数据集上使用本方法得到的三维重建效果图；上图为fr3/structure_notexture_near序列，是带有结构特征的弱纹理近景场景；下图为fr3/structure_notexture_far序列，是带有结构特征的弱纹理远景场景；
图8是对microsoftkinectforwindows采集的场景三维重建的效果比较图，从上到下依次为zhou方法的重建结果，运行时间：6000songpus；elasticfusion的重建结果，运行时间：135songpus；本发明的重建结果，运行时间：51soncpus。
结合本文中所公开的实施例描述的方法或算法的步骤可以用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(ram)、内存、只读存储器(rom)、电可编程rom、电可擦除可编程rom、寄存器、硬盘、可移动磁盘、cd-rom、或技术领域内所公知的任意其它形式的存储介质中。
为了清楚地说明电子硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以电子硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。本领域技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。
术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法所固有的要素。
至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。