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

厲夫兵，蘇永琪，陳文劍

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

0 引 言

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

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

1 KD-Tree的創(chuàng)建

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

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

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

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

(1)

式中：centeri表示某一維度上的每個面元中心點坐標的值，mean表示某一維度上所有面元中心點坐標的平均值，F(xiàn)aceNumberInNode表示當前節(jié)點內三角面元的數(shù)量。

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

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

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

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

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

2 KD-Tree的遍歷

2.1 KD-Tree節(jié)點線索的生成

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

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

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

對比之下，傳統(tǒng)方法建立線索時6號節(jié)點的x-指向的是2號節(jié)點，改進后6號節(jié)點的x-指向的是9號節(jié)點，改進后的方法得到了更優(yōu)的線索指向。

2.2 KD-Tree的遍歷過程

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

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

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

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

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

3 射線追蹤并行計算

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

3.1 CUDA并行編程模型

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

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

圖3 CUDA線程層次結構

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

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

3.2 射線追蹤并行化過程

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

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

圖4 射線追蹤并行化實現(xiàn)方案

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

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

4 仿真結果與討論

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

為對算法有效性進行驗證，本文模擬射線在T型腔內進行彈射，并對其軌跡進行追蹤。腔體及發(fā)射源如圖5所示，圖5(a)為通過ANSYS Mechanical APDL 16.0對腔體表面進行三角網(wǎng)格剖分，剖分后三角面元總數(shù)為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的虛擬散點面，射線源點與散點面上的散點構成射線集合。射線可以通過出射面射出。實際應用中，射線一般是電磁波、聲波等，多次反射后會有能量衰減，因此設置射線在腔體內的彈射次數(shù)最多不超過30次。

圖5 腔體及發(fā)射源

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

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

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

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

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

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

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

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

5 結束語

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