翻譯版權所有，轉載請註明。

​

八叉樹由於劃分規則且與座標軸平行，是很常見的空間劃分技術。八叉樹是射線追蹤加速普遍採用的方法，射線穿過八叉樹時，可只針對八叉樹包含的節點進行求交。與簡單的

O(MxN)演算法相比，這種射線-物體求交測試次數大大減少。理論上說，這是加速射線追蹤非常有效的方法。但實踐中，對於有大量單元的複雜場景八叉樹顯然有用，但對於較簡單的場景，它開銷過大。

       

八叉樹遍歷的方法可以分為自底向上和自頂向下兩類。自底向上只針對葉子節點，而自頂向下方法從根節點開始遍歷至葉節點。一般來說，自底向上一類的演算法主要是3D DDA演算法（實際上3維中的情形是Bresenham畫線演算法），原因在於DDA在2維光柵影象和繪圖儀上是非常快的畫線演算法。主要的缺點在於它需要確定面鄰居，且一般要求八叉樹是固定深度的（這樣就能保證面的鄰居數目一定）。

3D DDA演算法建立時採用了陣列實現，這樣所有型別的鄰居（面、邊、頂點）隱含著是確定的。但樹結構中確定鄰居並不簡單，並且在動態場景中需要不斷重新計算鄰居。

自頂向下方法利用資料結構本身的特點使問題簡化了。主要的優點之一在於它不用考慮八叉樹多長時間重建一次，而且也不用考慮八叉樹是否固定深度。缺點在於遞迴造成的效能和記憶體開銷。實際上，八叉樹遞迴效能也很好，原因就在於其分叉數很大，也就是不需要很大的劃分深度就可以得到合理的體素解析度。實際的射線追蹤場景中自頂向下和自底向上劃分方法的速度差只在較深劃分的時候才會變大。這裡我介紹自頂向下方法，它和

Arvo在第一期Ray Tracing news中給出的演算法及Agate的HERO演算法緊密相關。

1

射線-節點相交

簡化起見，我們從二維情形開始（四叉樹），擴充套件到三維情形。

射線定義為向量對（o,d）。o為原點，d為歸一化方向向量。任何位於射線上的2d

點可用引數方程[x(t),y(t)]表示：

X(t) = o.x + t * dx, y(t) = o.y + t*dy,t>=0

對於四叉樹中的節點q，我們僅需知道其最小和最大頂點，因為它是軸對齊的長方形。左上角用(x0,y0),右下角用（x1,y1）表示。形式上，q表示為下列點的集合：

X0(q)<=x<=x1(q),y0(q)<=y<=y1(q)

函式x0,x1,y0,y1給出了四叉樹的範圍。實現中，我們僅儲存範圍，而不是計算出範圍。根據上述定義，若射線和節點相交當且僅當存在一個正的t值，使得：

x 0(q)<=x(t)<=x1(q),y0(q)<=y(t)<=y1(q)更一般的實現中，首先我們將射線延長至節點的邊界，即取t值使x=x0(q)，然後x=x1(q)，然後是y0和y1。設上述t為tx0，tx1等等，則對於節點q和射線r=(o,d)，有：

tx0(q,r) = (x0(q) – o.x)/d.x

tx1(q,r) = (x1(q) – o.x)/d.x

ty0(q,r) = (y0(q) – o.y)/d.y

ty1(q,r) = (y1(q) – o.y)/d.y

因此，若射線與節點相交，當且僅當存在某個t>=0，使得：

tx0(q,r)<=t<=t< ty1(q,r)。

換句話說，有[tx0(q,r), tx1(q,r)）和[ty0(q,r), ty1(q,r)）兩個區間，若兩個區間重疊，且有正的t值，則相交。上式可以進一步簡化為：

tmin = max(tx0(q,r),ty0(q,r))

tmax =min(tx1(q,r),ty1(q,r))

if(tmin，則相交。

實際應用中，為避免處理負數，我們可以直接將其用0代替。這樣如果存在一個區間，則必定是正的區間。如果tmin = tmax，則交點在反方向（如果已經將負值用0代替，即

tmin = tmax = 0），或者與節點相交於一點（tmin = tmax > 0）。如果存在一個穿過頂點的物體，穿過一個節點的頂點並非不可能。這是實際有可能發生的特殊情形，但極為罕見。因此不值得對其進行處理。實際上我們是這樣確定某個節點的交點的。如果節點不是葉節點，繼續向下遍歷檢索其子節點。儘管需要對4個子節點求交，但只需要對可能相交節點的四個子節點求交。

注意，如果射線與座標軸平行（僅有1個非零分量）或位於座標平面內（有兩個非零分量），可將對應軸的初始t值和終止t值儲存為0，這樣它就可以總是從[tmin, tmax]區間中剔除，除非正方向沒有交點。

2

目標

最終會得到一個射線相交的葉子節點列表。可以從根節點開始測試是否相交，如果相交，根節點是連結串列的頭節點。設根節點非葉節點，則可以測試根節點的四個子節點。這樣可以生成四個子連結串列，替換主連結串列中的根節點。這樣替換得到的連結串列中節點的子節點再次替換其各自的主節點。雖然可以採用非遞迴的方式實現演算法，但遞迴更容易理解。

與實際的實現更接近的是，判斷節點是否和射線相交，如果相交，對子節點以特殊的順序進行相交判斷。如果子節點為葉節點，將其新增至連結串列。最後所有與射線相交的葉節點被加入連結串列，然後可以對連結串列進行任意操作。這時無需再訪問八叉樹，因為連結串列包含了對節點的指標及相應的引數。

將該演算法從四叉樹擴充套件至八叉樹，只需增加一個z分量即可。因此對於八叉樹，需要計算

tz0(q,r)和tz1(q,r)等等。

3 深度優先搜尋

當沿八叉樹向下搜尋時，空間劃分的特性可以利用。一個節點的子節點包含於節點之內。新劃分採用中間截面，即子節點也同樣包含在範圍[x0,x1)，[y0,y1)，[z0,z1)內。此外，此外唯一需要考慮的三元數為(x1/2，y1/2，z1/2)。因此，無需計算節點及其子節點的所有t值（56個），只需要提前算好上述範圍的t值（12個），就可直接在後續判斷中使用。特別是，某個中間平面正好是兩個平面的中間位置，因此其t值也是對應t值的一半。因此我們無需顯式計算中間平面t值，而只需將父節點的t值求平均。也就是說我們只需顯式計算根節點的t0和t1，然後求平均就可以得到子節點的t值。

4 基本演算法的虛擬碼

 Base algorithm Pseudocode :

assume all negative t-values are stored as zero-valued...

void ray_step( octree *o, ray *r ) {     compute  tx0;     compute  tx1;     compute  ty0;     compute  ty1;     compute  tz0;     compute  tz1;

    proc_subtree(tx0,ty0,tz0,tx1,ty1,tz1, o->root); }

void proc_subtree( float tx0, float ty0, float tz0,                               float tx1, float ty1, float tz1,                               node *n ) {     if ( !(Max(tx0,ty0,tz0) < Min(tx1,ty1,tz1)) )         return;

    if (n->isLeaf)     {         addtoList(n);         return;     }

    float txM = 0.5 * (tx0 + tx1);     float tyM = 0.5 * (ty0 + ty1);     float tzM = 0.5 * (tz0 + tz1);

    // Note, this is based on the assumption that the children are ordered in a particular     // manner.  Different octree libraries will have to adjust.     proc_subtree(tx0,ty0,tz0,txM,tyM,tzM,n->child[0]);     proc_subtree(tx0,ty0,tzM,txM,tyM,tz1,n->child[1]);     proc_subtree(tx0,tyM,tz0,txM,ty1,tzM,n->child[2]);     proc_subtree(tx0,tyM,tzM,txM,ty1,tz1,n->child[3]);     proc_subtree(txM,ty0,tz0,tx1,tyM,tzM,n->child[4]);     proc_subtree(txM,ty0,tzM,tx1,tyM,tz1,n->child[5]);     proc_subtree(txM,tyM,tz0,tx1,ty1,tzM,n->child[6]);     proc_subtree(txM,txM,tzM,tx1,ty1,tz1,n->child[7]); } 

5 未排序連結串列

雖然上述演算法可以實現，但需要指出其存在一個效率較低的部分，即生成的連結串列一般都不是有序的。雖然連結串列包含所有訪問過的葉節點，但該連結串列並不是從最近到最遠進行排列的（也就是t值從最小到最大排列）。

有幾種方法可以解決這個問題。最簡單的辦法是在生成連結串列後對其按照t值進行排序（因此可以生成一個自射線原點的按照t值排序的有序連結串列）。這要求在連結串列結構中包含一個t值，可以用tmin或者tmax，用哪一個關係不大。最複雜的是修改addtoList函式，使它將所選的節點根據合適的t值插入連結串列中。實際上我們不需要在連結串列節點中加入t值，只需要在proc_subtree函式中求得t值時呼叫addtoList函式。問題是這涉及到對每個節點都要進行連結串列搜尋和插入操作，當射線與越來越多的節點相交時，會產生很大的開銷。更復雜的方法是對子節點的相交測試進行排序，使得到的連結串列是有序的。這需要求取射線的進入和退出平面，並根據此更改搜尋順序。雖然這個很可能實現，但他會帶來很大的難以預測的分支判斷開銷。在樹的平均劃分深度比較大的場景中，很可能連結串列的後排序會比預排序開銷更大（射線相交的葉節點會非常多）。在平均劃分較少的場合，連結串列會足夠短，這樣後排序的開銷會較小。

6 確定射線進入和退出表面

只要找到了入射平面，剩下的問題就是求解與入射平面中4個子節點哪一個相交。但是這裡有個例外，因為前面假設了所有的射線分量都是正的，因此子節點7不可能是第一個入射的節點。為確定與哪個節點首先相交，需要利用中點平面。其實很簡單，已知tmin代表射線的入射點，因此只需判斷中點平面在tmin之前還是之後被穿過。

如果入射XY平面，則首先相交的子節點是0-2-4-6.

如果入射YZ平面，則首先相交的子節點是0-1-2-3.

如果入射XZ平面，則首先相交的子節點是0-1-4-5.

這裡的便利之處在於上述四種可能性只與兩個二進位制位有關。XY平面情形下節點編號與1-2位有關，YZ情形下節點編號與0-1位有關，XZ情形下節點編號與0-2位有關。並且，三種情形都有與節點0相交的可能。這樣，每種情形下判斷分支數目自然減少了，原因是可以將節點0作為初始值，然後採用邏輯or設定正確的結果。

Entry Plane Conditions Bit to Set

XY txM < tz0  tyM < tz0 2  1

YZ tyM < tx0  tzM < tx0 1  0

XZ txM < ty0  tzM < ty0 2  0

注意兩個條件都需要計算檢驗，如果都成立，則需要將對應位執行邏輯or操作。

 

確定後續節點

 

確定初始相交的子節點後，還需要確定後續的相交子節點。此處採用順序確定的方法，而不是一下建立整個子節點清單然後進行處理。如果已知當前位於哪個節點，以及該節點的出射平面，則其後續相交節點顯而易見。注意子節點的出射t值等於是中點平面的t值，但應將其作為子節點的t1值。如此可以確定子節點的出射平面，且已知當前位於哪個子節點，可以確定下一個鄰居節點（鄰居節點在離開整個根節點時為空）。

Current Sub-node Exit Plane  YZ Exit Plane  XZ Exit Plane  XY

0 4 2 1

1 5 3 Exit

2 6 Exit 3

3 7 Exit Exit

4 Exit 6 5

5 Exit 7 Exit

6 Exit Exit 7

7 Exit Exit Exit

簡而言之，在proc_subtree函式中，不需要對所有子節點遞迴呼叫proc_subtree函式，只需設定一個引數記錄當前所在的子節點。該currentChild的初始值為前述首次相交子節點編號。然後我們遞迴呼叫currentChild，順序確定下一個相交節點，對後者進行處理。這個迴圈直到currentChild的值為Exit時結束。

此後我們就可以進行精確的交點遍歷，但對於映象的射線來說編號是錯誤的。比如，射線的x方向分量為負，射線順序經過子節點4-6-2。但映象的射線會被認為順序經過子節點0-2-6。需要做的是將序號進行變化，令射線認為其經過0-2-6，但實際上訪問子節點4-6-2。為此，需要建立一個變換函式對序號進行變換。比如，若ray->d.x為負，節點的順序不是01234567，而是45670123.若ray->d.y為負，則為23016745.若d.x和d.y均為負，則為67452301.

幸運的是，此處還有竅門可供利用。可以看出每個座標軸有對應的位。這說明，若射線方向與座標軸方向相反，只需要將其對應位翻轉。將對應位翻轉，只需要使用xor算符。由此可得變換函式f(i)為：

f(i) = i^a

a=4sx + 2sy + sz，其中s#當某座標軸的射線方向為負時=1，為正時=0.

請注意該變換是一種欺騙演算法，使程式認為射線的各方向分量為正。我們並沒有改變順序查詢節點函式中子節點的引數。We simply index a different node andpretend that it is the same one.

比如，在原偽碼處理節點2時，呼叫proc_subtree(tx0,tym,tz0,txm,ty1,tzm,n->child[2])。在改進版中，呼叫proc_subtree(tx0,tym,tz0,txm,ty1,tzm,n->child[2^a])，此處假設a值已確定。注意演算法永遠不會訪問這些節點，它只是將其新增到相交搜尋列表。因此才能騙過程式，因為程式永遠不會真的檢查child[2]是否正確。

上述採用的子節點順序的後果是需要將你的八叉樹編號調整至與本文一致。目前這種順序編碼方法的便利性被我們充分採用，但還有其他幾種排列方式也有相同的性質（儘管位對應順序不同）。

Pseudocode with ordered searching :

byte a;

// In practice, it may be worth passing in the ray by value or passing in a copy of the ray // because of the fact the ray_step() function is destructive to the ray data. void ray_step( octree *o, ray *r ) {     a = 0;     if (r->d.x < 0)     {         r->o.x = o->size - r->o.x;         r->d.x = -(r->d.x);         a |= 4;     }     if (r->d.y < 0)     {         r->o.y = o->size - r->o.y;         r->d.y = -(r->d.y);         a |= 2;     }     if (r->d.z < 0)     {         r->o.z = o->size - r->o.z;         r->d.z = -(r->d.z);         a |= 1;     }

    compute  tx0;     compute  tx1;     compute  ty0;     compute  ty1;     compute  tz0;     compute  tz1;

    float tmin = Max(tx0,ty0,tz0);     float tmax = Min(tx1,ty1,tz1);

    if ( (tmin < tmax) && (tmax > 0.0f) )         proc_subtree(tx0,ty0,tz0,tx1,ty1,tz1, o->root); }

void proc_subtree( float tx0, float ty0, float tz0,                               float tx1, float ty1, float tz1,                               node *n ) {     int currNode;

    if ( (tx1 <= 0.0f ) || (ty1 <= 0.0f) || (tz1< = 0.0f) )         return;

    if (n->isLeaf)     {         addtoList(n);         return;     }

    float txM = 0.5 * (tx0 + tx1);     float tyM = 0.5 * (ty0 + ty1);     float tzM = 0.5 * (tz0 + tz1);

    // Determining the first node requires knowing which of the t0's is the largest...     // as well as comparing the tM's of the other axes against that largest t0.     // Hence, the function should only require the 3 t0-values and the 3 tM-values.     currNode = find_firstNode(tx0,ty0,tz0,txM,tyM,tzM);

    do {         // next_Node() takes the t1 values for a child (which may or may not have tM's of the parent)         // and determines the next node.  Rather than passing in the currNode value, we pass in possible values         // for the next node.  A value of 8 refers to an exit from the parent.         // While having more parameters does use more stack bandwidth, it allows for a smaller function         // with fewer branches and less redundant code.  The possibilities for the next node are passed in         // the same respective order as the t-values.  Hence if the first parameter is found as the greatest, the         // fourth parameter will be the return value.  If the 2nd parameter is the greatest, the 5th will be returned, etc.         switch(currNode) {         

case 0 : proc_subtree(tx0,ty0,tz0,txM,tyM,tzM,n->child[a]);                     

currNode = next_Node(txM,tyM,tzM,4,2,1);                     break; 

        case 1 : proc_subtree(tx0,ty0,tzM,txM,tyM,tz1,n->child[1^a]);                     currNode = next_Node(txM,tyM,tz1,5,3,8);                     break; 

        case 2 : proc_subtree(tx0,tyM,tz0,txM,ty1,tzM,n->child[2^a]);                     currNode = next_Node(txM,ty1,tzM,6,8,3);                     break;     

    case 3 : proc_subtree(tx0,tyM,tzM,txM,ty1,tz1,n->child[3^a]);                     currNode = next_Node(txM,ty1,tz1,7,8,8);                     break;       

  case 4 : proc_subtree(txM,ty0,tz0,tx1,tyM,tzM,n->child[4^a]);                     currNode = next_Node(tx1,tyM,tzM,8,6,5);                     break;     

    case 5 : proc_subtree(txM,ty0,tzM,tx1,tyM,tz1,n->child[5^a]);                     currNode = next_Node(tx1,tyM,tz1,8,7,8);                     break;      

   case 6 : proc_subtree(txM,tyM,tz0,tx1,ty1,tzM,n->child[6^a]);                     currNode = next_Node(tx1,ty1,tzM,8,8,7);                     break;       

  case 7 : proc_subtree(txM,txM,tzM,tx1,ty1,tz1,n->child[7]);               

      currNode = 8;                     break;         

}   

  } while (currNode < 8); }