基于CUDA的等值线云图并行生成算法

2022-12-06杜小甫刘辉林刘鹤丹

小型微型计算机系统 2022年12期

关键词：等值线数组等值

杜小甫，刘辉林，刘鹤丹

1(东北大学计算机科学与工程学院，沈阳 110819)

2(厦门大学嘉庚学院信息科学与技术学院，福建漳州 363105)

1 引言

1.1 等值线和等值线云图

等值线，是指标量场中具有相同标量值的各点的连线.等值线作为一种标量场可视化技术，应用非常广泛.根据标量场物理含义不同，常见的等值线包括等温线，等压线，等高线，等势线等.

将相邻的两条等值线之间的区域，根据一定的规则填充，则得到了等值线云图，这也是一种常见的标量场可视化技术.

等值线生成算法可以分为两大类：网格序列法和等值线追踪法.所谓网格序列法是指将整个物理场划分为若干网格单元，然后对每个网格单元求内部的等值线段，当所有的网格单元都求完后，各个独立的等值线段必然连接成完整的等值线.蒋瑜[1]等人提出了一种快速构建Delaunay三角网格算法，然后利用内插法求出每个网格单元内的等值线段，最后采用三次方Bezier曲线平滑生成等值线.这是一种典型的不规则网格内网格序列法算法，具有算法简单、速度快、对物理场表达能力强等优点.等值线追踪法是指从物理场中指定位置出发，沿着一条等值线不断追踪，直到等值线超出物理场边界或闭合.等值线追踪法一般更适用于规则网格，具有天然连续，利于曲线化等优点.Guo[2]等人提出利于反距离加权插值法对离散的气温数据插值得到规则矩形网格，再计算等值点、追踪等值线并对等值线进行平滑处理，这是一种典型的等值线追踪算法.冯建设[3]等人指出，由等值线填充产生云图最关键的环节是确定区域属性值和排序，不得不反复进行网格顶点与区域、区域与区域之间的包含关系判别，严重影响运行效率.他们提出基于有序标记，可以显著减少云图填充的时间复杂度.

无论哪种方法，等值线生成过程中都要处理很多网格单元和很多等值点连接成的等值线段，是一种典型的计算密集型算法.并且在等值线计算中，很多计算不存在先后依赖关系，因此非常适合做并行化处理.王志斌[4]等人讨论了一种CPU并行算法，使用插值算法计算每个像素的温度值，然后直接点绘制，最终得到等值线云图.徐涛[5]等人讨论了一种CPU并行算法，本质上是基于规则网格的等值线追踪算法，加速比2倍多.黄熙祥[6]提出了一种基于网格边标记的等值线并行生成算法，本质上也是一种基于规则网格的CPU并行等值线追踪算法，加速比最大达到3.852.

近年来，随着GPU算力不断提高，相关开发环境逐渐成熟，GPU并行计算的应用越来越广.在等值线和云图算法领域，GPU也被大量应用.陈家辉[7]提出了一种GPU并行算法，进行等值线追踪以及等值线段的合并，最后利用CUDA编程实现了上述算法.Dong[8]等人设计了一种基于CUDA的等值线跟踪算法，将等价点存储在纹理存储器中，并采用索引数组而不是链表结构.钱宸[9]等人提出一种“区间块搜索的等值线提取算法”，本质上是基于规则网格的等值线追踪算法.Cao[10]等人提出一种采用块搜索的等值线跟踪算法，减少了网格遍历次数，避免了对被排除单元的搜索.Dong[11]等人提出了一种不依赖于坐标，直接根据标量值进行等值线跟踪的算法.Song[12]等人针对USGS地图中部分等高线会被标记文字遮挡的问题，提出了一种GPU并行算法，可以去除遮挡，实现等高线重连.

1.2 CUDA编程架构

CUDA(Compute Unified Device Architecture)，是NVIDIA公司推出的通用GPU并行计算架构，非常适于进行并行程序开发，广泛应用在图像处理、深度学习、区块链等计算密集型应用场景中.目前CUDA最新版本为11，支持C/C++，FORTRAN，python，matlab语言.

NVIDIA显卡中包含很多内核，每个内核是一个独立的计算单元.为了更好的利用这些内核，将若干个内核组成一个流多处理器SM，然后由若干个SM组成完整的显卡.与这种三层硬件结构对应，CUDA框架中也将线程组织为3层.首先整个显卡对应着一个网格Grid，然后将网格划分为若干块Block，最后在每个Block中包含若干线程Thread.

一个CUDA程序包含运行在CPU端的host程序和运行在GPU端的device程序，每个device程序都被封装为一个核函数.每个核函数会被分配到所有的GPU内核上并行运行，提高程序速度.如果一个程序比较复杂，还可以将其按照先后流程划分为多个核函数，分别进行并行计算.

Liu[13]等人对航空图像处理中的梯度聚类算法进行了基于深度学习的并行研究，利用CUDA对航空图像数据进行大规模并行处理.Domkiv[14]等人提出了利用CUDA结合区块链技术来恢复加密PDF文件的密钥的并行算法.刘端阳[15]等人提出对KNN算法中距离计算、距离排序和决定分类标号三阶段采用CUDA并行化.曾浩洋[16]等人提出了一种基于天空区域分割的图像去雾算法，利用CUDA编程实现快速去雾.

综上所述，目前所有等值线GPU并行算法都是基于规则网格等值线追踪法的.对于不规则网格中网格序列法的GPU并行算法没有展开研究.

网格序列法不需要考虑同一条等值线上前后等值点之间的依赖关系，天然更适合进行并行化.相对于规则网格，不规则网格可以更好的贴合实际物理场的外形，被广泛应用在各类数据分析软件中.基于以上原因，本文提出了一种应用于不规则网格中的网格序列法等值线以及等值线云图的GPU并行生成算法，并利用CUDA编程实现了上述算法.不规则网格根据外形可以分为很多种，本文中只讨论最简单的三节点三角形网格.

2 基于网格序列法的等值线云图算法

2.1 等值线段生成算法

网格序列法计算等值线云图的过程可以分为两个步骤：1)生成等值线，2)生成云图.在第一个步骤中，最终得到很多等值线段，每条等值线段由两个等值点连接组成，而每个等值点指的是某一个网格单元的某条边与某条等值线的交点.

三角形网格中，某网格单元如果与某等值线相交，则必定有两条边与该等值线相交产生两个交点.将两个交点连接，就得到一条等值线段.某网格单元可能会与多条等值线相交，而这些具有不同标量值的等值线段必定不会相交.当所有网格单元都求出等值线段，相邻的网格单元中具有相同标量值的等值线段必然连接成完整的等值线.此时等值线是一条折线，如果有必要，可以通过Bezier曲线等方法对其平滑化.

上述过程中，一个重要的计算是求网格边与等值线的交点.如果某条网格边的两个顶点标量值和某条等值线的标量值存在一大一小关系，这条边上必定存在等值点.使用双线性插值法求交点坐标，公式如下：

(1)

其中(x0，y0)、(x1，y1)是当前网格边两个顶点的坐标，v0和v1是两个顶点的标量值，(x，y)是所求交点的坐标.

按照顺序对所有网格的所有边，与所有等值线一一判断是否相交，如果相交则利用公式(1)和公式(2)求出等值交点.则这些等值点必定成对出现，并且必定按照如下3个规则排序：

1)等值点按网格单元编号排列.

2)等值点按标量场值由小到大排列.

3)等值点按所在网格边编号由小到大排列.

2.2 等值线云图生成算法

在前步骤中，一个网格单元可能会被多条等值线段分割成多个独立的多边形，我们称之为等值多边形.三角形网格中，这些多边形有可能是三角形、四边形或五边形.如图1所示，是两个等值线实例.

图1 三角形网格中等值线示例图

根据上一小节总结的有序性，我们前期工作[17]中提出了一种基于 “等值多边形”的云图生成算法.该算法根据第一条等值线与网格单元3个顶点的大小关系，将所有网格单元的处理分为“两小一大”型和“一小两大”型.

1)“两小一大”型：某网格3个顶点中有一个的场值比第一条等值线的场值大，另外两个顶点的场值比第一条等值线的场值小，如图1(a)所示.这种情况下，第1条等值线会切下一块等值四边形，后续每条等值线都会切下一个等值四边形，最后一条等值线切下一个等值四边形后，还剩下最后一个等值三角形.

2)“一小两大”型：某网格3个顶点中有一个的场值比第一条等值线的场值小，另外两个顶点的场值比第一条等值线的场值大.如图1(b)所示.这种情况下，第1条等值线会切下一块等值三角形，后续的发展情况稍微复杂一些.首先可能会有一些等值线仍然小于剩下两个顶点，这些等值线会不断的切下一些等值四边形；然后有可能从某条等值线开始，其值大于剩下两个顶点中的某一个，仍然小于另一个，这条等值线会切下一个等值五边形；如果后续还有等值线与当前网格单元相交，则转化为“两小一大”型，不断切下等值四边形，直到最后剩下一个等值三角形.

3 算法并行化

3.1 等值线生成算法并行化

前文所述等值线生成算法的串行实现非常简单，只需要写一个两层循环.外层循环对所有网格单元遍历，里层循环对某个网格单元与所有等值线相交情况遍历.循环体内首先判断当前等值线与当前网格是否相交，如果相交，则根据公式(1)和公式(2)计算得到网格边上的等值点.这个算法的并行化，有几个关键点需要重点考虑.

1)线程划分.串行算法中循环往往是计算量最密集的地方，往往也是并行化的重点.对上述双层循环进行分析，我们发现不同网格单元之间的处理并不存在前后依赖关系；而同一网格内部，与不同等值线的相交处理，也没有前后依赖关系.根据“线程划分粒度越小越好”的原则，我们利用GPU的单个线程处理每个单元与每条等值线的相交.在每个线程内部，不存在任何循环，只需要判断一个网格单元和一条等值线是否相交，并进行等值点、等值线段的计算即可.因此，线程总数等于网格单元数量乘以等值线数量.

2)GPU内存设计.GPU并行程序中，必须尽量减少不同线程之间内存访问冲突的可能性.因为当不同线程不得不同时访问相同内存时，这些线程不能并行运行，只能串行运行.这将极大影响并行速度，因此应该让不同线程尽量访问不同的内存.

本算法中，输入数据有两组，分别是网格单元信息数据和等值线数据.输出数据则只有一组，即所有等值点数据，这个“等值点数组”的设计是重点.

我们定义一个double型一维数组ContourNodes来保存所有的等值点数据，其中每相邻4个double型数据保存1个等值点的x、y、z坐标和标量场值v.一条等值线如果和网格单元相交，必定产生两个等值点，所以在最坏情况下，等值点数组ContourNodes的长度为：MeshUnitNum×ContourNum×2×4，其中MeshUnitNum是网格单元的数量，ContourNum是等值线的数量.为了不造成内存访问冲突，我们按照最坏情况设计等值点数组：double ContourNodes[MeshUnitNum×ContourNum×2×4].在GPU所有内核并行计算时，每个内核都将计算结果保存到对应位置的元素中，这样就不会造成内存访问冲突.当然，实际运行中，很多网格单元和很多等值线之间不会相交，这样ContourNodes很多元素并没有使用，这会造成一定空间浪费.但是这样保证了并行化顺利进行，以空间换取时间效率是值得的.

在核函数中，我们首先取得当前线程的线程编号.然后我们根据如下两个公式计算得到网格单元编号和等值线编号.

meshID=TID/ContourNum

(3)

contourID=TID%ContourNum

(4)

其中TID为当前线程编号.如上所述，我们对每个网格单元与每条等值线的相交处理分配一个线程，那么用线程编号对等值线总数进行整除运算，得到的就是当前网格编号；用线程编号对等值线总数进行求余运算，得到的就是当前等值线的编号.

在计算出某条等值线和当前网格单元的两个交点后，第一个交点的xyz坐标存入ContourNodes数组中从(meshID×ContourNum+contourID)×8开始的连续3个元素中，将当前等值线标量值写入(meshID×ContourNum+contourID)×8+3元素中；第2个交点的xyz坐标存入ContourNodes数组中从(meshID×ContourNum+contourID)×8+4开始的连续3个元素中，将当前等值线标量值写入(meshID×ContourNum+contourID)×8+7元素中.等值线并行生成算法的核函数1的完整流程图如图2所示.

图2 等值线并行算法流程图

3.2 等值线云图生成算法并行化

根据前面等到的等值点数据计算云图数据的串行算法也非常简单，循环遍历所有网格单元，根据3个顶点标量值与第一条相交等值线的大小关系，将网格处理分为“两小一大”和“一小两大”两种情况，对两种情况分别处理就好.而在这个算法的并行化过程中，有几点问题需要注意.

1)线程划分.等值线云图处理中，必须按照等值线段的先后顺序生成若干等值多边形，同个网格单元中不同等值线段存在先后依赖.因此，线程只能划分为网格单元.一个线程中不可避免的会出现循环，来处理多条等值线段，生成多个等值多边形.

2)输入数据压缩.上一小节计算得到的ContourNodes数组中存在大量无效数据，如果直接把它复制到GPU内存中，必然会造成很大空间浪费.针对这个问题，我们提出了一种压缩方法.

首先，在定义ContourNodes数组时，将他所有元素初始化为DBL_MAX，这个值是C++语言中最大double型数值.

然后，经过等值线算法处理后，如果这个数组内某些元素值改变了，则保存的是真正有意义的等值点数据；如果某些元素值还保留DBL_MAX值，则意味着该网格单元和该等值线没有交点.我们将ContourNodes数组中所有DBL_MAX值压缩掉，压缩后的数组命名为CNsBeenCompressed.同时定义一个长度为2倍MeshUnitNum的int型数组，命名为IndexOfEachUnit.IndexOfEachUnit数组中第(i/2)个元素保存了第i个网格单元在CNsBeenCompressed数组中的起始索引位置；第(i/2+1)个元素保存了第i个网格单元中等值点的个数.根据这两个值，我们可以从CNsBeenCompressed数组中得到每一个网格单元中的所有等值点数据，而又节省了大量空间.

3)输出数据格式.为了避免线程访问冲突，应该将不同线程的计算结果(若干个等值多边形数据)，写入不同的内存空间中.在算法执行之前，不能提前知道每个网格单元将产生多少个等值多边形，也无法知道将得到的是三角形，四边形，还是五边形.只能按照最坏情况设计double型数组NephogramPolygons.最坏情况下每条等值线与每个网格单元都相交，将产生MeshUnitNum×(ContourNum+1)个等值多边形.而多边形最坏情况是一个五边形，需要保存5个顶点信息，每个顶点则需要保存xyz坐标和标量值v.因此数组NephogramPolygons的长度是MeshUnitNum×(ContourNum+1)×5×4.

云图算法执行后，NephogramPolygons数组同样可能存在大量的空间浪费，很有可能其中很多元素并没有保存真实有意义的数据.所以我们同样可以将NephogramPolygons数组所有元素初始化为DBL_MAX值，当整个网格处理完毕，将NephogramPolygons数组内容复制到CPU内存中，再对其内容进行过滤，得到所有真实等值多边形的数据.

如图3所示是等值线云图并行算法的核函数2的整体流程图.首先判断当前网格是否包含等值点：如果不包含则调用WholeUnitToOneTriangle函数将网格单元整体转换为一个等值三角形，并存入NephogramPolygons对应位置；如果包含则调用SplitUnitToPolygons函数.SplitUnitToPolygons函数将处理流程分为“两小一大”和“一小两大”两种情况，分别调用函数TwoSmallOneBig和OneSmallTwoBig.

图3 等值线云图并行算法流程图

如图4所示，是“两小一大”型的函数流程图.该函数比较简单，处理流程中不会出现分支.

图4 “两小一大”函数流程图

如图5所示，是OneSmallTwoBig函数流程图.OneSmallTwoBig函数中后续后续的发展情况稍微复杂一些，出现了几种不同的分支处理.

图5 “两大一小”函数流程图

3.3 完整的等值线云图并行算法

前面分别讨论了两个核函数，核函数1用于并行计算等值点对组成的等值线段，核函数2用于并行计算等值多边形.之所以不将两个过程合并在一起，是因为两个核函数中线程划分的粒度不同.核函数1中每个线程只处理一个网格单元和一条等值线的相交，而在核函数2中每个线程要处理一个网格单元和所有等值线的相交.另外两个核函数中对GPU内存的使用也有所不同，所以我们不得不将完整的并行算法拆分为两个步骤.

如图6所示，是完整的等值线云图并行算法示意图.CPU端主机程序启动后，第1步进行内存准备.在CPU内存和GPU显存中创建对应的3个数组，分别用于保存网格单元数据、等值线数据和等值点数据;第2步从本地磁盘文件中读入网格数据和等值线数据并存入CPU内存中的网格单元数组和等值线数组，然后拷贝到GPU内存对应数组中，作为输入数据;第3步，启动核函数1，计算得到所有网格单元中所有成对出现的等值点数据，保存到GPU内存中的等值点数组中，并将其拷贝回CPU内存中的等值点数组作为输出.到此为止，完成了第1阶段，完成求解等值点的操作;第4步，对得到的等值点数组进行数据压缩，并拷贝到GPU内存中创建对应数组，作为第2阶段的输入;第5步启动核函数2，利用网格单元数组、等值线数组和压缩后的等值点数组作为输入，计算得到所有等值多边形数据，并拷贝回CPU内存中存入等值多边形数组.至此，完成了等值线云图的所有计算过程，再利用OpenGL等技术完成绘制.

图6 等值线云图并行算法整体流程图

4 算法实现与分析

4.1 与其他云图生成算法比较分析

以上我们基于“分割穷举”思想设计了一种等值线云图的生成算法，适用于三角形不规则网格.该算法与其他现有算法相比，在速度上有较大提升，具体有两方面原因.

1)云图串行算法时间复杂度降低

与现有的其他等值线云图算法进行比较，我们可以分成两个步骤来讨论算法效率.首先第一步生成等值线.在这个步骤中，无论是基于等值线序列法还是基于网格序列法的算法，时间效率都是O(m×n)，本质上都需要对每个网格单元与每条等值线的相交情况进行一一判断和处理.在这个步骤中，我们的算法并没有速度优势.如果考虑到后期等值线由折线顺滑为曲线的过程，我们的网格序列法还要进行等值线段的排序，会产生额外的计算量.从这个角度看，我们的算法在生成等值线段这一步骤反而更慢.当然，在本文中，我们并没有必要将等值线曲线化，所以这个额外的计算量其实并没有产生.

第2个步骤是根据等值线数据填充生成云图，这个步骤现有的绝大多数算法都是基于等值线序列法[2，3]的.而在等值线序列法中，从等值线数据出发完成云图填充，需要进行额外的计算，例如需要完成等值区域的排序、判断包含关系、凹多边形区域处理等.例如在文献[3]中需要额外的计算量的时间复杂度为O(n2)，其中n是等值线条数.

而在我们提出的云图算法中，云图中等值多边形的计算是混合在等值线段计算之后同步完成的(当然也可以先计算等值线段，保存后再计算等值多边形).从时间复杂度上看，完整的等值线云图算法的全部时间复杂度就是O(m×n)，并没有增加.当然，从等值线段计算出等值多边形，还需要花费一些额外计算，但是并没有在数量级上增加算法的时间消耗.因此，如果考虑云图整体的算法效率，我们的云图串行算法比其他现有基于等值线数据填充思想的云图算法要快.

2)云图并行算法速度显著提升

我们的云图算法的第2个优势是，非常适合于进行并行化，这也是本算法速度显著提升的根本原因.

在我们的云图串行算法中，计算过程被分割到每个网格单元内部，在不同网格单元之间不存在任何先后依赖关系.这就允许云图串行算法并行化.而实际结果证明，这种并行化后的算法，时间效率提示了一个数量级.

4.2 云图生成算法实现与结果分析

我们利用C++语言实现了上述并行算法，开发环境为Win10+VS2019+CUDA11.使用普通笔记本电脑，配置为intel i7-10510U处理器，16G内存，NVIDIA GeForce MX250显卡，2G显存.显卡CUDA版本为10，拥有3个流多处理器，每个流多处理器上有128个流处理器内核，共384个内核.

为了验证算法效果，我们选取了两个典型温度场数据进行了等值线和云图的计算以及绘制试验.我们利用OpenGL实现了图像绘制.在OpenGL中，多边形的颜色填充可以采用两种不同的模式.第1种模式可以为整个多边形指定一种固定的颜色，然后用这个单一颜色填充，称为平面着色(FlatShading).第2种模式，可以给多边形每个顶点指定不同的颜色，然后用这些不同颜色自动实现过渡填充，这种方式称为光滑着色处理(Smooth Shading).OpenGL中利用函数“void glShadeModel(GLenum mode)”来实现两种模式的设置，当参数mode取值为GL_FLAT或GL_SMOOTH时，分别代表平面着色和光滑着色.

采用平面着色模式时，填充的等值线云图具有明确的条纹分层，我们称为条纹云图.采用光滑着色模式时，填充的等值线云图具有光滑的颜色过渡，效果非常像像素级云图，我们称为伪像素云图.如图7所示，是一组涡轮盘横截面的温度场效果图.从上到下分别为等值线效果图、条纹云图效果图、和伪像素云图效果图.

图7 等值线和云图效果图1

如图8所示，是一组异形垫片的温度场效果图.从上到下分别为等值线效果图、条纹云图效果图、和伪像素云图效果图.从图7和图8中，可以看出条纹云图和伪像素云图的区别很明显.也能看出等值线以折线形式实现，并没有光滑成曲线形式.

图8 等值线和云图效果图2

为了更好的观察并行算法的加速效果，我们利用Delaunay算法生成了20组三角网格，网格单元数量分别是10，16，24，34，40，72，98，188，256，382，480，922，1410，2114，2808，3700，4408，5266，6448，7623.然后针对这20个网格分别创建了一个温度场，在其中都设置20条等值线.最后利用上述思想分别编写了CPU串行算法和GPU并行算法，分别运行得到结果.

我们首先单独运行等值线算法，分别得到了20组串行算法耗时和并行算法耗时数据，如表1所示.再根据表1数据绘制了20组数据的加速比变化曲线，如图9所示.

图9 等值线并行算法加速比变化图

表1 串行、并行等值线算法耗时

从图9中可以看到，随着网格数量的增加，并行算法的优势很快体现出来.当网格数量达到200左右时，乘以等值线数量20，总共需要约200×20=4000个线程并行处理.这时基本达到了试验用笔记本的性能极限，后续网格数量继续增加，并行算法的加速比并没有继续提高.

然后使用同样的20组数据，我们完整运行了等值线云图算法，得到的结果如表2所示.再根据表2数据绘制了云图算法的20组数据加速比变化曲线，如图10所示.

表2 串行、并行云图算法耗时

对比图10与图9，可以发现有3处明显不同：

图10 云图并行算法加速比变化图

1)当网格数量较少时，云图并行算法加速比甚至小于1.这说明在网格数量较少时，并行算法并没有体现出优势，反而因为内存分配、数据拷贝、线程调度等处理的耗时，最终需要花费比串行算法给多的时间，才能完成同样的计算.

2)云图算法加速比递增的速度明显要小等值线算法的.这主要是因为等值线算法中，所需线程总数等于网格个数乘以等值线条数.而在云图算法中，线程总数就等于网格个数.因此随着网格的增加，所需的线程数量在等值线算法中增加的更快，也就能更充分的利用GPU线程多的优势，因此在等值线算法中，加速比更快增加到最优程度.

3)同样的计算机配置，云图算法最终达到的加速比只有10左右，明显比等值线并行算法最高加速比13要低.这主要是因为云图算法中每个线程需要进行的计算量比较大，导致最终达到的加速效果要差一些.