一种基于体素八叉树的碰撞算法研究

情况。碰撞检测方法主要有层次包围盒法[3]和空间分解法[4]。常见的包围盒算法有AABB包围盒法、Sphere球包围盒法、方向包围盒及K-DOP包围盒法等。包围盒算法的原理是根据物体形状，分别用立方体、圆等多边体将物体包围，并自动检测两个物体包围盒的距离，当距离值小于指定阈值时，则认定两个物体发生碰撞。国内外研究包围盒方法均有一定历史，国内最具代表性的是国防科大[5]提出的基于固定方向的凸包包围盒方法，该方法能适应软组织变形过程中拓扑结构的变化。2014年，Sulaiman等人根据AABB包围盒，提出一种在虚拟狭窄环境中进行阶段性碰撞检测的方法;2013年，Emre等人利用多线程提高碰撞检测效率;2014年，陈琳等人研究一种混合包围碰撞树算法，有效提高了基于三角面片的碰撞检测效率;2015年，彭晏飞[6]等人针对碰撞检测的实时性和逼真度较差的缺陷，提出一种新的混合碰撞检测算法，实验结果证明，与经典的Rapid算法相比，该算法有效地减少了碰撞检测的时间开销，提高了碰撞检测的实时性和真实感;2016年，孙劲光[7]提出一种基于形状分类的包围盒碰撞检测优化算法，实验结果表明，该算法缩短了相交测试的时间，提高了碰撞检测的效率。

空间分解法是为了保证虚拟空间中不规则物体碰撞的精度，原理是将物体所在空间分割成若干个小的立方体，根据物体是否占据同一个空间节点来判断是否发生碰撞。该方法是当前研究的热点。例如，Fan[8]采用线性八叉树结构实现空间碰撞检测;李山[9]等人提出了基于八叉树细分的并行碰撞检测算法;台湾大学的Mei[10]等人利用空间八叉树结构来检测虚拟工具和机器人的碰撞。近年来，有些学者利用硬件加速和并行技术提高算法效率，如浙江工业大学Mei K J等人[11]研究利用GPU加速粒子流体动力学模拟流血效果;青岛大学吴世宇等人利用GPU实现光栅化，有效提高模型的体素化和碰撞检测效率。空间分解法是一种常用的粗略碰撞测试方案，这将大大降低组合测试的数量，因此有利于提高碰撞效率。

总体来看，包围盒法检测规则几何体的实时性和精确性较好。对于复杂几何体，虽然可以利用改进的包围盒法来实现，但精度仍低于实时仿真的标准。空间分解法的优点是精确度高，适用于不规则物体，缺点是内存需求大，计算复杂，实时性差，因此对计算机硬件性能和算法有一定要求。而利用两种方法的混合碰撞算法[12]是提高仿真实时性和精确性的有效方法。

2 碰撞检测

本文提出一种基于哈希八叉树结构的多方向体素化算法。针对实时仿真中复杂软体模型数据量多、计算量大的特点，采用一种基于欧式距离场的多方向快速体素法，并利用哈希表存储八叉树节点及模型离散体素;采用一种结合包围盒和空间分解的碰撞检测算法，通过粗略碰撞和精确碰撞两阶段实现碰撞检测，同时实时反馈碰撞点位置和力。

2.1 八叉树结构和存储

八叉树是处理3D空间常见的方法和数据结构。使用八叉树存储体素可以有效减少体素的高内存需求，且更易理解离散化模型。一个优秀的八叉树数据结构能高效和方便地搜索相邻节点，因此可以有效提高对象内部的体素化和碰撞检测速度。

如图1所示，八叉树根节点是一个轴对齐的立方体包围盒，随后将其沿[x、y、z]轴划分为8个子立方体，这8个子立方体就是子节点，继续按上述方法对子节点进行递归划分，直至达到指定深度。在本文中，八叉树单元中的深度是体素的[2l]倍，[l]为八叉树单元的深度。

常用的传统八叉树数据结构如图2所示，原理如图3所示。*NextNode指向子节点的下一节点地址;*FirstChild由子节点指向下个叶子节点的首地址;*object为模型指针;center是节点中心点，是int整型变量;m_size是采样空间中的实际坐标，原本是double型变量，在这里为了方便计算和哈希存储，扩大1 000倍后存为整型量。

为了快速查询节点，本文采用基于索引的哈希算法来存储数据结构，如图4所示。哈希算法本质上是一种散列表，是一种以“空间换时间”的算法，是关键码值（Key value）而直接进行访问的数据结构（见图5）。也就是说，其通过把关键码值映射到表中一个位置来访问记录，以加快查找速度。这个映射函数叫做散列函数，存放记录的数组叫做散列表。而当使用哈希表进行查询时，就是再次使用哈希函数将key转换为对应的数组下标，并定位到该空间获取value，如此一来，就可以充分利用到数组的定位性能进行数据定位。图4中的数据结构省去两个单链表指针，只新增一个int型的value值，大大节约内存。

本文哈希函数使用乘法哈希：

[Hash（key）=floor（m∗（K∗A-floor（K∗A）））]            （1）

式中，[floor（）]表示小数向下取整;[K]取黄金比例0.618;[m]为210，哈希表使用STL容器中的hasp_map;[A]由center中的[x、y、z]确定，即体素点的相对空间位置，公式为：

[A=n2（z-1）+n（x-1）+y]                           （2）

式中：[n=8l]。

2.2 体素化

体素化可以理解为二维像素化到三维空间的推广，就是以三维采样的方式将模型离散化，以数个正立方体的形式来表现模型，当立方体越小时，模型精度越高，同时计算量也越大。

本文采用一种基于八叉树的体素化算法。该方法分为两个阶段：一是根据八叉树理论采用递归的方法细分空间，直到指定深度;二是提取软组织所占空间网格，即体素化。首先输入3DS格式的体网格模型，再归一化到指定空间，然后对空间进行八叉树细分，再根据欧式距离法判定模型表面体素点，并计算投影向量与三角面片法向量关系来标记模型体素点，最后将体素点镶嵌于八叉树节点中。

为了方便描述，将3D空间简化为2D空间。对于一个二维物体，由于已知模型边界信息，因此考虑从四个方向逐行扫描模型包围盒空间，扫描示意图如图6所示。在扫描过程中，将模型外部体素标记为0，当遇到第一个非0体素点时，则进入下一行，直至完成这一方向。在经过4个方向扫描后，可取得模型外部体素点的集合。对于整个包围盒数据，除去外部体素，剩余则为模型表面和内部体素，经过渲染则实现了物体的体素化。

另外，判断每一扫描行单个体素点是否为模型边界，需要进行基于欧式法的距离判定和三角形投影判定，具体步骤如下：第一步：预计算模型表面三角片所在平面、法向量等信息并存储;第二步：根据欧式距离法，计算体素点与三角片距离，当该距离小于指定阈值时，表示体素点在模型表面，进入第三步，否则循环判断下一个体素点;第三步：计算体素中心点在三角片上的投影向量，通过计算投影向量与三角片法向量的乘积，判断该点是否在三角片内，如果为真，表示该点是第一个模型边界点，循环进入下一行。

2.3 碰撞算法

本文的碰撞算法是结合基于空间分解法的粗略碰撞和基本图元法的精确碰撞[13]。粗略碰撞通过计算刚体与节点中心点距离判断刚体末端与软体是否在同一节点内;精确碰撞通过计算刚体末端与软体表面三角形单元距离判断是否产生碰撞。

2.3.1 粗略碰撞检测。该阶段使用已构造的空间Hash表，每个八叉树节点存在对应的Hash表中。由于八叉树节点数量远低于体素点和模型点的数量，因此，在粗略碰撞阶段，可以快速剔除远离潜在碰撞点的八叉樹节点，从而有效加快碰撞检测速度。

为了方便描述粗略碰撞的原理，这里用2D图表示3D场景（见图7）。设点[A（x0，y0，z0）]为虚拟工具末端初始点，点[B（x1，y1，z1）]为虚拟工具末端移动终点。图中黑色正方形为八叉树节点，边长[l=12d]，[d]为八叉树深度，点[O（x2，y2，z2）]为其中心点，灰色为模型体素点。当虚拟工具在点[A]时，与点[O]的距离为：

[distAO=（x0-x1）2+（y0-y1）2+（z0-z1）2]          （3）

此时，[distAO>22l]，则判定不存在粗略碰撞。

当虚拟工具末端点在[B]时，与点[O]距离为：

[distBO=（x2-x1）2+（y2-y1）2+（z2-z1）2]           （4）

此时，[distBO<22l]，则判定存在粗略碰撞，进入精确碰撞检测阶段。

2.3.2 精确碰撞检测。精确碰撞检测是当粗略碰撞检测发生后，如果虚拟工具末端点与软组织在同一节点内，则进行基本图元检测，称为点与三角形的相交测试，具体是判断虚拟工具与软组织表面是否产生碰撞。点与三角形测试使用较多的方法是未预处理三角形的计算质心法和经过预计算的三角形法两种。

本文使用另外一种方法：针对三角形[ABC]和一点[P]，先进行平移操作，使得[P]点位于原点，然后测试原点是否位于平移后的三角形内部。当且仅当三角形[PAB]、[PBC]、[PCA]同为顺时针或逆时针排列时，顶点[P]位于三角形[ABC]内部（见图8）。

经过平移后，顶点[P]即为原点[O]，当发生碰撞时，由于[O]在三角形内部，所以该方法等价于计算式（5）的三个公式是否指向同一方向，即测试[u·v]是否大于等于0和[u·w]是否大于等于0。该方法的代码如图9所示。

[u=B×Cv=C×Aw=A×B]                               （5）

由于叉积成本较大，根据拉格朗日恒等式记为5点积形式，该方法代码如图10所示。

3 实验结果

本文对上述算法进行实验验证，实验使用的PC机是Win7、4.0GB的64位操作系统，利用双线程实现基于CHAI3D和OPENGL开源库的力觉和视觉渲染。

针对不同模型进行不同分辨率下的体素化，得到仿真效果图如图11所示，体素化算法测试结果如表1所示。由图11可知，该算法可以实现正常的视觉效果;由表1可知，在128和256体素率下，模型均可以较好地实现体素化。

碰撞检测实验过程：用户操作力反馈设备末端执行器（见图12），位姿映射到虚拟场景中刚体模型位姿，控制刚体模型触碰软体模型并进行按压操作，界面实时显示刚体的相对位姿、反馈力及刷新频率。

图13为碰撞检测效果图，表2是碰撞检测结果。通过表2可知，不断移动虚拟工具，刷新频率平均在600Hz以上，可以满足视觉要求，检测时间在毫秒级，可以满足实时性要求，同时可以实时反馈碰撞点的位置和产生的力。

4 结语

在实时仿真中，体素化算法处理复杂曲面软体模型和高分辨率模型时预加载时间过长，后期需要进一步优化算法和利用硬件GPU提升算法性能。另外，在刚体模型移动过程中，需要遍历整个软体模型，耗时较長，可以利用领域查找等方法进一步提高实时性。

参考文献：

[1]王文杰，刘国兵，陶磊.一种基于多传感器的客车防碰撞预警系统[J].河南科技，2018（9）：115-117.

[2]沈华，陈金良，周志靖.混合空域中无人机飞行防相撞技术[J].指挥信息系统与技术，2016（6）：24-29.

[3]左军涛，朱恩成，黄四牛.基于GPU胡航迹快速计算方法[J].指挥信息系统与技术，2011（4）：55-59.

[4]王杰，田宏安.无人机融入非隔离空域感知与规避技术[J].指挥信息系统与技术，2017（1）：27-32.

[5]龚随，侯进，钟李涛.基于椭球拟合胡人体一服装碰撞检测方法[J].计算机应用2，2019（1）：274-283.

[6]彭晏飞，卢真真.基于空间剖分和包围盒的快速碰撞检测算法[J].计算机应用与软件，2015（8）：150-153，165.

[7]孙劲光，吴素红，周积林.基于形状分类的包围盒碰撞检测优化算法[J].计算机应用与软件，2016（2）：242-245.

[8] Fan W S， Wang B， Paul J C， et al. An octree-based proxy for collision detection in large-scale particle systems[J]. Science China-Information Sciences，2013（1）：1-10.

[9]李山，赵伟，李菲.一种基于八叉树与流水线技术的快速碰撞检测算法[J].计算机与现代化，2011（1）：20-24.

[10]张璐鹏.适用于柔性体切割仿真的八叉树体模型生成算法[D].青岛：青岛大学，2018.

[11] Mei K J， Rong S. Collision detection for virtual machine tools and virtual robot arms using the Shared Triangles Extended Octrees method[J]. International Journal of Computer Integrated Manufacturing，2016（4）：355-373.

[12]鲍义东，吴冬梅.粘弹性软组织建模和切割及虚拟仿真实验研究[D].哈尔滨：哈尔滨工业大学，2016.

[13]朱卓，李宏胜.虚拟手术中人体软组织超粘弹性建模及穿刺仿真[J].应用力学学报，2019（2）：304-309.