Алгоритмы пересечения Ray-Octree

Question

Я ищу хороший алгоритм пересечения лучей-октетов, который дает мне листья, проходящие через иранский путь. Я планирую реализовать его на процессоре, так как я пока не хочу погружаться в CUDA :)

На данный момент мой Voxel raycaster просто выполняет 3D DDA (версия Amanatides/Woo) на неиерархическом массив из вокселей XxYxZ. Вы можете себе представить, что это довольно дорого, когда есть много свободного пространства, как показано на следующем рисунке (ярче красный = больше работы :)):

Я уже понял, что есть два вида алгоритмов для этой задачи: снизу вверх, который работает с листьев вверх, и сверху вниз, что является основным поиском глубины поиска.

Я уже нашел алгоритм Revelles от 2000 года, который называется An efficient parametric algorithm for octree traversal, который выглядит интересным, но довольно старым. Это алгоритм сверху вниз.

Самый популярный подход снизу вверх, кажется, K. Sung, DDA Octree Traversal Algorithm for Ray Tracing, Eurographics'91, North Holland-Elsevier, ISBN 0444 89096 3, p. 73-85. Проблема в том, что большинство алгоритмов обхода DDA Octree ожидают, что octree будет иметь равную глубину, чего я не хочу - пустые поддеревья должны быть просто нулевым указателем или чем-то подобным.

В более поздней литературе о разреженных VOXEL Octrees мне удалось прочитать, (прежде всего Laine's work on SVO's, все они, как представляется, на основе какой-то GPU Реализуемый версии ДВР (Amanatides/Woo стиль).

Итак, вот мой вопрос:? кто-нибудь есть опыт, реализующие основные, без излишеств, лучевое октодерево алгоритм пересечения что бы вы порекомендовали

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->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); 
} 
} 

Answer 2

сверху вниз работает очень хорошо для меня, верхняя части octree может быть указателем, так что большие пустые подтомы не занимают память, а нижняя часть более эффективна для реализации без указателей ... Сложность времени для удара по стене - log2 (N) (это очевидно y лучший случай). Рекурсивная реализация довольно проста, поэтому проще оптимизировать код. Вся математика может быть эффективно реализована с помощью целых SSE-операций - она ​​занимает около x30 циклов ЦП для вычисления новых координат XYZ для каждого подтомного прыжка. Кстати, публичные версии октодерева обходов хороши только для образования - освоить действительно эффективную реализацию может легко занять несколько месяцев ...

Стефан