algorithm - Ray-algoritmos de intersección de octree

graphics traversal

Question

Estoy buscando un buen algoritmo de intersección de octárbol de rayos, que me dé las hojas a través de las cuales pasa el rayo de manera iterativa. Estoy planeando implementarlo en la CPU, ya que no quiero sumergirme en CUDA todavía :)

En este momento, mi Voxel raycaster solo hace DDA 3D (versión Amanatides / Woo) en una matriz no jerárquica de voxels XxYxZ. Puede imaginar que esto es bastante costoso cuando hay mucho espacio vacío, como se muestra en la siguiente imagen (rojo más brillante = más trabajo :)):

Ya me he dado cuenta de que hay dos tipos de algoritmos para esta tarea: de abajo hacia arriba , que funciona desde las hojas hacia arriba y de arriba hacia abajo , que es básicamente una búsqueda en profundidad.

Ya he encontrado el algoritmo de Revelles del 2000, llamado Un algoritmo paramétrico eficiente para el recorrido de octree , que parece interesante, pero es bastante antiguo. Este es un algoritmo de arriba hacia abajo.

El enfoque de abajo hacia arriba más popular parece ser K. Sung, un algoritmo de recorrido de octetos de DDA para el trazado de rayos, Eurographics''91, North Holland-Elsevier, ISBN 0444 89096 3, pág. 73-85. El problema es que la mayoría de los algoritmos de recorrido de Octree DDA esperan que el octree sea de igual profundidad, lo que no quiero. Los subárboles vacíos deben ser solo un puntero nulo o algo similar.

En la literatura más reciente sobre Sparse Voxel Octrees que he logrado leer (en particular, el trabajo de Laine en SVO , todos parecen estar basados ​​en algún tipo de versión de DDA implementada por GPU (estilo Amanatides / Woo).

Ahora, aquí está mi pregunta: ¿Alguien tiene alguna experiencia en la implementación de un algoritmo básico de intersección de rayos y octetos? ¿Qué recomendarías?

Answer 1

La parte superior funciona muy bien para mí; la parte superior de octree puede estar basada en punteros, por lo que los grandes sub-volúmenes vacíos no toman memoria; la parte inferior es más eficiente para implementar sin punteros ... La complejidad del tiempo para golpear la pared es log2 (N) (aparentemente es el mejor caso). La implementación recursiva es bastante simple, por lo que es más fácil optimizar el código. Todas las matemáticas se pueden implementar de manera efectiva a través de operaciones SSE de números enteros: se necesitan aproximadamente 30 ciclos de CPU para calcular las nuevas coordenadas XYZ para cada salto de subvolumen. Por cierto, las versiones públicas de los recorridos de octree son buenas solo para la educación: dominar la implementación verdaderamente efectiva puede llevar varios meses fácilmente ...

Stefan

Answer 2

Para el registro, esta es mi implementación del documento de Revelles que terminé usando:

#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->size[0] - ray.origin[0];
    ray.direction[0] = - ray.direction[0];
    a |= 4 ; //bitwise OR (latest bits are XYZ)
}
if(ray.direction[1] < 0){
    ray.origin[1] = octree->size[1] - ray.origin[1];
    ray.direction[1] = - ray.direction[1];
    a |= 2 ; 
}
if(ray.direction[2] < 0){
    ray.origin[2] = octree->size[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);
}
}