融合尺度与信心指数的三维表面重建算法

2024-02-21李亚兰陆汝华黄健全蒋纯志

计算机集成制造系统 2024年1期

李亚兰,陆汝华,黄健全,蒋纯志,李翔+

(1.湘南学院物理与电子电气工程学院微电子与光电子技术湖南省普通高等学校重点实验室,湖南郴州 423000;2.先进嵌入式计算技术与智能医疗系统湖南省工程研究中心,湖南郴州 423000;3.湘南学院计算机与人工智能学院,湖南郴州 423000)

0 引言

随着计算机技术的发展,三维表面重建在逆向工程[1]、工业设计制造与工业视觉测量[2-3]、历史遗迹恢复、旅游和影视娱乐[4]等众多领域受到越来越多的关注,成为了当前计算机领域研究的热点问题之一。早期基于体积的距离图像处理(Volumetric Range Image Processing,VRIP)算法[5]通过假设某一深度值下的众多深度观测值服从高斯分布,在已知众多观测值的条件下求取物体表面的最大似然解,该算法以增量的方式处理不断增加的深度图,对于深度图中误差进行平滑,不易重现表面精细几何结构。泊松表面重构(Poisson Surface Reconstruction,PSR)算法[6]直接处理带有方向的三维点云,通过隐函数来表示三维表面,使隐函数上三维采样点处的梯度场近似对应输入三维点云的方向场,并将该近似问题通过散度操作转换为空间泊松问题来处理,该算法重建的三维模型存在部分伪封闭表面,细小结构有平滑过度的现象。文献[7]结合三维表面三角面片周长和三维表面采样点与输入点之间的距离值对伪表面进行分割,还原了真实表面。屏蔽泊松表面重构(Screened Poisson Surface Reconstruction,SPSR)算法[8]对PSR算法[6]的泊松方程增加了位置约束项,鼓励重建出的表面穿过输入三维点的坐标位置,在一定程度上缓解了PSR算法[6]的过度平滑效应。光滑有向距离表面重建(Smooth Signed Distance surface reconstruction,SSD)算法[9]则在PSR算法[6]的基础上使用变分的方法通过有向符号场函数求解物体的光滑近似表面。文献[10]对PSR算法[6]、SPSR算法[8]和SSD算法[9]进行了系统的分析,并提出增加可视外壳约束条件来提高户外复杂场景重建质量的方法。文献[11]为每一个三维采样点引入了新的尺度信息,根据尺度信息自适应调整三维表面近似函数,以抑制点云中的噪声。三维点云作为三维表面重建的输入信号,包含了丰富的三维几何信息同时也存在各类噪声,充分挖掘有用信息,尽量抑制噪声影响是提高三维表面重建质量的主要途径。上述表面重建算法输入点云的信息挖掘利用尚存在扩展空间,本文充分融合点云信心指数与尺度信息,构建了三维表面基函数与权重函数,以八叉树结构组织三维点云并自适应调整三维采样密度,抑制噪声影响,还原精细结构,提高三维表面的重建质量。

1 算法原理

本文构造的三维表面隐函数主要包括基函数和权重函数两部分,联合点云尺度与信心指数自动调整基函数覆盖范围和权重值大小,并以八叉树数据组织形式自适应调整表面隐函数的采样密度。

1.1 表面隐函数构造

1.1.1 整体函数构造

在VRIP算法[5]的启发下构造了表面隐函数。VRIP算法[5]通过使式(1)成立,求得场景三维表面最大似然解。

(1)

其中:D(x)是物体某一三维点p上的累积有向符号函数,x为p点世界坐标,di(x)是从第i幅深度图中获取的p点有向符号函数,ωi(x)是从第i幅深度图获取的p点权重函数。

构造了以输入三维点云为对象的物体表面连续隐函数如式(2):

(2)

其中:F(x)是物体某一点p的隐函数,bi(x)是p点邻域点pi的基函数,wi(x)是p点邻域点pi的权重函数。

1.1.2 基函数构造

将点pi的世界坐标xiw转换到以点p为原心,p点法向方法为x轴方向的局部坐标系,得其局部坐标为:

xil=Ri(xiw-x)。

(3)

其中:坐标xiw、xil和x分别为pi点的世界坐标,局部坐标和p点的世界坐标。Ri为pi点从世界坐标变换到局部坐标的旋转变换矩阵。

pi点的基函数设为式(4):

(4)

其中:si和ci为pi点尺度和信心指数;xil、yil和zil为pi点局部坐标系下的3个坐标分量。基函数的选取既要充分覆盖pi点周围的影响区域,又不能将区域延展太大而过渡平滑数据特征。为此,选择延展性较好的高斯函数,并用pi点的信心指数ci和尺度si来控制高斯函数的方差,以自适应调整pi点周围基函数覆盖区域大小。信心指数高,尺度大则可以适当扩大基函数范围;信心指数低,尺度小则相应缩小基函数范围。

pi点在x轴方向的基函数为高斯函数的导数:

(5)

yoz平面与p点法向方向垂直,y轴和z轴可互换,为了保证yoz平面旋转不变,可设pi点在y轴方向和z轴方向的基函数为相同的形式,如式(6)和式(7):

(6)

(7)

1.1.3 权重函数构造

pi点权重函数为:

wi=wix(xil)wiy(yil)wiz(zil)。

(8)

其中wix、wiy和wiz为pi点在x、y和z轴方向权重。wix需采用非对称形式,如图1a所示,其中x轴方向为物体表面点法向方向,正半轴为表面外向,负半轴为表面内向。由于物体表面内的点是不可观测的,为了抑制表面内点的影响,表面内点权重值随内部深度增加快速下降,而物体表面外点权重相对取值较高。如图2所示为wix权重快速下降的两种情况,图2中箭头方向为物体表面p点法向方向,即为物体表面外向。在图2a中,由于遮挡的原因,pi位于p点所处表面内部,一方面,在重构p点周围表面时,pi点会导致该重构表面偏离真实表面分布情况;另一方面,由于遮挡,pi点所处位置的不确定性增高,从而降低p点重构表面精度及pi点权重wix。在图2b中,由于物体很薄,使得位于物体上、下两表面的p点、pi点距离很近。为了降低下表面pi点对上表面p点周围表面重建的影响,降低pi点权重wix。如图1a所示表面内、外点权重都随|xil|增大而减小,因为pi点离p越远,则pi点对p点的影响越小。当xil范围扩大到[-3sici,3sici]以外时,对应的基函数值已经非常小。因此,在[-3sici,3sici]以外将权重设为0。图1a所示的非对称函数是通过设置三阶贝塞尔函数的4个控制点得到的,左侧曲线4个控制点坐标为(0,1)、(0.55,0.96)、(0.29,0.05)、(1,0);右侧曲线的4个控制点坐标为(-1,0)、(-0.34,0.06)、(-0.37,0.3)、(0,1)。

图1 权重函数曲线

图2 两种wix取低值的情况

如图1b所示,因为yoz平面与法方向垂直,所以在yoz平面上采用对称的权重函数。由于yoz平面旋转不变,y轴和z轴可互换,y轴方向与z轴方向采用相同的权重函数曲线。图1b所示的对称函数是通过设置三阶贝塞尔函数的4个控制点得到的,4个控制点坐标为(0,1)、(0.36,0.97)、(0.57,0.03)、(1,0)。与wix一样,距离越远的点影响越小,信心指数越低的点重建精越低,这些点的权重也需相应降低,反之则权重值相应增加。

1.2 八叉树数据组织结构

三维点云数据是通过八叉树的数据结构结合点云尺度进行存储和采样的。

1.2.1 数据存储

三维点云的存储规则如图3和图4所示。首先,以第一个插入的三维点p1为根节点的中心rc,以p1点的尺度s1为边长l1建立根节点,如图3a中红色立方体所示。然后,从第2个三维点p2开始,按该点在根节点内、外两种情况插入八叉树中,若待插入点p2在根结点外,如图3a所示,且p2点x、y和z坐标值均大于p1点x、y和z坐标值,则以当前根节点为基准向p2点的方向(即x、y和z轴的正方向)生长形成新的根节点,如图3b所示。并使新根节点(包含红色立方体在内,体积为红色立方体的8倍)边长是旧根节点(红色立方体)的两倍。同时将旧根节点中心rc移至新根节点中心位置,即图3b中蓝色实心点位置。若生长一次还不能包含p2点,则重复上述生长过程直至新根结点包含p2点,如图3b所示。若待插入点已经在当前根结点内,则根据不等式(9)选择插入节点级数,即插入的第i级节点边长li要小于待插入点尺度si,大于等于0.5si,这样就可以使小尺度点插入细节点即高级数节点,大尺度点插入粗节点即低级数节点。由于相邻级数节点边长按一半递减规律划分,不等式(9)又可转换为不等式(10),即插入点尺度si大于插入级数节点边长li,小于等于低一级节点边长2li。

图3 待插入点在根节点外

图4 待插入点在根节点内

0.5si≤li

(9)

(10)

如图4a所示,待插入点p3已经在当前根节点内。则根据p3点的位置和尺度s3大小,将p3点放置在边长小于p3点尺度s3的叶节点(灰色小立方体)内,如图4b所示。插入操作如下:将p3点尺度s3与当前根节点的边长l1进行比较,由于s3l3,所以将p3点放在对应象限的第三级叶节点(灰色小立方体)中。而其他象限叶节点仍保留为第二级,这样根据插入点云尺度自适应对八叉树进行划分,节省了存储空间,并且为后续物体表面以自适应分辨率采样隐函数提供基础。按以上规律,将所有三维点插入八叉树数据结构中完成点云数据存储。

1.2.2 隐函数采样

(11)

图5 邻近点示意图

式中:d为v点与叶节点n中某一点的距离,ln为叶节点n的边长。

计算出表面隐函数采样值后,再利用任意八叉树的无约束等值面提取方法[12]提取隐函数等值面,即为物体三维表面。

2 实验结果与数据分析

先使用运动恢复结构系统(Structure From Motion,SFM)[13]和基于窗口匹配的多视图三维重建方法[14]获取场景稠密三维点云,再对同一稠密三维点云分别使用SSD[9]、PSR[6]、SPSR[8]和本算法4种不同的表面重建算法重建其三维表面。通过大量实验,并从定性和定量两个角度对4种表面重建算法的重建结果进行了对比和分析。以下实验均在ubuntu 20.04系统、8核intel i7-9700 CPU、32 G内存和2T机械硬盘的实验条件下得到。

2.1 实验结果定性评价及分析

如图6和图7所示为对BYLOW等[15]拍摄的DOOR场景和PARK_GATE场景表面重建结果对比图。考虑到4种表面重建算法顶点数相近的条件,DOOR场景前3种算法设置的深度值为10,而PARK_GATE场景前3种算法设置的深度值为12,本算法无需设置深度值,将根据点云位置与点云尺度大小自动设置八叉树局部深度。由图6a～图6d可见,4种算法均可完整地重建出整个DOOR场景。图6e～图6h对应于图6a～图6d正门上方红色框内局部区域的放大细节。由图6h可见,本算法重建的动物面部周围细小结构更加清晰。图6i～图6l对应图6a～图6d左门柱上黄色框内局部区域的放大细节。由图6l可见,门柱上花瓣花芯以及下部凹凸起伏的细小雕刻更清晰完整。如图7a～图d所示为4种算法重建的PARK_GATE场景正面。由图7d可见,本算法重建的地面、墙面和植物更完整清晰。图7e～图7h对应于图7a～图7d正门右方红色大框内局部区域放大细节。由图7h可见,本算法重建的墙面、文字标牌细节更清晰,特别是人面部和头部细小结构清晰可见。图7i～图7l对应图7a～图7d正门右方黄色小框内局部区域的放大细节。由图7l可见,本算法对较难重建的花朵绿叶等复杂多变结构的重建结果相对更自然合理。

图6 DOOR场景表面重建结果对比

图7 PARK GATE场景表面重建结果对比

2.2 实验结果定量评价及分析

2.2.1 评价原理

三维模型的重建精度与完整度是衡量三维重建算法的主要性能指标。重建精度越高,表明重建结果越逼真,重建完整度越高表明重建的场景越全面,重建精度与完整度越高表明重建算法性能越好。然而,重建精度与重建完整度通常相互制约,提高重建精度在一定程度上会降低重建完整度,反之亦然。实验根据丹麦科技大学标准三维重建数据库(Danmarks Tekniske Universitet multiple-view stereopsis data sets,DTU)[16]的定量评价机制对4种表面重建算法的重建结果进行了定量评价。其评价原理如图8所示,图中蓝色R为重建三维表面,蓝色点为重建表面上的点;绿色G为标准三维表面,绿色点为标准表面上的点。在计算精度时如图8a所示,对于R中的每一点,计算该点到G面的垂直距离,如点pr1到G的垂直距离为红色虚线所示,该线段越短,表明重建点越接近真实面,精度越高。在计算完整度时则反过来操作,如图8b所示,对于G中的每一点,计算该点到R面的垂直距离,如点pg1到R的垂直距离为红色虚线所示,该线段足够短则表明物体该位置被成功重建,该线段太长则表明该位置没有被成功重建,如pg2点对应的物体位置没有被重建出来。DTU库没有提供标准三维表面,但是提供了通过高精度激光测量系统得到标准三维点云作为参照对像(Ground Truth,GT),即图8中的绿点。DTU使用点到点的最近距离近似替代点到面的垂直距离,如图8a中,DTU找到G中离pr1点最近的点,并将该点与pr1点的距离即图8a中蓝色虚线近似的替代pr1点到G面的垂直距离即红色虚线。同样的方式在图8b中用绿色虚线替代红色虚线,当点云非常稠密时替换误差非常小,常常可以忽略。因此在计算精度时,对于R中的每一点,计算其与G中最近点的距离为该点的精度值,将R中所有点的精度值求平均即可得到R的平均精度。将所有点的精度值由小到大排序,处于90%位置的精度值被Middlebury标准数据库[17]定义为R的重建精度。在计算完整度时,对于G中的每一点,计算其与R中最近点的距离为该点的完整度值,将所有点的完整度值取平均即为R的平均完整度。Middlebury标准数据库[17]将R中所有完整度值大于1.25 mm的点数与R的所有点数的比值定义为R的重建完整度。由于点云的稀疏分布会影响两指标计算的可靠性,因此DTU对点云进行了均匀采样。由于标准点云的本身不完整可能会影响精度计算的可靠性,如图8a中,pr2点到G中最近点的距离可能因G面的不完整而变得较长,精度变低,因此DTU给出了标准点云MASK,在MASK以外的点云不参与评价。

图8 评价原理示意图(绿点:真值点云;蓝点:重建点云)

由于重建的三维表面与DTU库提供的标准点云不在同一坐标系下,先使用DTU库提供的摄像机矩阵将重建表面转换至标准点云坐标系下,再通过式(12)将重建表面进行微小的旋转、平移和缩放,使重建表面与标准点云对齐。

(12)

其中:M0,M1分别为对齐前后重建表面顶点坐标。Rq是由单位四元数(a1～a4)组成的旋转矩阵。a5为缩放系数,a6～a8为平移矢量。a1～a8的详细求解流程如图9所示。先将a1～a8的初始值设为[1 0 0 0 1 0 0 0],再由i、j控制内、外两层循环,内外循环对应的误差初值为调整前重建表面与标准点云间的平均距离。内循环一次,更新一个参数,并固定其余7个参数不变;外循环一次,8个参数均更新一次。内循环时,根据式(13)确定当前待确定参数ai的步长Δi,即当ai≠0时,步长为ai的5%,当ai=0时,步长为0.01。再分别在ai的基础上增、减一个步长,求得对应误差err1、err2。根据式(14),结合两次参数与误差值预测增长步长Δpre,求出对应误差err3。在3次误差中选取其中最小误差和对应参数ai的值,并将其更新内循环误差err和参数ai,如此循环8次完成8个参数的一次更新。再根据式(15)计算外循环误差Err与内循环误差err的总体相对误差der,若der小于阈值the,在本实验中the取0.001,则外循环次数加1,并更新外循环误差Err的值为此时内循环误差err的值;否则外循环次数清零,重新开始计数。若外循环连续50次使得总体相对误差der在0.001以内,则结束外循环,此时8个参数的值即为最终取值。求解过程中误差定义为重建表面顶点与标准点云之间的平均距离,即为图8a中的重建平均精度。对于重建表面中的一个顶点,求取其在标准点云中的最近点,并计算两点距离为该重建顶点的精度,将重建表面所有顶点的精度求平均值得重建表面平均精度,即为图9中的误差。在实际求解过程中,为了提高求解速度,没有将重建表面全部顶点参与计算,仅选择了重建表面有代表性的多个部位顶点参与计算。

图9 重建点云对齐参数求解流程图

(13)

Δpre=ai-(err1-err2)×Δi/2(err1+err2-2×err);

(14)

der=(Err-err)/Err。

(15)

如图10所示为DTU库SCAN07场景重建表面模型与标准点云模型的对齐过程。图10a为对齐前两模型位置分布,靠左较小的为重建表面,靠右较大的为标准点云。由图10b可见,重建表面从对齐前的靠左位置转换到了对齐后的靠右位置。由图10c可见,对齐前靠左位置的重建表面,对齐后转换到靠右的位置与标准点云重叠在一起。

图10 两模型对齐过程

2.2.2 评价结果与分析

如图11和图12所示为DTU库SCAN07和SCAN46场景的重建结果对比。两图中结果排列顺序为,从上到下依次为SSD[9]、PSR[6]、SPSR[8]、本算法重建表面和GT;从左到右,第一列为场景整体正视图,第二列为第一列中红框内细节放大,第三列为相同手工剪裁后模型结果。考虑到4种重建表面顶点数相近的条件,SSD[9]、PSR[6]、SPSR[8]重建SCAN07和SCAN46场景表面时的深度设置为10和12,本算法无需设置深度,八叉树将根据点云位置与尺度自动调整局部深度。由图11和图12第一列可见,4种表面重建算法均能较完整地重建两个场景。由图11第二列细节可见,本算法重建兔子的鼻子、嘴巴及周围起伏细节更清晰更接近GT。仔细观察图12第二列细节可见,本算法重建的墙面上微小起伏仍然可见,墙面到窗户的转折以及凸出镜面的细小窗框更加清晰,没有出现丢失部分窗框或边沿转折被平滑的情况,明显优于其他3个表面重建结果。

图11 SCAN07场景表面重建结果对比

图12 SCAN46场景表面重建结果对比

虽然DTU库给出了MASK用于屏蔽不能参与评价的位置,但实验发现仍有小部分GT空洞对应位置的重建表面参与了评价。因此,为使评价公平,手工剪裁GT中包含空洞的位置,并对所有重建表面进行相同的剪裁,如图11和图12第三列所示,再将剪裁后的结果进行定量评价。

如表1所示为SCAN07场景定量评价结果。从中可见SSD[9]、PSR[6]、SPSR[8]和本算法重建表面顶点数量接近,分为1 022 218个、1 004 398个、1 016 716个和1 012 061个。在4种重建算法中,本算法重建表面平均精度与平均完整度均略高于其他3种算法。其中,平均精度为0.21 mm,略优于PSR[6]和SPSR[8]重建表面平均精度0.22 mm,与SSD[9]重建表面平均精度0.21 mm持平,但本算法重建表面精度0.41 mm略高于SSD[9]重建表面精度0.42 mm。本算法重建表面平均完整度为0.34 mm,优于SSD[9]、PSR[6]、SPSR[8]重建表面平均完整度0.35 mm、0.36 mm、0.36 mm。由于获取点云时在深度融合环节产生了冗余信息导致点云数量较大,有2 066 050个空间点,从而使得点云的平均完整度0.26 mm高于4种重建表面平均完整度0.35 mm,0.36 mm,0.36 mm,0.34 mm。又因为精度与完整度普遍存在相互约束,所以点云平均精度0.30 mm低于4种重建表面平均精度0.21 mm,0.22 mm,0.22 mm,0.21 mm。

表1 SCAN07场景重建结果定量评价与对比

如表2所示为SCAN46场景定量评价结果。由表可见SSD[9]、PSR[6]、SPSR[8]和本算法重建表面顶点数量相差不大,分别为1 569 791个、1 597 947个、1 638 734个和1 699 518个,点云点数最多为6 160 510个。点云平均完整度完整度相对最高,为0.23 mm99.36%;平均精度精度相对最低,为0.58 mm1.06 mm。4种重建法中,本算法重建表面平均精度与平均完整度均为最优。其中,平均精度0.35 mm明显优于SSD[9]、PSR[6]和SPSR[8]重建表面平均精度0.41 mm、0.47 mm和0.41 mm。本算法重建表面平均完整度0.39 mm明显优于SSD[9]、PSR[6]重建表面平均完整度0.45 mm、0.51 mm,略优于SPSR[8]重建表面平均完整度0.40 mm。

表2 SCAN46场景重建结果定量评价与对比

总体而言,SCAN46场景凹凸起伏大,细小结构多;SCAN07场景相对而言凹凸起伏小,细小结构少。在重建两个场景时本重建算法平均精度与平均完整度均处于相对靠前的位置。由于表面隐函数根据点云尺度大小和信心指数自动调整点云权重和基函数影响范围,在重建凹凸起伏较大的结构时抑制因遮挡等原因造成的噪声影响更加明显;由于八叉树根据尺度大小自适应调整局部深度与采样密度,较好的保留了表面细小结构;因此,相比SCAN07场景,本算法重建SCAN46场景的平均精度与完整度领先优势更加明显。由于SSD[9]存在平滑效应,SPSR[8]存在抑制平滑效应,因此在重建SCAN07场景时,SSD[9]重建平均精度与平均完整度优于SPSR;而在重建SCAN46场景时,SPSR[8]重建平均精度与平均完整度优于SSD[9]。