光线求交加速算法：kd-树

空间二分树,即Binary space partitioning (BSP)树利用分割平面自适应地细分空间。 BSP树以包围整个场景的边界框开始。如果框中的图元数量大于某个阈值，则该框将被平面分成两个区域。然后，将图元与它们重叠的区域相关联。如果图元同时位于两个区域，则图元将与两个区域都关联(这与BVH不同，在BVH中，每个图元在拆分后仅分配给两个子组之一)。

拆分过程将递归地继续进行，直到树中的每个叶区域包含足够少的图元，或者直到达到最大深度为止。因为可以将拆分平面放置在整体边界内的任意位置，并且可以将3D空间的不同部分细化到不同的程度，所以BSP树可以轻松处理几何图形的不均匀分布。

BSP树的两个变体是kd树和八叉树。 kd树仅将分割平面限制为垂直于坐标轴之一，这使得树的遍历和构造都更加有效，但在空间划分方式上有些灵活性。八叉树在每个步骤中使用三个垂直于轴的平面将框同时划分为八个区域(通常是通过在每个方向上划分范围的中心)。这里，我们只介绍kd树。我们给出kd树类定义：

class KdTreeAccel : public Aggregate {public:KdTreeAccel(std::vector<std::shared_ptr<Primitive>> p,int isectCost = 80, int traversalCost = 1,Float emptyBonus = 0.5, int maxPrims = 1, int maxDepth = -1);Bounds3f WorldBound() const { return bounds; }~KdTreeAccel();bool Intersect(const Ray &ray, SurfaceInteraction *isect) const;bool IntersectP(const Ray &ray) const;private:void buildTree(int nodeNum, const Bounds3f &bounds,const std::vector<Bounds3f> &primBounds, int *primNums,int nprims, int depth,const std::unique_ptr<BoundEdge[]> edges[3], int *prims0,int *prims1, int badRefines = 0);const int isectCost, traversalCost, maxPrims;const Float emptyBonus;std::vector<std::shared_ptr<Primitive>> primitives;std::vector<int> primitiveIndices;KdAccelNode *nodes;int nAllocedNodes, nextFreeNode;Bounds3f bounds;
};

kd树构建的过程如下图：

通过沿三个坐标轴之一递归拆分场景的边界框来构建kd树。在这里，第一个分割沿x轴。它的放置位置使三角形在右侧区域恰好是单独的，其余图元在左侧结束。然后使用与轴对齐的分割平面将左侧区域再细化几次。细化标准的细节(在每个步骤中使用哪个轴分割空间，在平面上沿该轴放置在哪个位置以及细化终止的点)在实践中都会严重影响树的性能。

树形式

kd树是二叉树，其中每个内部节点始终都具有两个子节点，并且树的叶子存储与其重叠的图元。每个内部节点必须提供以下三部分信息：

分割轴：以x，y或z轴中的哪一个作为分割轴。
分割位置：分割平面沿轴的位置。根据分割轴和分割位置，例如x=5，则表示取x轴为分割轴，分割平面(分割位置)为x=5。
子节点：有关如何到达其下方的两个子节点的信息。

每个叶节点仅需要记录哪些图元与它重叠。

确保所有内部节点和大部分叶节点仅使用8字节的内存是一件非常重要的事，因为这样做可以确保8个节点适合64字节的缓存行。因为树中通常有许多节点，并且由于每个射线通常访问许多节点，所以最小化节点表示的大小会大大提高缓存性能。当我们将大小减小到8个字节时，相比16字节的节点表示形式速度提高了近20％。

叶子和内部节点均由以下KdAccelNode结构表示。每个联合成员之后的注释指示将特定字段用于内部节点，叶节点还是两者都使用：

struct KdAccelNode {void InitLeaf(int *primNums, int np, std::vector<int> *primitiveIndices);void InitInterior(int axis, int ac, Float s) {split = s;flags = axis;aboveChild |= (ac << 2);}//一些获取属性的方法Float SplitPos() const { return split; }int nPrimitives() const { return nPrims >> 2; }int SplitAxis() const { return flags & 3; }bool IsLeaf() const { return (flags & 3) == 3; }int AboveChild() const { return aboveChild >> 2; }union {float split;                 // 内部结点int onePrimitive;            // 叶子结点int primitiveIndicesOffset;  // 叶子结点};private:union {int flags;       // 两者都使用int nPrims;      // 叶子结点int aboveChild;  // 内部结点};
};

KdAccelNode::flags变量的两个低阶位用于区分不同轴(x，y和z)拆分的内部节点(其中这些位分别保存值0、1和2)和叶节点(位保存值3)。将叶子节点存储在8个字节中相对容易：由于KdAccelNode::flags的低2位用于指示这是一片叶子，因此KdAccelNode::nPrims的高30位可用于记录有多少个图元重叠。然后，如果仅单个图元与KdAccelNode叶子重叠，则onePrimitive用来索引KdTreeAccel::primitives数组中的图元。如果多个图元重叠，则它们的索引将存储在KdTreeAccel::primitiveIndices的一段中。叶子的第一个索引的偏移量存储在KdAccelNode :: primitiveIndicesOffset中，其余索引紧密排列跟随其后。

叶子节点的标志和nPrims共享相同的存储。我们需要注意在初始化另一个数据时不要破坏其中一个数据。此外，图元的数量在存储之前必须向左移动两位，以便KdAccelNode :: flags的低两位都可以设置为1，以指示这是叶节点：

void KdAccelNode::InitLeaf(int *primNums, int np,std::vector<int> *primitiveIndices) {flags = 3;nPrims |= (np << 2);// 在叶子结点中保存图元的索引(即在KdTreeAccel::primitives数组中的位置)if (np == 0)//内部结点图元数量总是为0onePrimitive = 0;else if (np == 1)onePrimitive = primNums[0];else {primitiveIndicesOffset = primitiveIndices->size();for (int i = 0; i < np; ++i) primitiveIndices->push_back(primNums[i]);}
}

对于具有零个或一个重叠图元的叶节点，直接使用KdAccelNode :: onePrimitive记录索引(内部结点为0)。对于多个图元重叠的情况，可以用外部originalIndices数组的指针来存储。

将内部节点减小到8个字节叶比较容易。Float存储沿所选拆分轴的位置，节点在该拆分轴上拆分空间，并且如前所述，KdAccelNode :: flags的最低2位用于记录节点沿哪个轴拆分。剩下的就是存储足够的信息，以便在遍历树时找到该节点的两个子节点。

我们不存储两个指针或偏移量，而是仅存储一个子指针的方式设置节点：所有节点都分配在单个连续的内存块中，而内部节点的子节点负责拆分平面下方的空间始终存储在其父项随后的数组位置。另一个代表分裂平面上方空间的子元素将终止于数组中的其他位置，我们使用单个整数偏移量KdAccelNode :: aboveChild将其位置存储在节点数组中。此表示类似于BVH篇中BVH节点的表示。

考虑到所有这些约定，用于初始化内部节点的代码与InitLeaf()方法一样，首先将值分配给flags并计算索引移位后的值，之后将移位后的值和aboveChild求或运算，以免破坏标志中存储的位:

void InitInterior(int axis, int ac, Float s) {split = s;flags = axis;aboveChild |= (ac << 2);}

构建树

kd树使用递归自顶向下算法构建。在每个步骤中，我们都有一个与轴对齐的空间区域，以及一组与该区域重叠的图元。将该区域分为两个子区域并转换为内部节点，或者使用重叠基元创建叶节点，从而终止递归。

如有关KdAccelNodes的讨论中所述，所有树节点都存储在连续的数组中。KdTreeAccel::nextFreeNode记录此数组中下一个可用节点，而KdTreeAccel::nAllocedNodes记录已分配的总数。通过将它们都设置为0且在启动时不分配任何节点，此处的实现可确保在初始化树的第一个节点时立即进行分配。

如果未将最大树深提供给构造函数，则还需要确定最大树深。尽管树的构建过程通常会在合理的深度自然终止，但是限制最大深度非常重要，这样，用于树的内存量就不会在病理情况下不受限制地增长。我们发现值8 + 1.3log(N)为各种场景提供了合理的最大深度，构建kd树代码如下：

KdTreeAccel::KdTreeAccel(std::vector<std::shared_ptr<Primitive>> p,int isectCost, int traversalCost, Float emptyBonus,int maxPrims, int maxDepth): isectCost(isectCost),traversalCost(traversalCost),maxPrims(maxPrims),emptyBonus(emptyBonus),primitives(std::move(p)) {//构建kd树nextFreeNode = nAllocedNodes = 0;if (maxDepth <= 0)maxDepth = std::round(8 + 1.3f * Log2Int(int64_t(primitives.size())));// +计算所有图元的边界框// +为kd-树构建分配工作内存// +初始化primNums以进行kd-树构建// +递归构建kd-树
}

由于图元的边界框会反复使用，因此我们首先用一个vector保存所有的图元边界框：

// =计算所有图元的边界框
std::vector<Bounds3f> primBounds;
for (const std::shared_ptr<Primitive> &prim : primitives) {Bounds3f b = prim->WorldBound();bounds = Union(bounds, b);//bounds为KdTreeAccel的一个成员变量，保存整个图元所在空间的边界框primBounds.push_back(b);
}

我们需要分配一些内存用以kd-树的构建，我们首先介绍边界框边缘结构体BoundEdge：

enum class EdgeType { Start, End };//起始边缘和末尾边缘
struct BoundEdge {BoundEdge() { }BoundEdge(Float t, int primNum, bool starting): t(t), primNum(primNum) {type = starting ? EdgeType::Start : EdgeType::End;//光线方向反向则交换首尾边缘 }float t;//沿轴分割位置int primNum;//图元索引EdgeType type;
};

如我们沿着x轴作为分割轴，如下图：

则A图元边界框会有 $a_{0}$ 和 $a_{1}$ 两个边缘。因此我们可以使用BoundEdge来保存所有图元边界框边缘信息。因此我们分配内存代码如下：

// =为kd-树构建分配工作内存
std::unique_ptr<BoundEdge[]> edges[3];
for (int i = 0; i < 3; ++i)//为x,y,z都分配空间edges[i].reset(new BoundEdge[2 * primitives.size()]);//每个框两个面
std::unique_ptr<int[]> prims0(new int[primitives.size()]);
std::unique_ptr<int[]> prims1(new int[(maxDepth + 1) * primitives.size()]);

prims0和prims1保存图元索引。prims0分配的内存用以处理临近当前结点的子结点，prims1分配的内存用以处理第二个子结点。由于在深度优先遍历的过程中临近的子结点会在下一次遍历中直接使用，因此只需要分配primitives.size()的大小(每次递归时可以更新内容)。而当前结点在深度优先遍历时，第二个子结点会在何时遍历我们并不能预测(我们可能会在很后面才处理根附近结点的第二个子结点)，所以分配大小为(maxDepth + 1) * primitives.size()，分别处理不同层结点的第二个子结点。过程如下图所示：

树构造函数的参数之一是图元索引数组，表明哪些图元与当前节点重叠。因为所有图元都与根节点重叠(当递归开始时)，所以我们初始化该数组的值从零到Primitives.size()-1:

// =初始化primNums以进行kd-树构建
std::unique_ptr<int[]> primNums(new int[primitives.size()]);
for (size_t i = 0; i < primitives.size(); ++i)primNums[i] = i;

为每个树节点调用KdTreeAccel :: buildTree()。它负责确定该节点是内部节点还是叶节点，并适当地更新数据结构:

// =递归构建kd-树
buildTree(0, bounds, primBounds, primNums.get(), primitives.size(), maxDepth, edges, prims0.get(), prims1.get());

KdTreeAccel::buildTree()的主要参数是KdAccelNodes数组中要用于创建的节点偏移量nodeNum；提供节点覆盖空间区域的边界框nodeBounds；以及与之重叠的图元的索引primNums。其余参数将在后面进行描述：

void KdTreeAccel::buildTree(int nodeNum, const Bounds3f &nodeBounds,const std::vector<Bounds3f> &allPrimBounds, int *primNums,int nPrimitives, int depth,const std::unique_ptr<BoundEdge[]> edges[3], int *prims0, int *prims1, int badRefines) {// +从nodes数组中获取下一个空闲结点// +如果达到递归边界(该区域中的图元数量足够少，或者已经达到最大深度)，创建叶子结点// +初始化内部结点，继续递归
}

如果所有分配的节点都已用完，则将使用两倍的条目重新分配节点内存，并复制旧值。第一次调用KdTreeAccel :: buildTree()时，KdTreeAccel::nAllocedNodes为0，并分配了512个树节点的初始块。

// =从nodes数组中获取下一个空闲结点
if (nextFreeNode == nAllocedNodes) {int nNewAllocNodes = std::max(2 * nAllocedNodes, 512);KdAccelNode *n = AllocAligned<KdAccelNode>(nNewAllocNodes);if (nAllocedNodes > 0) {memcpy(n, nodes, nAllocedNodes * sizeof(KdAccelNode));FreeAligned(nodes);}nodes = n;nAllocedNodes = nNewAllocNodes;
}
++nextFreeNode;

如果该区域中的图元数量足够少，或者已经达到最大深度，则会创建一个叶节点(停止递归)。深度参数以树的最大深度开始，并在每个级别递减:

// =如果达到递归边界(该区域中的图元数量足够少，或者已经达到最大深度)，创建叶子结点
if (nPrimitives <= maxPrims || depth == 0) {nodes[nodeNum].InitLeaf(primNums, nPrimitives, &primitiveIndices);return;
}

如果这是一个内部节点，则必须选择一个分割平面，相对于该平面对图元进行拆分，然后递归：

// =初始化内部结点，继续递归// +为内部结点选择分割轴和位置// +如果没有找到分割点，则创建叶子结点// +对图元进行划分// +递归初始化子结点

我们的实现使用SAH选择一个拆分。SAH适用于kd树和BVH。在此，为节点中的一系列候选分割平面计算估算的成本，并选择给出最低成本的分割。

在这里的实现中，用户可以设置求交成本 $t_{isect}$ 和遍历成本 $t_{trav}$ 。它们的默认值分别是80和1。最终，这两个值的比率决定了树构建算法的行为。与用于BVH构造的值相比，这些值之间更大的比率反映了以下事实：访问kd-tree节点比访问BVH节点便宜。

对用于BVH树的SAH的一种修改是，对于kd树，应该优先选择那些子结点没有图元相互重叠的拆分方案，因为穿过这些区域的光线可以立即前进到下一个没有任何射线图元相交测试的结点。因此，未分割区域(直接将该结点作为叶子结点)和分割区域的修订成本分别为:

如果分割后两个区域之一完全为空，则 $b_{e}$ 为0到1之间的值。否则当两个区域都有图元则 $b_{e}$ 为0。

给定一种计算成本模型概率的方法，唯一要解决的问题是如何生成候选分割位置以及如何有效地计算每个候选的成本。可以证明，该模型的最低成本是选择与图元边界框之一的面重合的拆分来实现，无需在中间位置考虑拆分。这时我们可以使用上面介绍的BoundEdge结构，该结构保存的边缘信息可以用来计算最低成本。

在确定了创建叶子的估计成本之后，KdTreeAccel :: buildTree()选择一个轴尝试进行分割，并为每个候选分割计算成本。 bestAxis和bestOffset记录了迄今为止成本最低(bestCost)的轴和边界框边缘索引。 invTotalSA初始化为节点表面积的倒数。当计算光线通过每个候选子节点的概率时，将使用其值：

// =为内部结点选择分割轴和位置
int bestAxis = -1, bestOffset = -1;
float bestCost = Infinity;
float oldCost = isectCost * Float(nPrimitives);//全部用叶子结点代替时的花费
float totalSA = nodeBounds.SurfaceArea();
float invTotalSA = 1 / totalSA;
Vector3f d = nodeBounds.pMax - nodeBounds.pMin;int axis = nodeBounds.MaximumExtent();// 选择x,y,z中范围最大作为划分轴
int retries = 0;
retrySplit:
// +为相应轴初始化边缘
// +计算相应轴的所有分割点的成本，以找出最小成本分割轴

该方法首先尝试找选择x,y,z中范围最大作为划分轴。如果成功，则此选择有助于提供趋于呈正方形的空间区域。如果未能沿该轴找到良好的分割，它将返回并依次尝试其他分割。

寻找好的分割点过程：首先，使用重叠图元的边界框初始化轴的edges数组。然后，将数组沿轴从低到高排序，以便可以从头到尾扫过框的边缘：

// =为相应轴初始化边缘
for (int i = 0; i < nPrimitives; ++i) {int pn = primNums[i];const Bounds3f &bounds = allPrimBounds[pn];edges[axis][2 * i] =     BoundEdge(bounds.pMin[axis], pn, true);edges[axis][2 * i + 1] = BoundEdge(bounds.pMax[axis], pn, false);
}
// +为相应轴排序

利用sort函数给相应轴排序：

//为相应轴排序
std::sort(&edges[axis][0], &edges[axis][2*nPrimitives],[](const BoundEdge &e0, const BoundEdge &e1) -> bool {if (e0.t == e1.t)return (int)e0.type < (int)e1.type;else return e0.t < e1.t; });

得到边缘的排序数组后，我们希望针对每个边缘的分割快速计算出成本函数。穿过每个子节点的射线的概率很容易使用它们的表面积来计算，并且在拆分的每一侧上图元的数量由变量nBelow和nAbove跟踪。我们保持它们的值更新，以便如果我们选择在edgeT处进行特定遍历循环拆分，则nBelow将给出最终在拆分平面以下的图元数量，而nAbove将给出在拆分平面以上的图元数量。

在第一个边缘，根据定义，所有图元必须在该边缘上方，因此nAbove初始化为nPrimitives，nBelow设置为0。当循环在边界框范围的末端考虑拆分时，需要将nAbove减1，因为该框之前位于分割平面上方，所以如果在该点进行分割，则该框将不再位于分割平面上方。同样，在计算拆分成本之后，如果拆分候选对象位于边界框范围的开始位置，则该框将位于所有后续拆分的下侧。过程如下图(红线为分割的边缘)：

循环主体开始和结尾处的测试会更新这两种情况的原始计数。代码如下：

// =计算相应轴的所有分割点的成本，以找出最小成本分割轴
int nBelow = 0, nAbove = nPrimitives;
for (int i = 0; i < 2 * nPrimitives; ++i) {if (edges[axis][i].type == EdgeType::End) --nAbove;float edgeT = edges[axis][i].t;if (edgeT > nodeBounds.pMin[axis] &&edgeT < nodeBounds.pMax[axis]) {// +计算分割线为第i条边缘时的成本}if (edges[axis][i].type == EdgeType::Start) ++nBelow;
}

我们接下来利用上文提到的成本计算公式写出代码：

// =计算分割线为第i条边缘时的成本
// 计算分割点在edgeT时两个子区域的表面积
int otherAxis0 = (axis + 1) % 3, otherAxis1 = (axis + 2) % 3;
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;//概率pb
float pAbove = aboveSA * invTotalSA;//概率pa
float eb = (nAbove == 0 || nBelow == 0) ? emptyBonus : 0;
//成本计算
float cost = traversalCost + isectCost * (1 - eb) * (pBelow * nBelow + pAbove * nAbove);// 更新最低成本
if (cost < bestCost) {bestCost = cost;bestAxis = axis;bestOffset = i;
}

但是如果沿轴计算始终无法找到分割点，如下图：

我们设定一个重置次数retries，如果其他两个轴都找不到分割点(此时重置次数为2)，则直接创建叶子结点，保存所有图元。

最佳拆分的成本也可能仍高于完全不拆分节点的成本(oldCost)。如果情况更糟(bestCost > 4 * oldCost)并且图元较少，则会立即创建一个叶子节点。否则，badRefines会跟踪到目前为止在树的当前节点上方已进行了多少次较差的拆分。我们允许有一些稍微差的拆分，因为以后的拆分可能会找到更好的拆分(只包含很少图元)。直到badRefines值为3，我们也立即创建叶子结点。代码如下：

// =如果没有找到分割点，则创建叶子结点
if (bestAxis == -1 && retries < 2) {++retries;axis = (axis + 1) % 3;goto retrySplit;
}
if (bestCost > oldCost) ++badRefines;
if ((bestCost > 4 * oldCost && nPrimitives < 16) || bestAxis == -1 || badRefines == 3) {nodes[nodeNum].InitLeaf(primNums, nPrimitives, &primitiveIndices);return;
}

选择拆分位置后，可以使用边界框边缘将图元分类为拆分的上方，下方或两侧，方法与跟踪早期代码中的nBelow和nAbove相同。注意，在下面的循环中跳过了数组中的bestOffset条目。这样就不会将其边界框边缘用于拆分的图元不正确地分类为在拆分的两侧。

// =对图元进行划分(如果图元穿过划分线，则该图元会分配到两个区域)
int n0 = 0, n1 = 0;
for (int i = 0; i < bestOffset; ++i)if (edges[bestAxis][i].type == EdgeType::Start)prims0[n0++] = edges[bestAxis][i].primNum;
for (int i = bestOffset + 1; i < 2 * nPrimitives; ++i)if (edges[bestAxis][i].type == EdgeType::End)prims1[n1++] = edges[bestAxis][i].primNum;

回想一下，kd-树节点数组中此节点“下一个”子节点的节点号是当前节点号加1。从树的那一侧返回递归后，将nextFreeNode偏移量用于“上方”第二个子结点。唯一重要的细节是，prims0内存直接传递给两个孩子，以供重用。而prims1指针则先向前推进nPrimitives。这是必需的，因为当前对KdTreeAccel :: buildTree()的调用总是先递归处理第一个子结点，第二个子结点总是在左分支全部处理完成后才开始处理，因此之前结点的第二个子结点都会保存在指针起点为prims1+k*nPrimitives的内存块中，其中k为结点的深度。(可以根据之前的分配prims0和prims1内存时的图进行理解)

// =递归初始化子结点
float tSplit = edges[bestAxis][bestOffset].t;
Bounds3f bounds0 = nodeBounds, bounds1 = nodeBounds;
bounds0.pMax[bestAxis] = bounds1.pMin[bestAxis] = tSplit;//计算两个子结点的边界
buildTree(nodeNum + 1, bounds0, allPrimBounds, prims0, n0,depth - 1, edges, prims0, prims1 + nPrimitives, badRefines);//递归构建第一个子结点
int aboveChild = nextFreeNode;
nodes[nodeNum].InitInterior(bestAxis, aboveChild, tSplit);
buildTree(aboveChild, bounds1, allPrimBounds, prims1, n1, depth - 1, edges, prims0, prims1 + nPrimitives, badRefines);//递归构建第二个子结点

遍历树

下图显示了光线穿过树的基本过程：

将光线与树的整体边界相交会得到初始的tMin和tMax值，在图中用点标记。如果射线错过了整个基本图元边界，则此方法可以立即返回false。否则，它将开始从根开始下降到树中。在每个内部节点处，它确定射线首先进入两个子节点中的哪个子节点，并按顺序处理两个子节点。当光线从树上离开或找到最接近的相交点时，遍历结束。以上图为例：(a)光线和根结点判断相交成功。(b)把tSplit作为子结点的分割面，进行深度优先遍历。(c)近结点未与光线相交，因此判断光线是否和远结点相交。(d)光线与该结点相交，然后继续用深度优先遍历该结点。

遍历树的过程基本都在光线求交中，代码如下：

bool KdTreeAccel::Intersect(const Ray &ray, SurfaceInteraction *isect) const {// +计算kd树整体边界框范围内射线的初始参数范围// +用ray遍历kd树的准备// +根据光线的顺序遍历结点
}

该算法首先找到射线与树重叠的总参数范围[tmin，tmax]，如果没有重叠则立即退出：

// =计算kd树整体边界框范围内射线的初始参数范围
float tMin, tMax;
if (!bounds.IntersectP(ray, &tMin, &tMax)) return false;

KdToDo结构的数组用于记录尚未为射线处理的节点。它是有序的，因此数组中最后一个元素是应该考虑的下一个节点。此数组中所需的最大条目数是kd树的最大深度，实际上，以下使用的数组大小就足够了。

struct KdToDo {const KdAccelNode *node;Float tMin, tMax;
};// =用ray遍历kd树的准备
vector3f invDir(1 / ray.d.x, 1 / ray.d.y, 1 / ray.d.z);
constexpr int maxTodo = 64;
KdToDo todo[maxTodo];
int todoPos = 0;

每次遍历循环都处理单个叶节点或内部节点。 tMin和tMax值将始终保持射线与当前节点重叠的参数范围：

// =根据光线的顺序遍历结点
bool hit = false;
const KdAccelNode *node = &nodes[0];
while (node != nullptr) {// +获取当前结点内最近的交点，进入下个结点时退出遍历(下个结点tmin更新，将大于ray.tMax)。 if (!node->IsLeaf()) {// +处理内部结点} else {// +检查叶节点内的交点// +从todo数组中获取下一个要处理的节点}
}
return hit;

如果光线已经和图元找到了交点(ray.tmax会更新为交点值)，我们还需要在当前结点的边界框tmin内继续遍历，直到确定没有更近的交点：

// =获取当前结点内最近的交点，进入下个结点时退出遍历(下个结点tmin更新，将大于ray.tMax)。
if (ray.tMax < tMin) break;

对于内部树节点，首先要做的是使光线与节点的分割平面相交。给定相交点，我们可以确定需要处理一个还是两个子节点，以及光线以什么顺序穿过它们：

// =处理内部结点// +沿着光线计算到分割平面交点的距离// +获取子结点顺序// +进入下一个子结点，如果光线经过两个子结点，则将第二个子结点入队

到分裂平面距离的计算方法与计算射线与边界框的交点所用的方法相同。每次循环时，我们使用预先计算的invDir值保存一个除数：

// =沿着光线计算到分割平面交点的距离
int axis = node->SplitAxis();
Float tPlane = (node->SplitPos() - ray.o[axis]) * invDir[axis];

现在，必须确定射线遇到子节点的顺序，以便沿着射线以从前到后的顺序遍历树。下图显示了两种可能的几何形状：

射线的原点相对于分离平面的位置足以区分这两种情况。还有一种极端的情况，光线的原点能位于分裂平面上(ray.o[axis] == node->SplitPos())，这需要使用光线方向来确定子结点遍历的顺序。

// =获取子结点顺序
const KdAccelNode *firstChild, *secondChild;//根据射线起点位置确定先遍历哪个子结点
int belowFirst = (ray.o[axis] <  node->SplitPos()) ||(ray.o[axis] == node->SplitPos() && ray.d[axis] <= 0);
if (belowFirst) {firstChild = node + 1;secondChild = &nodes[node->AboveChild()];
} else {firstChild = &nodes[node->AboveChild()];secondChild = node + 1;
}

可能不必同时处理该节点的两个子节点。下显示了一些情况：

其中光线仅通过其中一个子代。(a)顶端光线与分裂平面tSplit相交，超出了光线的最大位置tMax，因此不会进入远子结点。底射线背向分裂平面，由负tSplit值表示。(b)射线在tMin值之前与平面tSplit相交，表明近的子结点不需要处理。代码表示：

// =进入下一个子结点，如果光线经过两个子结点，则将第二个子结点入队
if (tPlane > tMax || tPlane <= 0)//上图(a)的情况node = firstChild;
else if (tPlane < tMin)//上图(b)情况node = secondChild;
else {//两个子结点都要处理//入队第二个子结点todo[todoPos].node = secondChild;todo[todoPos].tMin = tPlane;todo[todoPos].tMax = tMax;++todoPos;node = firstChild;tMax = tPlane;
}

处理完成内部结点的情况，接下来考虑当前节点是叶子结点的情况，则针对叶中的图元执行相交测试：

// =检查叶节点内的交点
int nPrimitives = node->nPrimitives();
if (nPrimitives == 1) {const std::shared_ptr<Primitive> &p = primitives[node->onePrimitive]; if (p->Intersect(ray, isect)) hit = true;
} else {for (int i = 0; i < nPrimitives; ++i) {int index = primitiveIndices[node->primitiveIndicesOffset + i];const std::shared_ptr<Primitive> &p = primitives[index];if (p->Intersect(ray, isect)) hit = true;}
}

在叶节点上进行求交测试之后，要处理的下一个节点将从todo数组中获取。如果没有更多的节点，则射线已经穿过树而没有击中任何物体：

// =从todo数组中获取下一个要处理的节点
if (todoPos > 0) {--todoPos;node = todo[todoPos].node;tMin = todo[todoPos].tMin;tMax = todo[todoPos].tMax;
}
elsebreak;

kd-树到这里就完成了。虽然SVH是从图元进行分割，BSP是从空间进行分割，但二者的树形还是非常类似的。

光线求交加速算法：kd-树相关推荐

光线求交加速算法：边界体积层次结构(Bounding Volume Hierarchies)1-BVH引入
光线求交加速算法:边界体积层次结构(Bounding Volume Hierarchies)1 BVH引入光线和物体求交的加速算法中,最常见的是物体(图元)细分和空间细分.边界体积层次结构(BVH) ...
光线求交加速算法：边界体积层次结构(Bounding Volume Hierarchies)3-LBVH(Linear Bounding Volume Hierarchies)
光线求交加速算法:边界体积层次结构(Bounding Volume Hierarchies)3 尽管使用表面积启发式方法(SAH)构建边界体积层次结构会产生很好的结果,但是该方法确实存在两个缺点:首先 ...
光线求交加速算法：边界体积层次结构(Bounding Volume Hierarchies)2-表面积启发式法(The Surface Area Heuristic)
光线求交加速算法:边界体积层次结构(Bounding Volume Hierarchies)2 上篇的两种图元分区方法(Middle,EqualCounts)对于某些图元分布可以很好地工作,但是在实践 ...
用python手写KNN算法+kd树及其BBF优化（原理与实现）（下篇）
用python手写KNN算法+kd树及其BBF优化(原理与实现)(下篇) 接上一篇用python手写KNN算法+kd树及其BBF优化(原理与实现)(上篇) 我们使用training2和test2两个数 ...
K近邻快速算法 -- KD树、BBF改进算法
K近邻算法即是查找与当前点(向量)距离最近的K个点(向量),距离计算一般用欧氏距离. 最简单的方法就是穷举法:计算每个向量与当前向量的欧氏距离,选取最小的K个为所求.但这种方法计算量太大,无法应对大样 ...
数据结构和算法——kd树
一.K-近邻算法 K-近邻算法是一种典型的无参监督学习算法,对于一个监督学习任务来说,其 m m个训练样本为: {(X(1),y(1)),(X(2),y(2)),⋯,(X(m),y(m))} \lef ...
k近邻算法——kd树
一.kd树概念 kd树(K-Dimensional Tree)是一种对K维空间中的实例点进行存储以便对其进行快速检索的树形数据结构. kd树是二叉树,表示对K维空间的一个划分(partitio ...
多方隐私集合求交高效算法解读
前言之前我们介绍过两方的「隐私集合求交算法」,可以应用到计算广告的实际效果,寻找联系人,联邦学习的特征对齐等场景,例如:在新的APP上找到共同的微信好友.开会时找到所有参会者共同的空闲时间等,但是这 ...
基于RSA盲签名的隐私集合求交PSI算法
记录一种PSI算法. 在<Practical Private Set Intersection Protocols with Linear Computational and Bandwidth ...

光线求交加速算法：kd-树

光线求交加速算法：kd-树

树形式

构建树

遍历树

光线求交加速算法：kd-树相关推荐

最新文章

热门文章