范良，張曉蓉，吳亞東，陳呈，王昉

1. 西南科技大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院，四川 綿陽 621010；2. 四川輕化工大學(xué)計算機(jī)科學(xué)與工程學(xué)院，四川 自貢 643002；3. 空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621000；4. 中國空氣動力研究與發(fā)展中心計算空氣動力研究所，四川 綿陽 621000

1 引言

體可視化是科學(xué)可視化中的重要研究方向，被廣泛應(yīng)用于醫(yī)學(xué)影像、計算流體力學(xué)（computational fluid dynamics，CFD）、氣象氣候?qū)W等領(lǐng)域，體繪制算法是體可視化領(lǐng)域的重要研究內(nèi)容[1]。體繪制直接從三維的體數(shù)據(jù)中生成半透明的二維屏幕圖像，空間表現(xiàn)能力強(qiáng)，因此，也被稱為直接體繪制?？茖W(xué)家可以通過體繪制結(jié)果看穿數(shù)據(jù)體，探索數(shù)據(jù)內(nèi)部結(jié)構(gòu)，獲取有用信息。與面繪制相比，體繪制結(jié)果多了一個維度的信息，因此，可視化結(jié)果不僅更加精確，而且更能體現(xiàn)數(shù)據(jù)本身的特點(diǎn)[2]。

體數(shù)據(jù)是體繪制算法的輸入，分為結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格[3]。其中，結(jié)構(gòu)網(wǎng)格具有規(guī)則的幾何結(jié)構(gòu)和拓?fù)浣Y(jié)構(gòu)，相反，非結(jié)構(gòu)網(wǎng)格的幾何結(jié)構(gòu)和拓?fù)浣Y(jié)構(gòu)不規(guī)則，網(wǎng)格單元可以是任意形狀。在三維數(shù)據(jù)場變化劇烈的區(qū)域可以用體積較小的網(wǎng)格單元表示，變化比較平滑的區(qū)域則使用體積較大的網(wǎng)格單元表示，對于同一個數(shù)據(jù)場，可以用更少的非結(jié)構(gòu)網(wǎng)格單元進(jìn)行描述，有效節(jié)省了存儲空間[4]。因此，非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)在科學(xué)計算與仿真領(lǐng)域的應(yīng)用更加廣泛。然而，由于非結(jié)構(gòu)網(wǎng)格的幾何結(jié)構(gòu)和拓?fù)浣Y(jié)構(gòu)不規(guī)則，相對于結(jié)構(gòu)網(wǎng)格，非結(jié)構(gòu)網(wǎng)格的體繪制算法設(shè)計比較困難，涉及復(fù)雜的計算與龐大的存儲量。目前，非結(jié)構(gòu)網(wǎng)格體繪制算法可以通過基于圖像空間、基于對象空間和兩者混合的技術(shù)實(shí)現(xiàn)[2]。光線投射算法[5]是十分具代表性的基于圖像空間實(shí)現(xiàn)的體繪制技術(shù)。然而，對于非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)而言，無法根據(jù)幾何結(jié)構(gòu)和拓?fù)浣Y(jié)構(gòu)計算光線離開當(dāng)前網(wǎng)格單元后即將進(jìn)入的下一個網(wǎng)格單元，采用局部或全局搜索下一個網(wǎng)格單元需要極大的開銷。Weiler M等人[6]將沿著光線的采樣點(diǎn)限制在單元面上，同時，采用預(yù)積分技術(shù)，將預(yù)積分的顏色和透明度值存入一張紋理圖，并且將光線進(jìn)入和離開單元點(diǎn)的標(biāo)量值以及光線在網(wǎng)格單元內(nèi)穿過的長度作為參數(shù)查詢紋理圖。此外，他們還精心設(shè)計了紋理數(shù)據(jù)結(jié)構(gòu)，將體數(shù)據(jù)以紋理的形式載入顯存，通過硬件技術(shù)加速算法。投影四面體算法[7]是十分著名的基于對象空間實(shí)現(xiàn)的體繪制技術(shù)，首先，計算所有網(wǎng)格單元的可見性順序；然后，根據(jù)可見性順序?qū)⑺拿骟w單元在圖像平面上的投影拆分為一系列三角形；最后，將三角形的頂點(diǎn)數(shù)據(jù)載入顯存，通過GPU光柵化并合成最終結(jié)果。此外，Wylie B等人[8]利用可編程的頂點(diǎn)著色器，設(shè)計了基于硬件加速的投影四面體算法，不需要額外的數(shù)據(jù)結(jié)構(gòu)與數(shù)據(jù)預(yù)處理。HAVS（hardware-assisted visibility svrting）[9]算法是最具代表性的基于圖像/對象空間混合技術(shù)實(shí)現(xiàn)的非結(jié)構(gòu)網(wǎng)格體繪制算法，在對象空間，對網(wǎng)格單元執(zhí)行局部排序，并生成一個近似有序的網(wǎng)格單元片段列表；在圖像空間，使用固定深度的排序網(wǎng)絡(luò)對片段列表執(zhí)行升序排序。其中，對象空間計算在CPU上進(jìn)行，圖像空間計算在GPU上進(jìn)行。HAVS算法在簡化CPU的處理運(yùn)算的同時，將大量排序計算交給GPU執(zhí)行，進(jìn)一步提升了算法性能。

由于各種科學(xué)計算與模擬仿真的精度越來越高，非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)規(guī)模呈爆炸式增長。體繪制包含數(shù)據(jù)重構(gòu)、傳遞函數(shù)映射和繪制合成，涉及大量的浮點(diǎn)運(yùn)算，串行的算法難以滿足需求。隨著技術(shù)發(fā)展，大量的并行計算軟件和硬件出現(xiàn)，這為科學(xué)研究提供了更加強(qiáng)勁的計算能力，基于進(jìn)程并行的體繪制算法研究成為大規(guī)模非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)體可視化設(shè)計的趨勢[10]。Fogal T等人[11]實(shí)現(xiàn)了基于集群的并行體可視化算法，使用KD樹組織數(shù)據(jù)，進(jìn)程之間采用信息傳遞接口（massage passing interface，MPI）通信，可高效繪制大規(guī)模標(biāo)量場數(shù)據(jù)集。Chen J H等人[12]設(shè)計了分布式的并行體繪制算法，并對大規(guī)模燃燒模擬數(shù)據(jù)的多個變量同時進(jìn)行體繪制，允許用戶實(shí)時交互體繪制結(jié)果。本文在現(xiàn)有研究成果的基礎(chǔ)上，研究并實(shí)現(xiàn)了一種基于sort-last[13]架構(gòu)的非結(jié)構(gòu)網(wǎng)格并行體繪制算法，并集成到國產(chǎn)自主可控科學(xué)可視化軟件拓視（TopViz）中，主要貢獻(xiàn)如下。

● 實(shí)現(xiàn)了基于KD樹的并行體數(shù)據(jù)分割算法和可見性排序算法。首先，使用預(yù)處理程序確保每個進(jìn)程讀取的網(wǎng)格數(shù)量大致相等；其次，以每個進(jìn)程上的數(shù)據(jù)為輸入，構(gòu)造局部KD樹；然后，將所有進(jìn)程上的局部KD樹合成為全局KD樹；最后，通過相機(jī)位置和分割平面計算全局KD樹葉子節(jié)點(diǎn)的可見性順序。

● 實(shí)現(xiàn)了基于二叉樹的并行圖像合成算法。根據(jù)可見性順序?qū)λ羞M(jìn)程上的圖片兩兩配對并合成中間結(jié)果，不斷重復(fù)配對與合成的過程，直到合成最終結(jié)果。

● 實(shí)現(xiàn)了兩層LOD交互模型。首先，同步所有進(jìn)程的數(shù)據(jù)邊界信息；然后，在主進(jìn)程上渲染數(shù)據(jù)邊框模型。當(dāng)處于交互狀態(tài)時，渲染邊框模型；停止交互時，渲染體繪制模型。

● 設(shè)計實(shí)驗(yàn)驗(yàn)證，并分析了算法性能。實(shí)驗(yàn)結(jié)果顯示，本文的并行體數(shù)據(jù)分割算法具有良好的性能和負(fù)載均衡能力；并行體繪制算法可以被很好地應(yīng)用于大規(guī)模非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)體可視化，且交互時延在毫秒級別，滿足實(shí)時交互的需求。

2 基于進(jìn)程并行的非結(jié)構(gòu)網(wǎng)格體繪制算法

非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)的體繪制算法涉及龐大的計算量與存儲量，串行算法效率較低，難以滿足大規(guī)模非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)體可視化實(shí)時交互的需求。因此，基于進(jìn)程并行的體繪制算法是大規(guī)模非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)體可視化十分有效的解決方案。本節(jié)將詳細(xì)描述基于進(jìn)程并行的非結(jié)構(gòu)網(wǎng)格體繪制算法流程，算法包含數(shù)據(jù)分割、繪制和圖片合成3個階段，如圖1所示。在數(shù)據(jù)分割階段，首先，每個進(jìn)程讀取數(shù)量大致相等的網(wǎng)格單元；然后，將每個進(jìn)程上的體數(shù)據(jù)作為輸入，構(gòu)造局部KD樹；最后，通過MPI將所有進(jìn)程上的局部KD樹合成為全局KD樹。在繪制階段，每個進(jìn)程使用獨(dú)立的體繪制管線計算圖像。在圖片合成階段，按照可見性順序配對所有進(jìn)程上的圖片，并合成中間結(jié)果，不斷重復(fù)配對和合成過程，直到合成最終結(jié)果。

圖1 基于進(jìn)程并行的非結(jié)構(gòu)網(wǎng)格體繪制算法流程

2.1 數(shù)據(jù)分割

并行體繪制算法的第一步是將大規(guī)模的網(wǎng)格數(shù)據(jù)拆分為若干個規(guī)模較小的網(wǎng)格。由于繪制過程的性能取決于繪制最慢的進(jìn)程效率，因此，負(fù)載均衡是數(shù)據(jù)分割算法設(shè)計首要考慮的問題。本節(jié)將詳細(xì)描述基于KD樹的并行體數(shù)據(jù)分割算法的設(shè)計，為了確保負(fù)載均衡，最終分割結(jié)果中KD樹葉子節(jié)點(diǎn)之間的數(shù)據(jù)量應(yīng)基本相同。

2.1.1 局部KD樹構(gòu)造

KD樹是一種分割K維空間的高效數(shù)據(jù)結(jié)構(gòu)，屬于二叉空間分割樹的特殊情況，被廣泛應(yīng)用于各種數(shù)據(jù)分割和搜索程序中[14]。其中，KD樹的葉子節(jié)點(diǎn)表示K維空間中的點(diǎn)，非葉子節(jié)點(diǎn)表示將K維空間一分為二的超平面或超立方體。如果要將KD樹算法應(yīng)用于體數(shù)據(jù)分割，則需要賦予節(jié)點(diǎn)新的意義，其中，葉子節(jié)點(diǎn)表示非結(jié)構(gòu)網(wǎng)格單元集合，非葉子節(jié)點(diǎn)則表示將空間一分為二的超平面。體數(shù)據(jù)是三維的，因此，使用與坐標(biāo)軸平行的分割平面可以最大限度地簡化數(shù)據(jù)分割的運(yùn)算過程。

算法1給出了KD樹的構(gòu)造過程。對于一個非葉子節(jié)點(diǎn)，首先，計算所有網(wǎng)格單元的質(zhì)心坐標(biāo)；其次，計算分割平面的方向，為了確保與分割平面相交的網(wǎng)格單元最少和避免數(shù)據(jù)分割后被拉長，選擇數(shù)據(jù)邊界跨度最長的坐標(biāo)軸方向作為分割平面方向；然后，計算所有網(wǎng)格單元沿著分割平面方向的質(zhì)心坐標(biāo)中位數(shù)，為了確保左右子樹的網(wǎng)格單元數(shù)量大致相等，將質(zhì)心坐標(biāo)的中位數(shù)作為分割平面的位置坐標(biāo)；最后，沿著分割平面將數(shù)據(jù)一分為二，對于與分割平面相交的網(wǎng)格單元，復(fù)制一份并劃分到分割平面兩側(cè)。KD樹構(gòu)造是一個遞歸的過程，不斷拆分非葉子節(jié)點(diǎn)，直到生成指定數(shù)目的葉子節(jié)點(diǎn)。

算法1KD樹構(gòu)造

輸入：非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)

輸出：KD樹根節(jié)點(diǎn)指針

步驟1：計算所有網(wǎng)格單元的質(zhì)心坐標(biāo)；

步驟2：根據(jù)網(wǎng)格數(shù)據(jù)邊界值計算分割平面方向；

步驟3：計算質(zhì)心坐標(biāo)的中位數(shù)，將中位數(shù)質(zhì)心坐標(biāo)作為分割面的位置坐標(biāo)；

步驟4：沿著分割面將網(wǎng)格數(shù)據(jù)一分為二；

步驟5：遞歸執(zhí)行步驟2～4，直到生成指定數(shù)目的葉子節(jié)點(diǎn)。

通常，一個大規(guī)模的網(wǎng)格數(shù)據(jù)在超級計算機(jī)上使用多個并行程序模擬計算生成，以多塊（multi-block）數(shù)據(jù)的形式寫入磁盤，最好也用并行讀取的方法處理。雖然單塊數(shù)據(jù)的規(guī)模不大，但是單塊數(shù)據(jù)之間的網(wǎng)格規(guī)模可能差距很大。因此，本文首先使用一個預(yù)處理程序統(tǒng)計大規(guī)模網(wǎng)格數(shù)據(jù)中每一塊數(shù)據(jù)的網(wǎng)格單元數(shù)量，并保存到一個列表中；然后，根據(jù)列表網(wǎng)格數(shù)量信息計算每個進(jìn)程應(yīng)讀取的單塊數(shù)據(jù)編號，確保每個進(jìn)程讀取的網(wǎng)格單元數(shù)量大致相等；最后，每個進(jìn)程根據(jù)單塊數(shù)據(jù)的編號并行讀取體數(shù)據(jù)。

在每個進(jìn)程上使用KD樹算法分割體數(shù)據(jù)，這個過程被稱為局部KD樹的構(gòu)造。預(yù)處理過程確保了每個進(jìn)程上讀取的網(wǎng)格單元數(shù)量大致相等，因此，局部KD樹構(gòu)造算法是負(fù)載均衡的。每個進(jìn)程上的局部KD樹具有相同的非葉子節(jié)點(diǎn)，即分割平面相同，所有的樹節(jié)點(diǎn)同步生成。因此，首先構(gòu)造一個非葉子節(jié)點(diǎn)，同步所有進(jìn)程上該節(jié)點(diǎn)對應(yīng)的體數(shù)據(jù)邊界信息，根據(jù)邊界范圍選擇分割平面的方向；然后，并行計算所有進(jìn)程上該節(jié)點(diǎn)對應(yīng)的體數(shù)據(jù)網(wǎng)格單元沿著分割平面方向質(zhì)心坐標(biāo)的中位數(shù)[15]，將中位數(shù)的坐標(biāo)作為分割平面的位置坐標(biāo)；最后，沿著分割平面將該節(jié)點(diǎn)一分為二。遞歸構(gòu)造每一個非葉子節(jié)點(diǎn)，直到生成指定數(shù)目的葉子節(jié)點(diǎn)。為了最大化地確保KD樹平衡，以及便于后續(xù)設(shè)計并行圖像合成算法，本文構(gòu)造的KD樹均為滿二叉樹，即葉子節(jié)點(diǎn)數(shù)量為2的正整數(shù)冪。如果葉子節(jié)點(diǎn)數(shù)量為N，則KD樹的高度為lbN。

2.1.2 全局KD樹合成

全局KD樹與局部KD樹具有相同的非葉子節(jié)點(diǎn)，不同的是，每個進(jìn)程只存儲全局KD樹中一個葉子節(jié)點(diǎn)的數(shù)據(jù)，這些葉子節(jié)點(diǎn)的數(shù)據(jù)為最終分割結(jié)果。進(jìn)程之間的局部KD樹對應(yīng)的葉子節(jié)點(diǎn)數(shù)據(jù)交換過程被稱為全局KD樹的合成。數(shù)據(jù)分割需要保證每個進(jìn)程上都有數(shù)據(jù)，因此，進(jìn)程數(shù)量應(yīng)該等于KD樹的葉子節(jié)點(diǎn)數(shù)量，即2的正整數(shù)冪。

對于局部KD樹，只有其中一個葉子節(jié)點(diǎn)的數(shù)據(jù)屬于當(dāng)前進(jìn)程，因此，需要進(jìn)行N對N的數(shù)據(jù)交換操作。首先，給所有局部KD樹的葉子節(jié)點(diǎn)編號（從0開始）；然后，通過MPI，根據(jù)進(jìn)程秩將屬于其他進(jìn)程的葉子節(jié)點(diǎn)數(shù)據(jù)發(fā)送至對應(yīng)的進(jìn)程；最后，將每個進(jìn)程上的數(shù)據(jù)合并為完整的網(wǎng)格。

2.1.3 基于KD樹的可見性排序算法

對于與分割平面相交的網(wǎng)格單元（也被稱為ghost單元[2]），本文采取的處理方式是復(fù)制一份，然后將其劃分到分割平面兩側(cè)。該處理方式可以確保分割平面兩側(cè)的網(wǎng)格數(shù)據(jù)具有嚴(yán)格的可見性順序，不存在相互遮擋的情況。其中，可見性順序[16]用于描述空間中兩個對象的遮擋情況，在給定視線方向的情況下，如果對象a遮擋了對象b，則可見性順序Oa＜Ob；反之，則可見性順序Oa＞Ob。對于KD樹的非葉子節(jié)點(diǎn)的左右孩子節(jié)點(diǎn)，以分割平面為界，定義與相機(jī)位置處于同一側(cè)的孩子節(jié)點(diǎn)為近點(diǎn)，另一側(cè)為遠(yuǎn)點(diǎn)。近點(diǎn)的可見性順序明顯小于遠(yuǎn)點(diǎn)的可見性順序，因此，通過相機(jī)位置和分割平面可以計算KD樹中所有葉子節(jié)點(diǎn)的可見性順序。

算法2給出了基于KD樹的可見性順序計算過程。將全局KD樹作為輸入，從根節(jié)點(diǎn)開始，根據(jù)相機(jī)位置確定近點(diǎn)與遠(yuǎn)點(diǎn)；按照“先近點(diǎn)，后遠(yuǎn)點(diǎn)”的順序計算每個非葉子節(jié)點(diǎn)的可見性順序；如果遇到葉子節(jié)點(diǎn)，則將葉子節(jié)點(diǎn)對應(yīng)的進(jìn)程秩添加到可見性順序列表。遞歸處理每個節(jié)點(diǎn)，直到所有葉子節(jié)點(diǎn)的可見性順序計算完畢。

算法2可見性順序計算

輸入：全局KD樹

輸出：可見性順序列表list

步驟1：如果當(dāng)前節(jié)點(diǎn)為葉子節(jié)點(diǎn)，將該節(jié)點(diǎn)的進(jìn)程秩添加到list；

步驟2：如果當(dāng)前節(jié)點(diǎn)為非葉子節(jié)點(diǎn)，根據(jù)相機(jī)位置確定近點(diǎn)與遠(yuǎn)點(diǎn)：

· 計算近點(diǎn)可見性順序；

· 計算遠(yuǎn)點(diǎn)可見性順序；

步驟3：遞歸執(zhí)行步驟1～2，直到所有葉子節(jié)點(diǎn)的可見性順序計算完畢。

2.2 繪制

并行體繪制算法的第二個階段（繪制）以第2.1節(jié)中KD樹算法的分割結(jié)果為輸入，在每個進(jìn)程上執(zhí)行獨(dú)立的可視化管線，從而計算體繪制圖像。本文采用體繪制管線技術(shù)，如圖2所示，整個體繪制管線分為數(shù)據(jù)處理、體繪制算法、體繪制參數(shù)和渲染4個部分，其中，管線中每個階段的輸出作為下一個階段的輸入。在數(shù)據(jù)處理部分，使用Reader加載KD樹算法的分割結(jié)果，如果需要對體數(shù)據(jù)進(jìn)行裁剪和四面體化等操作，可以通過Filter完成；在體繪制算法部分，使用高效的線程并行投影四面體算法[17]計算圖元信息；在體繪制參數(shù)部分，建立體數(shù)據(jù)標(biāo)量值與顏色/透明度之間的映射關(guān)系；在最后的渲染部分，使用OpenGL API函數(shù)將體繪制圖元數(shù)據(jù)與傳遞函數(shù)參數(shù)載入顯存，通過GPU光柵化并合成最終的體繪制圖像。體繪制管線技術(shù)具有高度的可擴(kuò)展性與靈活性，管線每個部分只提供通用的接口，實(shí)現(xiàn)體繪制管線則需要提供所有通用接口的具體實(shí)現(xiàn)。如果要研究不同的體繪制算法性能，只需要提供新算法的Mapper即可。

圖2 體繪制管線

2.3 圖片合成

如圖3所示，基于sort-last架構(gòu)的并行體繪制算法最后需要將所有進(jìn)程上的體繪制圖片合成最終結(jié)果。由于體繪制結(jié)果是半透明的二維圖片，因此，需要根據(jù)可見性順序按照從前向后的順序進(jìn)行合成?？梢允褂胦ver運(yùn)算將兩張體繪制圖片合成一張圖片，其中，over運(yùn)算不滿足交換律，但是滿足結(jié)合律[18]。因此，在不改變合成順序的前提下，可以對所有進(jìn)程上的體繪制圖片任意分組，采取分組合成策略。over運(yùn)算的性質(zhì)為并行圖像合成算法設(shè)計提供了條件[15]。

圖3 圖片合成

本文采用基于二叉樹的圖像合成算法，首先，根據(jù)第2.1.3節(jié)中的方法計算所有進(jìn)程體繪制圖片的可見性順序列表；然后，根據(jù)可見性順序列表對體繪制圖片進(jìn)行兩兩配對；最后，通過MPI將圖片發(fā)送至其配對進(jìn)程，并使用over運(yùn)算合成中間結(jié)果。不斷重復(fù)配對與合成的過程，直到所有中間結(jié)果都被合成為最終圖片。與串行的圖像合成算法相比，雖然基于樹的合成算法增加了通信次數(shù)和合成次數(shù)，但是它將時間重疊利用，減少了資源閑置，算法并行粒度更高，效率優(yōu)勢更明顯。

以4個進(jìn)程為例，二叉樹圖像合成算法流程如圖4所示。假設(shè)4個進(jìn)程上的體繪制圖片可見性順序?yàn)镺0＜O2＜O1＜O3，第一次合成，2號進(jìn)程和3號進(jìn)程分別將圖片發(fā)送至0號進(jìn)程和1號進(jìn)程，在0號進(jìn)程和1號進(jìn)程上合成中間結(jié)果；第二次合成，1號進(jìn)程將中間結(jié)果發(fā)送至0號進(jìn)程，在0號進(jìn)程上合成最終結(jié)果。

圖4 二叉樹圖像合成算法過程（4個進(jìn)程）

2.4 兩層LOD模型

體繪制算法涉及大量的距離運(yùn)算、求交運(yùn)算和插值運(yùn)算，同時，圖片合成與進(jìn)程之間的通信開銷不可忽略。對于大規(guī)模非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)體可視化，難以滿足實(shí)時交互的需求，因此，本文采用層次細(xì)節(jié)（level of detail，LOD）技術(shù)優(yōu)化用戶交互體驗(yàn)。如圖5所示，一共設(shè)計了兩層LOD模型，首先，所有從進(jìn)程將數(shù)據(jù)邊界信息發(fā)送至主進(jìn)程（0號進(jìn)程），在主進(jìn)程上同步邊界信息，并渲染一個剛好能包圍住體數(shù)據(jù)的邊框模型。當(dāng)處于交互狀態(tài)（按下鼠標(biāo)）時，渲染邊框模型；當(dāng)交互狀態(tài)結(jié)束（釋放鼠標(biāo)）時，開始渲染體繪制模型；當(dāng)體繪制模型渲染完成后，替換邊框模型。使用兩層LOD模型，用戶在交互時可以根據(jù)邊框模型判斷當(dāng)前體數(shù)據(jù)的位置和角度，同時，邊框模型的渲染在主進(jìn)程上進(jìn)行，只需要極小的開銷，并且不需要圖片合成，因此，可以滿足實(shí)時交互的需求。

圖5 兩層LOD模型

3 實(shí)驗(yàn)分析

為了研究算法的性能，本文設(shè)計了節(jié)點(diǎn)內(nèi)并行實(shí)驗(yàn)和跨節(jié)點(diǎn)并行實(shí)驗(yàn)。節(jié)點(diǎn)內(nèi)并行實(shí)驗(yàn)配置為高性能Intel服務(wù)器，其中，CPU為4個20核心的Intel(R) Xeon(R) E5-2698 v4＠ 2.20 GHz，內(nèi)存大小為500 GB?？绻?jié)點(diǎn)并行實(shí)驗(yàn)配置為國產(chǎn)高性能服務(wù)器，其中，CPU為飛騰處理器，內(nèi)存大小為120 GB。所有的算法均使用C/C++語言實(shí)現(xiàn)，并行框架采用MPI，編譯器的版本為gcc/g++ 8.3，操作系統(tǒng)為64位Ubuntu 16.04。所有的算法執(zhí)行時間通過C語言API函數(shù)ftime()獲取。同時，為了避免偶然性，所有的結(jié)果均通過多次重復(fù)實(shí)驗(yàn)計算平均值得出。

3.1 數(shù)據(jù)分割算法性能

為了研究數(shù)據(jù)分割算法性能，采用節(jié)點(diǎn)內(nèi)并行實(shí)驗(yàn)配置，從數(shù)據(jù)分割開銷、負(fù)載均衡測試和VTK（visualization toolkit）算法對比3個方面設(shè)計實(shí)驗(yàn)。首先，在不同進(jìn)程數(shù)量下研究數(shù)據(jù)分割后的網(wǎng)格數(shù)量增加情況；然后，對數(shù)據(jù)分割算法的負(fù)載均衡能力進(jìn)行測試；最后，對比了本文算法與VTK算法的性能。

3.1.1 數(shù)據(jù)分割開銷

在構(gòu)造局部KD樹時，與分割平面相交的單元會被復(fù)制一份，并被劃分到分割平面的兩側(cè)。因此，數(shù)據(jù)分割后，網(wǎng)格單元總數(shù)會增加。表1給出了使用不同數(shù)量的進(jìn)程分割數(shù)據(jù)后的網(wǎng)格規(guī)模。與預(yù)期相符，數(shù)據(jù)分割后，網(wǎng)格規(guī)模增大；并且，隨著進(jìn)程數(shù)量增加，更多的網(wǎng)格單元與分割平面相交，產(chǎn)生了更多的ghost單元，因此，網(wǎng)格規(guī)模隨之增大。ghost單元是數(shù)據(jù)分割開銷的重要部分，直接影響后續(xù)并行體繪制算法的效率。

表1 網(wǎng)格規(guī)模

3.1.2 負(fù)載均衡測試

一個好的數(shù)據(jù)分割算法應(yīng)該具有良好的負(fù)載均衡能力。本節(jié)將最多網(wǎng)格進(jìn)程與最少網(wǎng)格進(jìn)程的網(wǎng)格單元數(shù)量差值作為衡量標(biāo)準(zhǔn)，該值越小，說明算法負(fù)載均衡能力越好。表2給出了不同進(jìn)程數(shù)量下的網(wǎng)格單元數(shù)量差值，結(jié)果顯示，所有的差值均在4 000以內(nèi)，全局KD樹葉子節(jié)點(diǎn)之間的網(wǎng)格單元數(shù)量大致相等，本文數(shù)據(jù)分割算法的負(fù)載均衡能力較好。

表2 負(fù)載均衡測試結(jié)果

3.1.3 VTK算法對比

為了研究本文中的數(shù)據(jù)分割算法性能，將開源可視化框架VTK中的D3（vtkDistributedDataFilter）類作為對比，在不同進(jìn)程數(shù)量下進(jìn)行實(shí)驗(yàn)，結(jié)果如圖6所示。對于VTK算法，當(dāng)進(jìn)程數(shù)量在2個和8個之間時，算法執(zhí)行時間有所下降；但是，當(dāng)進(jìn)程數(shù)量超過8個時，算法執(zhí)行時間反而隨著進(jìn)程數(shù)量的增加而增加。本文算法則是隨著進(jìn)程數(shù)量的增加，執(zhí)行時間降低。出現(xiàn)該結(jié)果的原因是，VTK算法以塊數(shù)據(jù)為單位并行讀取數(shù)據(jù)，為了確保每個進(jìn)程上的網(wǎng)格數(shù)量大致相等，需要在數(shù)據(jù)讀取完畢后進(jìn)行一次同步操作；而本文算法使用預(yù)處理程序統(tǒng)計每個數(shù)據(jù)塊的網(wǎng)格數(shù)量，確保了每個進(jìn)程上讀取的網(wǎng)格數(shù)量大致相等，不需要同步。隨著進(jìn)程數(shù)量增加，VTK算法的同步操作需要更大的通信開銷，算法執(zhí)行時間因此增加；與VTK算法相比，本文算法少一次數(shù)據(jù)同步開銷，隨著進(jìn)程數(shù)量增加，局部KD樹構(gòu)造時間下降，算法執(zhí)行時間隨之下降，性能優(yōu)勢更加明顯。

圖6 數(shù)據(jù)分割算法性能對比

3.2 并行體繪制算法性能

為了研究本文的并行體繪制算法性能，首先，使用節(jié)點(diǎn)內(nèi)并行實(shí)驗(yàn)配置，測試算法在不同進(jìn)程數(shù)量下的算法執(zhí)行時間和加速比；然后，使用跨節(jié)點(diǎn)并行實(shí)驗(yàn)配置，測試算法在大規(guī)模網(wǎng)格體繪制下的性能，并結(jié)合實(shí)驗(yàn)數(shù)據(jù)分析了算法的效率和瓶頸。

3.2.1 算法性能

由于本文使用的實(shí)驗(yàn)數(shù)據(jù)均是曲面網(wǎng)格，經(jīng)過四面體化處理后，網(wǎng)格規(guī)模將會大幅度增加，分別為DATASET1（113 920 117）、DATASET2（46 412 858）、DATASET3（25 984 608）。所有的并行體繪制圖片分辨率均為1 000×800。表3給出了本文算法在不同進(jìn)程數(shù)量下的交互時間和繪制時間，結(jié)果顯示，隨著進(jìn)程數(shù)量增加，算法執(zhí)行時間縮短，使用64個進(jìn)程并行繪制時獲得最佳性能。由于使用了兩層LOD交互模型，交互開銷非常低，均為毫秒級別，并且基本不受進(jìn)程數(shù)量影響，可以滿足實(shí)時交互需求。其中，在64個進(jìn)程下，飛行器數(shù)據(jù)DATASET3周圍的流場體繪制結(jié)果如圖7所示。

圖7 飛行器DATASET3周圍的流場體繪制結(jié)果（64個進(jìn)程）

表3 并行體繪制算法性能

3.2.2 加速比

為了進(jìn)一步研究算法性能和瓶頸，統(tǒng)計了本文算法相對于串行算法的加速比。如圖8所示，當(dāng)進(jìn)程數(shù)量在8個以內(nèi)時，增加進(jìn)程數(shù)量，算法性能提升不明顯；當(dāng)進(jìn)程數(shù)量大于8個時，算法加速比開始明顯提升；當(dāng)進(jìn)程數(shù)量為64個時，獲得最大平均加速比，為37.6。這是因?yàn)?，進(jìn)程數(shù)量在8個以內(nèi)時，單個進(jìn)程需要處理的網(wǎng)格單元數(shù)量較大，繪制階段占據(jù)整個算法執(zhí)行時間的絕大部分；進(jìn)程數(shù)量大于8個后，當(dāng)進(jìn)程數(shù)量翻番時，單個進(jìn)程需要處理的網(wǎng)格數(shù)量減半，繪制階段占據(jù)整個算法執(zhí)行時間的比例明顯下降，加速比隨之明顯提升。

圖8 并行體繪制算法加速比

3.2.3 跨節(jié)點(diǎn)性能

由于單個計算節(jié)點(diǎn)性能有限，為了研究大規(guī)模網(wǎng)格體繪制性能，本節(jié)實(shí)驗(yàn)采用了32個計算節(jié)點(diǎn)。其中，使用的3個非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)規(guī)模為10億級別，分別為yA31（12.1億）、DATASET4（29.5億）、DATASET5（46.5億），所有的并行體繪制圖片分辨率均為1 000×800。表4給出了實(shí)驗(yàn)結(jié)果，其中，“-”表示進(jìn)程數(shù)量太少，不足以計算體繪制結(jié)果。根據(jù)實(shí)驗(yàn)結(jié)果可知，在不同進(jìn)程數(shù)量下，交互開銷基本一致，并且都在毫秒級別，滿足了實(shí)時交互需求。其中，行星yA31撞擊地球海面模擬數(shù)據(jù)的并行體繪制結(jié)果如圖9所示。

圖9 行星yA31撞擊地球海面模擬數(shù)據(jù)的并行體繪制結(jié)果

表4 并行體繪制算法性能（32個節(jié)點(diǎn)）

4 結(jié)束語

本文提出了一種基于sort-last架構(gòu)的非結(jié)構(gòu)網(wǎng)格并行體繪制算法。整個算法分為數(shù)據(jù)分割、繪制和圖片合成3個階段。在數(shù)據(jù)分割階段，設(shè)計了基于KD樹的并行體數(shù)據(jù)分割算法，首先，使用預(yù)處理程序確保每個進(jìn)程讀取的網(wǎng)格數(shù)量大致相等；然后，在每個進(jìn)程上同步構(gòu)造局部KD樹；最后，使用MPI API函數(shù)進(jìn)行數(shù)據(jù)交換并合成全局KD樹。此外，還根據(jù)相機(jī)位置和分割平面計算了全局KD樹葉子節(jié)點(diǎn)的可見性順序。在繪制階段，在每個進(jìn)程上建立可視化管線，使用高效的線程并行投影四面體算法計算體繪制圖片。在圖片合成階段，設(shè)計了基于二叉樹的圖像合成算法，根據(jù)可見性順序?qū)w繪制圖片分組合成。為了提升用戶的交互體驗(yàn)，本文還設(shè)計了兩層LOD模型，處于交互狀態(tài)時，渲染開銷較低的邊框模型；停止交互后，再渲染開銷較大的體繪制模型。

在Intel服務(wù)器和國產(chǎn)服務(wù)器上 測試了本文的數(shù)據(jù)分割算法和并行體繪制算法。對于數(shù)據(jù)分割算法，首先，研究了不同進(jìn)程數(shù)量下的網(wǎng)格規(guī)模增加情況；然后，測試了數(shù)據(jù)分割算法的負(fù)載均衡能力；最后，對比了本文算法與VTK算法的性能。實(shí)驗(yàn)結(jié)果顯示，本文算法在不同進(jìn)程數(shù)量下可保持良好的負(fù)載均衡能力；同時，與VTK算法相比，本文算法更具性能優(yōu)勢。對于并行體繪制算法，設(shè)計了節(jié)點(diǎn)內(nèi)并行實(shí)驗(yàn)和跨節(jié)點(diǎn)并行實(shí)驗(yàn)，首先，測試了在節(jié)點(diǎn)內(nèi)并行實(shí)驗(yàn)配置下算法的執(zhí)行時間，并計算了本文算法相對于串行體繪制算法的加速比；然后，測試了在32個計算節(jié)點(diǎn)并行實(shí)驗(yàn)配置下算法的執(zhí)行時間。實(shí)驗(yàn)結(jié)果顯示，本文算法在64個進(jìn)程下可以獲得的最大平均加速比為37.6，由于采用了兩層LOD模型，交互時間均控制在毫秒級別，本文算法可以很好地應(yīng)用于大規(guī)模非結(jié)構(gòu)網(wǎng)格體可視化，滿足實(shí)時交互的需求。

本文實(shí)現(xiàn)了基于sort-last架構(gòu)的并行體可視化原型系統(tǒng)，但是，還存在一些需要繼續(xù)優(yōu)化的地方。在數(shù)據(jù)分割算法中，并行計算分割平面的方向和位置是一個耗時的過程，下一階段的研究工作可以考慮保存這些KD樹分割平面信息，避免同一個數(shù)據(jù)被再次分割時重復(fù)計算分割平面。在圖片合成算法中，使用over運(yùn)算逐像素串行地合成兩張圖片，下一階段的研究工作可以考慮使用線程并行的方法優(yōu)化over運(yùn)算。在用戶交互體驗(yàn)優(yōu)化方面，由于體繪制模型渲染需要極大的開銷，下一階段的研究工作可以考慮過濾無效的鼠標(biāo)事件，避免造成無效的體繪制更新。