一種基于sort-last架構(gòu)的非結(jié)構(gòu)網(wǎng)格并行體可視化算法

范良，張曉蓉，吳亞?wèn)|，陳呈，王昉

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

通常，一個(gè)大規(guī)模的網(wǎng)格數(shù)據(jù)在超級(jí)計(jì)算機(jī)上使用多個(gè)并行程序模擬計(jì)算生成，以多塊（multi-block）數(shù)據(jù)的形式寫(xiě)入磁盤(pán)，最好也用并行讀取的方法處理。雖然單塊數(shù)據(jù)的規(guī)模不大，但是單塊數(shù)據(jù)之間的網(wǎng)格規(guī)?？赡懿罹嗪艽蟆Ｒ虼?，本文首先使用一個(gè)預(yù)處理程序統(tǒng)計(jì)大規(guī)模網(wǎng)格數(shù)據(jù)中每一塊數(shù)據(jù)的網(wǎng)格單元數(shù)量，并保存到一個(gè)列表中；然后，根據(jù)列表網(wǎng)格數(shù)量信息計(jì)算每個(gè)進(jìn)程應(yīng)讀取的單塊數(shù)據(jù)編號(hào)，確保每個(gè)進(jìn)程讀取的網(wǎng)格單元數(shù)量大致相等；最后，每個(gè)進(jìn)程根據(jù)單塊數(shù)據(jù)的編號(hào)并行讀取體數(shù)據(jù)。

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

2.1.2 全局KD樹(shù)合成

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

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

2.1.3 基于KD樹(shù)的可見(jiàn)性排序算法

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

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

算法2可見(jiàn)性順序計(jì)算

輸入：全局KD樹(shù)

輸出：可見(jiàn)性順序列表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)：

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

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

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

2.2 繪制

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

圖2 體繪制管線(xiàn)

2.3 圖片合成

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

圖3 圖片合成

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

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

圖4 二叉樹(shù)圖像合成算法過(guò)程（4個(gè)進(jìn)程）

2.4 兩層LOD模型

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

圖5 兩層LOD模型

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

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

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

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

3.1.1 數(shù)據(jù)分割開(kāi)銷(xiāo)

在構(gòu)造局部KD樹(shù)時(shí)，與分割平面相交的單元會(huì)被復(fù)制一份，并被劃分到分割平面的兩側(cè)。因此，數(shù)據(jù)分割后，網(wǎng)格單元總數(shù)會(huì)增加。表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ù)分割開(kāi)銷(xiāo)的重要部分，直接影響后續(xù)并行體繪制算法的效率。

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

3.1.2 負(fù)載均衡測(cè)試

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

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

3.1.3 VTK算法對(duì)比

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

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

3.2 并行體繪制算法性能

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

3.2.1 算法性能

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

圖7 飛行器DATASET3周?chē)牧鲌?chǎng)體繪制結(jié)果（64個(gè)進(jìn)程）

表3 并行體繪制算法性能

3.2.2 加速比

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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

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

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