基于線索KD-Tree的射線追蹤并行計算

厲夫兵，蘇永琪，陳文劍

(1.北京信息科技大學 信息與通信工程學院，北京 100101； 2.哈爾濱工程大學 水聲工程學院，黑龍江 哈爾濱 150001)

0 引 言

射線追蹤法[1]在通信與計算機領域具有很廣泛的應用。通信領域中，利用射線追蹤法可以分析目標的散射特性。計算機領域中，通過射線追蹤可以使二維平面圖像呈現更真實的三維效果，常用于圖像渲染?？臻g剖分技術和并行計算可以運用于射線追蹤計算過程中，解決射線追蹤效率低下的問題。1980年J. Bentley提出KD-Tree空間剖分技術，已被運用于射線追蹤靜態場景的計算中，旨在將射線追蹤計算過程的時間復雜度降低[2，3]。此外，隨著GPU架構技術的不斷發展，GPU逐漸展現出了強大的并行計算能力。NVIDIA公司提出的統一計算設備體系結構(CUDA)[4]作為開發環境，可以將大量數據進行合理分配后，交由GPU進行并行計算。射線追蹤過程中會面對大量射線數據，因此GPU并行計算可以有效提高射線追蹤計算效率。

本文針對傳統射線追蹤計算效率低下的問題，首先對目標場景進行三角剖分建模，其次采用KD-Tree數據結構對場景空間進行合理劃分，并且通過添加線索指針的方式實現KD-Tree的無堆棧遍歷，以此避免大部分無效的求交運算[5]，最后在CUDA平臺下設計射線追蹤在CPU-GPU異構計算的實現方案，合理分配計算資源，將射線追蹤過程并行化[6]。通過KD-Tree與并行計算結合，改善射線追蹤的計算時間。

1 KD-Tree的創建

KD-Tree是二叉空間剖分樹的一種，它是利用平面的定策性質，使用一個平面將目標模型的空間劃分為兩部分，之后再遞歸地對已經劃分好的空間進行分割，直到達到提前設定的條件[7]。

目前構建KD-Tree的一種方法是通過表面積啟發式(surface area heuristic，SAH)[8]選取合適的剖分面，遞歸將目標場景剖分為若干個軸對齊長方體包圍盒。但當包圍盒內三角面元數量較多時使用SAH算法計算代價值會出現誤差，因此本文結合目標空間每一維度的方差與維度上面元坐標的中心點來確定剖分的維度與剖分面。創建KD-Tree的具體步驟如下：

(1)基于仿真模型上的三角面元和節點信息構建包圍盒，作為KD-Tree的根節點。

(2)選取合適的分割軸與剖分平面，切割根包圍盒生成左右子節點。為選擇合適的分割軸，需要計算節點內各個面元的中心點坐標分別在x，y，z維上的方差。任一維度方差的計算公式如式(1)所示

(1)

式中：centeri表示某一維度上的每個面元中心點坐標的值，mean表示某一維度上所有面元中心點坐標的平均值，FaceNumberInNode表示當前節點內三角面元的數量。

在確定分割軸后，可以對分割維度上三角面元的中心點進行排序，選取中間值作為剖分面。

(3)對于新生成的節點重復步驟(2)的內容，直到節點內的三角面元數量達到提前設定的閾值，節點不可再分割。

自此，KD-Tree構建完成，剖分過程中所形成的節點是KD-Tree的中間節點，剖分結束后的不可再分割節點是KD-Tree的葉子節點。

腔體模型的剖分示例和分割模型生成的KD-Tree如圖1所示。根包圍盒為N0，面L0將N0切割為N1，N2兩個節點，面L1，L2分別將N1，N2切割為N3，N4，N5，N6這4個節點，為節省內存開銷，在將節點分割后將模型的面元信息與點信息都保存至生成的子節點中，因此最終模型的信息都會保存在葉子節點中[9]。由于模型剖分的密度不同，每次切割并不能保證左右子節點中的面元數量相同，因此每個葉子節點中的三角面元的數量會隨著剖分的不同而不同。對于完全處于包圍盒內的面元，如圖1(a)中的2號與4號面元，它們分別只屬于N4，N5節點；而對于屬于分割平面上的面元，例如1號面元既屬于N3節點又屬于N4節點，3號面元既屬于N5節點又屬于N6節點。

圖1 模型分割及二叉樹生成

2 KD-Tree的遍歷

2.1 KD-Tree節點線索的生成

為實現KD-Tree的無堆棧遍歷，本文通過線索指針的方式建立節點之間的索引。Popov等提出了一種KD-Tree的線索創建算法，該方法現常被用于KD-Tree的無堆棧遍歷[10]，文獻[11]根據其算法思想進行了一定程度的簡化與改進，提出了一個節點間線索建立算法?！熬€索”實際上是KD-Tree節點之間的指向關系。線索建立的過程分為“繼承”和“更新”，“繼承”是指節點首先要繼承其父節點的所有線索信息，“更新”是指更新節點分割軸方向上的指向關系。平面包圍盒剖分結果及射線過程如圖2(a)所示，下面通過圖2(a)介紹傳統與改進后建立線索的方式。以3號與6號節點為例，通過傳統方式建立線索時，3號節點首先繼承其父節點1號節點的線索信息，即x軸正向索引(以下簡稱x+)、x軸負向索引(以下簡稱x-)、y軸正向索引(以下簡稱y+)、y軸負向索引(以下簡稱y-)均為NULL，其次由于1號節點的分割軸為x軸，因此根據包圍盒的位置關系將3號節點的x-更新為2號節點，至此3號節點線索信息更新完畢；6號節點首先繼承其父節點3號節點的線索信息，即x-為2號節點，x+、y+、y-均為NULL；其次6號節點的父節點分割軸為y軸，所以將6號節點的y-更新為7號節點，6號節點的線索信息也更新完畢。

圖2 平面模型中射線追蹤路徑

傳統方式建立線索是在分割形成節點的同時就將節點線索隨之建立，當其相鄰節點再度進行分割后，分割信息無法及時反饋至先前節點中，因此無法達到最優的線索指向關系。改進后的方法優化步驟主要在“更新”的過程中，每個節點在建立線索時需要了解其相鄰節點的完整分割信息，因此建立線索的過程需要在建樹完成后進行。通過改進方法建立線索時，3號節點依舊要先“繼承”1號節點的所有線索關系，在“更新”過程中，需要判斷3號節點的兄弟節點(2號節點)的分割軸(x軸)與其父節點(1號節點)的分割軸(x軸)是否相同，若相同則將3號節點的x-更新為2號節點的右孩子節點(5號節點)，隨后繼續更新。更新完成的條件一是直到3號節點的x-為葉子節點，二是x-指向節點的分割軸與3號節點父節點(1號節點)的分割軸不再相同。由于5號節點的分割軸仍與1號節點相同，因此最終3號節點的x-更新為9號節點；6號節點直接繼承3號節點索引，x-為9號節點，其余方向索引均為NULL，在更新過程中，由于6號節點的父節點(3號節點)的分割軸為y軸，因此將6號節點y-更新為7號節點，更新結束。

對比之下，傳統方法建立線索時6號節點的x-指向的是2號節點，改進后6號節點的x-指向的是9號節點，改進后的方法得到了更優的線索指向。

2.2 KD-Tree的遍歷過程

射線追蹤的過程實際上是射線遍歷線索KD-Tree尋找葉子節點并與節點內三角面元進行相交檢測的過程。圖2(a)演示了一條射線在平面模型中反射的路徑示意：

(1)射線通過入射位置定位到6號葉子節點，再通過Moller-Trumbore方法與節點內部的三角面元遍歷求交，求出交點后更新射線為反射射線[12]；

(2)反射射線確定出射位置后查找線索確定下一個待遍歷節點為9號節點，遍歷節點內面元求交，繼續更新反射射線，并通過出射位置查找對應線索指向節點為8號節點，射線在8號節點內繼續遍歷；

(3)在8號節點內求得交點與反射射線后，射線進入4號節點中遍歷后并未與節點中的三角面元發生相交，并且出射位置指向的線索為NULL，該條射線射出模型，遍歷結束。

因此，單根射線遍歷過程就是射線與葉子節點內面元求交后，通過線索信息查找到下一個待遍歷的節點繼續求交，直到射線在某節點出射位置的線索指向NULL時，遍歷結束。射線在KD-Tree中的遍歷過程如圖2(b)所示，圖2(b)中虛線(實線)表示傳統(改進)線索建立方式下射線遍歷線索KD-Tree的過程，由圖可見，改進后明顯減少了部分中間節點的遍歷，節省了開銷。

3 射線追蹤并行計算

傳統串行計算中，未進行追蹤的單根射線需要等待在其之前的射線追蹤完畢后，方可執行計算，而每根射線在與三角面元求交時都會產生巨大的時間開銷，使得串行計算效率低下。單根射線在追蹤過程中都是相互獨立的，因此很適合采用GPU對射線追蹤計算過程進行并行加速，提高計算效率。

3.1 CUDA并行編程模型

CUDA是一種CPU與GPU相結合的異構運算平臺，它將CUDA C語言作為編程語言，并由CPU負責處理邏輯關系復雜的事務，GPU負責處理需要高度計算的事務。相應的，CPU一般被稱作主機端(Host)，GPU一般被稱作設備端(Device)[13]。在CUDA編程中，運行在GPU端的程序需要通過核函數(Kernel)“標識”出來。

核函數的啟動需要調度網格(Grid)、線程塊(Block)及線程(Thread)，這三者最高可組織為三維形式。三者二維類型的層次結構如圖3所示：一個Grid中包含多個Block，一個Block中包含多個Thread[14]。其中Thread是執行計算的單位，因此并行計算的過程實際上是每個Thread各自執行Kernel的過程。

圖3 CUDA線程層次結構

由于計算量的不同，在進行并行計算時需要選擇合適的計算資源，即選取合適數量的Thread進行計算，程序層面上，通過線程塊數量(BlockperGrid)指定執行計算的線程塊數量，通過塊內線程數量(ThreadperBlock)指定每個線程塊中執行計算的線程數量，分配的線程總數量通過線程塊數量與塊內線程數量相乘得到。本文為進行并行計算資源的合理分配，設置線程總數量與射線的數量保持接近，由于一般情況下BlockperGrid和ThreadperBlock都設置為32的倍數，為同時滿足這兩個條件，一些情況下射線數量會略大于線程總數量，這時需要少數線程執行兩根射線的計算過程。

CUDA架構可以調用GPU內存資源對數據進行處理，GPU內存包括全局內存、共享內存、常量內存、紋理內存等，各個存儲器的位置與特性不都相同。由于共享內存讀寫速度快，但空間較小，而常量內存與紋理內存只能進行讀操作，因此本文將數據全部存至設備端全局內存中，每個線程直接在全局內存中進行讀寫操作。

3.2 射線追蹤并行化過程

串行計算的全部過程在CPU上執行，每條射線依次遍歷線索KD-Tree計算其在場景中的傳播路徑。將串行計算程序改為并行計算程序一般有兩種方法，一是整個過程直接由GPU完成，二是將需要高度計算的部分由GPU完成，其它的部分仍由CPU完成[15]。

在射線追蹤過程中，每條射線在場景中進行彈射時是相互獨立的，因此將射線追蹤的計算過程進行并行化的核心思想是將單根射線分配至GPU的單個核心進行射線的路徑計算，對射線追蹤進行并行化實現的方案如圖4所示，并行化過程中，構建線索KD-Tree的任務邏輯比較復雜，因此由CPU完成；而射線遍歷線索KD-Tree的過程存在大量射線與節點內三角面元的求交計算，因此需要由GPU并行完成。

圖4 射線追蹤并行化實現方案

由于CPU與GPU無法互相直接讀取存儲在對方內存里的數據，因此需要進行數據傳輸以完成計算。需要傳輸的數據分為兩部分：一部分是射線集合，一部分是線索KD-Tree節點。射線集合是預設條件，在CPU中靜態存儲，因此可以通過cudaMemcpy函數直接由CPU端拷貝至GPU端；本文在第1節末描述了構成模型的三角面元都將分別存儲在各個葉子節點中，且剖分后每個葉子節點中面元數量并不相同，因此葉子節點指針需要動態分配內存，這造成KD-Tree無法直接通過cudaMemcpy函數進行傳輸，本文通過計算建樹完成后每個葉子節點的內存信息，并在GPU端開辟相同大小的內存將其逐一拷貝，以實現葉子節點內面元數據的傳輸。而對于根節點和中間節點，由于其內的面元數據已被釋放，因此仍按照靜態存儲方式進行拷貝。

數據傳輸完成后，由一個線程(Thread)負責一條射線的追蹤，多條射線并行進行追蹤，得到每條射線的彈射次數及每次反射時發生相交的三角面元編號，最終將計算結果回傳至CPU。

4 仿真結果與討論

本文實驗數據通過一臺計算機仿真得到，該計算機操作系統為Windows 10，CPU型號為Intel(R)Core(TM)i5-10210U CPU @ 1.60 GHz，GPU型號為NVIDIA GeForce MX350，編譯環境為Microsoft Visual Studio 2019(CPU端)，CUDA toolkit 11.4(GPU端)。

為對算法有效性進行驗證，本文模擬射線在T型腔內進行彈射，并對其軌跡進行追蹤。腔體及發射源如圖5所示，圖5(a)為通過ANSYS Mechanical APDL 16.0對腔體表面進行三角網格剖分，剖分后三角面元總數為9478，圖5(b)為腔體模型在y軸方向上的長度為4.5 m，x軸、z軸方向上的長度分別為通過1 m、0.5 m，腔體的3個開口面邊長均為0.5 m。x軸方向上的腔體開口0.05 m處設置0.5 m×0.5 m的虛擬散點面，射線源點與散點面上的散點構成射線集合。射線可以通過出射面射出。實際應用中，射線一般是電磁波、聲波等，多次反射后會有能量衰減，因此設置射線在腔體內的彈射次數最多不超過30次。

圖5 腔體及發射源

表1對比了在射線數目為1500，KD-Tree樹高為6，節點總數為63、其中葉子節點總數為32時，線索建立方式改進前后遍歷節點數量的變化及串行、并行計算所耗費的時間，由表1可以看出，改進后的線索建立方式比改進前遍歷節點數量減少了4832，減少了射線查找KD-Tree尋找葉子節點過程的開銷。從運行時間上看，無論是串行計算還是并行計算，改進后算法的程序執行效率比改進前都有一定程度的提高。

表1 不同線索建立方式計算時間

實際情況中，運用射線追蹤算法對目標進行特性分析時，為了得到精確的計算結果，往往需要發射大量射線進行追蹤，因此本文討論了射線數量成倍增長的情況下并行計算相比于串行計算的效率提升。當KD-Tree樹高為6，線程塊與塊內線程數同為128，射線數目變化的情況下，CPU、GPU端的計算時間見表2。從表2中的時間及加速比可以得知，射線數目數量級增長時，加速效果越來越好，當射線數目達到十萬數量級時，加速比達到20倍左右。實際場景中的射線數目往往更加巨大，因此，通過GPU并行計算得到的效率提升會越來越明顯。

表3給出了在射線數目為15 000，線程塊與塊內線程數同為128時，KD-Tree構建不同時的程序運行時間，從表中可以看出，對于串行計算來說，樹高越高程序運行速核度越快，這是由于CPU適合處理KD-Tree這類邏輯結構較強的事務，葉子節點數量越多，內部三角面元數量越少，射線遍歷節點的速度越快；對于并行計算來說，GPU核心處理復雜邏輯事務的能力相比CPU較差，分支過多時，運行時間部分花費在遍歷樹尋找節點的過程中，但分支過少時又會有大量的無效求交計算，因此在樹高太高、太低時都會對程序的執行效率造成影響。

表2 射線數目變化對運行時間的影響

表3 KD-Tree層高對運行時間的影響

函數參數配置對程序運行效率也有一定影響，合理的參數配置更容易使得資源得到合理利用。核函數參數配置主要體現在線程塊數量(BlockperGrid)以及塊內線程數量(ThreadperBlock)，表4是在射線數目為15 000、KD-Tree層高為6時，BlockperGrid和ThreadperBlock不同的情況下程序的運行時間。由表可知，當BlockperGrid和ThreadperBlock都設置為128時程序運行效果最好。其中ThreadperBlock對程序運行時間影響更大，這是由于線程是執行計算的單位，程序會將一些寄存器分配給每個線程，寄存器的讀寫速度比全局內存快，如果ThreadperBlock設置太大容易造成可分配寄存器數量不足，線程將從全局內存讀取數據，全局內存的訪問延遲高于寄存器因此會造成運行時間增多；而ThreadperBlock設置太小又會造成一些寄存器資源浪費，同樣會對運行時間造成影響。

表4 核函數參數配置對運行時間的影響

5 結束語

本文首先介紹了一種KD-Tree線索建立方式，并進行了有效性驗證，其次針對傳統的基于線索KD-Tree的射線追蹤串行計算進行并行化實現，并對改變前后的運行時間進行了對比，實驗數據表明并行算法使得計算效率得到提高，計算時間大大減少，此外還探討了KD-Tree的構建質量以及核函數參數的設置對于程序執行效率的影響。在以后的研究中，還需要考慮以下問題，GPU全局內存的讀寫速度相比于其它內存較慢，可以將模型數據在不同內存區進行合理的規劃分配，進一步提高計算效率。