摘要
针对大规模粒子高表现可视化需求，提出基于包络面重构的大规模粒子并行绘制算法。该算法以连续曲面的形式表示,绘制大规模粒子的团簇表面及其物理量分布。对算法进行了分布式并行化，从而可以通过大规模并行来处理亿以上规模的粒子数据。在算法实现上，还解决了并行计算时的块间裂缝问题，并提出了快速查找邻域粒子的方法，同时，基于可见性对粒子数据进行剔除，提高了绘制效率。由此，可以通过带光照效果的光滑曲面来高表现展示大规模粒子数据中的团簇结构及其物理量分布。实验结果表明，该算法在512核上可在5 s内完成上亿粒子的绘制，并行效率可达60%。该算法已成功应用到大规模并行非平衡分子动力学模拟等实际模拟应用中。
Abstract
A parallel rendering algorithm based on wrapping surface reconstruction was proposed for large-scale particles in distributed environments so as to visualize the particles in high quality. In the algorithm, particle clusters were represented and then rendered in the form of a series of continuous surfaces, where the distribution of the physical variable was also shown. The algorithm was parallelized in distributed environments, thus more than a hundred million particles can be visualized using a lot of processing cores. In terms of algorithm implementation, the issue of inter-block cracks during parallel computation was be solved, and the method for rapidly finding adjacent particles was presented. Meanwhile, based on visibility culling, the particle data was filtered and thus the rendering efficiency was improved. As a result, smooth surfaces with lighting can be used to expressively exhibit inner structures and physical variable distributions of particle clusters for large-scale particles. Experiment results demonstrate that using the proposed algorithm, the rendering of more than 100 million particles is realized in 5 seconds on 512 processing cores with about 60% parallel efficiency. The proposed algorithm has been successfully applied to practical simulation applications such as massively parallel non-equilibrium molecular dynamics simulations.
粒子模拟出现在材料、工程力学、电磁环境、激光聚变、裂变能源等多个应用领域中，是科学计算的重要形式之一。其中一类粒子模拟方法是对连续物质空间抽取样本，利用一系列样本粒子来模拟物质的运动变化规律，然后通过统计或积分分析物质的宏观性质。这类模拟方法的优点，即非常适合模拟物质的畸变、断裂、破碎等大变形问题，常用方法包括分子动力学、光滑粒子流体动力学、物质点法等。在这类模拟方法中，人们不断增加粒子的抽样密度以实现更高的计算精度，由此输出急剧膨胀的粒子数据集，大规模粒子数据给可视化系统带来了极大的挑战。
在目前的可视化系统中，常常用球体表示一个个粒子，并以离散化球面的方式送入OpenGL绘制。根据离散化精度的不同，一个球面可能对应生成10~1 000个OpenGL绘制图元。对于数千万甚至数亿粒子，绘制中OpenGL需要处理的图元量可达到数亿至数百亿以上，这对目前的并行绘制系统来说是难以承受的。
在面对大规模粒子时，离散球面绘制形式不仅交互速度慢，而且大量小球极度混杂、相互遮挡，难以高质量展现模拟结果中的三维结构，不能满足材料破坏、力学毁伤等研究中的可视化需求。此外，粒子模拟的原连续物质虽然可能畸变、断裂甚至破碎，但在空间中仍然呈现宏观上的复杂形状与结构特征，也需要不同于全离散化小球的可视化表现形式。因此，目前亟须发展针对大规模粒子的高表现可视化算法，以清晰展现粒子模拟中的三维结构特征及其演化发展规律。
一个主要的对策，就是对粒子团簇的外表面进行曲面重构，通过绘制外表面的方法来表达粒子团簇的空间结构。这类方法以基于距离场的团簇表面重构方法为代表[1-3]，目前已可生成光顺的团簇表面几何，可视化表现力较高，然而其算法实现或者是串行的，或者是在一个CPU的多核或者单GPU上基于共享内存并行的，限于内存容量，难以处理大规模粒子数据。
单机/单结点或者单GPU的处理能力是有限的，难以适应规模快速增长的粒子数据可视化需求，此外，在部分国产高性能计算机上也没有配置GPU，因而，面向超大规模粒子数据，发展可扩展的分布式并行粒子绘制算法是很有必要的。为此，本文提出了基于包络面重构的大规模粒子并行绘制算法。该算法以一系列连续曲面的形式表示大规模粒子数据中的一个个团簇，并在团簇表面上绘制出粒子数据的物理量值分布。本文对算法进行了分布式并行化，从而可以通过大规模并行来处理上亿甚至更大规模的粒子数据。同时，本文从多个方面对算法进行了优化，以提升大规模粒子的绘制处理效率。
1 相关工作
最常用的粒子绘制方法，就是将粒子看作三维点球，将其表面离散成三角形网格，绘制出真实感光照效果。该方法适合于绘制少量的粒子数据，有时也用于粗略绘制出大规模粒子的空间分布形态。针对大规模粒子绘制效率低的问题，Wang等[4]提出了一个基于不可见剔除的筛选绘制方法，可剔除大量的被遮挡粒子，显著提升绘制效率。另外，研究者也在不断提升粒子的绘制效果，例如可基于OSPRay的光线跟踪方法[5]，绘制出具有真实感光影效果的粒子场景。然而，绘制大规模粒子时，视野中粒子极多，每个粒子占有屏幕区域极小，不仅难以呈现高级光影效果，其绘制也是极其费时的。
用离散球面表示一个个粒子，在处理大规模粒子时，其绘制效果难以自然形成团簇整体上的光滑表面效果。因而，引入体绘制方法来融合绘制粒子数据，形成比较连续的空间形态。抛雪球法将每个粒子看成具有深度、颜色但无几何连接的点元，投影到屏幕空间进行融合与插补[6]。Wang等[7]改进了该方法，在反向过程中优化梯度设计，正向过程中规则化点元分布，取得了很好的点云表面绘制效果。Fraedrich等[8]则利用插值核函数从粒子集合重构数据场，进而实施光线投射体绘制，可有效展现粒子集合中的三维体特征。在此基础上，Hochstetter等[9]进一步通过双向自适应采样，提升了光线投射体绘制的速度，以快速绘制动态变化的粒子集合。然而，体绘制生成的团簇边界常常显得比较模糊，而且与视点相关，不利于对团簇结构做进一步分析。
近年来，基于团簇外表面重构的粒子可视化方法得到了更多的关注。这类方法可生成连续的团簇外表面几何，也可导入高级绘制工具，由此来高表现展示粒子团簇的内部结构和空间分布。在这类方法中，生成表面的光顺效果、对复杂拓扑结构的适应性以及算法处理的速度是主要考虑的问题。
第一种方法是液滴表面生成方法[10]，该方法将粒子看成一个个的液滴球体，当粒子靠近时液滴自然融合，由此逐步形成越来越复杂的液体形态。该方法对少量粒子效果很好，但在处理大规模粒子表示的平整曲面时，例如一块平面、锥面或球面等，由于每个粒子附近都有明显的小球体隆起感，因而生成曲面光顺性差，视觉上负面效应明显。虽然Müller等[11]通过密度平均方法降低了这种效应，但无法从根本上消除平整表面上的隆起感问题。
第二种方法是对粒子数据重构标量场，通过提取隐式曲面的方法来获取粒子团簇包络面。这里，标量场重构计算常常是影响表面质量和计算时间的关键因素。Zhu和Bridson[1]将带符号距离场从单粒子球推广到了大量粒子的团簇，通过从带符号距离场中提取等值面的方法来构造、生成粒子团簇的外表面，但该方法用于凹的粒子团簇则容易出错。Solenthaler等[2]通过在距离场中加入校正因子来处理这个错误，可较好地适应粒子团簇的复杂结构，生成质量较高的团簇外表面。
针对标量场重构，Akinci等[3]基于OpenMP对重构算法进行了并行化，取得了较高的曲面重构速度。基于三层自适应格式，Akinci等[12]又对算法进行了优化，在高曲率区域用高分辨率网格保持精度，在低曲率区域则用低分辨率网格减少计算量。基于GPU加速，Huang等[13]和Du等[14]对曲面重构算法做了进一步的扩展与应用。Wu等[15]则基于两层网格优化了曲面重构算法，粗层使用哈希方法索引曲面粒子，细层实施标量场重构，由此加快了重构速度。近年来，Yang和Gao[16]对曲面重构算法的全过程进行了并行化与优化，重构速度提升了一个量级。然而，上述实现都是在一个CPU的多核或者单GPU上基于共享内存完成的，由于资源有限，难以处理大规模粒子数据。
在标量场重构中，Yu和Turk[17]还引入了各向异性核函数，其方向基于主成分分析计算，但该方法需要邻域粒子搜索，因而非常耗时，而且仍会出现凹区域出错问题。Wang等[18]进一步将基于各向异性核的曲面重构引入粒子多相流模拟，以抽取由水平集定义的流体表面和多相界面，得到了高质量曲面重构效果。在各向异性核曲面重构中，Chen等[19]提出了检测外围曲面顶点的方法，可只在曲面外围区域构造标量场，由此节省了一半的计算时间。
进一步工作还集中于边界粒子优化处理方面，用于提高曲面重构的精度和效率，包括基于空间哈希网格的重构策略[20]、边界粒子自适应重采样技术[21]、粒子预处理缺陷修复技术[22]、基于区间分析的边界粒子检测方法[23]、流体界面的层次边界检测方法[24]、改进的基于邻域粒子分布各向异性粒子变换方法[25]、基于离散指示函数的流体粒子表面重建方法[26]等。
2 算法流程
针对大规模粒子高表现可视化需求，本文提出了基于包络面重构的大规模粒子并行绘制算法。该算法对三维空间离散点集生成拓扑二维的连续曲面，对多个点集团簇则对应生成一系列连续曲面。首先，大规模粒子数据是分块存储在并行文件系统上的，本文基于可见性剔除只读入可见区域中的各个粒子数据块（patch）。为了避免计算区域交叠造成距离场计算的不一致性，本文将计算区域重新划分给各个处理器核，并重新分发粒子数据，各处理器核收集落于增广区域之内的所有粒子。进而，各处理器核对所负责的计算区域构造均匀网格分块，并根据局域粒子位置关系计算一个带符号距离场。计算中，为了快速查询局域粒子，本文让粒子与网格关联并有序组织。为了节省计算开销，将靠近目标曲面的粒子判定为曲面粒子，曲面粒子附近的网格标记为计算有效网格，后续只对有效网格计算距离场并抽取等值面。然后，利用移动立方体（marching cubes，MC）方法[27]从距离场网格中抽取值为0的隐式曲面，即为所得曲面。最后，基于就近法则，将粒子上的物理量赋值给所得曲面，即可基于云图方法完成后续的连续曲面光照绘制。该算法的流程如图1所示。下文将具体介绍该算法中的主要关键技术。
图1基于包络面重构的大规模粒子并行绘制算法流程
Fig.1Flow chart of the parallel rendering algorithm based on wrapping surface for large-scale particles in distributed environments
3 粒子分布式载入
3.1 粒子可见性剔除
本文算法的输入是大规模粒子模拟输出的结果数据，这些数据常常是分块存储于并行文件系统中，分块数目可以达到上万个甚至更多。本文首先基于可见性对粒子数据进行筛选剔除，利用OpenGL的视景锥对数据场进行粗筛，将落在视景锥之外因而对用户不可见的数据块剔除（见图2）。OpenGL使用视景投影变换将三维空间中的几何元素映射到屏幕像素区域，只有位于视景锥之内的几何元素才可能最终出现在屏幕上。那些整个位于视景锥之外的数据块，包括由此重构出来的包络面特征数据，对观察者来说是不可见的，因而，在读数据时无须读取它们，在绘制时更不需要处理它们。
图2基于可见性剔除粒子数据
Fig.2Culling particle data based on visibility
具体地，本文基于合约机制[28]实现大规模粒子数据的可见性剔除。首先，从元数据信息获取数据块序列的空间包围盒，构建一棵空间区间树，并根据当前绘制窗口的视点参数构建视景锥。其次，从根节点开始，递归处理空间区间树：当某节点的包围盒全部在视景锥之外因而对视景锥不可见时，则该节点之下的所有叶子（数据块）均不可见，停止遍历其子节点；否则，该节点的包围盒与视景锥有交，继续对该节点的子节点进行可见性判断，如此递归下去，直至到达叶子。最终得到所有与视景锥有交的数据块列表，将它们的读入需求写入合约，即可在可视化流程的下行执行阶段实际只读入这些数据块。
3.2 粒子重新分发
在分布式并行环境下，本文基于“数据并行”原则，将经过可见性剔除之后的数据块均衡分配给各个处理器核，一般采用连续分配策略。如果数据块数目为M，处理器核数目为N，则连续分配策略将把数据块列表划分为N段，依次把数目为M/N的数据块分配给处理器核P0，P1，···，PN-1。然后，各个处理器核按照所分配的数据块号，实际读入各个数据块，并送入后续处理流程。
在后续处理中，各个处理器核对读入的数据块构造带符号距离场φ（r）并分块抽取等值面φ（r）=0。为了保证所得等值面不重不漏，分布式计算中每个隐函数数据场（带符号距离场）分块应是空间中不重叠，为了算法高效，本文要求各个分块的空间包围盒不重叠。然而，从文件系统读入的粒子初始分块不能保证这一点，不同分块的粒子子集可能纠缠交错，致使其空间包围盒重叠，见图3。因而，本文需要对粒子数据进行重新划分，并分发给每个处理器核。
图3初始数据中两个分块（不同颜色标出）
Fig.3Two initial data patches (in different colors)
具体地，首先并行收集全体粒子的实际空间包围盒，将它递归划分为多个长方体子区域，指定为各核的计算区域。然后，将该子区域在各个维度的两侧各扩大R（粒子影响半径）的距离，并对各自增广区域收集落于其内的所有粒子，见图4。这样，后续计算各个子区域带符号距离场时，其局域的所有粒子皆在本地，由此保证区域边界附近的计算结果在所有处理器核上是一致的。这样，分布式计算的数据场以及后面的等值面都能保证分块附近的一致性。
图4子区域粒子收集二维示意图
Fig.4Sketch map for collecting particles in an extended region
粒子数据重新分发的优点在于，即使初始数据块的数目少于处理器核数，或者甚至只有一个数据块，都可以重新划分子区域，并均衡地分配给各个处理器核，同时，由于对子区域进行了扩展，使得各个子区域的距离场重构与等值面抽取都可以独立进行，无须通信中间结果，也没有块间裂缝问题。
4 粒子有效管理
4.1 粒子有序组织
在距离场计算中，对距离场每个节点都有寻找附近粒子的需要。若每次都遍历所有粒子，判断并找寻距离小于阈值的粒子，则这样的计算复杂度就太高了。因而，本文先对各个粒子，计算其所在单元的编号；然后依据单元编号，对所有粒子进行排序（使用快排算法，复杂度为O（nlgn））；接着，对每个单元，本文记录其内粒子的最小序号k=kmin，若其粒子的最大序号为kmax，则下一个单元记录的最小序号是kmax+1，即使它是一个空单元（不包含粒子）。由此，对任何一个单元i，立即知道其内的粒子序号为{k|ki≤k＜ki+1}，其中ki和ki+1分别是单元i和下一个单元i+1中记录的粒子最小序号。
粒子有序组织之后，在具体计算一个网格节点的带符号距离场标量值时，本文会根据影响半径R确定附近的相关单元（对均匀网格而言，复杂度为O（1）），然后根据单元访问其内的所有粒子，由此可以找到所有的局域粒子，避免了每次都对全部粒子遍历找距离满足阈值者。
4.2 标记曲面粒子
为计算区域构造一个光滑颜色场[11]，即
(1)
其中，r为位置向量，rj为位置r影响区域内的所有粒子，W（·，h）为影响半径为h的一个核函数，mj和ρj分别是粒子的质量和密度。该颜色场在粒子团簇边界处变化剧烈，因而其梯度较大，也就是说，当一个粒子处的梯度大于某个阈值，即
(2)
则它一定位于边界附近，称为曲面粒子。此处，本文取W（s，h）=g（s/h），其中
(3)
同时，令所有粒子都有相同的质量和密度，所以不妨设mj=1，ρj=1。于是，梯度计算公式为
(4)
根据这些公式，将曲面粒子标记出来，然后将曲面粒子附近的网格判定为计算有效网格（见图5），也就是说，后续只对曲面粒子附近的网格节点计算上述距离场，也只在曲面粒子附近进行等值面抽取计算，由此减少计算开销。下文中，梯度阈值l使用经验值0.125，用户也可根据粒子的疏密程度进行交互调节。
图5只重构曲面粒子附近的距离场
Fig.5Reconstruct the distance field around surface particles only
5 粒子团簇包络面重构
5.1 带符号距离场计算
引入带符号距离场标量场[1-2]，它将单个粒子球面的距离场推广到大量粒子的集合，由此即可用值为0的隐式曲面描述粒子团簇的包络面。带符号距离场标量场的重构质量是决定团簇表面生成质量的关键。具体地，首先在计算区域布置一套均匀网格，然后在每一个节点位置计算带符号距离场值。带符号距离场在一个网格节点r位置的值按以下公式计算：
(5)
其中，r为平均粒子半径或者光滑粒子的平均距离，f是一个校正因子，则
(6)
其中，，thigh=3.5，tlow=0.4，EVmax为的最大特征值，是局域粒子的加权平均，即
(7)
其中，R为带符号距离场计算中的粒子影响半径，g（s）为加权平均的核函数。即
(8)
为一个3×3的矩阵，它的一个元素为
其中，； r=（r1， r2，r3）T； rj=（r1j，r2j，r3j）T，p（s）=-6×（1-s2）2/R2；。
5.2 粒子量赋值曲面
当带符号距离场计算完毕后，即可用MC方法抽取值为0的等值面，也即对大规模输入粒子各个团簇外表面构造连续曲面。MC方法为经典的等值面抽取方法，在此无须赘述。
以上计算只是得到了粒子团簇包络面的几何信息。显然，原来粒子一般是带有物理量的，表征了粒子所在位置的物理属性，这一属性有时也需要表现在绘制曲面上。为此，本文对所得曲面的每一个顶点赋上物理值。具体地，本文以近似广度优化顺序扫描顶点周围的网格单元，找到其中的粒子，并比较距离最小者，由此快速得到局域最近粒子，将该粒子的物理量赋给该曲面顶点。注意，在之前粒子数据核间分发时，这些物理量也需一并包含。
至此，针对大规模粒子数据，既构造了几何包络面，又赋上了物理量，因而可以基于云图方式绘制该粒子数据，并通过光滑曲面光照效果来提升其表现力。
6 实验结果
将本文算法应用到大规模并行非平衡分子动力学模拟，该模拟研究分析衰减冲击载荷下金属熔化破碎机理。模拟输出的数据包含约600万粒子，在浪潮集群一个结点24个处理器核上完成绘制，如图6所示。由于物理空间是一个非常细长的区域，此处只给出局部放大图，以便更好观察绘制效果。本文对两种方法进行了比较，一是传统的粒子球绘制方法，一是本文基于包络面绘制方法。从图中可以看出，基于包络面绘制方法可以清晰展现熔化金属破碎后形成的各种团簇外表面形态，而且，通过透明化绘制（见图7）可以更好地观察熔化金属由连续体断裂破碎至液态颗粒时三维空间复杂分布。相比较而言，粒子球绘制方法则因大量小球极度簇拥，团簇表面的方向杂乱无章，因而绘制效果显得暗淡无光，无法自然呈现清晰光洁的团簇表面，又因不支持透明绘制而无法看清被遮挡的内部结构。
图6粒子球绘制与基于包络面绘制效果比较
Fig.6Comparison of the rendering results based on little spheres and wrapping surfaces
图7基于包络面绘制下以透明形式展示的清晰三维结构
Fig.7The clear three-dimensional structures displayed in a transparent form based on wrapping surfaces rendering
为了测试大规模粒子绘制性能，本文将该算法用于大体系动态破坏过程分子动力学实际模拟结果的可视化，该模拟输出结果中粒子数目达到1.2亿，包含103个时间步，时变数据总量达到736 GB。并行性能测试分别使用了32、64、128、256和512个处理器核，测试结果见图8、图9和图10。本文测试了粒子包络面并行绘制算法的时间开销以及采取可见性剔除方法后的时间开销，图8展示了并行实现中可见性剔除带来的性能提升效果，从图中可以看出，在使用不同数目处理器核的情况下，可见性剔除均可以提升9倍以上的绘制速度。在512核上，可见性剔除带来的性能提升，将粒子绘制时间降低到了5 s以内。同时，本文也测试了算法的并行可扩展性，如图9和图10所示。当处理器核数不断增加时，绘制时间随之快速下降，如图9所示，本文相应计算了算法的并行效率，如图10所示，从图中看出，本文算法的并行效率在60%以上，这显示了算法良好的并行可扩展性。
图8并行实现中可见性剔除带来的性能提升
Fig.8Rendering performance enhanced based on visibility culling in the parallel testing
图9并行绘制算法的时间开销
Fig.9Time costs of the proposed distributed parallel rendering algorithm
图10并行绘制算法的并行效率
Fig.10Parallel efficiencies of the proposed distributed parallel rendering algorithm
尽管如此，本文算法仍然存在一定的不足。从图10可以看出，当处理器核数不断增加的时候，并行效率虽然在60%之上，但其数值还是在不断地下降着。并行效率下降的原因，除了读数据文件、消息传递接口标准（message passing interface，MPI）聚合通信等的扩展性较差，还在于增广区域带来的计算量增加。当处理器核数增加时，各核处理的子区域持续变小，但仍然要在各个维度两侧增加R的扩展区域，这使得总计算量实际在不断增加，影响了算法的可扩展性。
扩展区域的负面影响在图8所示的可见性剔除效果中也有所体现。图8显示，可见性剔除带来的性能优化效果随核数增加而下降。这是由于使用了可见性剔除后，各核上处理子区域变小了，但由于始终带着宽度R的扩展区域，使得可见性剔除前后计算量的比值打了折扣。随着处理器核的增加，这个折扣后的比值越来越小。例如，理想情况下，可见性剔除前后一个核上子区域分别设为D×D×D和d×d×d，D＞d，实际计算量为（D+2R）3和（d+2R）3，若核数增加8倍，对应各维度子区域减半，则实际计算量分别变为（D/2+2R）3和（d/2+2R）3，显然
(9)
所以
(10)
若核数继续增加，则该比值持续变小。计算量比值的减小，一定程度上造成了加速比的减小。
7 结论与展望
针对大规模粒子高表现可视化需求，本文提出了基于包络面重构的大规模粒子并行绘制算法，该算法从带符号距离场中提取，绘制粒子团簇的外表面。在分布式并行实现时，该算法首先基于可见性对粒子数据进行剔除，只处理可见区域内的重构计算，由此节省了大量的内存和时间开销；其次，对增广区域完成粒子收集与距离场计算，避免包络面并行生成时的块间裂缝问题。进一步，各核在计算带符号距离场时，令粒子与网格关联并有序组织，以快速查询邻域粒子。在抽取曲面几何后，基于就近法则，将粒子上的物理量赋值给所得曲面。由此，即可基于曲面光照效果来高表现展示粒子团簇的空间结构及其物理量值分布。
该算法已成功应用到分子动力学等大规模实际模拟，对上亿粒子进行高表现可视化，可有效表现清晰光洁的粒子团簇表面，并通过透明化显示出内部/背后被遮挡的结构，绘制效果显著好于传统的粒子球绘制，有力支撑了相关应用模拟结果的后处理可视分析。通过高效并行实现策略，该算法取得了60%的并行效率，在512核上实现了上亿粒子5 s以内完成绘制。
本文算法仍存在一定的不足，未来可能在以下几个方面开展工作，以进一步提高算法的执行速度与并行可扩展性：
1）动态负载平衡策略。针对团簇特征的空间任意性，探讨曲面粒子和有效网格核间分布的不均衡问题，设计任务再划分与重分配策略，平衡计算和通信开销，提升不均衡分布特征下的并行绘制效率。
2）稀疏区域剔除。并行分配任务将粒子数据包围盒划分为子区域时，剔除掉空白或稀疏粒子区域，对于其中的孤立粒子，直接生成包络球面，由此减小计算开销。
3）混合并行实现。本文算法目前是进程并行，可考虑结点间进程、结点内线程的混合并行实现机制，由此降低子区域的划分个数，从而减小MPI聚合通信的开销以及子区域扩展部分的计算量。
4）粒子数据下采样。针对超大规模粒子数据，可考虑进行下采样处理，基于采样后的代表性粒子来重构距离场，由此降低算法的处理开销。
