八叉树搜索-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); }