张璐鹏　贾世宇

摘要： 针对任何三角几何模型的体素化需求，本文基于距离场和八叉树结构的特性，提出将任意三角几何模型离散为体素模型的算法。提取任意三角几何模型的无符号距离场的偏移表面，自动删除内部表面，然后从偏移表面计算符号距离场，并调整阈值，得到与原模型接近的重构表面，最后对重构表面进行体素化，得到适用于切割仿真的体素模型，为验证算法的有效性，本文对基于空间分割的体素化算法和本文提出的算法进行对比测试。测试结果表明，该算法对于包含自相交和非连通的输入三角几何的体素化效果更好，更为贴近原模型，且时间和内存成本较低，效果较好。该算法为仿真柔性体切割的体素模型提供一种有效方案。

关键词： 距离场； 八叉树； 非流形； 自相交； 三角几何； 体素化

中图分类号： TP391.41文献标识码： A

近年来，隨着可编程图形处理器和计算机图形学的发展，医学仿真在医学领域的应用越来越广泛。医学仿真中的手术仿真通过精确的人体组织器官模型，真实模拟组织器官在手术器械的外力交互作用下变形乃至被切割等操作过程，而为实现人体器官柔性体变形和切割的真实感，模型需要被体素化后才能够实现相应的数学计算进行仿真。目前，用于生成体素模型的方法有很多。Huang J等人[1]提出利用模型法线到表面函数作为距离标准，将模型离散为体素模型的方法；S.Oomes等人[2]采用尺度空间理论（Scalespace Theory）转换模型为反锯齿体素数据；S.Laine[3]提出基于拓扑结构的体素化方法。以上研究适用于绝大部分几何模型体素化的方法，虽然这些方法精度和速度较高，但体素化应用对象只用于流形几何，普适性较差。针对人体器官模型无定向、非流形或包含自相交的三角几何网格的特点，本文提出基于距离场和表面重构算法填补模型表面孔洞，界定模型在数学意义上的内外之分，然后对模型进行体素化处理，从而产生适用于切割仿真的体素模型。

1相关工作

1.1距离场

距离场的概念由Rosenfeld和Pfaltz在1966年与数字图像处理相关的论文中提出[4]。距离场存在有无符号区分，符号指的是某个点位于模型内部还是外部。只有当几何模型为密封的流形几何体时，才能计算出其符号距离场，而当几何模型为非流形时，便无法明确符号的定义，也就无法从数学意义上界定模型的内外。

A.Fuhrmann等人[5]提出通过分析三角网格法线方向产生的棱柱来快速计算局部距离场的方法；Rong G等人[6]提出在常数时间内计算输入网格近似距离场的jump flooding算法；S.P.Mauch[7]根据扫描转换解决静态HamiltonJacobi方程，进而计算精确距离场的CSC算法等。上述算法全部建立在封闭流形几何的基础上，无法计算非流形几何体的距离场。由于没有明确定义划分非流形几何体的内部、外部和边界，所以在计算距离场时，通常通过从三角网格多个方向进行投影，生成二进制体数据[8]，如果三角网格包含孔洞或自交，将会体现在扫描转换后的结果中。虽然可以通过对每个体素的扫描转换，并使用奇偶规则判定重构合理的体数据[9]，但仍然存在几何体结构化的异常值处理困难的问题。

1.2体素化

1987年，A.Kaufman[10]首次提出体素化概念，根据三维模型的网格数据，转换为包含模型表面信息和内部属性的体数据集。近年来，相关学者们提出了大量的体素化算法，J.Huang等人[10]提出的保证拓扑一致性的体素化方法；Fang S等人[11]利用硬件加速方法实现体素化算法，但这些算法通常需要大量内存，并且输出一个非缩放的三维体素网格；M.Schwarz等人[12]提出一种基于八叉树的稀疏固体体素化算法，虽然内存要求大大降低，但是该方法仍受限于可用内存。

随着可编程图形处理器的快速发展，研究者们利用计算机图形处理器（graphics processing unit，GPU）加速改进体素化方法。C.Crassin等人[13]使用GPU光栅单元进行稀疏体素化；M.Ptzold等人[14]提出一种非网格化（gridfree）、非核心（outofcore）的GPU体素化方法。这些方法虽然利用GPU并行特性进行算法改进，但是仍然存在算法判断过程复杂，内存消耗较大，对硬件要求较高的缺陷。

2八叉树体模型生成算法

针对包含自相交、孔洞、无定向、不连通及重复的几何体在内的任何三角几何模型的体素化需求，本文提出一种基于八叉树的体模型生成算法。该方法分为三个阶段：一是计算得到输入三角网格的符号距离场，并重构模型表面网格；二是利用保证拓扑结构的Marching Cubes算法修复重构的表面网格，得到封闭的流形表面网格；三是利用八叉树对模型空间进行划分和扫描，区分模型内外部数据，进而获取体数据，实现模型体素化。

2.1计算距离场

距离场是一个标量场，表示场内每个点到给定对象上最近点的距离，该距离有符号，表示该点位于模型的内部或外部。定义为

dist（x，Ω）=inf‖p-x‖， （x∈R3，p∈Ω）（1）

式中，dist（x，Ω）表示从点x到模型Ω的无符号距离场；p是属于模型Ω表面网格上的点；inf‖p-x‖表示点p到x的距离下界；R3表示距离场内的点的集合。符号距离场是在原有无符号距离场的基础上，使用正负号表示点x位于三角几何模型Ω的内部或外部。通常情况下，判定曲面内部点的符号为负，曲面外部点的符号为正，若用σ表示模型Ω的内部距离场的点集，则符号距离场φx的定义为

2.1.1提取无符号距离场

计算无符号距离场前进行空间划分[16]，首先将输入三角网格用平行于坐标轴的包围盒包含，随后根据距离场偏移距离σ扩展已有包围盒，从而生成采样空间。根据所定义的划分分辨率h，将整个采样空间均匀划分，模型包围盒用固定值均匀划分。均匀空间划分如图1所示。

依次计算每个采样点到包围盒立方体的最小距离和最大距离，取最大距离中的最小值作为采样点到三角网格的最小距离的参考值，随后遍历所有包围盒立方体，比较参考值和采样点到包围盒立方体的最小距离。如果当前立方体的最小距离大于参考值，则忽略当前立方体；若小于参考值，则计算采样点到当前立方体中所有三角面的距离，并更新最小值，最终得到无符号距离场。以无符号距离场作为输入，使用Marching Cubes算法产生一个流形的偏移曲面，偏移曲面定义为

Sσ={x∈R3|dU（x，Ω）=σ}（3）

式中，dU（x，Ω）表示从采样点x到输入几何Ω的无符号距离场。生成的偏移曲面可能存在嵌套，其定义为

Sσ=（∏iEiσ）∏（∏iIiσ）（4）

式中，Eiσ代表外部组件；Iiσ表示内部组件；∏表示不相交的集合。使用UnionFind算法检测所有连通分支，去除所有內部组件，得到外部偏移曲面。

2.1.2计算符号距离场

通过计算得到一个流形的偏移曲面Eσ和输入几何的无符号距离场，通过使用伪法线测试[15]计算偏移曲面Eσ的符号距离场。由于蛮力执行伪法线测试运行效率较低，所以本文利用文献[16]证明所得的两条引理。

引理1如果x为偏移曲面Eσ外的一点，则符号距离场dS定义为

dS（x，Eσ）=dU（x，Ω）-σ≥0（5）

引理2如果点x满足dU（x，Ω）<σ，则点x一定位于偏移曲面Eσ的内部。

引理1表明在计算偏移曲面Eσ的外部符号距离场时，再次使用之前所计算的符号距离场，极大程度地缩减计算量。但是对于内部而言，当点x在内部无限接近于偏移曲面Eσ时，引理1的等式不成立，所以之前计算的无符号距离场不能应用于内部符号距离场的计算。针对这一问题，由引理2可知，网格上所有满足dU（x，Ω）<σ的点均无需进行伪法线测试。

首先对偏移曲面Eσ进行均匀的空间分割及边界检测 ，标记所有和偏移曲面相交的立方体，然后再次计算偏移曲面的无符号距离场，并执行S遍历，S遍历后可得到曲面网格点的符号，最后在计算出Eσ的符号距离场基础上，将其偏移-σ产生相对于输入三角网格的符号距离场。

2.2修复重构表面

为实现任意几何体的体素化，得到输入三角网格的符号距离场后，执行Marching Cubes算法可以得到输入几何体的重构表面Fσ。当输入几何趋于封闭流形的非流形几何时，重构表面Fσ是一个封闭的流形几何，可以直接进行体素化。但是，并非所有的非流形几何体都趋近于流形，当一个非流形几何由许多不连通的三角网格组成时，初次重构的表面效果不足以体素化。偏移曲面图如图2所示。

为解决上述问题，本文在已有重构表面的基础上进行修复。已知重构表面Fσ是以符号距离场作为等值面执行Marching Cubes所得，当输入三角网格本身为不连通的情况时，由于计算所得为精确符号距离场，所以距离场也是不连通的，导致重构表面效果较差。此时，精确的符号距离场定义为

dU（x，Ω）=σ， σ→0（6）

当符号距离场能够连通并覆盖整个输入几何体时，Marching Cubes的重构表面才能达到足以进行体素化效果，所以通过调整偏移量，让符号距离场覆盖整个输入几何体，即

dU（x，Ω）=σ+δ（7）

式中，δ作为调整的偏移量由用户输入。调整距离场后，再次利用Marching Cubes生成封闭的流形重构表面。

2.3体素化

现将封闭的流形重构表面体素化。针对重构表面再次构成包围盒，并在顶点（xmin，ymin，zmin）固定与欧式空间方向相同的坐标系，此时包围盒空间定义为体素空间V3，每个体素由一组唯一整数序列（i，j，k）确定，且体素空间内也存在“6邻接”、“18邻接”和“26邻接”三种拓扑关系[16]。体素空间确定后，基于八叉树原理进行空间划分[17]，一般八叉树划分方法仅对灰节点进行划分，本文则对3种节点均进行划分操作，生成均匀的体素网格，且体素标志位满足如下关系

f（i，j，k）=-1， 体素位于模型边界外0， 体素与模型边界相交1， 体素位于模型边界内（8）

通过空间细分获取模型的表面数据，随后利用逐行扫描根据模型表面数据快速区分模型内外部数据，从而获得体数据集，实现模型体素化。

3实验结果与分析

为验证算法的有效性，本文对基于空间分割的体素化算法和本文提出的算法进行对比测试。基于空间分割的体素化算法实质是基于MIN提出的使用空间分割快速计算距离场算法的基础上重构表面并体素化的方法，由于此方法对趋于流形几何的非流形几何体的体素化效果较好，因此用于作对比测试。斯坦福兔子模型、骷髅模型和棕榈植物模型的无符号距离场和符号距离场的测试结果如表1所示。

所有符号距离场均使用八核计算，3种模型均包含自相交、孔洞和非连通特性的三角几何。体素化效果对比如图3所示。本算法合理解决了具有精细和超薄特性的非流形几何体的符号距离场计算问题。图3b和图3c展示斯坦福兔子模型分别使用基于空间分割的体素化算法和本文算法的体素化效果图，可见两种算法在对趋近于封闭流形的几何模型进行体素化时，所产生的效果并无较大差别；图3e～i展示骷髅模型和棕榈树模型的体素化效果对比，经过表面修复后，本文算法对于包含自相交和非连通的输入三角几何的体素化效果更好，体素化后的模型更为贴近原模型，而基于空间分割的体素化算法的体素化模型则明显存在失真问题，无法精确还原模型。

4结束语

基于距离场和八叉树结构的特性，本文提出将任意三角几何模型离散为体素模型的算法，解决了非流形几何难以体素化的问题，采用表面修复，改善了精确距离场对重构“三角形汤”几何偏移曲面不连通的缺陷。实验结果表明，对非流形几何进行体素化时，相对基于空间分割的体素化算法，本算法能够实现更好的效果，考虑到三角几何结构的不确定性问题，可以适用于任意三角几何的体素化，对通用体素化算法的研究具有一定参考价值。本文采用Marching Cubes算法进行表面重构，虽然保证了拓扑结构的一致性，但仍存在等值面精度较低的问题，导致修复效果并不完美。下一步研究可致力于改善算法效率和提升算法精度，利用GPU加速运算、改善Marching Cubes算法或者选用更为高效的修复算法，使体素化的效果更为完美。

参考文献：

[1]Huang J， Yagel R， Kurzion Y. An Accurate Method To Voxelize Polygonal Meshes[C]∥IEEE Symposium on Volume Visualization. Chapel Hill， North Carolina， USA： IEEE， 1998： 119126.

[2]Oomes S， Snoeren P， Dijkstra T. 3D Shape Representation： Transforming Polygons into Voxels[J]. International Conference on Scalespace Theory in Computer Vision， 1997， 1252： 349352.

[3]Laine S. A Topological Approach to Voxelization[J]. Computer Graphics Forum， 2013， 32（4）： 7786.

[4]Rosenfeld A， Pfaltz J. Sequential Operations in Digital Picture Processing[J]. Journal of the ACM， 1966， 13（4）： 471494.

[5]Fuhrmann A， Sobotka G. Distance Fields for Rapid Collision Detection in Physically Based Modeling[C]∥Proceedings of Graphi Con 2003. Moscow， Russia： IEEE， 2003： 5865.

[6]Rong G， Tan T S. Variants of Jump Flooding Algorithm for Computing Discrete Voronoi Diagrams[C]∥4th International Symposium on. Voronoi Diagrams in Science and Engineering. Clamorgan， UK： IEEE， 2007： 176181.

[7]Mauch S P. Efficient Algorithms for Solving Static HamiltonJacobi Equations[M]. USA： Mauch Sean Patrick， 2003.

[8]Spillmann J， Wagner M， Teschner M. Robust Tetrahedral Meshing of Triangle Soups[C]∥Vision， Modeling， Visualization （VMV）. Oxford， UK： IEEE， 2006： 916.

[9]Kaufman A， Shimony E. 3D ScanConversion Algorithms for VoxelBased Graphics[C]∥Proceedings of the 1986 Workshop on Interactive 3D graphics. Chapel Hill， North Carolina， USA： IEEE， 1987： 4575.

[10]Huang J， Yagel R， Filippov V， et al. An Accurate Method for Voxelizing Polygon Meshes[C]∥IEEE Symposium on Volume Visualization. Chapel Hill， North Carolina， USA： IEEE， 1998： 119126.

[11]Fang S， Chen H. Hardware Accelerated Voxelization[J]. Computers & Graphics， 2000， 24（3）： 433442.

[12]Schwarz M， Seidel H P. Fast Parallel Surface and Solid Voxelization on GPUs[J]. Acm Transactions on Graphics， 2010， 29（6）： 110.

[13]Crassin C， Neyret F， Sainz M， et al. Interactive Indirect Illumination Using Voxel Cone Tracing[J]. Symposium on Interactive 3d Graphics & Game， 2011， 30（7）： 19211930.

[14]Ptzold M， Kolb A. GridFree OutofCore Voxelization to Sparse Voxel Octrees on GPU[C]∥Proceedings of the 7th Conference on HighPerformance Graphics. Dublin， Ireland： ACM， 2015： 95103.

[15]Nooruddin F S， Turk G. Simplification and Repair of Polygonal Models Using Volumetric Techniques[J]. IEEE Transactions on Visualization and Computer Graphics， 2003， 9（2）： 191205.

[16]Xu H， Barbi J. Signed Distance Fields for Polygon Soup Meshes[C]∥Proceedings of Graphics Interface 2014. Processing Society， Montreal， Q C， Canada： ACM. 2014： 3541.

[17]Min D， Jia S， Zhang X. Multithreaded Fast Distance Field Computation Using Spatial Partition[J]. Science & Technology Vision， 2015（8）： 153154.

[18]Brentzen J A， Aanaes H. Signed Distance Computation Using the Angle Weighted Pseudonormal[J]. IEEE Transactions on Visualization and Computer Graphics， 2005， 11（3）： 243253.

[19]Meagher D. Geometric Modeling Using Octree Encoding[J]. Computer Graphics and Image Processing， 1982， 19（2）： 129147.

[20]吳晓军， 刘伟军， 王天然. 基于八叉树的三维网格模型体素化方法[J]. 图学学报， 2005， 26（4）： 17.