一種基于sort-last架構的非結構網格并行體可視化算法

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

1. 西南科技大學計算機科學與技術學院，四川 綿陽 621010；2. 四川輕化工大學計算機科學與工程學院，四川 自貢 643002；3. 空氣動力學國家重點實驗室，四川 綿陽 621000；4. 中國空氣動力研究與發展中心計算空氣動力研究所，四川 綿陽 621000

1 引言

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

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

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

● 實現了基于KD樹的并行體數據分割算法和可見性排序算法。首先，使用預處理程序確保每個進程讀取的網格數量大致相等；其次，以每個進程上的數據為輸入，構造局部KD樹；然后，將所有進程上的局部KD樹合成為全局KD樹；最后，通過相機位置和分割平面計算全局KD樹葉子節點的可見性順序。

● 實現了基于二叉樹的并行圖像合成算法。根據可見性順序對所有進程上的圖片兩兩配對并合成中間結果，不斷重復配對與合成的過程，直到合成最終結果。

● 實現了兩層LOD交互模型。首先，同步所有進程的數據邊界信息；然后，在主進程上渲染數據邊框模型。當處于交互狀態時，渲染邊框模型；停止交互時，渲染體繪制模型。

● 設計實驗驗證，并分析了算法性能。實驗結果顯示，本文的并行體數據分割算法具有良好的性能和負載均衡能力；并行體繪制算法可以被很好地應用于大規模非結構網格數據體可視化，且交互時延在毫秒級別，滿足實時交互的需求。

2 基于進程并行的非結構網格體繪制算法

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

圖1 基于進程并行的非結構網格體繪制算法流程

2.1 數據分割

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

2.1.1 局部KD樹構造

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

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

算法1KD樹構造

輸入：非結構網格數據

輸出：KD樹根節點指針

步驟1：計算所有網格單元的質心坐標；

步驟2：根據網格數據邊界值計算分割平面方向；

步驟3：計算質心坐標的中位數，將中位數質心坐標作為分割面的位置坐標；

步驟4：沿著分割面將網格數據一分為二；

步驟5：遞歸執行步驟2～4，直到生成指定數目的葉子節點。

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

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

2.1.2 全局KD樹合成

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

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

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

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

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

算法2可見性順序計算

輸入：全局KD樹

輸出：可見性順序列表list

步驟1：如果當前節點為葉子節點，將該節點的進程秩添加到list；

步驟2：如果當前節點為非葉子節點，根據相機位置確定近點與遠點：

· 計算近點可見性順序；

· 計算遠點可見性順序；

步驟3：遞歸執行步驟1～2，直到所有葉子節點的可見性順序計算完畢。

2.2 繪制

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

圖2 體繪制管線

2.3 圖片合成

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

圖3 圖片合成

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

以4個進程為例，二叉樹圖像合成算法流程如圖4所示。假設4個進程上的體繪制圖片可見性順序為O0＜O2＜O1＜O3，第一次合成，2號進程和3號進程分別將圖片發送至0號進程和1號進程，在0號進程和1號進程上合成中間結果；第二次合成，1號進程將中間結果發送至0號進程，在0號進程上合成最終結果。

圖4 二叉樹圖像合成算法過程（4個進程）

2.4 兩層LOD模型

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

圖5 兩層LOD模型

3 實驗分析

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

3.1 數據分割算法性能

為了研究數據分割算法性能，采用節點內并行實驗配置，從數據分割開銷、負載均衡測試和VTK（visualization toolkit）算法對比3個方面設計實驗。首先，在不同進程數量下研究數據分割后的網格數量增加情況；然后，對數據分割算法的負載均衡能力進行測試；最后，對比了本文算法與VTK算法的性能。

3.1.1 數據分割開銷

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

表1 網格規模

3.1.2 負載均衡測試

一個好的數據分割算法應該具有良好的負載均衡能力。本節將最多網格進程與最少網格進程的網格單元數量差值作為衡量標準，該值越小，說明算法負載均衡能力越好。表2給出了不同進程數量下的網格單元數量差值，結果顯示，所有的差值均在4 000以內，全局KD樹葉子節點之間的網格單元數量大致相等，本文數據分割算法的負載均衡能力較好。

表2 負載均衡測試結果

3.1.3 VTK算法對比

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

圖6 數據分割算法性能對比

3.2 并行體繪制算法性能

為了研究本文的并行體繪制算法性能，首先，使用節點內并行實驗配置，測試算法在不同進程數量下的算法執行時間和加速比；然后，使用跨節點并行實驗配置，測試算法在大規模網格體繪制下的性能，并結合實驗數據分析了算法的效率和瓶頸。

3.2.1 算法性能

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

圖7 飛行器DATASET3周圍的流場體繪制結果（64個進程）

表3 并行體繪制算法性能

3.2.2 加速比

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

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

3.2.3 跨節點性能

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

圖9 行星yA31撞擊地球海面模擬數據的并行體繪制結果

表4 并行體繪制算法性能（32個節點）

4 結束語

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

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

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