基于散乱点云的快速体积计算法

胡晓彤，陶森柏

(天津科技大学计算机科学与信息工程学院，天津 300222)

基于散乱点云的快速体积计算法

胡晓彤，陶森柏

(天津科技大学计算机科学与信息工程学院，天津 300222)

三维可视化体积计算基本上都是先由散乱点云构建出表面网格模型，然后基于网格模型计算体积，存在计算量大、速度慢的缺点．针对此问题提出一种快速体积计算法，首先使用改进的增量式 Delaunay三角剖分对散乱点云进行四面体剖分；然后利用 K近邻计算散乱点的拟合曲面和最小生成树，得到各点的法向量；由各点法向量剔除体外四面体；最后计算各四面体体积之和从而得到总体积．实验表明，该算法不仅保证了计算准确度，而且较传统算法大大提高了效率．

散乱点云；四面体剖分；Delaunay三角剖分；法向量；K近邻

Abstract：Visual volume calculation in 3D space basically is based on mesh model nowadays which is constructed from scattered point cloud. It exposes inferiority like huge calculation and low speed when volume is needed only. According to that，an algorithm of rapid volume calculation was proposed. First，the convex hull of point cloud was subdivided into tetrahedron with improved incremental Delaunay triangulation. Second，fitting quadric surface and MST of points with KNN were calculated to get normal vectors. Third，those tetrahedron in vitro were removed by using normal vectors. Finally，all tetrahedron’s volume were added up to get object’s volume. Experimental results show that this method can improve efficiency greatly compared with common methods and keep precision of volume calculation.

Keywords：scattered point cloud；tetrahedron subdivision；Delaunay triangulation；normal vector；KNN

科学计算可视化(visualization in scientific computing)是 20世纪 80年代后期提出并发展起来的一个新领域．它指的是运用计算机图形及图像技术，将科学计算过程及结果的数据转化为屏幕上的图形及图像并进行交互处理的理论、方法和技术．对于可视化出来的图形进行交互测量在医学、地质勘探、天体物理等领域都有重要意义．

要实现可视化交互测量，先要获得目标物的三维点云，这通常基于机器视觉或设备采集实现．基于机器视觉即由多个摄像机获得图像序列，通过图像匹配，计算目标物的三维坐标；基于设备采集指由三坐标仪或激光扫描仪等设备采集散乱点云．然后进行测量，现行方法基本上都是基于不规则网格模型的计算，其中三角网格模型使用最为广泛．任意网格模型体积计算的常见方法有：坐标法[1]，即由封闭空间的表面基元的坐标行列式计算体积；切片法[2]，用一组互相平行的平面剖切模型，由剖面面积和平面间距计算体积；投影法[3]，计算三角网格模型的三角面片及其在投影平面上的投影所围成的凸五面体的带符号体积，所有带符号体积的代数和即为总体积．

这些方法都需要先重建出表面网格模型，然后在此基础上计算体积．虽然体积计算算法本身的精度和效率都较高，但是整个过程大部分时间耗费在了重建表面网格模型上．在单纯需要得到目标体积的情况下，把时间浪费在重建表面网格模型上是不需要的．因此，本文提出了一种基于不附带任何拓扑信息的散乱点云的四面体剖分体积计算法，该算法不需要重建表面网格模型，可直接由点云得到体积，因而提升了体积计算的效率．

1 四面体剖分体积算法

1.1 基于凸包分割的四面体剖分

四面体剖分本质上就是将点集S⊆R3的凸包分割成若干以 S中的点为顶点的四面体．传统的四面体剖分算法在浮点数系统中由于退化现象而不够健壮，所以本文采用基于增量式 Delaunay三角剖分的改进算法．Edelsbrunner和Seidel于1986年证明：任何维度为d的有限点集S⊆Rd，其 Delaunay三角剖分都可以从的凸包中获得[4]．

1.1.1 基本概念

(1)凸包：是指包含点集S⊆R3中任意两点连成的线段的最小凸多面体．

(2)Delaunay规则：三维空间中的四面体的外接球内不包含其他的散乱点．

(3)退化现象：在浮点数系统中，因计算机数值精度有限而导致拓扑错误．

1.1.2 算法

(1)三维 Delaunay三角剖分是由二维 Delaunay三角剖分扩展而来．首先构建一个初始四面体，形成初始化四面体网格．

(2)将散乱点插入当前四面体网格中，对于输入点P，使用随机行走方法来寻找包含P的四面体．先指定一个四面体T，如果 P位于该四面体内，则完成行走．如果不在四面体内，则随机指定一个三角面E，如果E所在的平面将T和P分割开(即T和P在平面的两边)，下一个访问的四面体就是共享E的邻近四面体；否则，就按预定的顺序遍历其他的面，直到找到分割开T和P的面．

(3)找到包含P的四面体，则分割该四面体为4个小四面体．

(4)如果P位于当前四面体网格外，则选择网格的一个可见面(即P在面的一侧)，连接P与该三角面的3个顶点构成新的四面体加入到四面体网格中，选择可见面时，尽量避免使新生成的四面体是狭长的.

(5)重复步骤(2)—(4)，直到所有散乱点都被插入四面体网格．

(6)验证Delaunay三角剖分的有效性．首先检查Delaunay三角剖分数据结构的连贯性，即四面体的邻接关系．然后验证各四面体的方向和由 Delaunay三角剖分获得的凸包的正确性．

该算法的改进主要在于步骤(2)的行走方法．传统算法使用的是线性行走方法．对于输入点 P，指定一个四面体，如果 P位于该四面体内，则完成行走；如果不在该四面体内，就构建一条射线，原点是当前四面体的体内一点，标记为 C，它的方向为 C→P，定位到与该射线相交的一个四面体的面，与当前四面体相邻并共用这个面的四面体就是包含 P的下一个候选四面体．这个方法虽然快捷，但有个缺陷，即当射线穿过四面体的顶点或边时与射线相交的下一个四面体不是相邻的四面体，导致 P的定位错误，而随机行走方法有效地避免了该现象[5]．

1.1.3 四面体剖分实例

测试用的部分散乱点云实例见图 1，点云数据来自 CGAL标准图库[6]．图 2为由图 1(b)散乱点云生成的物体凸包以及四面体网格．计算所有四面体的体积总和就可以得到物体凸包的体积，但是，凸包中还存在大量体外四面体，需要对体外四面体进行剔除，才能够实现准确的体积计算．

图1 测试用散乱点云Fig.1 Scattered point cloud for test

图2 四面体剖分结果Fig.2 Results of tetrahedron subdivision

上述算法在最坏情况下的时间复杂度为O(n2)，在输入点集为一般物体表面点集的情况下，算法复杂度接近O(n lg n)．使用 C++语言实现以上算法进行四面体剖分的耗时数据见表 1，运行环境是 Windows XP,32位系统，编译环境是 Visual Studio 2008，编译器为VC++9.0，硬件为Core2 Duo E7400 2.80 GHz，2.0,GB DDR2 800 MHz RAM．

表1 四面体剖分耗时Tab.1 Time-consuming of tetrahedron subdivision

1.2 法向量的计算

由于四面体剖分实质上是对凸包的分割，因此部分四面体是体外的冗余数据，需要剔除．由于散乱点云本身不包含任何拓扑信息，直接确定体外四面体比较困难，因此本文借由计算散乱点—四面体的顶点的法向量来描述四面体所处位置，从而剔除体外四面体．上述计算主要分为两步：基于拟合曲面计算法向量和法向量一致化．

1.2.1 基于拟合曲面计算法向量

为了计算各点的法向量，用二次曲面来拟合点和它的K近邻(指在欧几里德距离下，包含n个点的数据集S⊆Rd中，找到一个点的 K个最近点)．本文使用最小二乘法计算二次曲面参数，二次曲面方程可表示为

对于给定散乱点(xi，yi，zi)i＝1，2，…，N，使总误差 Q最小．

根据式(3)可得到方程的参数，得到二次曲面方程后就可以通过偏微分来计算该点的法向量[7]．

1.2.2 法向量一致化

因为初始法向量的指向内外不一，得到初始法向量后需要一致化调整，这个问题可以模型化为图的优化．这个图的节点由各散乱点构成，如果两个节点足够近的话则构建一条边，边的权值用 ni·nj表示，这样法向量一致化就可归结为图的总权值最大化问题，这里的主要问题是选取连接图中的哪对节点．由于物体表面可以认为是单连通元，所以图应该是连通的．最小生成树就是一个连接邻近点的简单连通图，因此基于散乱点的欧几里得距离建立最小生成树．由于初始最小生成树边的密度不能满足本文要求，因此需要添加边，如果两个节点中的一个点是另一个点的 K近邻则连接两点，由此生成的图称为黎曼图(Riemannian graph)；然后，再基于该图生成最小生成树，选择一个种子法向量，按深度优先准则在最小生成树中传播调整法向量方向[8]．调整规则就是：对于最小生成树中的两个近邻点，如果相应法向量的点积 ni·nj＜0，则 ni或 nj应当被反向．

这一步的时间复杂度为O(n2)，近邻个数 K的选择是影响法向量精度的一个重要因素，由于对象三维拓扑的复杂度不同，可将 K设置成可调整的变量．根据经验，在大部分表面比较平滑的情况下 K可以取偏小的值；如果表面纹理复杂，K可以取偏大的值．

图 3是奶牛点云的法向量效果图(经过多次测试，这里将 K设置为 18，效果较好)．可以看到，在散乱点比较密集的头部和局部平滑的躯干部获得的法向量准确度比较高，而散乱点比较稀疏的腿部和尾部得到的法向量则准确度相对低一些，部分点很难获得指向正确的法向量．

图3 奶牛点云的最终法向量Fig.3 Ultimate normal vectors of milk cow point cloud

1.3 体外四面体的剔除

获得点云的法向量以后，体外和体内四面体可通过法向量的指向判断：体内四面体顶点的法向量指向都背离四面体，而体外四面体顶点的法向量则指向四面体内部或附近区域．实验表明，在获得正确法向量的前提下，利用顶点的法向量是否和四面体外接球相交来剔除体外四面体可以取得理想的结果．具体算法是：遍历四面体，计算每个四面体的外接球，然后从四面体的每个顶点沿着法向量方向引出一条射线，如果4个顶点的射线都与外接球有交点(不包含射线原点)，可以认为该四面体为体外四面体，不予考虑；否则，计算该四面体的体积并加入总体积．如图 4(b)所示为图 4(a)中四面体的外接球和 4个顶点的沿法向量指向的射线，射线与外接球有 4个交点，就可以判断其为体外四面体．

图4 四面体剔除原理图Fig.4 Schematic of tetrahedron removing

这一步的时间复杂度为 O(n)，经大量程序验证，在能得到基本正确的法向量的前提下，本算法能剔除绝大部分体外四面体，证明了上述方法的可行性．

1.4 体积的计算

经过以上步骤，最后计算物体体积就变得非常简单，只需要根据式(4)计算形体内各四面体的体积之和，就可以得到物体的体积

2 实 验

基于上述方法，对图 1中散乱点云进行了测试，并与传统算法比较．由于传统算法的用时主要在构建表面网格模型过程，所以比较时忽略了传统算法的后续体积计算用时．本文选择被广泛使用的 Poisson重建法[9]重建表面网格模型，其将表面重建问题表示为 Poisson方程的解，通过估计模型的指示函数和提取等值面，对表面重建一个无缝的三角逼近．这样做有诸多优点，很多对隐式表面进行拟合的方法把数据分割到不同区域进行局部拟合，然后使用合成函数进一步合成这些局部拟合结果，相比而言，Poisson重建是一个全局方法，不借助于启发式的分割或合并，这样它能创建较光滑的表面，稳健地近似含有噪声的数据.

先使用本文算法计算物体体积，再使用 Poisson方法重建网格模型，比较两者耗时，然后比较散乱点云提供者给出的体积信息与本文计算结果，如表2所示．对于 3个实例的散乱点云，本算法的总计算时间仅为Poisson重建时间的20.9%～40.2%，准确率达到94.93%～98.07%．

表2 算法验证与比较Tab.2 Algorithm verifying and comparison

3 结 语

针对传统算法先构建表面网格模型再计算物体体积的流程，本文提出一种基于散乱点云直接计算物体体积的方法．本算法通过基于凸包分割的四面体剖分、法向量的计算、体外四面体的剔除和体积的计算等步骤，在保证准确度的前提下，大幅度提高了效率，使实时的快速体积计算成为可能．

算法的准确率和效率还有待提高，考虑在复杂点云的法向量计算及体外四面体剔除等方面作深入研究，进一步提高算法的准确率和效率．

［1］ 韦进. 三维空间任意多面体体积的一种坐标计算法[J]. 湖州师专学报：自然科学，1997，19(5)：67–73.

［2］ 庞明勇，戴文俊. 基于体积分布特征匹配的三维实体网格模型检索[J]. 系统仿真学报，2007，19(1)：30–34.

［3］ 王泉德. 任意三角网格模型体积的快速精确计算方法[J]. 计算机工程与应用，2009，45(18)：32–34，58.

［4］ Edelsbrunner H，Seidel R. Voronoi diagrams and arrangements[J]. Discrete & Computational Geometry，1986(1)：25–44.

［5］ Devillers O，Pion S，Teillaud M. Walking in a triangulartion[J]. International Journal of Foundations on Computer Science，2002(13)：181–199.

［6］ European Research Programs. CGAL[EB/OL]. [2010–03–20]. http：//www. cgal. org/download. html.

［7］ 魏永超，苏显渝. 基于方向角的散乱点云三角剖分算法[J]. 四川大学学报：工程科学，2009，41(4)：202–207.

［8］ Hoppe H，DeRose T，Duchamp T，et al. Surface reconstruction from unorganized points[J]. Computer Graphics，1992，71(8)：71–77.

［9］ Kazhdan M，Bolitho M，Hoppe H. Poisson surface reconstruction[C]// Polthier K，Sheffer A. Symposium on Geometry Processing. Switzerland：The Eurographics Association，2006：61–70.

Algorithm of Rapid Volume Calculation Based on Scattered Point Cloud

HU Xiao-tong，TAO Sen-bai
(College of Computer Science and Information Engineering，Tianjin University of Science & Technology，Tianjin 300222，China)

TP391

A

1672-6510(2011)01-0067-05

2010–10–18；

2010–12–03

胡晓彤（1971—），男，北京人，副教授，huxt@tust.edu.cn.