我的kd树遍历代码哪里错了？答案

【问题标题】：Where is my kd tree traversal code wrong?我的kd树遍历代码哪里错了？
【发布时间】：2014-10-11 14:22:22
【问题描述】：

我正在优化我的 C++ 光线追踪器。我正在通过 kdtrees 跟踪单条光线。到目前为止，我使用的是 Havran 的递归算法“B”，这对于 OOP 来说似乎很古老且过于夸张。我的新代码尽可能短（希望更容易被编译器优化并且更容易维护）：

    struct StackElement{
        KDTreeNode<PT>* node;
        float tmax;
        array<float, 3> origin;
    };

    //initializing explicit stack
    stack<StackElement> mystack;

    //initialize local variables
    KDTreeNode<PT>* node = tree.root;
    array<float, 3> origin {ray.origin[0], ray.origin[1], ray.origin[2]};
    const array<float, 3> direction {ray.direction[0], ray.direction[1], ray.direction[2]};
    const array<float, 3> invDirection {1.0f / ray.direction[0], 1.0f / ray.direction[1], 1.0f / ray.direction[2]};
    float tmax = numeric_limits<float>::max();
    float tClosestIntersection = numeric_limits<float>::max();
    bool notFullyTraversed = true;

    while(notFullyTraversed) {
        if (node->isLeaf()) {

            //test all primitives inside the leaf
            for (auto p : node->primitives()) {
                p->intersect(ray, tClosestIntersection, intersection, tmax);
            }

            //test if leaf + empty stack => return
            if (nodeStack.empty()) {
                notFullyTraversed = false;
            } else {
                //pop all stack
                origin = mystack.top().origin;
                tmax = mystack.top().tmax;
                node = mystack.top().node;
                mystack.pop();
            }
        } else {
            //get axis of node and its split plane
            const int axis = node->axis();
            const float plane = node->splitposition();

            //test if ray is not parallel to plane
            if ((fabs(direction[axis]) > EPSILON)) {
                const float t = (plane - origin[axis]) * invDirection[axis];

                //case of the ray intersecting the plane, then test both childs
                if (0.0f < t && t < tmax) {
                    //traverse near first, then far. Set tmax = t for near
                    tmax = t;
                    //push only far child onto stack
                    mystack.push({
                                 (origin[axis] > plane )  ? node->leftChild() : node->rightChild() ,
                                 tmax - t,
                                 {origin[0] + direction[0] * t, origin[1] + direction[1] * t, origin[2] + direction[2] * t}
                                 });
                }
            }
            //in every case: traverse near child first
            node = (origin[axis] > plane) ? node->rightChild() : node->leftChild();

        }
    }

    return intersection.found;

遍历远子的频率不够。我在哪里错过了相关案例？

【问题讨论】：

首先，不要过分努力地试图通过编写您认为“更快”的代码来“击败编译器”以获得速度。由于别名问题，您可能最终让编译器的优化器不接触代码。其次，您是否使用调试器来查看问题所在以及您的代码在哪里走错了路？最后，So far I was using Havran's recursive algorithm 'B', which seems antique and overblown for OOP 那么你有没有使用这种递归算法的程序，如果有，它可以工作吗？如果是这样，请将其与您的新程序并排运行，以查看您的新程序出错的地方。
当您使用调试器时，您的代码的哪一部分导致了问题？您在发布之前使用过调试器吗？您是否尝试过更改已处理的栅格数量以帮助找出出错的地方？
另外：使用单个堆栈并推送结构，而不是尝试保持三个堆栈同步。会给你更好的地方。
@Anteru：很棒的提示。对于调试：在发布之前我没有尝试过，因为它看起来像是大海捞针。更简单的模型没有问题，而且它需要的不仅仅是少量的光线。我的调试技能不是那么好，但我尝试了最后一个小时左右的时间来跟踪它。问题是：我知道存在差异。我推远节点的频率不够高，所以我肯定经常错过这种情况。顺便说一句：图片不正确。我在我的 kd 树构造中发现了一个错误，因此它只在一个轴上分裂。这就是为什么它看起来像茶壶的一部分。
它不应该是大海捞针。第一个任务是找到它失败的最简单、最小的输入案例。然后，添加大量调试输出（或在调试器中逐步执行），以便您将进度与手动计算进行比较。这应该可以让您非常快速地确定行为偏离预期的第一个点，然后您可以从那里向后工作。

标签： c++ c++11 raytracing tree-traversal kdtree

【解决方案1】：

一个问题很小（原始的，错误的代码）：

      //traverse near first, then far. Set tmax = t for near
      tmax = t;
      //push only far child onto stack
      mystack.push({ ... , tmax - t,  ...    });

它总是将 0.0f 作为远节点的退出距离推入堆栈，这意味着交叉点不接受正 t。

交换两条线可以解决这个问题。

我的恢复堆栈跟踪/决策仍然不同（Havran 需要多约 25% 的迭代），输出图片具有 99.5% 的相同像素。那是浮点舍入问题的一部分，但它仍然没有回答这个问题：简化的实现不能识别什么情况？或者这个版本中哪些操作（数值上）不够稳定？

【讨论】：