基于CUDA的等值線云圖并行生成算法

杜小甫，劉輝林，劉鶴丹

1(東北大學 計算機科學與工程學院，沈陽 110819)

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

1 引 言

1.1 等值線和等值線云圖

等值線，是指標量場中具有相同標量值的各點的連線.等值線作為一種標量場可視化技術，應用非常廣泛.根據(jù)標量場物理含義不同，常見的等值線包括等溫線，等壓線，等高線，等勢線等.

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

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

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

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

1.2 CUDA編程架構

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

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

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

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

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

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

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

2.1 等值線段生成算法

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

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

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

(1)

(1)

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

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

1)等值點按網(wǎng)格單元編號排列.

2)等值點按標量場值由小到大排列.

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

2.2 等值線云圖生成算法

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

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

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

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

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

3 算法并行化

3.1 等值線生成算法并行化

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

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

2)GPU內(nèi)存設計.GPU并行程序中，必須盡量減少不同線程之間內(nèi)存訪問沖突的可能性.因為當不同線程不得不同時訪問相同內(nèi)存時，這些線程不能并行運行，只能串行運行.這將極大影響并行速度，因此應該讓不同線程盡量訪問不同的內(nèi)存.

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

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

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

meshID=TID/ContourNum

(3)

contourID=TID%ContourNum

(4)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1)云圖串行算法時間復雜度降低

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

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

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

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

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

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

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

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

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

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

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

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

圖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.然后針對這20個網(wǎng)格分別創(chuàng)建了一個溫度場，在其中都設置20條等值線.最后利用上述思想分別編寫了CPU串行算法和GPU并行算法，分別運行得到結果.

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

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

表1 串行、并行等值線算法耗時

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

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

表2 串行、并行云圖算法耗時

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

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

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

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

3)同樣的計算機配置，云圖算法最終達到的加速比只有10左右，明顯比等值線并行算法最高加速比13要低.這主要是因為云圖算法中每個線程需要進行的計算量比較大，導致最終達到的加速效果要差一些.

5 結 論

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

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