PBRT阅读：第四章 体素和求交加速 第4.4节

4.4 KD-树加速器

二元空间分割 (Binary space partitioning, BSP) 树自适应地把空间分割成大小不一的区域。跟均匀网格相比，BSP对几何体分布不均匀的场景而言是更有效率的数据结构。BSP树的创建起始于一个包含整个场景的包围盒。如果盒子里体素的个数超过某个阀值，就要用一个平面把包围盒一分两半。体素跟它所重叠的那半空间联系在一起，如果体素跟两个空间都重叠，则跟两个空间都要联系起来。这个分割过程要递归地进行下去，直到每个叶子区域含有足够少的体素，或者递归的深度已经达到给定的最大值。因为分割平面可以在任意位置进行分割，并且可以对三维空间的任意部分做任意程度的分割，所以BSP树可以很容易地处理不均匀分布的几何体。

有两种常用的BSP树：kd-树和八叉树。kd-树要求分割平面必须垂直于某个坐标轴，这使得树的创建和遍历都很有效，但也牺牲了某些关于分割的灵活性。八叉树用三个和轴垂直的平面把空间分割成八个区域（通常在空间的中心处分割）。在这一节里，我们介绍KdTreeAccel类，它实现了光线求交进行加速的kd-树。

<KdTreeAccel Declarations> =
    class KdTreeAccel : public Aggregate {
    public:
        <KdTreeAccel Public Methods>
    private:
        <KdTreeAccel Private Data>
    };

KdTreeAccel的构造器的参数包括所要存放的体素，还有一些控制树的创建的一些参数。这些参数被存放在成员变量中以备后用。为了简单起见，KdTreeAccel要求所有体素是可求交的。因此，在创建树之前，构造器把所有不可求交的体素转化为可求交的。

<KdTreeAccel Method Definitions> =
    KdTreeAccel:: KdTreeAccel(const vector<Reference<Primitive>> &p,
                    int icost, int tcost,
                    float ebonus, int maxp, int maxDepth)
        : isectCost(icost), traversalCost(tcost), maxPrims(maxp), emptyBonus(ebonus) {
        vector< Reference<Primitive>> prims;
        for(u_int i = 0; i < p.size(); ++i)
            p->FullyRefine(prims);
        <Initialize mailboxes for KdTreeAccel>
        <Build kd-tree for accelerator>
    }

<KdTreeAccel Private Data> =
    int isectCost, traversalCost, maxPrims;
    float emptyBonus;

跟GridAccel 一样，kd-树也使用了信箱技术来避免重复的求交计算。实际上，它用到了同样的MailboxPrim结构。
<Initialize mailboxes for KdTreeAccel> =
    curMailBoxId = 0;
    nMailboxes = prims.size();
    mailboxPrims = (MailboxPrim *)AllocAligned(nMailboxes * sizeof(MailboxPrim));
    for(u_int i = 0; i < nMailboxes; ++i)
        new (&mailboxPrims)MailboxPrim(prims);

<KdTreeAccel Private Data> +=
    u_int nMailboxes;
    MailboxPrim *mailboxPrims;
    mutable int curMailboxId;

4.4.1 树的表示

kd-树是一种二叉树，它的每个内部结点都有两个子结点，每个叶子结点存放所重叠的体素。每个内部结点必须存放下列三类信息：
·    分割轴：该结点用哪一个轴（x,y或z轴）来分割的。
·    分割位置： 在分割轴上的分割平面的位置。
·    子结点：用来找到两个子结点的信息。

每个叶子结点只存放跟该区域重叠的体素。

我们要做些稍微复杂的工作，使得所有的内部结点和部分的叶子结点都用8个字节（假定浮点数和指针都是4个字节），这样就可以保证4个结点刚好占满32字节的缓存线(cache line)。由于树上有大量的结点，每条光线要用到大量的结点，减少结点的内存使用可以极大地提高缓存效率。我们最初用16个字节存放结点，改为8个字节后，有20%的速度提升。叶子结点和内部结点都用下面的KdAccelNode结构。注意每个union成员的注释说明了这个成员是用在内部结点，还是叶子结点，还是两者都用。

<KdAccelNode Declarations> =
    struct KdAccelNode {
        <KdAccelNode Methods>
        union {
            u_int flags;     // Both
            float split;    //Integer
            u_int nPrims;    // Leaf
        };
        union {
            u_int aboveChild;        //Interior
            MailboxPrim *onePrimitive;    //Leaf
            MailboxPrim **primitives;     //Leaf
        };
    };

KdAccelNode::flags的最低两位用来区别x、y 或z分割的内部结点(分别用值0,1,2)和叶子结点(用值3)。

相对来说，把叶子结点放入8个字节的内存比较容易：因为KdAccelNode::flags的低两位用于表示它是叶子结点， KdAccelNode::nPrims的高30位用来表示有多少体素跟它重叠。和GridAccel一样，如果只有一个体素跟它重叠，那么它的MailboxPrim指针就直接存放在KdAccelNode::onePrimitive。如果有多个体素重叠，则要动态地申请存放这些指针的数组，而KdAccelNode::primitives指向这个数组。

很容易对叶子结点初始化：在把体素的个数存放到结点之前，要把它左移两位，并把KdAccelNode:flags的低两位设置为3，用来说明这是个叶子结点。

<KdAccelNode Methods> =
    void initLeaf(int *primNums, int np, MailboxPrim *mailboxPrims,
            MemoryArea &arena) {
        nPrims = np <<2;
        flags |= 3;
        <Store MailboxPrim *s for leaf node>
    }

因为我们有了KdAccelNode::onePrimitive，对于那些有0个或1个体素的叶子结点，没有必要申请内存。如果有多个体素重叠，调用者要传入一个MemoryArena（见第A.2.4节），用来申请MailboxPrim的指针数组。MemoryArena类可以帮助减少对空间的浪费，并把这些数组放在一起，以提高缓存效率。

<Store MailboxPrim *s for leaf node> =
    if (np == 0)
        onePrimitive = NULL;
    else if (np == 1)
        onePrimitive = &mailboxPrims[primNums[0]);
    else {
        primitives = (MailboxPrim **)arena.Alloc(np * sizeof(MailboxPrim *));
        for(int i = 0; i < np; i++)
            primitive = &mailboxPrims[primNums];
    }

用8个字节表示内部结点需要多一些工作。前面解释过，KdAccelNode::flags的最低2位记录了分割轴。然而，分割位置KdAccelNode::split作为一个浮点数，也跟KdAccelNode::flags一起共用同一个内存地址。这似乎不可能做到--我们不能告诉编译器只用KdAccelNode::split的高30位当作一个浮点数使用。

但是，只要我们在设置KdAccelNode::split之后再设置KdAccelNode::flag的低2位，就可以达到上述的效果。这个技术得利于IEEE浮点数的内存布局：KdAccelNode::flag所用到的低2位只是占了浮点数的两个最低有效位，对浮点数值的影响很小。

虽然这个技巧相当出格儿，但为了用8字节存储树结点而得到的性能提升，还是很值得的。另外，我们用了几个KdAccelNode函数来隐藏所有这些复杂性，其它部分的实现并不受其影响。

我们不需要额外的内存存放内部结点指向两个子结点的指针，所有结点被存放在一个连续的内存区域中，代表分割面下面的区域的那个子结点被放置在紧靠其父结点的位置（这也提高了缓存效率，因为至少有一个子结点跟其父结点相邻）。而另一个代表分割面上面区域的那个子结点，将被放在数组的另一个位置，KdAccelNode::aboveChild指向这个位置。

有了上述的约定以后，我们就可以初始化内部结点了。在把分割轴写入KdAccelNode::flags之前，我们要设置分割位置（否则就会覆盖掉分割位置了）。

<KdAccelNode Methods> +=
    void initInterior(int axis, float s) {
        split = s;
        flags &= ~3;
        flags |= axis;
    }

最后，我们提供一些辅助函数来获得这些值，同时也隐藏了实现的复杂细节。

<KdAccelNode Methods> +=
    float SplitPos() const { return split;}
    int nPrimitives() const { return nPrims >> 2;}
    int SplitAxis() const { return flags &3;}
    bool IsLeaf() const { return (flags & 3) == 3; }

4.4.2 树的创建

kd-树是用自顶向下的递归算法创建的。每一步开始时，有一个轴对齐的空间区域和一组跟它重叠的体素。这个区域要么被分割成两个子区域并变成一个内部结点，要么创建一个叶子结点，结束递归过程。

在讨论KdAccelNode的时候，我们提到所有树结点都被存放在一个内存连续的数组中。KdTreeAccel::nextFreeNode记录了这个数组中下一个可用的结点。KdTreeAccel::nAllocedNodes记录了已经分配了内存的结点个数。我们把它们初始化为0，这里的实现保证当初始化第一个树结点时，立即分配内存。

如果用户没有提供给构造器树的最大深度，我们需要确定一个缺省值。虽然构造树的递归过程在正常情况下可以在某个深度上自然结束，但是设置最大深度值仍然很重要，这可以防止在有些病态情况下耗费大量的内存。我们发现对大多数的场景来讲，8+1.3log(N)是一个很合理的最大深度。

<Build kd-tree for accelerator> =
    nextFreeNode = nAllocedNodes = 0;
    if(maxDepth <= 0)
        maxDepth = Round2Int( 8 + 1.3f * log2Int(float(prims.size())));
    <Compute bounds for kd-tree construction>
    <Allocate working memory for kd-tree construction>
    <Initialize primNums for kd-tree construction>
    <Start recursive construction of kd-tree>
    <Free working memory for kd-tree construction>

<KdTreeAccel Private Data> +=
    KdAccelNode *nodes;
    int nAllocedNodes, nextFreeNode;

因为构造例程要重复使用体素的包围盒，所以在构造树之前，我们把它们放在一个vector中，这样就省去了对那些有可能很慢的Primitive::WorldBound()的重复调用。

<Compute bounds for kd-tree construction> =
    vector<BBox> primBounds;
    primBounds.reserve(prims.size());
    for (u_int i = 0; i < prims.size(); ++i) {
        BBox b = prim->WorldBound();
        bounds = Union(bounds, b);
        primBounds.push_back(b);
    }

<KdTreeAccel Private Data> +=
    BBox bounds;

树构造器的参数之一是一个体素索引值数组，用来指定那些和当前结点重叠的体素。因为所有体素都和根结点重叠，所以在递归过程开始之前，我们用0到prims.size()-1来初始化一个数组，并用此数组作为递归的开始。

<Initialize primNums for kd-tree construction> =
    int *primNums = new int[prims.size()];
    for (u_int i = 0; i < prims.size(); ++i) 
        primNums = I;

每个结点都要调用KdTreeAccel::buildTree()函数，这个函数负责确定该结点是内部结点还是叶子结点，并修改相关的数据结构。
最后三个参数edges,prims0和prims1，是指向片段<Allocate working memory for kd-tree construction>中数据的指针，这个片段在后面几页有介绍。

<Start recursive construction of kd-tree> =
    buildTree(0, bounds, primBounds, primNums,
            prims.size(), maxDepth, edges,
            prims0,prims1);

KdTreeAccel::buildTree()函数的主要参数是结点的数组偏置值nodeNum；给出结点覆盖区域的包围盒nodeBounds；跟它重叠的体素的索引值数组primNums。其它的参数会在以后加以描述。

<KdTreeAccel Method Definitions> +=
    void KdTreeAccel::buildTree(int nodeNum, const BBox &nodeBounds,
            const vector<BBox> &allPrimBounds, int *primNum,
            int nPrims, int depth, BoundEdge *edges[3],
            int *prims0, int *prims1, int badRefines) {
        <Get next free node from nodes array>
        <Initialize leaf node if termination criteria met>
        <Initialize interior node and continue recursion>
    }

如果所分配的结点都用完了，就会申请分配两倍的结点空间，并把旧值拷贝到新内存区域。当KdTreeAccel::buildTree()被第一次调用时，KdTreeAccel::nAllocNodes是0，所以要申请第一块结点内存空间。

<Get next free node from nodes array> =
    if (nextFreeNode == nAllocedNodes) {
        int nAlloc = max(2 * nAllocedNodes, 512);
        KdAccelNode *n = (KdAccelNode *)AllocAligned(nAlloc*sizeof(KdAccelNode));
        if(nAllocedNodes > 0) {
            memcpy(n, nodes, nAllocedNodes * sizeof(KdAccelNode));
            FreeAligned(nodes);
        }
        nodes = n;
        nAllocedNodes = nAlloc;
    }
    ++nextFreeNode;

如果区域内体素的个数已经足够少了，或者已经达到深度最大值，就要停止递归并创建一个叶子结点。参数depth最初被设为树的最大深度值，并在每一层递归中减1。

<Initialize leaf node if termination criteria met> =
    if(nPrims <= maxPrims || depth == 0) {
        nodes[nodeNum].initLeaf(primNums, nPrims, mailboxPrims, arena);
        return;
    }

如前所述，KdAccelNode::initLeaf用一个MemoryArena为可变长度的体素数组申请内存空间。因为arena是一个成员变量，当KdTreeAccel对象被撤销时，它所分配的空间将被完全释放。

<KdTreeAccel Private Data> +=
    MemoryArena arena;

如果这是一个内部结点，就必须选择一个分割平面，并将体素按照这个平面分类，再进行下一步的递归。

<Initialize interior node and continue recursion> =
    <Choose split axis position for interior node>
    <Create leaf if no good splits were found>
    <Classify primitives with respect to split>
    <Recursively initialize children nodes>

我们选择分割平面的依据是一个成本估算模型，它可以估算出进行光线求交测试的成本，其中包括遍历树结点所花费的时间和光线-体素求交测试所花费的时间。其目的就是使总成本最小化；我们要实现一个使结点上的成本最小化的贪心算法。结点上的几个候选分割平面上的成本都要估算到，然后选择成本最小的那个平面。

这个成本估算模型的原理很简单：如果我们把当前结点当做叶子结点，那么任何穿过这个区域的光线跟其中所重叠的体素都要进行求交测试，其成本是：
    
其中N是区域中体素的个数，t_i(i)是光线和第i个体素求交所花费的时间。
如果我们分割当前结点，那么光线在结点上的成本是：
        
其中t_t是光线遍历结点并决定要穿过哪个子结点的时间，P_B和P_A是光线分别穿过分割平面上、下区域的概率，b_i和a_i分别是分割平面上、下区域内的体素的索引值。N_B和N_A分别是分割平面上、下区域内的体素的索引值。分割的选择直接影响了分割面两边的概率值和体素个数。

在我们的实现中，我们做个简化约定：假设所有体素的求交测试时间都相同，这个假设跟实际情况相差并不太远，任何由它引起的误差并不对加速器的性能构成太大的影响。另一种方法是在Primitive类中加上估算求交时间的函数，用来返回一个求交运算所需CPU周期的估算值。求交成本t_i和遍历成本t_t可以由用户来设置，他们的缺省值分别是80和1。实际上，这两者的比率决定了算法的运行结果。

最后一点：我们倾向于选择那些能够使其中一个子结点没有体素的分割，因为光线可以立即跳过这个子结点转向下一个结点，而不做任何求交运算。这样，对于未分割的区域和分割的区域，它们的成本分别是：

            t_i N
和
    t_t + (1 – b_e) (P_BN_Bt_i + P_AN_At_i)

其中b_e是一个“奖励”值，如果两个子区域非空，它的值是0；否则，它是处于0到1之间的值。

概率PB和PA可以用几何概率的理论求得。如果凸体(convext volume) A被另一个凸体B所包含，那么一条随机光线穿过B同时也穿过A的条件概率是它们表面积S_A和S_B之比：

        p(A | B) = S_A/S_B
如图，穿过区域B和C的概率分别是S_B/S_A和S_C/S_A。
    
有了计算成本估算模型中的概率的计算方法后，下一个问题是如何生成分割候选位置，并对于每个候选位置，如何有效地计算其成本。可以看出这个估算模型的最小成本所对应的分割位置，是和某个体素的包围盒的某一面是重合的，所以就没有必要在中间位置上开刀了。（为了说服你自己，估算一下在中间位置分割时的成本吧）。在这里，我们要考虑区域中所有包围盒的面对三个坐标轴（其中一个或多个轴的）所产生的分割位置。

为了降低检查所有这些候选位置的成本，需要很仔细地设计算法的构造。为了计算每个候选的成本，我们要扫描所有包围盒到每个轴上的投影，并记录下成本最小的那个位置。每个包围盒在轴上的投影是两条（垂直于轴的）边，它们都用一个BoundEdge结构的实例表示。这个结构记录了边在轴上的位置，并标明它是始边(START)还是尾(END)边，还有它所关联的体素。见图，图中a0,a1,b0,b1,b2,c1,c2都要用一个BoundEdge实例表示。


<KdAccelNode Declarations> +=
    struct BoundEdge {
        <BoundEdge Public Methods>
        float t;
        int primNum;
        enum { START, END} type;
    };

<BoundEdge Public Methods> =
    BoundEdge(float tt, int pn, bool starting) {
        t = tt;
        primNum = pm;
        type = starting? START : END;
    }

为了计算任何结点上的成本，最多只需要2 * prims.size()个BoundEdge。所有这些边的内存只需做一次性申请，然后可以被新创建的结点再利用。片段<Free working memory for kd-tree construction>是用来释放这个内存空间的，这里略去。

<Allocate working memory for kd-tree construction>=
    BoundEdge *edge[3];
    for(int i = 0; i < 3; ++i)
        edge = new BoundEdge[2 * prims.size()];

计算出该结点作为叶子结点的成本后，KdTreeAccel::buildTree()选择一个分割轴，并计算每个候选分割位置的成本。bestAxis和bestOffset分别记录目前已经找到的最佳轴和最佳分割位置。invTotalSA被初始化为结点的表面积，它的值要用于计算光线穿过候选子结点的概率。

<Choose split axis position for interior node> =
    int bestAxis = -1, bestOffset = -1;
    float bestCost = INFINITY;
    float oldCost = isectCost * float(nPrims);
    Vector d = nodeBounds.pMax – nodeBounds.pMin;
    float totalSA = 2.f * (d.x * d.y + d.x * d.z + d.y * d.z);
    float invTotalSA = 1.f/totalSA;
    <Choose which axis to split along>
    int retries = 0;
    retrySplit:
    <Initialize edges for axis>
    <Compute cost of all split for axis to find best>

这个算法首先试着选择空间范围最大的那个轴，这样有助于形成趋于正方体的子区域。当然这是个依靠直觉的选择。如果在这个轴上不能找到好的分割，还要继续选择其它轴。

<Choose which axis to split along> =
    int axis;
    if(d.x > d.y && d.x > d.z) axis = 0;
    else axis = (d.y > d.z) > ? 1 : 2;

首先，我们用所有跟结点重叠的体素的包围盒来初始化关于该轴的BoundEdge数组。然后对这个数组从低到高进行排序，这样我们就可以从头到尾对这些边扫描，来寻找最佳选择。

<Initialize edges for axis> =
    for(int i = 0; i < nPrims; ++i) {
        int pn = primNums;
        const BBox &bbox = allPrimBounds[pn];
        edges[axis][2*i] = BoundEdge(bbox.pMin[axis], pn, true);
        edges[axis][2*i+1] = BoundEdge(bbox.pMax[axis], pn, false);
    }

    sort(&edges[axis][0], &edges[axis][2*nPrims]);

C++标准排序例程sort()需要用户提供一个比较函数。我们可以用BoundEdge::t进行比较，如果t值相同，我们用它的类型进行比较（很显然START < END）。

<BoundEdge Public Methods> +=
    bool operator < (const BoundEdge &e) const {
        if (t == e.t) 
            return (int)type < (int)e.type;
        else
            return t < e.t;
    }

边数组排序好后，我们就希望快速地计算出每个分割位置的成本。利用子结点的表面积，我们可以很容易地计算出光线穿过它们的概率。分割面两边的体素个数分别用变量nBelow和nAbove记录。在每次循环中，我们要更新这两个变量，使得nBelow的值是在分割平面之下的体素个数，nAbove的值是分割平面之上的体素个数。

对于第一条边，所有体素都在该边之上，所以nAbove等于nPrims，nBelow等于0。每当循环体考虑在一个包围盒的尾部进行分割时，nAbove需要减1，因为该包围盒已经被列为当前分割平面之上的了，所以当分割完成后，该平面就不能在下一个分割平面之上了。类似地，当计算完分割成本后，如果分割位置在包围盒范围的起始位置，那么这个包围盒就会在所有后继分割平面之下。所以，在循环体的开始和结束要更新这两个变量。

<Compute cost of all splits for axis to find best> =
    int nBelow = 0, nAbove = nPrims;
    for (int i = 0; i < 2 * nPrims; ++i) {
        if (edges[axis].type == BoundEdge::END)  --nAbove;
        float edget = edges[axis].t;
        if (edget > nodeBounds.pMin[axis] && 
            edget < nodeBounds.pMax[axis]) {
            <Compute cost for split at ith edge>
        }
        if (edges[axis].type == BoundEdge::START) ++nBelow;
    }

有了这些信息后，我们可以计算出一个分割所对应的成本。belowSA和aboveSA是两个候选子包围盒的表面积。给定一个轴的编号，我们可以利用数组otherAxis来快速地找出其它两个轴的编号（而无需使用条件语句）。

<Compute cost for split at ith edge> =
    int otherAxis[3][2] = { {1, 2}, {0, 2}, {0, 1}};
    int otherAxis0 = otherAxis[axis][0];
    int otherAxis1= otherAxis[axis][1];
    float belowSA = 2 * (d[otherAxis0] * d[otherAxis1]) +
                (edget - nodeBounds.pMin[axis]) *
                d[otherAxis0] + d[otherAxis1]));
    float aboveSA = 2 * (d[otherAxis0] * d[otherAxis1]) +
                (nodeBounds.pMax[axis] - edget) *
                d[otherAxis0] + d[otherAxis1]));
    float pBelow = belowSA * invTotalSA;
    float pAbove= aboveSA * invTotalSA;
    float eb = (nAbove == 0 || nBelow == 0) ? emptyBonus : 0.f;
    float cost = traversalCost + isectCost * (1.f - eb) * (pBelow *nBelow + pAbove *nAbove);
    <Update best split if this is the lowest cost so far>

如果计算出的成本是目前最好的，就要记录下分割位置的细节。

<Update best split if this is the lowest cost so far> =
    if (cost < bestCost) {
        bestCost = cost;
        bestAxis = axis;
        bestOffset = i;
    }

在前面的测试中，也有可能找不到分割点（下图显示了这样的例子：点虚线表示的多个包围盒都包围这一个树结点，没有令两个子区域含有更少体素的分割位置）。当这种情况发生时，我们要试其它两个轴，如果所有的轴都试过仍然没有找到分割位置，只好放弃，做一个叶子结点完事。
    
另外，也有可能最好的分割所对应的成本仍然比不分割时的成本大。如果这个成本确实很差，并且没有太多的体素，我们就立刻做个叶子结点。否则，我们用badRefines变量来记录当前坏分割的次数。我们容许几次坏分割，因为后面有可能出现不错的分割结果。

<Create leaf if no good splits were found > =
    if(bestAxis == -1 && retry < 2) {
        ++retries;
        axis = (axis + 1) % 3;
        goto retrySplit;
    }
    if(bestCost > oldCost)  ++badRefines;
    if((bestCost > 4.f * oldCost && nPrims < 16) || bestAxis == -1 || badRefines == 3) {
        nodes[nodeNum].initLeaf(primNums, nPrims, mailboxPrims, arena);
        return;
    }

有了分割位置后，就可以用包围盒边的数组对体素进行分类，即确定体素是在分割平面之上，之下，或者两边都在。其方法跟前面记录nBelow和nAbove变量的方法类似。注意我们跳过了bestOffset这个数组项，这是为了避免该体素同时被分类到两个子结点中。

<Classify primitives with respect to split> =
    int n0, n1 = 0;
    for(int i = 0; i < bestOffset; ++i)
        if(edges[bestAxis].type == BoundEdge::START)
            prims0[n0++] = edges[bestAxis].primNum;
    for(int i = bestOffset+1; i < 2 * nPrims; ++i)
        if(edges[bestAxis].type == BoundEdge::END)
            prims1[n1++] = edges[bestAxis].primNum;

回忆一下：在kd-树结点数组中，当前结点的“下面”的子结点是当前结点编号加1。这个子结点的递归结束后，nextFreeNode偏置值就要为“上面”的子结点所用。另一个重要细节是，prims0的内存可以直接传给两个子结点并再被利用，而prims1指针要向后移动nPrims个位置才可以传给两个子结点。其原因是prims1的值要做为参数被用到第二次递归调用中，如果它的指针不后移，内容就会被第一次递归覆盖掉了，第二次递归就会产生错误。然而，没有必要保留edges和prims0的内容。

<Recursively initialize children nodes> =
    float tsplit = edges[bestAxis][bestOffset].t;
    nodes[nodeNum].initInterior(bestAxis, tsplit);
    BBox bounds0 = nodeBounds, bounds1 = nodeBounds;
    bounds0.pMax[bestAxis] = bounds1.pMin[bestAxis] = tsplit;
    buildTree(nodeNum+1, bounds0,
        allPrimBounds, prims0, n0, depth-1, edges,
        prims0, prims1 + nPrims, badRefines);
    nodes[nodeNum].aboveChild = nextFreeNode;
    buildTree(nodes[nodeNum].aboveChild, bounds1,
        allPrimBounds, prims1, n1, depth-1, edges,
        prims0, prims1 + nPrims, badRefines);

可以看出，为了应付最坏的情况，prims1所需内存比prims0要大得多：

<Allocate working memory for kd-tree construction> +=

    int *prims0 = new int[prims.size()];
    int *prims1 = new int[(maxDepth + 1) * prims.size()];

4.4.3 遍历

如图，图中显示了光线遍历树的基本过程。首先，tmin和tmax被初始化为光线跟树的包围盒的交点的参数距离。跟网格加速器一样，如果光线和场景的包围盒没有交点，这个函数就立即返回false。否则，就要从树的根结点开始，遍历这个树。在每个内部结点，算法确定光线先进入两个子结点的哪一个，并安装先后顺序遍历子结点。当光线走出树，或者找到了最近交点，遍历就结束。

<KdTreeAccel Method Definitions> +=
    bool KdTreeAccel::Intersect(const Ray &ray, Intersection *isect) const {
        <Compute initial parametric range of ray inside kd-tree extent>
        <Prepare to traverse kd-tree for ray>
        <Traverse kd-tree node in order for ray>
    }

算法先寻找光线和包围盒相交的参数范围[tmin, tmax]，如果没有交点则立即退出。

<Compute initial parametric range of ray inside kd-tree extent> =
    float tmin, tmax;
    if(!bounds.IntersectP(ray, &tmin, &tmax))
        return false;

在树遍历开始之前，先要为光线赋予一个信箱编号，并计算出光线方向各个分量的倒数，这样就可以在循环体中避免大量的除法运算。KdToDo结构数组用来存放还没有处理的结点，它的最后一项应该是下一个要考虑的结点。数组的大小是kd-树的最大深度。

<Prepare to traverse kd-tree for ray> =
    int rayId = curMailboxId++;
    Vector invDir(1.f/ray.d.x, 1.f/ray.d.y, 1.f/ray.d.z);
    #define MAX_TODO 64
    KdToDo todo[MAX_TODO];
    int todoPos = 0;

<KdTreeAccel Declarions> +=
    struct KdToDo {
        const KdAccelNode *node;
        float tmin, tmax;
    };

遍历过程对nodes数组进行循环，每次循环处理一个叶子结点或内部结点。

<Traverse kd-tree nodes in order for ray> =
    bool hit = false;
    const KdAccelNode *node = &nodes[0];
    while(node != NULL) {
        <Bail out if we found a hit closer than the current node>
        if (!node ->IsLeaf()) {
            <Process kd-tree interior node>
        }
        else {
            <Check for intersections inside leaf node>
            <Grab next node to process from todo list>
        }
    }
    return hit;

有可能存在这样的情况，在以前的求交测试中找到了覆盖多个结点的体素与光线的交点。如果第一次发现这样的交点，并且交点在当前结点之外，就必须继续遍历树，直到所遇到结点的tmin比交点还要远。只有这时，我们才能确定没有比这个交点更近的其它的交点了。

<Bail out if we found a hit closer than the current node> =
    if (ray.maxt <tmin) break;

对于内部结点，要做的第一件事就是求得光线和分割平面的交点，并确定是否处理其中一个结点，还是两个结点都要处理，并且确定处理两个子结点的先后次序。

<Process kd-tree interior node > =
    <Compute parametric distance along ray to split plane>
    <Get node children pointers for ray>
    <Advance to next child node, possibly enqueue other child>

光线和分割平面的求交方法跟光线和包围盒交点的求法类似。

<Compute parametric distance along ray to split plane> =
    int axis = node->SplitAxis();
    float tplane = (node->SplitPos() – ray.o[axis]) * invDir[axis];

现在我们要确定光线进入子结点的次序，使得树的遍历是沿着光线从前向后的顺序进行的。很显然，只需比较光线原点和分割面位置，就可以确定先后次序了。（如图）。

<Get node children pointers for ray> =
    const KdAccelNode *firstChild, *secondChild;
    int belowFirst = ray.o[axis] <= node->SplitPos();
    if(belowFirst) {
        firstChild = node + 1;
        secondChild = &nodes[node->aboveChild];
    }
    else {
        firstChild = &nodes[node->aboveChild];
secondChild = node + 1;
    }

有时候没有必要对两个子结点都进行处理。如图，光线只穿过其中一个子结点。下面代码中的第一个if语句对应图中(a)的情形：如果光线背对远结点(far node)，或者并不穿过远结点(tsplit > tmax)，那么我们只需处理近结点(near node)。下面代码中的第二个if语句对应图中(b)的情形：如果光线不穿过近结点(near node)，就不必处理它。第二个else语句对两个子结点都处理：先处理近结点，并把远结点放入todo数组中。

    <Advance to next child node, possibly enqueue other child> =
        if(tplane > tmax || tplane < 0)
            node = firstChild;
        else if (tplane < tmin)
            node = secondChild;
        else {
            <Enqueue secondChild in todo list>
            node = firstChild;
            tmax = tplane;
        }
<Enqueue secondChild in todo list> =
    todo[todoPos].node = secondChild;
    todo[todoPos].tmin = tplane;
    todo[todoPos].tmax = tmax;
    ++todoPos;

如果当前结点是叶子结点，我们就可以对其中的体素进行求交测试了，其中也用到了信箱检查技术，从而避免了多余的求交测试。

<Check for intersections inside leaf node> =
    u_int nPrimitives = node->nPrimitives();
    if(nPrimitives == 1) {
        MailboxPrim *mp = node->onePrimitive;
        <Check one primitive inside leaf node>
    }
    else {
        MailboxPrim **prims = node->primitives;
        for(u_int i = 0; i < nPrimitives; ++i) {
            MailboxPrim *mp = prims;
            <Check one primitive inside leaf node>
        }
    }

处理单个的体素很简单：只需做些信箱检查工作，并调用体素的求交函数：

<Check one primitive inside leaf node> =
    if(mp->lastMailboxId != rayId) {
        mp->lastMailboxId = rayId;
        if(mp->primitive->Intersect(ray, isect))
            hit = true;
    }

做完叶子结点上的求交测试以后，要从todo数组中取下一个要处理的结点。如果数组是空的，光线就穿过了整个树，并没有发现交点。

<Grab next node to process from todo list> =
    if(todoPos > 0) {
        --todoPos;
        node = todo[todoPos].node;
        tmin = todo[todoPos].tmin;
        tmax = todo[todoPos].tmax;
    }
    else
        break;

跟GridAccel一样，KdTreeAccel也有一个IntersectP()函数，它跟Intersect()极为相似，不同的是，只要发现交点就立即返回，并不寻找最近的交点。这里也不讨论它的实现了。

<KdTreeAccel Public Methods> =
    bool IntersectP(const Ray &ray) const;