Ray / Octree traversal: Parametric algorithm implementation

This is my implementation of a Ray/Octree traversal algorithm using an adapted version of the Revelles 2000 paper An efficient parametric algorithm for octree traversal. I’ve made it compliant with IEEE double precision standard. This implementation is minimal and very bare-bones, but should be a good starting point if you want to play with Voxel raycasting.

If you plan on using this code for a project, please include a reference to this website, and my name. The details of the Creative Commons license I use can be found here.

// Ray-Octree intersection based on
// 
// An Efficient Parametric Algorithm for Octree Traversal (2000)
// by J. Revelles , C. Ureña , M. Lastra
//
// Implementation by Jeroen Baert - www.forceflow.be - [email protected]
//
// Licensed under a Creative Commons Attribute-NonCommercial Sharealike 3.0 Unported license
// which can be found at: http://creativecommons.org/licenses/by-nc-sa/3.0/
//
// TL:DR? You must attribute me, cannot use this for commercial purposes, 
// and if you share this, it should be under the same or a similar license!
//
// Thanks,
//
// Jeroen
//
 
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 (currNodesize[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);
	}
}

Comments are closed.