雷-八叉树相交算法

Question

我正在寻找一种好的光线-八叉树相交算法，它可以以迭代的方式给出光线通过的叶子节点。目前，我的体素光线投射器只是在非分层的XxYxZ体素数组上执行三维DDA (Amanatides/Woo版本)。当空间中有很多空白区域时，这很耗费时间，如下图所示（更亮的红色=更多的工作:)）：

我已经发现了两种解决此问题的算法：自底向上，从叶子节点向上工作，和自顶向下，基本上是深度优先搜索。

我已经找到了Revelles于2000年发布的算法，称为An efficient parametric algorithm for octree traversal ，看起来很有趣，但已经过时了。这是一种自顶向下的算法。

最流行的自底向上方法似乎是K. Sung，一个DDA八叉树遍历算法，用于光线跟踪，Eurographics'91，North Holland-Elsevier，ISBN 0444 89096 3，p. 73-85。问题是，大多数DDA八叉树遍历算法都希望八叉树具有相等的深度，而我不想这样做-空子树应该只是一个空指针或类似的东西。

在我已经阅读过的关于Sparse Voxel Octrees的更近期的文献中（尤其是Laine在SVO方面的工作），它们似乎都基于一些GPU实现版本的DDA (Amanatides/Woo样式)。

现在，我有一个问题：有没有人有实现基本的、不带花哨的光线八叉树相交算法的经验？你会推荐什么呢？

- Jeroen Baert

---

那是什么，可怕的兔子食人魔？ :) - RBarryYoung

这是斯坦福大学的犰狳模型，可以在此处获取：http://graphics.stanford.edu/data/3Dscanrep/。 - Jeroen Baert

谷歌Ingo Wald（他早期的作品有惊人的基础）。我现在有点懒，找不到正确的论文。此外，在搜索和思考中删除八叉树 - 正确的术语是KD-Tree。编辑链接：http://web.archive.org/web/20091211153343/http://graphics.cs.uni-sb.de/~wald/Publications/EG2001_IRCRT/InteractiveRenderingWithCoherentRayTracing.pdf - starmole

嘿，你可能会对这个Area51提案感兴趣。像这样的问题在那里都是受欢迎的。 - Daniel Node.js

Answer 1

顺便说一下，这是我最终采用的Revelles论文实现：

#include "octree_traversal.h"

using namespace std;

unsigned char a; // because an unsigned char is 8 bits

int first_node(double tx0, double ty0, double tz0, double txm, double tym, double tzm){
unsigned char answer = 0;   // initialize to 00000000
// select the entry plane and set bits
if(tx0 > ty0){
    if(tx0 > tz0){ // PLANE YZ
        if(tym < tx0) answer|=2;    // set bit at position 1
        if(tzm < tx0) answer|=1;    // set bit at position 0
        return (int) answer;
    }
}
else {
    if(ty0 > tz0){ // PLANE XZ
        if(txm < ty0) answer|=4;    // set bit at position 2
        if(tzm < ty0) answer|=1;    // set bit at position 0
        return (int) answer;
    }
}
// PLANE XY
if(txm < tz0) answer|=4;    // set bit at position 2
if(tym < tz0) answer|=2;    // set bit at position 1
return (int) answer;
}

int new_node(double txm, int x, double tym, int y, double tzm, int z){
if(txm < tym){
    if(txm < tzm){return x;}  // YZ plane
}
else{
    if(tym < tzm){return y;} // XZ plane
}
return z; // XY plane;
}

void proc_subtree (double tx0, double ty0, double tz0, double tx1, double ty1, double tz1, Node* node){
float txm, tym, tzm;
int currNode;

if(tx1 < 0 || ty1 < 0 || tz1 < 0) return;
if(node->terminal){
    cout << "Reached leaf node " << node->debug_ID << endl;
    return;
}
else{ cout << "Reached node " << node->debug_ID << endl;}

txm = 0.5*(tx0 + tx1);
tym = 0.5*(ty0 + ty1);
tzm = 0.5*(tz0 + tz1);

currNode = first_node(tx0,ty0,tz0,txm,tym,tzm);
do{
    switch (currNode)
    {
    case 0: { 
        proc_subtree(tx0,ty0,tz0,txm,tym,tzm,node->children[a]);
        currNode = new_node(txm,4,tym,2,tzm,1);
        break;}
    case 1: { 
        proc_subtree(tx0,ty0,tzm,txm,tym,tz1,node->children[1^a]);
        currNode = new_node(txm,5,tym,3,tz1,8);
        break;}
    case 2: { 
        proc_subtree(tx0,tym,tz0,txm,ty1,tzm,node->children[2^a]);
        currNode = new_node(txm,6,ty1,8,tzm,3);
        break;}
    case 3: { 
        proc_subtree(tx0,tym,tzm,txm,ty1,tz1,node->children[3^a]);
        currNode = new_node(txm,7,ty1,8,tz1,8);
        break;}
    case 4: { 
        proc_subtree(txm,ty0,tz0,tx1,tym,tzm,node->children[4^a]);
        currNode = new_node(tx1,8,tym,6,tzm,5);
        break;}
    case 5: { 
        proc_subtree(txm,ty0,tzm,tx1,tym,tz1,node->children[5^a]);
        currNode = new_node(tx1,8,tym,7,tz1,8);
        break;}
    case 6: { 
        proc_subtree(txm,tym,tz0,tx1,ty1,tzm,node->children[6^a]);
        currNode = new_node(tx1,8,ty1,8,tzm,7);
        break;}
    case 7: { 
        proc_subtree(txm,tym,tzm,tx1,ty1,tz1,node->children[7^a]);
        currNode = 8;
        break;}
    }
} while (currNode<8);
}

void ray_octree_traversal(Octree* octree, Ray ray){
a = 0;

// fixes for rays with negative direction
if(ray.direction[0] < 0){
    ray.origin[0] = octree->center[0] * 2 - ray.origin[0];//camera origin fix
    ray.direction[0] = - ray.direction[0];
    a |= 4 ; //bitwise OR (latest bits are XYZ)
}
if(ray.direction[1] < 0){
    ray.origin[1] = octree->center[1] * 2 - ray.origin[1];
    ray.direction[1] = - ray.direction[1];
    a |= 2 ; 
}
if(ray.direction[2] < 0){
    ray.origin[2] = octree->center[2] * 2 - ray.origin[2];
    ray.direction[2] = - ray.direction[2];
    a |= 1 ; 
}

double divx = 1 / ray.direction[0]; // IEEE stability fix
double divy = 1 / ray.direction[1];
double divz = 1 / ray.direction[2];

double tx0 = (octree->min[0] - ray.origin[0]) * divx;
double tx1 = (octree->max[0] - ray.origin[0]) * divx;
double ty0 = (octree->min[1] - ray.origin[1]) * divy;
double ty1 = (octree->max[1] - ray.origin[1]) * divy;
double tz0 = (octree->min[2] - ray.origin[2]) * divz;
double tz1 = (octree->max[2] - ray.origin[2]) * divz;

if( max(max(tx0,ty0),tz0) < min(min(tx1,ty1),tz1) ){
    proc_subtree(tx0,ty0,tz0,tx1,ty1,tz1,octree->root);
}
}

- Jeroen Baert

---

在这个解决方案中，X/Z轴似乎被交换了（第0位是Z，第2位是X）。有特殊的原因吗？ - paniq

那是一段时间以前的事了，但我可能为了OpenGL兼容性而这样做过。 - Jeroen Baert

为什么要进行这样的变换？我不理解这个。镜像的起点不应该是这样的，为什么使用 size 而不是 octree->center，我认为使用 center 来计算镜像的起点是正确的，就像这样：ray.origin[0] = octree->center[0] * 2 - ray.origin[0]; - benlong

可能需要检查并加入一个小的差值。 - Jeroen Baert

如果我只想知道现有八叉树中最近的占用节点，我该如何使用这个算法？ - Exia

Answer 2

自顶向下的方法对我非常有效；八叉树的上部分可以基于指针，因此大的空子卷不会占用内存；而下部分更适合实现无指针...撞墙的时间复杂度是log2(N)（显然是最佳情况）。递归实现相当简单，因此更容易优化代码。所有数学计算都可以通过整数SSE操作有效地实现 - 计算每个子卷跳跃的新XYZ坐标大约需要x30个CPU周期。顺便说一下，公共版本的八叉树遍历仅适用于教育 - 要真正掌握有效的实现可能需要几个月的时间...

斯特凡

- Stefan

---

我不明白为什么它应该是O(log2(N))（除非您使用某些奇怪的N定义）。此外，您建议的与Revelles等人的论文相同，而OP已经提到了这一点。 - Mikola

你说的是哪个自顶向下的算法？是参数化的那个吗？ - Jeroen Baert