局部相关点平衡的三维激光点云骨架提取

周明全，李纯辉，王丽清，张雨禾，耿国华

（西北大学 信息科学与技术学院，陕西 西安 710127）

1 引言

通过激光扫描设备扫描获取物体的点云模型，具有简单便捷、准确率高等优点［1］，成为目前最常用的数字化手段之一［2］。由一组骨架点和骨架点之间的拓扑关系［2-4］构成的点云模型骨架线能够有效描述点云模型的整体形状和拓扑结构［4-5］。因此，骨架线的计算、提取和匹配等成为计算机图形学［6-7］、计算机视觉［8-9］和机器智能［10］等领域的研究热点，并被广泛应用于形状分析［11］与重建［12］、路径规划［13-14］、动作和手势识别［15-17］、模型检索［18-20］以及植物分析［21-22］等任务中。基于骨架分析的相关研究都需要以骨架提取为基础，因此骨架线的提取至关重要。常见的骨架类型按照空间维度可以分为2D骨架和3D骨架两类［12］，其中2D骨架是由一组位于模型各分支的中间曲 线 构 成 的［4，6，23］，而3D骨架［24］则是由 中 心曲面和中心曲线共同构成的。相较于3D骨架，2D骨架线具有信息压缩程度高、拓扑关系明确等优势，因此其应用也更为广泛［9］。

目前已有诸多研究人员开展了点云模型2D骨架线提取相关工作，具体可以分为：（1）基于质心Voronoi细 分（Centroidal Voronoi Tessellation，CVT）的方法［25-26］。该方法是骨架提取中应用最多的策略。例如：Lu等［26］通过变分特征将CVT概念推广到图中，将CVT应用在骨架提取的任务中，在提取骨架的同时还实现了对模型的形状分割，但该方法时间复杂度较高，针对大规模点云模型处理时间较长。（2）基于平均曲率流（Mean Curvature Flow，MCF）的方法［27-29］，该方法能提取到高质量的骨架。例如：Chuang等［29］发表了一种基于有限元跟踪的快速MCF的方法，随 后Tagliasacchi等［27］将MCF应 用 在 骨 架 提 取的任务中，通过分析流动的微分特征，发现MCF局部增加了形状各向异性，利用其面积最小化的特性将曲率流推向极值，从而折叠输入网格几何体并获得骨架结构。（3）基于拉普拉斯收缩的算法。Cao等［4］提出了一种基于拉普拉斯收缩的曲线骨架提取算法（Laplacian-Based Contraction，LBC）。但是这些方法只能应用在经过降噪、去除异常值等操作处理之后的模型上，难以应用于使用激光扫描设备所获取的点云模型。（4）基于深度学习的方法。Atienza等［30］基于Pyramid UNets提出了PSPU-SkelNet，适用于从二维形状点云模型中提取曲线骨架任务。Lin等［24］基于中轴变换（Medial Axis Transform，MAT）提出了Point2Skeleton（P2S），用于从三维点云模型中提取二维骨架，从而指导表面重构和分割。Dou等［31］基于覆盖轴（Coverage Axis）提出一种深度学习的提取骨架的算法，其结果接近基于MAT的骨架提取算法。

适用于激光扫描点云模型的曲线骨架提取方法，大多都基于几何学相关原理［2-7］。Tagliasacchi等［2］基于旋转对称轴（Rotational Symmetry Axis，ROSA）的概念，构造了一种可以针对缺失大量数据的不完整点云模型的曲线骨架提取方法。Huang等［6］将L1中轴（L1-Medial）应用于曲线骨架提取任务中，将L1中轴局部化自适应提取不同位置的中值点作为骨架点。Qin等［7］通过构建基于质量驱动拓扑感知的最优传输方案提取粗糙曲线骨架（MdCS），然后使用能量正则化完成几何优化，最终获得平滑曲线骨架。ROSA虽然能够应对大面积缺失数据的原始扫描点云模型，但是旋转对称轴的计算需要无噪声扰动的点云，因此针对噪声点云效果较差。L1中轴骨架和MdCS虽然能应对存在噪声、异常值和缺失数据的原始扫描点云模型，但是L1中轴骨架不能处理多分支结构的模型，MdCS拓扑结构准确度不高，且计算所需时间资源消耗巨大。

针对上述工作中存在的问题，受文献［32-34］的启发，本文提出一种快速提取激光扫描点云骨架线的方法。该方法将初始骨架点定义为三维模型中具有对称法向的最近相关点的中点，通过中点状态的平衡化迭代得到最终的骨架点；然后使用主成分分析法（Principal Component Analysis，PCA）分析骨架点的方向一致性，从而得到骨架曲线段；最后将使用广度优先搜索（Breadth First Search，BFS）合并后的曲线段利用拉普拉斯平滑进行优化并得到最终曲线骨架。本文实验将此算法应用在通过激光扫描设备得到的点云模型上，验证了算法的有效性、鲁棒性和高效性。将提取的激光扫描点云的骨架线应用在点云模型的动作变换任务中，展示了骨架线的实际作用。

2 基本原理

本文的骨架提取方法主要可以分为两步：提取骨架点和生成骨架线。首先，通过选择具有对称法向和最近距离的对称点的中点，获取初始骨架点，如图1（a）所示；应用射线法去除模型外的初始骨架点，如图1（b）所示；通过迭代来平衡骨架点并筛选权重大的骨架点组成最终骨架点，如图1（c）所示。然后，基于PCA判断骨架点的方向一致性并重建骨架曲线，采用广度优先遍历算法合并骨架曲线，效果如图1（d）所示；最后将平滑后的骨架曲线相连得到最终的曲线骨架，如图1（e）所示；基于骨架的模型动作变换实际效果如图1（f）所示。

图1 本文曲线骨架提取算法的过程展示Fig.1 Pipeline of proposed method

2.1 提取骨架点

骨架点可以看作模型的部分采样点投影到模型内部得到的局部中心点［6］，也就是说，初始的骨架点是模型的点集通过采样处理得到的子集；另外，直接对采样点进行计算，还可以有效提高骨架线提取的效率，避免遍历点云上的所有点数据。除了基于效率的考虑，对模型下采样还能减少模型本身的一些不利因素对最终骨架的影响。例如，当模型不均匀时，通过均匀采样便可以避免这一缺陷；当遇到噪声时，采样也会减少异常值噪声对骨架提取任务的干扰。因此，本文算法首先进行随机采样，使得采样点能均匀的覆盖模型，从而避免骨架线变形。具体而言，对于一个离 散 点云模型P={pi}i∈I⊂R3，其中I是 点 云P的点数，将采样点组成的集合记作Ps，从P中使用随机算法获得点的索引，本文实验中采样点的规模为0.1×I。

2.1.1 法向估计和重定向

本文中初始骨架点需要基于法向求解，所以首先需要估计采样点的法向［35-36］。本文使用PCA对Ps中的点进行法向估计，即：ps及其在P上的邻域点共同构成点集Pt，对Pt进行特征分解得到特征值和特征向量，将最小特征值对应的特征向量pn t作为ps的法向。具体如下，利用式（1）计算邻域点集Pt的中心点pm：

此时ps和法向pnt应该满足pnt使 得Pt中 所 有的邻域点在方向pnt上的投影分布最密集：

式中：k=20，将pt-pm记作则有

由于上述过程计算的法向具有二义性［37］，如图2（a）所示，而这种二义性会影响骨架点的求解过程，所以还需要对法向的二义性进行消除，即法向的重定向。具体而言，对于两个距离很近的采样点pi，pj，如果他们的法向pni，pnj同向，则满足≈1。本文使用深度优先算法（Depth First Search，DFS）来快速完成重定向。具体而言，选择距离Pt的中心点PsO最远的采样点pi，如果‖(pi+pni)-PsO‖＞‖pi-PsO‖，则pni=-pni，反 之则不改变pni。随后使用K最邻近算法（K-Nearest Neighbor，KNN）可以直接得到pi的距离最近的一个Ps中的点pj，如果pni⋅pnj/(|pni|×|pnj|)≤-0.5，则pnj=-pnj，并将pj标记为已经访问的状态，此时令pi=pj。通过KNN算法寻找一个没有被标记为已访问状态的采样点作为pj，重复上述过程，直到Pt中所有的点都被标记为访问的状态。此时Pt中所有的点都将指向模型的同侧，消除了法向的二义性，如图2（b）所示。用该重定向方法重定向之后的法向，都指向模型的内部。

图2 法向二义性示意图Fig.2 Schematic diagram of normal vector ambiguity

2.1.2 骨架点的求解

要进行骨架点识别，本节首先做以下定义：

定义一：对于采样点集合Ps中的任意一点pi，存在pj∈P，j∈S{i}是和pi具有相反法向的点中距离pi最近的点，将pj和pi定义为一组相关点。

定义二：一组相关点的中点是初始骨架点。

对于Ps中的任意一点pi，从Ps中识别出一组法向与pni相反的点，在这里法向pni和pnj使用一个带阈值的点积来判断。在得到的这一组法向中，进一步找到特定点pj，使得pj到pi的欧氏距离最小，随后计算出pj和pi的中点作为初始骨架点。当模型的分支之间的距离较小时，识别到的相关点可能并不在一个分支上，导致提取到无效的冗余骨架点，如图1（a）中手臂和身体之间产生的无效骨架点，此时需要通过算法将不在模型中的相关点之间的中点，也就是无效的骨架点去除掉。具体而言，本文使用射线法来去除无效骨架点，从而得到有效骨架点。射线法具体是指从任意一点出发，从该点射出的射线，如果与模型边界的交点数量为奇数，则该点在模型内部；反之，如果交点为偶数，则该点在模型的外部，如图3所示（彩图见期刊电子版）。

图3 （a）中所需判断的点（被标记为红色）在模型外部，黄色射线与模型的交点（被标记为绿色）有两处，是偶数。（b）中所需判断的点（被标记为红色）处于模型内部，黄色射线与模型的交点（被标记为绿色）有一处，为奇数Fig.3 Figure（a）shows that when the point（marked red）to be judged is outside the model，the intersection of the yellow ray and the model（marked green）has two even numbers.The point（marked red）to be judged in the figure（b）is inside the model，and there is an odd number at the intersection of the yellow ray and the model（marked green）

在将无效的骨架点去除之后，效果如图1（b）所示。当模型足够简单时，例如圆柱模型，这些点可以直接形成模型P的骨架，但是当模型不是简单的基本几何形状时，还需要使初始骨架点进一步收敛到最终骨架点的位置。由于此时的有效骨架点qi∈Q可以看作是模型上的采样点pi经过位置变换得到的，受文献［38］的启发，在垂直于骨架点法向qni的平面上，在每个有效的骨架点qi处均匀地发射出m条射线，其中第j条和第j+m/2条射线的方向是相反的。从骨架点qi出发，沿着每条射线方向迭代地向前搜索，直到靠近模型上的任意一点，将qi到该点的距离记作dj，将qi到模型上该点的方向标记为向量vqji。本文实验中，m的数量为16条，被分为8对，每对由两个方向相反的向量组成。对于点集Q中的任意一点qi，定义qi的能量Eqi为：

由上述定义可得到骨架点Q的平衡状态如式（4）所示：

其中：NQ是Q中点的数量，r是模型的最小边界球的半径。根据上述定义，可知qi和qi+1之间的关系如式（5）所示：

其中，L=0.01是步长参数。当式（3）中能量梯度为0时，意味着任意一个qi的E=0，也就是qi的位置处于模型横截面的中心位置，如图4（c）所示。根据上述过程，便可以完成骨架点的提取。在获取初始骨架点下个位置的同时也可以得到下次移动的方向，可以加速初始骨架点平衡化的过程。但是此时的骨架点数量还是远远大于最终用来生成骨架点的数量，所以还需要进一步从骨架点中筛选合适的骨架点。这里使用二次邻域中权值比较的方式来完成骨架点的选择。具体而言，对于一个骨架点qi，其权值被定义为骨架点qi的rknn邻域内点的数量N。对于Q中任意一个骨架点qi，将qi的rknn邻域中的点记作Q r，则对于Q r中所有的骨架点而言，应该保留的qr应该满足式（6）：

图4 （a）展示的是得到垂直于qni的平面的方法，pi+1是pi距离最近的采样点，（b）中骨架点qi处于不平衡的位置，（c）中骨架点qi处于平衡位置Fig.4（a）shows the method of obtaining the plane perpendicular to qni，and pi+1 is the nearest sampling point from pi，in（b），the skeleton point qi is in an unbalanced position，and in（c），the skeleton point qi is in an balanced position

2.2 骨架生成

上述迭代过程能够得到一组骨架点Q={qi}，i∈NQ，它是输入点云模型P的各个分支的中间地带，仅表示了模型P骨架的位置信息，而完整的骨架线还包含了骨架点之间及骨架分支之间的拓扑连接信息。所以本节计算骨架点Q的拓扑连接信息。

首先，使用PCA来检查qi和它的邻域点集是否满足方向一致性的条件，随后将满足要求的邻域点集生成模型曲线骨架拓扑关系的子集，进一步通过广度优先搜索算法将存在连接关系的子集合并得到模型的骨架曲线分支，再通过拉普拉斯平滑骨架曲线分支，最后将平滑后的骨架曲线分支组织成完整的曲线骨架。

PCA常被用来检测骨架点是否达到收敛要求［6，23］，即：对于骨架点qi，计算qi及其k邻域内点集的协方差矩阵的特征值和特征向量，如式（7）所示：

其中：所有的特征值满足λ0i≤λ1i≤λ2i都是实值，对应的特征向量{v0i，v1i，v2i}形成一个正交框架，即qi和其邻域点的主成分。定义qi和其邻域点的方向一致性为σi，则有如式（8）的关系：

σi的实际含义是当σi越接近1，即λ0i，λ1i和λ2i相比越小，更多的点分布在沿着分支对齐的位置，也就是方向一致性越高。实践中，设置k=3，对于任意qi而言，如果σi＞0.9，便可以认为qi和它的两个邻域点已经满足骨架点的收敛要 求，可 以 认 为 对 于 点qi和qi+1，qi+1和qi+2而言，存在拓扑关系。随后遍历所有骨架点，便可以得到达到要求的所有的两个点之间的拓扑关系，记 作R={(rmn)}m∈Q，n∈Q，m≠n，rmn=(qm，qn)，如图5（a）所示。

图5 R和RL的区别，（a）是R，（b）是根据R通过广度优先搜索得到RLFig.5 Difference between R and RL.（a）is R，and（b）is RLobtained through breadth first search according to R

得到R之后，使用广度优先搜索，增长(qm，qn)的关系，得到RL，如图5（b）所示。增长后的效果如图1（c）所示。随后对于RL中任意rnm=(qm，qm+1，…，qn-1，qn)两端的端点qm和qn，对qm寻找最近骨架点qm+1，构成方向向量vm+1m=qm+1-qm，同理得qn的方向向量vn n-1=qn

qn-1，使它们在此方向上迭代寻找合适的骨架点相连即可得到完整的骨架。由于在此过程中，在尺度变化剧烈时，可能会导致骨架曲线的变化不够平滑。为了减少这种剧烈变化带来的质量上的影响，采用拉普拉斯收缩以获得平滑的骨架曲线。对于拉普拉斯平滑策略点qm的定义如式（9）所示：

其中：i=m+1，m+2，…，n-1。

经过平滑后，便可以得到完整的曲线骨架，如图1（d）所示。

3 实验以及结果分析

本文实验所用硬件和实验环境包括：CPU是AMD的 锐 龙7代3700X，内 存 为32 GB，GPU为NVIDIA GeForce RTX 2060 SUPER，Visual Studio 2017 Professional，PCL（Point Cloud Library）1.8.0，开发语言为C++。为了测试和评估算法的实际性能和效果，本文实验采用了公开数据集Visionair Repository中的多个三维模型，以及部分利用FARO Focus Premium激光扫描仪激光扫描设备扫描获取的点云模型作为测试数据，测试本文方法及对比方法的结果。

3.1 本文方法实验结果

图6所示为本文算法应用于鹿、手掌、人、雕塑模型的实验结果。图6（a）中鹿模型的难点在于，鹿的脖子部位存在轴向偏转，骨架点位置易偏移，导致骨架线失真；然而，从图6中的结果可以看出本文算法受到的扰动较小，原因是虽然存在轴向偏转，但轴向偏转并没有导致模型脖子部位的形状产生较大变化，所以对骨架生成没有产生干扰。图6（b）中手掌模型的难点在于多个距离很近的手指会导致骨架点出现在手指之间的位置，产生错误的骨架点和骨架线。在本文中，为了避免这种现象的产生，采用了射线法来去除掉出现在模型之外的错误位置的骨架点，最终结果证实了此算法的有效性。图6（c）中人模型对于骨架提取算法的挑战在于可能会在手臂和身体之间产生错误的分支；图6（d）中雕塑模型的难点在于存在复杂的拓扑结构，然而本文算法都能有效应对这些难点，提取出完整平滑的曲线骨架。

图6 本文算法提取的模型的曲线骨架（右侧）和原模型（左侧）Fig.6 Curve skeleton of some models（right）and the display of the original model（left）extracted by this algorithm

为了进一步验证本文算法的有效性和优势，本文将其与目前3种前沿方法进行比较，即：L1中轴骨架算法、MdCS算法和Point2Skeleton算法。Point2Skeleton所使用的模型已按照其要求通过均匀采样，将点云下采样到2 000点后作为其输入。图7（a）为4种算法应用于稀疏点云狗模型上的结果。由于狗模型本身的点数较少，狗模型的腹部存在一些点的空白，导致使用点的坐标作为输入的方法提取到曲线骨架的难度更高。L1中轴骨架算法由于内切球的半径增长，导致无法提取到其他腿部的骨架。MdCS算法在狗模型的身体上提取到了错误的拓扑结构。Point2 Skeleton算法将狗的身体提取为二维骨架，虽然存在拓扑连接错误和狗模型的右前腿的骨架缺失，但整体效果好于L1中轴骨架算法。然而，本文算法能够得到较为满意的效果，这是因为虽然模型上点缺失的部分难以估计法向，但本文采用了迭代提取骨架点的方法，降低了这种稀疏采样对骨架线提取的影响。图7（b）所示为马模型的骨架线提取结果。马模型在骨架线提取任务中的难点在于，四肢和身躯之间存在巨大的尺度变化，这种尺度变化会影响基于几何性质收缩或者内切球增长的算法的有效性。从图中可以看出，L1中轴骨架算法继续受到内切球半径增长导致分支骨架缺失，MdCS算法即使存在拓扑检查的过程，但依然不能应对马腿和身体之间的巨大尺度变化。Point2Skeleton算法在马头部位与真实情况相差较大，在身体部位提取的二维骨架失真严重。而本文算法能够得到合理的拓扑关系并生成完整的骨架线。图7（c）所示为人模型，人腿部和腹部的连接点的选择是一个难点，导致L1中轴骨架算法选择的连接点不在模型的内部，MdCS算法在腹部依然存在拓扑结构错误，Point2Skeleton算法则未提取到模型的左手。而本文算法在腹部未选择合并，导致腿部骨架和身体骨架合二为一。

图7 各种方法提取的曲线骨架效果对比Fig.7 Effect comparison of curve skeleton extracted by various methods

3.2 鲁棒性实验

为了验证本文算法的鲁棒性，本节针对添加高斯噪声的模型开展实验。本节实验中，定义高斯噪声的标准差为2×dmin。其中dmin是模型中点之间最小距离。结果如图8所示。MdCS算法在应对有噪声的点云模型时效果好于没有噪声的模型，原因在于，添加噪声后，模型点的坐标发生了变化，导致算法在拓扑检查时使用了较大的半径，减少了拓扑错误。由于MAT易受噪声影响，所以Point2Skeleton算法在应对噪声时的表现变差，产生错误拓扑连接的可能性相较于不受噪声干扰时有所提高。受到噪声干扰时，Point2Skeleton生成了更少的拓扑连接，由此虽然拓扑连接错误的数量下降了，但是骨架的质量却没有提高。

图8 输入带噪声的模型时各种方法提取的曲线骨架效果对比Fig.8 Effect comparison of curve skeleton extracted by various methods when inputting the model with noise

由于较低程度的高斯噪声对PCA中最小特征值的影响有限，所以对基于PCA的法向计算的影响也非常有限。并且中点状态的平衡化迭代可以在得到居中性更好地骨架点的同时，降低噪声的干扰，所以本文算法受到噪声影响较小。

为了定量地分析几种算法的实际效果，统计了几种方法在应对模型无噪声和有噪声的情况下的骨架点数量、拓扑连接的数量、交叉点数量和拓扑连接错误的数量。其中，拓扑连接错误的数量等于骨架缺少的分支数量和去掉的拓扑连接数量的和，去掉的拓扑连接数量为去掉拓扑能使拓扑连接正确所需去除的数量。如表1所示。从表1可以看出，L1中轴骨架线提取算法提取的骨架点数量较多，拓扑连接错误的数量也较多，骨架分支缺失的情况较多。MdCS算法由于提取的骨架线可能存在环结构，所以拓扑连接错误较多，而且可能存在骨架点和拓扑连接在模型外的情况。Point2Skeleton具有固定的输出，即：100个骨架点和骨架点的拓扑关系。由于其输出的是二维骨架，二维骨架的连接状态较多，所以拓扑连接关系出现错误的可能性较大，这也反应在了图8中。本文算法提取的骨架线和骨架点数量中等，交叉点较少，拓扑连接数量错误较少，并且受到噪声影响时，提取的骨架点数量和拓扑连接数量变化不大，表现出了对噪声的鲁棒性。

表1 几种方法的骨架点数量、拓扑连接数量、交叉点数量、拓扑连接错误数量的统计Tab.1 Statistics of the number of skeleton points，topological connections，junction points and topological connection errors of several methods

3.3 算法效率

利用相同的软硬件设备运行L1中轴骨架线提取算法、MdCS算法、Point2Skeleton算法和本文算法，并对计算耗时进行统计，结果如表2所示。从表2中可以看到，本文算法的时间损耗小于其他三种算法，即：略小于Point2Skeleton算法，相较于L1中轴骨架算法速度提升了10倍以上，比MdCS算法速度提升了百倍以上。由于Point2Skeleton是基于深度学习的算法，所以在计算骨架时相较于传统算法更有优势，而该方法实际耗时表现略差与本文方法的原因是，Point2Skeleton在加载深度学习模型和保存运算结果步骤中耗时较多。

表2 各算法的运行时间Tab.2 Running time of different algorithms（s）

本文算法速度较快的原因是，使用了降采样的方式降低了参与运算的点的数量。骨架点迭代时，由于骨架点平衡状态判断的同时，也会在常数时间内获取骨架点迭代的方向，这种设计加速了骨架点的平衡移动过程。

3.4 应用

骨架线常被应用于动作捕捉识别［15-17］，基于骨架线可以改变模型的姿态和动作，对于点云重建、数据生成以及形状变化具有重要作用。直接通过三维图形编辑软件制作动画的效率瓶颈是导入的模型需要手动绑定系统内预置的骨架，包括骨架的三维坐标和拓扑连接关系。此过程不仅耗费时间，并且因为系统预置骨架的种类有限，也会影响软件的通用性。通过骨架线提取算法得到骨架线后，三维图形图像编辑软件可以自动获取骨架点的三维坐标、构建拓扑连接，不再需要手动调整参数和绑定骨架设定，可以提高动画制作的效率。

图9（a）中展示的是，使用FARO Focus Premium激光扫描仪得到玩具马（右上黄棕色）的三维点云模型（右下）；通过本文的算法提取到计算扫描点云的骨架线，如图9（b）所示；将骨架线和点云模型导入三维图形编辑软件Blender中，可以轻松实现对于模型动作的改变，图9（c）中展示的是对马的激光扫描点云基于骨架线的动作变形，将马的右前腿从弯曲动作改变为向前翘起。图1（e）和1（f）展示了基于骨架线将人模型的右脚和右手臂抬起后的效果（彩图见期刊电子版）。

图9 基于本文骨架线指导的激光点云的动作变换Fig.9 Action transformation of lidar point cloud based on the guidance of curve skeleton

4 结论

本文提出了一种适用于激光扫描点云的骨架线的方法，并在此骨架基础上实现了对模型的变形。通过将骨架点定义为具有对称法向且距离最近相关点的中点，通过迭代使中点平衡化得到骨架点。然后采取主成分分析法寻找满足方向一致性要求的曲线段，使用广度优先搜索算法合并各曲线段，最后将拉普拉斯平滑后的曲线段相连得到完整的骨架线。通过在多个点云模型上开展对比实验，验证了本文算法的有效性、鲁棒性和高效性，将原本处理由33 041个点组成的点云所需86.110 s缩短到4.356 s。将提取的骨架应用在激光雷达点云的动作变换任务中，展示骨架线的实际作用。