基于距离场的网格模型骨架提取

(浙江工业大学 计算机科学与技术学院，浙江 杭州 310023)

3D模型在计算机辅助设计、计算机图形学和科学可视化等领域中扮演重要角色。作为重要研究对象，3D模型被广泛应用，包括三维模型特征提取[1]、孔洞修复[2]、模型重建[3]和遥感等。传统计算机图形学中的三维形状表示方法一般表示的是模型表面信息，并不描述模型内部结构和模型拓扑信息[4-5]。而模型的三维表示是“昂贵的”，许多应用需要模型的可替代的紧凑表示[6]。有这样的一种“线状”或“棒状”表示方法，称之为“骨架表示”或“曲线骨架”[7]。Blum[8]以中轴的形式建立了骨架化的基础，该骨架由具有较低维度的对称平面/轴组成。理论上讲，Blum的骨架或中轴是使用草火传播过程定义的。假设物体是干草材质，在物体边界用火点燃，火焰以均匀的速度在物体内部传播。骨架被定义为两条独立火线相遇点的集合。

模型骨架提取方法可以分为以下几类：1) 拓扑细化法(Topology thinning)，又称为模拟烧草模型法。此类算法从边界开始，反复迭代地逐层剥离离散后的模型形成骨架。绝大部分细化方法依靠简单点作用于离散空间(体素)，去除之后不会改变形状的拓扑结构的点(体素)称之为简单点[9]。为提供拓扑保护，文献[10]提出了两种分别用于提取“表面骨架”和“曲线骨架”的8 次子迭代算法，其中每次子迭代只能去除某种特定类型的边界点。文献[11]提出了一种拓扑保护的三维骨架化方法，该方法主要计算体素标记为D6距离模型的面骨架或曲线骨架。考虑笛卡尔超立方网格，文献[12]提供了一个数学框架，其中给出了迭代细化过程去除点的显式布尔条件。2) 基于距离场的方法。首先，需要估计一个距离图，然后视中轴为距离图的脊，通过脊跟踪算法在此距离图中抽取得到中轴。物体内部每点的距离值DT被定义为此点离模型边界点的最短距离[13]。由于骨架对模型整体而言是处于中心位置应具有较大的DT值，因此DT能作为骨架提取的一种重要参考指标。Bitter等[14]通过梯度搜寻方法检测非均匀梯度的相邻部分，并标记这些点为候选体素。Bouix等[15]用散度计算选取优先函数并由用户自定义阈值来进行体素的简单点移除。3) 基于拓扑与几何分析的方法，通过构造模型的Voronoi图和Reeb图来得到骨架。Voronoi图将空间划分成子区域，Voronoi图的内部边和面可用于提取形状的中轴面(骨架)的近似值。笔者提出了一种基于距离场的骨架提取方法，通过寻找距离域的“脊线” 从而找到骨架。笔者主要的贡献在于提出了最大内切球算法来对候选体素进行约束，能很好地去除冗余体素点。笔者提供了几个模型的实验结果并将之与经典骨架提取方法进行比较，实验结果表明：笔者方法能得到一个质量较好的三维模型骨架。

1 算法流程

笔者提出了基于距离场的骨架提取方法，通过包围盒技术确定包裹模型的最小立方体空间。根据所求得的包围盒将模型体素化，从而将模型划分为一个个立方块空间。再计算每块体素的距离值，找出单位空间内距离值最大的体素作为骨架点的潜在体素。比较过后，所得的候选体素并非全部属于骨架点所在体素，提出最大内切球逼近算法来对体素做出约束，去除冗余的体素，从而找出最终的骨架点。最后，根据所求得的体素，按照一定规则有序连接形成骨架。具体算法流程如图1所示。

图1 算法流程Fig.1 Pipeline of skeleton extraction method

如图1所示，笔者方法主要步骤如下：

1) 模型内部体素化。此过程主要包括包围盒计算和模型外部体素去除。由于笔者方法的核心思想是基于距离场来找到处于中心的那条“脊线”。因此，模型外部体素可用信息不大。根据需求，对于包围盒的选取笔者选择的是AABB盒，主要是其操作简单且存储空间小。在判断体素是否外部体素时，充分利用了模型边界点的法向量指向模型外部这一信息，根据向量积的正负，成功地去除掉外部体素。

2) 单位空间内最大距离值体素的选取及最大内切球约束。模型内部体素化后，从边界到中心，体素的距离值是递增的。因此，位于模型中心的骨架的距离值一定是最大的，所要找的就是这类体素。因此，事先设定一单位空间，然后找出该空间内的最大距离值的体素，比较完所有体素后得到初始筛选的候选体素。然而，这些体素并非全部都属于骨架。因此，利用最大内切球再对这些候选体素进行约束，对于每个候选体素找出其最大内切球，将此最大内切球所包裹的所有体素再进行距离值比较，选择其中距离值最大的体素与最大内切球球心距离值(半径)作比较。若该体素距离值大于该最大内切球半径，则将此体素作为该局部的骨架点。若该体素距离值小于该最大内切球半径，则将此最大内切球的球心作为该局部的体素点。

3) 骨架点连接形成骨架。此部分主要是考虑骨架点连接的规则，根据已求得的骨架点形成骨架。必须保证骨架点连接不能重复、乱序和错漏。

2 最大内切球逼近

三维模型骨架就是模型M中所有最大内切球球心的集合[6]。笔者提出的最大内切球逼近算法主要思想：从模型内部任一点出发，计算其与表面最近点之间的矢量并向模型内部延伸来调整点的位置，从而将初始点位置从低距离值往高距离值的场内移动，经多次调整内部点位置便可以找到三维模型的最大内切球及其球心。

设三维网格模型M(x,y,z)={vk,k=1,2,…,N}，其中(x,y,z)是顶点v的坐标，N是顶点数目，则距离值函数为

(1)

式中v为距离点p最近的模型表面点。图2为Human模型局部放大距离场灰度编码图，颜色越浅代表距离值越大(模型中心为灰白色，模型边界为灰黑色)。

图2 距离场灰度编码Fig.2 Distance field color coding

定义函数Extend_Len(p′,p″,l)来调整点的位置。假设给定点p′和点p″，沿矢量p′p″延长距离l得到点p‴，则

p‴(x‴,y‴,z‴)=Extend_Len(p′,p″,l)

(2)

结合公式(1，2)可得

(3)

最大内切球逼近算法流程如下：

Step1最近点计算。针对模型内部任一点p1，遍历模型表面点，找到与其距离最近的模型表面点v1，并计算p1的距离值以及矢量v1p1。

Step2反向延伸调整位置，找到距离新点的模型最近点并计算两点间的矢量。沿矢量v1p1方向向p1的距离值等值线内部延长距离l2，其中l2=β×L(0<β<1.0)，L为模型网格平均边长。由式(2)确定模型内新点p2=Extend_Len(v1,p1,l2)，并找到距离p2最近的模型边界点v2，计算p2的距离值以及矢量v2p2。

Step3重复执行Step 2，计算pk的距离值以及矢量vkpk，直到‖vk+1pk+1‖<‖vkpk‖为止，即当pk+1的距离值小于pk的距离值，此时pk的距离值达到最大。则以pk为球心，‖vkpk‖为半径可以得到模型最大内切球，其球心为模型骨架点。

图3为horse模型的最大内切球逼近效果图。图3(a)为模型内部初始点位置，图3(b)为其逼近的最大内切球效果图，其中β取0.1，耗时52 ms。

图3 最大内切球逼近Fig.3 Maximum inscribed ball approximation

3 三维模型体素化

体素化是将模型的几何形式表示转换成一种最接近该物体的离散表示形式，产生体数据集。表示模型的空间体素跟表示图像的二维像素比较相似，只不过从二维的面扩展到三维的立方体单元，而且基于体素的三维模型技术被广泛应用。笔者算法主要思路就是操作离散化模型，找到属于骨架的体素点。

3.1 三维模型体素化

基于几种包围盒技术理论，在选择包围体类型时充分考虑各种包围体类型的优缺点以及笔者方法的需求性，最终选择轴对齐包围盒(AABB)，主要是因为其构造比较简单，存储空间小。

设三维网格模型M(x,y,z)={vk,k=1,2,…,N}，其中(x,y,z)是顶点v坐标，N是顶点数目。遍历这N个顶点，找到Xmin,Xmax,Ymin,Ymax,Zmin,Zmax这6 个值，从而确定包围盒的8 个顶点。依次为：P1(Xmin,Ymax,Zmin)，P2(Xmin,Ymin,Zmin)，P3(Xmin,Ymin,Zmax)，P4(Xmin,Ymax,Zmax)，P5(Xmax,Ymax,Zmin)，P6(Xmax,Ymin,Zmin)，P7(Xmax,Ymin,Zmax)，P8(Xmax,Ymax,Zmax)。图4为Human模型AABB包围盒，其中图4(a)为Human模型初始图，图4(b)是其包围盒生成图。

图4 包围盒生成Fig.4 Bounding box generation

3.2 模型内部体素提取

离散化操作中，选取精度为d(2 个单位的平均网格边长度)来将包围盒所包裹的空间划分成一个个体素块。根据得到的包围盒信息，全部体素数量应为Num=Nx×Ny×Nz。其中，Nx=(Xmax-Xmin)/d，Ny=(Ymax-Ymin)/d，Nz=(Zmax-Zmin)/d。

这些体素主要分为3 部分：模型外部体素、模型边界体素和模型内部体素。实验中，实验外部体素通过判定条件剔除，模型边界体素也不可能是骨架所在的体素，真正需要的是模型内部体素。图5(a)提供了人体模型体素化后(已剔除外部体素)的效果图，为更清楚观察体素化后的内部情况，图5(b)提供了身躯的局部放大图。其中，精度为2 个单位的平均网格边长度。为清楚观察每块体素，放置半径为0.5 个平均网格边长度单位的小球在每块体素的中心。外部体素判定是通过向量乘积判定的，如图6所示。对于任一点p(x,y)，计算其与模型表面最近的点，可以利用模型表面点法向量始终指向模型外部这一信息来判定。图6(a)中，模型内部点p(x,y)，p1(x1,y1)是与之最近的模型表面点，由法向量指向外部可以得到：pp1·n1>0，其中n1是p1(x1,y1)的法向量。图6(b)中，对于模型外部点p(x,y)，pp2·n2<0，n2是p2(x2,y2)的法向量。

图5 模型体素化Fig.5 Model voxelization

图6 点的判定Fig.6 Point determination

4 三维模型骨架提取

三维模型骨架提取主要包括对体素化后的模型进行距离值局部最大值求解以及最大内切球约束。距离值局部最大值求解是用一个小的立方块空间去遍历整个模型内部，不能存在错漏，找出该空间内的距离值最大的体素。该筛选后的体素包含处于骨架上的体素点以及需要去除的冗余体素点。最大内切球约束是对筛选后的体素进行最大内切球逼近，找出最大内切球所包围的距离值最大的体素，再将该最大值与该最大内切球球心距离值(半径)作比较。若该体素的距离值较大，则将该体素点作为骨架点。若最大内切球半径大，则将最大内切球球心作为骨架点。经过1～2 次的最大内切球约束，能较好地去出冗余体素。最后再按照一定规则将所求得的骨架点进行有序连接形成骨架。

4.1 距离值局部最大求解

对于体素化模型，只有模型中心位置的体素属于骨架点，本部分主要考虑如何剔除大量不属于骨架位置的体素。实验过程中，考虑使用边长为l的立方块遍历模型内部空间，找出每个空间内距离值最大的体素块，从而完成模型体素快的初步剔除。为提升效率，对每一体素块的数据结构定义为

Voxel=(x,y,z,distance,flag)

式中：x,y,z为体素的位置信息；distance为体素的距离值信息；flag为体素的标志位，表示体素是否被比较，体素标志位设置初始值flag=0。为将笔者方法阐述清楚，先在二维情况下对笔者方法作个引入，如图7所示。

图7 骨架提取原理图Fig.7 Skeleton extraction schematic

如图7(a)所示，对一已体素化的形状，图7(a)中用点来表示体素，该点位于体素中心。图7(a)中所有点初始设置为灰色，对应着三维情况下标志位初始值flag=0。图7(a)中方框对应着三维情况下的立方块，所要比较的就是方框内所包裹的所有体素的距离值，距离值最大的体素标记为深黑色，相应地三维空间内的flag=1，其余的比较过的体素标记为淡灰色，对应的flag=-1。因此，对于所有体素点进行局部距离值最大求解时，会对标志位flag先进行判定，如此可大大提高效率。图7(b)为局部距离值最大求解的结果，所选中的初始骨架点已标记为深黑色。

图8为人体模型距离值局部最大求解得到的候选体素,采用的立方块边长为l=3.25×d。从图8可以看出：距离值局部最大求解已经从人体模型中剔除了大量的体素块。但骨架点主要集中于人体模型中心位置，因此需要再次对这些体素块进行选择。主要考虑使用最大内切球约束的方法去挑选属于骨架的体素。

图8 候选体素Fig.8 Candidate voxel

4.2 最大内切球约束

此部分是对距离值局部最大求解得到的候选体素筛选做详细说明。如图7(b)所示，并不是所有深黑色点均为骨架点，需要进行选择，笔者主要是通过最大内切球的方法进行选择。如图7(c)所示，对深黑色体素点进行最大内切圆逼近，找出该最大内切圆内所包裹的体素点，并将其中距离值最大的体素点的距离值与最大内切圆圆心距离值(半径)进行比较。若体素的距离值相比最大内切圆半径大，则选择该体素为骨架点。反之，则令该最大内切圆圆心作为骨架点。对应三维中，是对所有标志位flag=1的体素进行最大内切球逼近，对于该球所包裹的体素点进行比较。选中为骨架点的体素点，令其标志位flag=2。对未选中的体素点，令其标志位flag=-1。最终标志位flag=2的体素即为所选中的骨架点，如图7(d)所示。

图9为对图8中候选体素进行最大内切球约束效果图。图9(a)是经过1 次最大内切球约束的结果，但仍然有冗余点的存在，需要进一步的选择，如人体模型的胸腔部分。图9(b)为2 次约束的效果图。因为最大内切球逼近能保证所找的骨架点较好地居于骨架正中心，因此图9(b)中所选择的骨架点有较好的中心性。

图9 最大内切球约束Fig.9 Maximum inscribed ball constraint

4.3 骨架生成

本部分主要介绍利用图9所选的骨架点进行骨架生成。如上文所述，所得到的骨架点标志位flag=2。基于距离优先的思想，骨架生成过程如下：

1) 设定2 存储序列A1,A2，用A1存储所有flag=2的骨架点。

2) 从A1中任意选择并剔除一个骨架点，令其flag=3然后置于A2。

3) 遍历A2，计算每个骨架点与A1的所有骨架点之间的距离长度并进行排序，选择最小的那对骨架点。将该对骨架点连接，并将此flag=2的骨架点从A1剔除且令flag=3，然后置于A2。

4) 多次重复步骤3)，直至A1为空。

在实验过程中对步骤3)的连接进行了限制，要求选择的一对骨架点连线不能超出模型外部(通过判断连线四分点是否都处于模型内部)，否则舍弃此种连接，选择其他相离最近的一对骨架点。

需要说明的是，经过距离值局部最大求解和最大内切球约束后，不会出现骨架点距离较近的情况。对于分支较好的模型，当分支骨架点flag=3，也能较好地参与步骤3)的选择连接，因为步骤4)的停止条件是A1为空，即不存在flag=2的骨架点，因此不会出现漏连的情况。对于拓扑结构复杂的模型，处理的效果还不足，将来需要对连接规则进一步的优化。

5 实验结果与分析

本节算法已经在Microsoft visual studio 2005开发环境下实现，利用OSG进行模型的渲染。程序的运行环境为IntelCore2 Duo E7500，2.93 GHz CPU，2 GB内存。为证明笔者方法的实用性，本部分展示了几个模型的骨架提取效果图以及与其他骨架提取方法比较。

图10是从图9骨架点所生成的人体模型骨架，顶点数12 500，面片数目25 000。耗时7.6 s。从图10可以看出：笔者方法所提取的骨架基本上能获取到人体模型的拓扑结构，且能较好地居于模型正中心，因为所有被选取的骨架点都经过了最大内切球的约束，能很好地保证骨架的中心性。

图10 人体模型骨架Fig.10 Human body model skeleton

图11(a)是Torso模型初始图，顶点数目99 447，面片198 384。图11(b)是文献[16]提取效果图，图11(c)是笔者方法提取效果图。从图11可以看出：笔者方法提取的结果图与文献[16]提取的结果图在效果上相近，都能大致反映原始模型的拓扑结构。文献[16]提取的结果图在骨架平滑性上处理的不够，但相比笔者方法，其提取的骨架更居于模型中心。提取Torso模型耗时44 s，主要是因为在体素化模型并去除模型外部体素过程中耗费了大量时间。距离值局部最大求解过程中，笔者方法采用的立方块边长l=2.5×d个平均网格边长度，体素化精度为d=8 个平均网格边长度。

图11 骨架对比Fig.11 Skeleton contrast

图12为笔者方法对horse模型提取效果图，顶点数目19 851，面片数目39 698。图12(a)为初始模型，图12(b，c)分别为Wu等[17]和Lien等[18]提取的骨架效果图，图12(d)为笔者方法效果图。采用的立方块边长l=4×d个平均网格边长度，体素化精度为d=2 个平均网格边长度，消耗时间为7.8 s。从图12可以看出：提取的效果与文献[17]和文献[18]相近，笔者方法提取的骨架基本上能获取到模型的拓扑结构。但笔者方法提取的骨架更接近模型中心位置，比如horse模型腿部位置，而文献[17]和文献[18]做的相对差一些。

图12 骨架对比Fig.12 Skeleton contrast

图13为Liu等[16]、Tagliasacchi等[19]和笔者算法对Homer模型的骨架提取效果图。Homer模型顶点数目5 103，面片数目10 202，笔者采用的立方块边长l=3.0×d个平均网格边长度，体素化精度为d=1.6 个平均网格边长度，耗时2.8 s。从图13可以看出：这3 种方法提取的模型骨架都能较好地体现输入模型的拓扑结构，相比Tagliasacchi等[19]的效果，笔者方法提取的骨架更为完整，图13(c)模型脚部骨架有一定的缺失。Liu等[16]提取的骨架虽然结构完整，但其平滑性不如笔者方法。但Liu等[16]、Tagliasacchi等[19]的鲁棒性较强，能处理许多复杂模型，笔者在这一点上做的还不够，将来需要改进。

图13 骨架对比Fig.13 Skeleton contrast

6 结 论

提出了一种基于距离场的全自动骨架生成方法。首先计算输入的三维模型的包围盒，并根据包围盒尺寸将该空间离散化，分散成一个个立方块空间；之后，通过计算剔除模型外部体素，保留模型内部体素并计算体素的距离信息，从而完成模型的体素化；然后，对存储着距离信息的模型内部体素进行距离值局部最大求解，完成对体素的初步选取；接着，利用最大内切球逼近算法对保留下来的候选体素进行最大内切球约束得到骨架点；最后，根据选择的骨架点，按照一定的规则进行连接形成骨架。笔者方法最大的创新之处在于用最大内切球逼近算法约束候选体素，能有效地去除冗余体素并对候选体素的位置进行一定的矫正。由于在模型体素化过程中要耗费大量时间，对于顶点数较多的模型消耗的时间会相对长些，如图11所示的Torso模型。未来方向将考虑如何提升体素化效率以及处理更多复杂模型。