杜小甫，劉輝林，劉鶴丹

1(東北大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院，沈陽(yáng) 110819)

2(廈門大學(xué) 嘉庚學(xué)院 信息科學(xué)與技術(shù)學(xué)院，福建 漳州 363105)

1 引 言

1.1 等值線和等值線云圖

等值線，是指標(biāo)量場(chǎng)中具有相同標(biāo)量值的各點(diǎn)的連線.等值線作為一種標(biāo)量場(chǎng)可視化技術(shù)，應(yīng)用非常廣泛.根據(jù)標(biāo)量場(chǎng)物理含義不同，常見(jiàn)的等值線包括等溫線，等壓線，等高線，等勢(shì)線等.

將相鄰的兩條等值線之間的區(qū)域，根據(jù)一定的規(guī)則填充，則得到了等值線云圖，這也是一種常見(jiàn)的標(biāo)量場(chǎng)可視化技術(shù).

等值線生成算法可以分為兩大類：網(wǎng)格序列法和等值線追蹤法.所謂網(wǎng)格序列法是指將整個(gè)物理場(chǎng)劃分為若干網(wǎng)格單元，然后對(duì)每個(gè)網(wǎng)格單元求內(nèi)部的等值線段，當(dāng)所有的網(wǎng)格單元都求完后，各個(gè)獨(dú)立的等值線段必然連接成完整的等值線.蔣瑜[1]等人提出了一種快速構(gòu)建Delaunay三角網(wǎng)格算法，然后利用內(nèi)插法求出每個(gè)網(wǎng)格單元內(nèi)的等值線段，最后采用三次方Bezier曲線平滑生成等值線.這是一種典型的不規(guī)則網(wǎng)格內(nèi)網(wǎng)格序列法算法，具有算法簡(jiǎn)單、速度快、對(duì)物理場(chǎng)表達(dá)能力強(qiáng)等優(yōu)點(diǎn).等值線追蹤法是指從物理場(chǎng)中指定位置出發(fā)，沿著一條等值線不斷追蹤，直到等值線超出物理場(chǎng)邊界或閉合.等值線追蹤法一般更適用于規(guī)則網(wǎng)格，具有天然連續(xù)，利于曲線化等優(yōu)點(diǎn).Guo[2]等人提出利于反距離加權(quán)插值法對(duì)離散的氣溫?cái)?shù)據(jù)插值得到規(guī)則矩形網(wǎng)格，再計(jì)算等值點(diǎn)、追蹤等值線并對(duì)等值線進(jìn)行平滑處理，這是一種典型的等值線追蹤算法.馮建設(shè)[3]等人指出，由等值線填充產(chǎn)生云圖最關(guān)鍵的環(huán)節(jié)是確定區(qū)域?qū)傩灾岛团判颍坏貌环磸?fù)進(jìn)行網(wǎng)格頂點(diǎn)與區(qū)域、區(qū)域與區(qū)域之間的包含關(guān)系判別，嚴(yán)重影響運(yùn)行效率.他們提出基于有序標(biāo)記，可以顯著減少云圖填充的時(shí)間復(fù)雜度.

無(wú)論哪種方法，等值線生成過(guò)程中都要處理很多網(wǎng)格單元和很多等值點(diǎn)連接成的等值線段，是一種典型的計(jì)算密集型算法.并且在等值線計(jì)算中，很多計(jì)算不存在先后依賴關(guān)系，因此非常適合做并行化處理.王志斌[4]等人討論了一種CPU并行算法，使用插值算法計(jì)算每個(gè)像素的溫度值，然后直接點(diǎn)繪制，最終得到等值線云圖.徐濤[5]等人討論了一種CPU并行算法，本質(zhì)上是基于規(guī)則網(wǎng)格的等值線追蹤算法，加速比2倍多.黃熙祥[6]提出了一種基于網(wǎng)格邊標(biāo)記的等值線并行生成算法，本質(zhì)上也是一種基于規(guī)則網(wǎng)格的CPU并行等值線追蹤算法，加速比最大達(dá)到3.852.

近年來(lái)，隨著GPU算力不斷提高，相關(guān)開(kāi)發(fā)環(huán)境逐漸成熟，GPU并行計(jì)算的應(yīng)用越來(lái)越廣.在等值線和云圖算法領(lǐng)域，GPU也被大量應(yīng)用.陳家輝[7]提出了一種GPU并行算法，進(jìn)行等值線追蹤以及等值線段的合并，最后利用CUDA編程實(shí)現(xiàn)了上述算法.Dong[8]等人設(shè)計(jì)了一種基于CUDA的等值線跟蹤算法，將等價(jià)點(diǎn)存儲(chǔ)在紋理存儲(chǔ)器中，并采用索引數(shù)組而不是鏈表結(jié)構(gòu).錢宸[9]等人提出一種“區(qū)間塊搜索的等值線提取算法”，本質(zhì)上是基于規(guī)則網(wǎng)格的等值線追蹤算法.Cao[10]等人提出一種采用塊搜索的等值線跟蹤算法，減少了網(wǎng)格遍歷次數(shù)，避免了對(duì)被排除單元的搜索.Dong[11]等人提出了一種不依賴于坐標(biāo)，直接根據(jù)標(biāo)量值進(jìn)行等值線跟蹤的算法.Song[12]等人針對(duì)USGS地圖中部分等高線會(huì)被標(biāo)記文字遮擋的問(wèn)題，提出了一種GPU并行算法，可以去除遮擋，實(shí)現(xiàn)等高線重連.

1.2 CUDA編程架構(gòu)

CUDA(Compute Unified Device Architecture)，是NVIDIA公司推出的通用GPU并行計(jì)算架構(gòu)，非常適于進(jìn)行并行程序開(kāi)發(fā)，廣泛應(yīng)用在圖像處理、深度學(xué)習(xí)、區(qū)塊鏈等計(jì)算密集型應(yīng)用場(chǎng)景中.目前CUDA最新版本為11，支持C/C++，F(xiàn)ORTRAN，python，matlab語(yǔ)言.

NVIDIA顯卡中包含很多內(nèi)核，每個(gè)內(nèi)核是一個(gè)獨(dú)立的計(jì)算單元.為了更好的利用這些內(nèi)核，將若干個(gè)內(nèi)核組成一個(gè)流多處理器SM，然后由若干個(gè)SM組成完整的顯卡.與這種三層硬件結(jié)構(gòu)對(duì)應(yīng)，CUDA框架中也將線程組織為3層.首先整個(gè)顯卡對(duì)應(yīng)著一個(gè)網(wǎng)格Grid，然后將網(wǎng)格劃分為若干塊Block，最后在每個(gè)Block中包含若干線程Thread.

一個(gè)CUDA程序包含運(yùn)行在CPU端的host程序和運(yùn)行在GPU端的device程序，每個(gè)device程序都被封裝為一個(gè)核函數(shù).每個(gè)核函數(shù)會(huì)被分配到所有的GPU內(nèi)核上并行運(yùn)行，提高程序速度.如果一個(gè)程序比較復(fù)雜，還可以將其按照先后流程劃分為多個(gè)核函數(shù)，分別進(jìn)行并行計(jì)算.

Liu[13]等人對(duì)航空?qǐng)D像處理中的梯度聚類算法進(jìn)行了基于深度學(xué)習(xí)的并行研究，利用CUDA對(duì)航空?qǐng)D像數(shù)據(jù)進(jìn)行大規(guī)模并行處理.Domkiv[14]等人提出了利用CUDA結(jié)合區(qū)塊鏈技術(shù)來(lái)恢復(fù)加密PDF文件的密鑰的并行算法.劉端陽(yáng)[15]等人提出對(duì)KNN算法中距離計(jì)算、距離排序和決定分類標(biāo)號(hào)三階段采用CUDA并行化.曾浩洋[16]等人提出了一種基于天空區(qū)域分割的圖像去霧算法，利用CUDA編程實(shí)現(xiàn)快速去霧.

綜上所述，目前所有等值線GPU并行算法都是基于規(guī)則網(wǎng)格等值線追蹤法的.對(duì)于不規(guī)則網(wǎng)格中網(wǎng)格序列法的GPU并行算法沒(méi)有展開(kāi)研究.

網(wǎng)格序列法不需要考慮同一條等值線上前后等值點(diǎn)之間的依賴關(guān)系，天然更適合進(jìn)行并行化.相對(duì)于規(guī)則網(wǎng)格，不規(guī)則網(wǎng)格可以更好的貼合實(shí)際物理場(chǎng)的外形，被廣泛應(yīng)用在各類數(shù)據(jù)分析軟件中.基于以上原因，本文提出了一種應(yīng)用于不規(guī)則網(wǎng)格中的網(wǎng)格序列法等值線以及等值線云圖的GPU并行生成算法，并利用CUDA編程實(shí)現(xiàn)了上述算法.不規(guī)則網(wǎng)格根據(jù)外形可以分為很多種，本文中只討論最簡(jiǎn)單的三節(jié)點(diǎn)三角形網(wǎng)格.

2 基于網(wǎng)格序列法的等值線云圖算法

2.1 等值線段生成算法

網(wǎng)格序列法計(jì)算等值線云圖的過(guò)程可以分為兩個(gè)步驟：1)生成等值線，2)生成云圖.在第一個(gè)步驟中，最終得到很多等值線段，每條等值線段由兩個(gè)等值點(diǎn)連接組成，而每個(gè)等值點(diǎn)指的是某一個(gè)網(wǎng)格單元的某條邊與某條等值線的交點(diǎn).

三角形網(wǎng)格中，某網(wǎng)格單元如果與某等值線相交，則必定有兩條邊與該等值線相交產(chǎn)生兩個(gè)交點(diǎn).將兩個(gè)交點(diǎn)連接，就得到一條等值線段.某網(wǎng)格單元可能會(huì)與多條等值線相交，而這些具有不同標(biāo)量值的等值線段必定不會(huì)相交.當(dāng)所有網(wǎng)格單元都求出等值線段，相鄰的網(wǎng)格單元中具有相同標(biāo)量值的等值線段必然連接成完整的等值線.此時(shí)等值線是一條折線，如果有必要，可以通過(guò)Bezier曲線等方法對(duì)其平滑化.

上述過(guò)程中，一個(gè)重要的計(jì)算是求網(wǎng)格邊與等值線的交點(diǎn).如果某條網(wǎng)格邊的兩個(gè)頂點(diǎn)標(biāo)量值和某條等值線的標(biāo)量值存在一大一小關(guān)系，這條邊上必定存在等值點(diǎn).使用雙線性插值法求交點(diǎn)坐標(biāo)，公式如下：

(1)

(1)

其中(x0，y0)、(x1，y1)是當(dāng)前網(wǎng)格邊兩個(gè)頂點(diǎn)的坐標(biāo)，v0和v1是兩個(gè)頂點(diǎn)的標(biāo)量值，(x，y)是所求交點(diǎn)的坐標(biāo).

按照順序?qū)λ芯W(wǎng)格的所有邊，與所有等值線一一判斷是否相交，如果相交則利用公式(1)和公式(2)求出等值交點(diǎn).則這些等值點(diǎn)必定成對(duì)出現(xiàn)，并且必定按照如下3個(gè)規(guī)則排序：

1)等值點(diǎn)按網(wǎng)格單元編號(hào)排列.

2)等值點(diǎn)按標(biāo)量場(chǎng)值由小到大排列.

3)等值點(diǎn)按所在網(wǎng)格邊編號(hào)由小到大排列.

2.2 等值線云圖生成算法

在前步驟中，一個(gè)網(wǎng)格單元可能會(huì)被多條等值線段分割成多個(gè)獨(dú)立的多邊形，我們稱之為等值多邊形.三角形網(wǎng)格中，這些多邊形有可能是三角形、四邊形或五邊形.如圖1所示，是兩個(gè)等值線實(shí)例.

圖1 三角形網(wǎng)格中等值線示例圖

根據(jù)上一小節(jié)總結(jié)的有序性，我們前期工作[17]中提出了一種基于 “等值多邊形”的云圖生成算法.該算法根據(jù)第一條等值線與網(wǎng)格單元3個(gè)頂點(diǎn)的大小關(guān)系，將所有網(wǎng)格單元的處理分為“兩小一大”型和“一小兩大”型.

1)“兩小一大”型：某網(wǎng)格3個(gè)頂點(diǎn)中有一個(gè)的場(chǎng)值比第一條等值線的場(chǎng)值大，另外兩個(gè)頂點(diǎn)的場(chǎng)值比第一條等值線的場(chǎng)值小，如圖1(a)所示.這種情況下，第1條等值線會(huì)切下一塊等值四邊形，后續(xù)每條等值線都會(huì)切下一個(gè)等值四邊形，最后一條等值線切下一個(gè)等值四邊形后，還剩下最后一個(gè)等值三角形.

2)“一小兩大”型：某網(wǎng)格3個(gè)頂點(diǎn)中有一個(gè)的場(chǎng)值比第一條等值線的場(chǎng)值小，另外兩個(gè)頂點(diǎn)的場(chǎng)值比第一條等值線的場(chǎng)值大.如圖1(b)所示.這種情況下，第1條等值線會(huì)切下一塊等值三角形，后續(xù)的發(fā)展情況稍微復(fù)雜一些.首先可能會(huì)有一些等值線仍然小于剩下兩個(gè)頂點(diǎn)，這些等值線會(huì)不斷的切下一些等值四邊形；然后有可能從某條等值線開(kāi)始，其值大于剩下兩個(gè)頂點(diǎn)中的某一個(gè)，仍然小于另一個(gè)，這條等值線會(huì)切下一個(gè)等值五邊形；如果后續(xù)還有等值線與當(dāng)前網(wǎng)格單元相交，則轉(zhuǎn)化為“兩小一大”型，不斷切下等值四邊形，直到最后剩下一個(gè)等值三角形.

3 算法并行化

3.1 等值線生成算法并行化

前文所述等值線生成算法的串行實(shí)現(xiàn)非常簡(jiǎn)單，只需要寫一個(gè)兩層循環(huán).外層循環(huán)對(duì)所有網(wǎng)格單元遍歷，里層循環(huán)對(duì)某個(gè)網(wǎng)格單元與所有等值線相交情況遍歷.循環(huán)體內(nèi)首先判斷當(dāng)前等值線與當(dāng)前網(wǎng)格是否相交，如果相交，則根據(jù)公式(1)和公式(2)計(jì)算得到網(wǎng)格邊上的等值點(diǎn).這個(gè)算法的并行化，有幾個(gè)關(guān)鍵點(diǎn)需要重點(diǎn)考慮.

1)線程劃分.串行算法中循環(huán)往往是計(jì)算量最密集的地方，往往也是并行化的重點(diǎn).對(duì)上述雙層循環(huán)進(jìn)行分析，我們發(fā)現(xiàn)不同網(wǎng)格單元之間的處理并不存在前后依賴關(guān)系；而同一網(wǎng)格內(nèi)部，與不同等值線的相交處理，也沒(méi)有前后依賴關(guān)系.根據(jù)“線程劃分粒度越小越好”的原則，我們利用GPU的單個(gè)線程處理每個(gè)單元與每條等值線的相交.在每個(gè)線程內(nèi)部，不存在任何循環(huán)，只需要判斷一個(gè)網(wǎng)格單元和一條等值線是否相交，并進(jìn)行等值點(diǎn)、等值線段的計(jì)算即可.因此，線程總數(shù)等于網(wǎng)格單元數(shù)量乘以等值線數(shù)量.

2)GPU內(nèi)存設(shè)計(jì).GPU并行程序中，必須盡量減少不同線程之間內(nèi)存訪問(wèn)沖突的可能性.因?yàn)楫?dāng)不同線程不得不同時(shí)訪問(wèn)相同內(nèi)存時(shí)，這些線程不能并行運(yùn)行，只能串行運(yùn)行.這將極大影響并行速度，因此應(yīng)該讓不同線程盡量訪問(wèn)不同的內(nèi)存.

本算法中，輸入數(shù)據(jù)有兩組，分別是網(wǎng)格單元信息數(shù)據(jù)和等值線數(shù)據(jù).輸出數(shù)據(jù)則只有一組，即所有等值點(diǎn)數(shù)據(jù)，這個(gè)“等值點(diǎn)數(shù)組”的設(shè)計(jì)是重點(diǎn).

我們定義一個(gè)double型一維數(shù)組ContourNodes來(lái)保存所有的等值點(diǎn)數(shù)據(jù)，其中每相鄰4個(gè)double型數(shù)據(jù)保存1個(gè)等值點(diǎn)的x、y、z坐標(biāo)和標(biāo)量場(chǎng)值v.一條等值線如果和網(wǎng)格單元相交，必定產(chǎn)生兩個(gè)等值點(diǎn)，所以在最壞情況下，等值點(diǎn)數(shù)組ContourNodes的長(zhǎng)度為：MeshUnitNum×ContourNum×2×4，其中MeshUnitNum是網(wǎng)格單元的數(shù)量，ContourNum是等值線的數(shù)量.為了不造成內(nèi)存訪問(wèn)沖突，我們按照最壞情況設(shè)計(jì)等值點(diǎn)數(shù)組：double ContourNodes[MeshUnitNum×ContourNum×2×4].在GPU所有內(nèi)核并行計(jì)算時(shí)，每個(gè)內(nèi)核都將計(jì)算結(jié)果保存到對(duì)應(yīng)位置的元素中，這樣就不會(huì)造成內(nèi)存訪問(wèn)沖突.當(dāng)然，實(shí)際運(yùn)行中，很多網(wǎng)格單元和很多等值線之間不會(huì)相交，這樣ContourNodes很多元素并沒(méi)有使用，這會(huì)造成一定空間浪費(fèi).但是這樣保證了并行化順利進(jìn)行，以空間換取時(shí)間效率是值得的.

在核函數(shù)中，我們首先取得當(dāng)前線程的線程編號(hào).然后我們根據(jù)如下兩個(gè)公式計(jì)算得到網(wǎng)格單元編號(hào)和等值線編號(hào).

meshID=TID/ContourNum

(3)

contourID=TID%ContourNum

(4)

其中TID為當(dāng)前線程編號(hào).如上所述，我們對(duì)每個(gè)網(wǎng)格單元與每條等值線的相交處理分配一個(gè)線程，那么用線程編號(hào)對(duì)等值線總數(shù)進(jìn)行整除運(yùn)算，得到的就是當(dāng)前網(wǎng)格編號(hào)；用線程編號(hào)對(duì)等值線總數(shù)進(jìn)行求余運(yùn)算，得到的就是當(dāng)前等值線的編號(hào).

在計(jì)算出某條等值線和當(dāng)前網(wǎng)格單元的兩個(gè)交點(diǎn)后，第一個(gè)交點(diǎn)的xyz坐標(biāo)存入ContourNodes數(shù)組中從(meshID×ContourNum+contourID)×8開(kāi)始的連續(xù)3個(gè)元素中，將當(dāng)前等值線標(biāo)量值寫入(meshID×ContourNum+contourID)×8+3元素中；第2個(gè)交點(diǎn)的xyz坐標(biāo)存入ContourNodes數(shù)組中從(meshID×ContourNum+contourID)×8+4開(kāi)始的連續(xù)3個(gè)元素中，將當(dāng)前等值線標(biāo)量值寫入(meshID×ContourNum+contourID)×8+7元素中.等值線并行生成算法的核函數(shù)1的完整流程圖如圖2所示.

圖2 等值線并行算法流程圖

3.2 等值線云圖生成算法并行化

根據(jù)前面等到的等值點(diǎn)數(shù)據(jù)計(jì)算云圖數(shù)據(jù)的串行算法也非常簡(jiǎn)單，循環(huán)遍歷所有網(wǎng)格單元，根據(jù)3個(gè)頂點(diǎn)標(biāo)量值與第一條相交等值線的大小關(guān)系，將網(wǎng)格處理分為“兩小一大”和“一小兩大”兩種情況，對(duì)兩種情況分別處理就好.而在這個(gè)算法的并行化過(guò)程中，有幾點(diǎn)問(wèn)題需要注意.

1)線程劃分.等值線云圖處理中，必須按照等值線段的先后順序生成若干等值多邊形，同個(gè)網(wǎng)格單元中不同等值線段存在先后依賴.因此，線程只能劃分為網(wǎng)格單元.一個(gè)線程中不可避免的會(huì)出現(xiàn)循環(huán)，來(lái)處理多條等值線段，生成多個(gè)等值多邊形.

2)輸入數(shù)據(jù)壓縮.上一小節(jié)計(jì)算得到的ContourNodes數(shù)組中存在大量無(wú)效數(shù)據(jù)，如果直接把它復(fù)制到GPU內(nèi)存中，必然會(huì)造成很大空間浪費(fèi).針對(duì)這個(gè)問(wèn)題，我們提出了一種壓縮方法.

首先，在定義ContourNodes數(shù)組時(shí)，將他所有元素初始化為DBL_MAX，這個(gè)值是C++語(yǔ)言中最大double型數(shù)值.

然后，經(jīng)過(guò)等值線算法處理后，如果這個(gè)數(shù)組內(nèi)某些元素值改變了，則保存的是真正有意義的等值點(diǎn)數(shù)據(jù)；如果某些元素值還保留DBL_MAX值，則意味著該網(wǎng)格單元和該等值線沒(méi)有交點(diǎn).我們將ContourNodes數(shù)組中所有DBL_MAX值壓縮掉，壓縮后的數(shù)組命名為CNsBeenCompressed.同時(shí)定義一個(gè)長(zhǎng)度為2倍MeshUnitNum的int型數(shù)組，命名為IndexOfEachUnit.IndexOfEachUnit數(shù)組中第(i/2)個(gè)元素保存了第i個(gè)網(wǎng)格單元在CNsBeenCompressed數(shù)組中的起始索引位置；第(i/2+1)個(gè)元素保存了第i個(gè)網(wǎng)格單元中等值點(diǎn)的個(gè)數(shù).根據(jù)這兩個(gè)值，我們可以從CNsBeenCompressed數(shù)組中得到每一個(gè)網(wǎng)格單元中的所有等值點(diǎn)數(shù)據(jù)，而又節(jié)省了大量空間.

3)輸出數(shù)據(jù)格式.為了避免線程訪問(wèn)沖突，應(yīng)該將不同線程的計(jì)算結(jié)果(若干個(gè)等值多邊形數(shù)據(jù))，寫入不同的內(nèi)存空間中.在算法執(zhí)行之前，不能提前知道每個(gè)網(wǎng)格單元將產(chǎn)生多少個(gè)等值多邊形，也無(wú)法知道將得到的是三角形，四邊形，還是五邊形.只能按照最壞情況設(shè)計(jì)double型數(shù)組NephogramPolygons.最壞情況下每條等值線與每個(gè)網(wǎng)格單元都相交，將產(chǎn)生MeshUnitNum×(ContourNum+1)個(gè)等值多邊形.而多邊形最壞情況是一個(gè)五邊形，需要保存5個(gè)頂點(diǎn)信息，每個(gè)頂點(diǎn)則需要保存xyz坐標(biāo)和標(biāo)量值v.因此數(shù)組NephogramPolygons的長(zhǎng)度是MeshUnitNum×(ContourNum+1)×5×4.

云圖算法執(zhí)行后，NephogramPolygons數(shù)組同樣可能存在大量的空間浪費(fèi)，很有可能其中很多元素并沒(méi)有保存真實(shí)有意義的數(shù)據(jù).所以我們同樣可以將NephogramPolygons數(shù)組所有元素初始化為DBL_MAX值，當(dāng)整個(gè)網(wǎng)格處理完畢，將NephogramPolygons數(shù)組內(nèi)容復(fù)制到CPU內(nèi)存中，再對(duì)其內(nèi)容進(jìn)行過(guò)濾，得到所有真實(shí)等值多邊形的數(shù)據(jù).

如圖3所示是等值線云圖并行算法的核函數(shù)2的整體流程圖.首先判斷當(dāng)前網(wǎng)格是否包含等值點(diǎn)：如果不包含則調(diào)用WholeUnitToOneTriangle函數(shù)將網(wǎng)格單元整體轉(zhuǎn)換為一個(gè)等值三角形，并存入NephogramPolygons對(duì)應(yīng)位置；如果包含則調(diào)用SplitUnitToPolygons函數(shù).SplitUnitToPolygons函數(shù)將處理流程分為“兩小一大”和“一小兩大”兩種情況，分別調(diào)用函數(shù)TwoSmallOneBig和OneSmallTwoBig.

圖3 等值線云圖并行算法流程圖

如圖4所示，是“兩小一大”型的函數(shù)流程圖.該函數(shù)比較簡(jiǎn)單，處理流程中不會(huì)出現(xiàn)分支.

圖4 “兩小一大”函數(shù)流程圖

如圖5所示，是OneSmallTwoBig函數(shù)流程圖.OneSmallTwoBig函數(shù)中后續(xù)后續(xù)的發(fā)展情況稍微復(fù)雜一些，出現(xiàn)了幾種不同的分支處理.

圖5 “兩大一小”函數(shù)流程圖

3.3 完整的等值線云圖并行算法

前面分別討論了兩個(gè)核函數(shù)，核函數(shù)1用于并行計(jì)算等值點(diǎn)對(duì)組成的等值線段，核函數(shù)2用于并行計(jì)算等值多邊形.之所以不將兩個(gè)過(guò)程合并在一起，是因?yàn)閮蓚€(gè)核函數(shù)中線程劃分的粒度不同.核函數(shù)1中每個(gè)線程只處理一個(gè)網(wǎng)格單元和一條等值線的相交，而在核函數(shù)2中每個(gè)線程要處理一個(gè)網(wǎng)格單元和所有等值線的相交.另外兩個(gè)核函數(shù)中對(duì)GPU內(nèi)存的使用也有所不同，所以我們不得不將完整的并行算法拆分為兩個(gè)步驟.

如圖6所示，是完整的等值線云圖并行算法示意圖.CPU端主機(jī)程序啟動(dòng)后，第1步進(jìn)行內(nèi)存準(zhǔn)備.在CPU內(nèi)存和GPU顯存中創(chuàng)建對(duì)應(yīng)的3個(gè)數(shù)組，分別用于保存網(wǎng)格單元數(shù)據(jù)、等值線數(shù)據(jù)和等值點(diǎn)數(shù)據(jù);第2步從本地磁盤文件中讀入網(wǎng)格數(shù)據(jù)和等值線數(shù)據(jù)并存入CPU內(nèi)存中的網(wǎng)格單元數(shù)組和等值線數(shù)組，然后拷貝到GPU內(nèi)存對(duì)應(yīng)數(shù)組中，作為輸入數(shù)據(jù);第3步，啟動(dòng)核函數(shù)1，計(jì)算得到所有網(wǎng)格單元中所有成對(duì)出現(xiàn)的等值點(diǎn)數(shù)據(jù)，保存到GPU內(nèi)存中的等值點(diǎn)數(shù)組中，并將其拷貝回CPU內(nèi)存中的等值點(diǎn)數(shù)組作為輸出.到此為止，完成了第1階段，完成求解等值點(diǎn)的操作;第4步，對(duì)得到的等值點(diǎn)數(shù)組進(jìn)行數(shù)據(jù)壓縮，并拷貝到GPU內(nèi)存中創(chuàng)建對(duì)應(yīng)數(shù)組，作為第2階段的輸入;第5步啟動(dòng)核函數(shù)2，利用網(wǎng)格單元數(shù)組、等值線數(shù)組和壓縮后的等值點(diǎn)數(shù)組作為輸入，計(jì)算得到所有等值多邊形數(shù)據(jù)，并拷貝回CPU內(nèi)存中存入等值多邊形數(shù)組.至此，完成了等值線云圖的所有計(jì)算過(guò)程，再利用OpenGL等技術(shù)完成繪制.

圖6 等值線云圖并行算法整體流程圖

4 算法實(shí)現(xiàn)與分析

4.1 與其他云圖生成算法比較分析

以上我們基于“分割窮舉”思想設(shè)計(jì)了一種等值線云圖的生成算法，適用于三角形不規(guī)則網(wǎng)格.該算法與其他現(xiàn)有算法相比，在速度上有較大提升，具體有兩方面原因.

1)云圖串行算法時(shí)間復(fù)雜度降低

與現(xiàn)有的其他等值線云圖算法進(jìn)行比較，我們可以分成兩個(gè)步驟來(lái)討論算法效率.首先第一步生成等值線.在這個(gè)步驟中，無(wú)論是基于等值線序列法還是基于網(wǎng)格序列法的算法，時(shí)間效率都是O(m×n)，本質(zhì)上都需要對(duì)每個(gè)網(wǎng)格單元與每條等值線的相交情況進(jìn)行一一判斷和處理.在這個(gè)步驟中，我們的算法并沒(méi)有速度優(yōu)勢(shì).如果考慮到后期等值線由折線順滑為曲線的過(guò)程，我們的網(wǎng)格序列法還要進(jìn)行等值線段的排序，會(huì)產(chǎn)生額外的計(jì)算量.從這個(gè)角度看，我們的算法在生成等值線段這一步驟反而更慢.當(dāng)然，在本文中，我們并沒(méi)有必要將等值線曲線化，所以這個(gè)額外的計(jì)算量其實(shí)并沒(méi)有產(chǎn)生.

第2個(gè)步驟是根據(jù)等值線數(shù)據(jù)填充生成云圖，這個(gè)步驟現(xiàn)有的絕大多數(shù)算法都是基于等值線序列法[2，3]的.而在等值線序列法中，從等值線數(shù)據(jù)出發(fā)完成云圖填充，需要進(jìn)行額外的計(jì)算，例如需要完成等值區(qū)域的排序、判斷包含關(guān)系、凹多邊形區(qū)域處理等.例如在文獻(xiàn)[3]中需要額外的計(jì)算量的時(shí)間復(fù)雜度為O(n2)，其中n是等值線條數(shù).

而在我們提出的云圖算法中，云圖中等值多邊形的計(jì)算是混合在等值線段計(jì)算之后同步完成的(當(dāng)然也可以先計(jì)算等值線段，保存后再計(jì)算等值多邊形).從時(shí)間復(fù)雜度上看，完整的等值線云圖算法的全部時(shí)間復(fù)雜度就是O(m×n)，并沒(méi)有增加.當(dāng)然，從等值線段計(jì)算出等值多邊形，還需要花費(fèi)一些額外計(jì)算，但是并沒(méi)有在數(shù)量級(jí)上增加算法的時(shí)間消耗.因此，如果考慮云圖整體的算法效率，我們的云圖串行算法比其他現(xiàn)有基于等值線數(shù)據(jù)填充思想的云圖算法要快.

2)云圖并行算法速度顯著提升

我們的云圖算法的第2個(gè)優(yōu)勢(shì)是，非常適合于進(jìn)行并行化，這也是本算法速度顯著提升的根本原因.

在我們的云圖串行算法中，計(jì)算過(guò)程被分割到每個(gè)網(wǎng)格單元內(nèi)部，在不同網(wǎng)格單元之間不存在任何先后依賴關(guān)系.這就允許云圖串行算法并行化.而實(shí)際結(jié)果證明，這種并行化后的算法，時(shí)間效率提示了一個(gè)數(shù)量級(jí).

4.2 云圖生成算法實(shí)現(xiàn)與結(jié)果分析

我們利用C++語(yǔ)言實(shí)現(xiàn)了上述并行算法，開(kāi)發(fā)環(huán)境為Win10+VS2019+CUDA11.使用普通筆記本電腦，配置為intel i7-10510U處理器，16G內(nèi)存，NVIDIA GeForce MX250顯卡，2G顯存.顯卡CUDA版本為10，擁有3個(gè)流多處理器，每個(gè)流多處理器上有128個(gè)流處理器內(nèi)核，共384個(gè)內(nèi)核.

為了驗(yàn)證算法效果，我們選取了兩個(gè)典型溫度場(chǎng)數(shù)據(jù)進(jìn)行了等值線和云圖的計(jì)算以及繪制試驗(yàn).我們利用OpenGL實(shí)現(xiàn)了圖像繪制.在OpenGL中，多邊形的顏色填充可以采用兩種不同的模式.第1種模式可以為整個(gè)多邊形指定一種固定的顏色，然后用這個(gè)單一顏色填充，稱為平面著色(FlatShading).第2種模式，可以給多邊形每個(gè)頂點(diǎn)指定不同的顏色，然后用這些不同顏色自動(dòng)實(shí)現(xiàn)過(guò)渡填充，這種方式稱為光滑著色處理(Smooth Shading).OpenGL中利用函數(shù)“void glShadeModel(GLenum mode)”來(lái)實(shí)現(xiàn)兩種模式的設(shè)置，當(dāng)參數(shù)mode取值為GL_FLAT或GL_SMOOTH時(shí)，分別代表平面著色和光滑著色.

采用平面著色模式時(shí)，填充的等值線云圖具有明確的條紋分層，我們稱為條紋云圖.采用光滑著色模式時(shí)，填充的等值線云圖具有光滑的顏色過(guò)渡，效果非常像像素級(jí)云圖，我們稱為偽像素云圖.如圖7所示，是一組渦輪盤橫截面的溫度場(chǎng)效果圖.從上到下分別為等值線效果圖、條紋云圖效果圖、和偽像素云圖效果圖.

圖7 等值線和云圖效果圖1

如圖8所示，是一組異形墊片的溫度場(chǎng)效果圖.從上到下分別為等值線效果圖、條紋云圖效果圖、和偽像素云圖效果圖.從圖7和圖8中，可以看出條紋云圖和偽像素云圖的區(qū)別很明顯.也能看出等值線以折線形式實(shí)現(xiàn)，并沒(méi)有光滑成曲線形式.

圖8 等值線和云圖效果圖2

為了更好的觀察并行算法的加速效果，我們利用Delaunay算法生成了20組三角網(wǎng)格，網(wǎng)格單元數(shù)量分別是10，16，24，34，40，72，98，188，256，382，480，922，1410，2114，2808，3700，4408，5266，6448，7623.然后針對(duì)這20個(gè)網(wǎng)格分別創(chuàng)建了一個(gè)溫度場(chǎng)，在其中都設(shè)置20條等值線.最后利用上述思想分別編寫了CPU串行算法和GPU并行算法，分別運(yùn)行得到結(jié)果.

我們首先單獨(dú)運(yùn)行等值線算法，分別得到了20組串行算法耗時(shí)和并行算法耗時(shí)數(shù)據(jù)，如表1所示.再根據(jù)表1數(shù)據(jù)繪制了20組數(shù)據(jù)的加速比變化曲線，如圖9所示.

圖9 等值線并行算法加速比變化圖

表1 串行、并行等值線算法耗時(shí)

從圖9中可以看到，隨著網(wǎng)格數(shù)量的增加，并行算法的優(yōu)勢(shì)很快體現(xiàn)出來(lái).當(dāng)網(wǎng)格數(shù)量達(dá)到200左右時(shí)，乘以等值線數(shù)量20，總共需要約200×20=4000個(gè)線程并行處理.這時(shí)基本達(dá)到了試驗(yàn)用筆記本的性能極限，后續(xù)網(wǎng)格數(shù)量繼續(xù)增加，并行算法的加速比并沒(méi)有繼續(xù)提高.

然后使用同樣的20組數(shù)據(jù)，我們完整運(yùn)行了等值線云圖算法，得到的結(jié)果如表2所示.再根據(jù)表2數(shù)據(jù)繪制了云圖算法的20組數(shù)據(jù)加速比變化曲線，如圖10所示.

表2 串行、并行云圖算法耗時(shí)

對(duì)比圖10與圖9，可以發(fā)現(xiàn)有3處明顯不同：

圖10 云圖并行算法加速比變化圖

1)當(dāng)網(wǎng)格數(shù)量較少時(shí)，云圖并行算法加速比甚至小于1.這說(shuō)明在網(wǎng)格數(shù)量較少時(shí)，并行算法并沒(méi)有體現(xiàn)出優(yōu)勢(shì)，反而因?yàn)閮?nèi)存分配、數(shù)據(jù)拷貝、線程調(diào)度等處理的耗時(shí)，最終需要花費(fèi)比串行算法給多的時(shí)間，才能完成同樣的計(jì)算.

2)云圖算法加速比遞增的速度明顯要小等值線算法的.這主要是因?yàn)榈戎稻€算法中，所需線程總數(shù)等于網(wǎng)格個(gè)數(shù)乘以等值線條數(shù).而在云圖算法中，線程總數(shù)就等于網(wǎng)格個(gè)數(shù).因此隨著網(wǎng)格的增加，所需的線程數(shù)量在等值線算法中增加的更快，也就能更充分的利用GPU線程多的優(yōu)勢(shì)，因此在等值線算法中，加速比更快增加到最優(yōu)程度.

3)同樣的計(jì)算機(jī)配置，云圖算法最終達(dá)到的加速比只有10左右，明顯比等值線并行算法最高加速比13要低.這主要是因?yàn)樵茍D算法中每個(gè)線程需要進(jìn)行的計(jì)算量比較大，導(dǎo)致最終達(dá)到的加速效果要差一些.

5 結(jié) 論

本文首先討論了不規(guī)則網(wǎng)格中，一種基于“網(wǎng)格序列法”和“分割窮舉”思想的等值線云圖生成算法.然后討論了這種算法的并行化策略，并基于CUDA編程實(shí)現(xiàn)了并行算法.最后利用20組數(shù)據(jù)進(jìn)行了試驗(yàn)，結(jié)果表明本文提出的云圖并行算法是簡(jiǎn)單高效的.

在后續(xù)工作中，我們將研究該算法在其他網(wǎng)格類型中的實(shí)現(xiàn)，例如在四邊形網(wǎng)格中實(shí)現(xiàn)云圖并行算法.另外，我們也會(huì)深入探討其他可視化方法的并行化可能，例如流線、橢球圖標(biāo)等算法的并行化.