八叉树搜索-1
翻譯版權所有,轉載請注明。
?
八叉樹由于劃分規則且與坐標軸平行,是很常見的空間劃分技術。八叉樹是射線追蹤加速普遍采用的方法,射線穿過八叉樹時,可只針對八叉樹包含的節點進行求交。與簡單的
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); }?
總結
- 上一篇: 中国PostgreSQL培训认证推出“浪
- 下一篇: 装备升级