1 /*

     2  * Pyrit Ray Tracer

     3  * file: octree.cc

     4  *

     5  * Radek Brich, 2006-2007

     6  */

     7 

     8 #include "octree.h"

     9 

    10 OctreeNode::~OctreeNode()

    11 {

    12 	if (shapes != NULL)

    13 		delete shapes;

    14 	else

    15 		delete[] children;

    16 }

    17 

    18 void OctreeNode::subdivide(BBox bbox, int maxdepth)

    19 {

    20 	if (maxdepth <= 0 || shapes->size() <= 4)

    21 		return;

    22 

    23 	// make children

    24 	children = new OctreeNode[8];

    25 

    26 	// evaluate centres for axes

    27 	const Float xsplit = (bbox.L.x + bbox.H.x)*0.5;

    28 	const Float ysplit = (bbox.L.y + bbox.H.y)*0.5;

    29 	const Float zsplit = (bbox.L.z + bbox.H.z)*0.5;

    30 

    31 	// set bounding boxes for children

    32 	BBox childbb[8] = {bbox, bbox, bbox, bbox, bbox, bbox, bbox, bbox};

    33 	for (int i = 0; i < 4; i++)

    34 	{

    35 		// this is little obfuscated, so on right are listed affected children

    36 		// the idea is to cut every axis once per child, making 8 combinations

    37 		childbb[i].H.x = xsplit;	// 0,1,2,3

    38 		childbb[i+4].L.x = xsplit;  // 4,5,6,7

    39 		childbb[i+(i>>1<<1)].H.y = ysplit;	// 0,1,4,5

    40 		childbb[i+(i>>1<<1)+2].L.y = ysplit;// 2,3,6,7

    41 		childbb[i<<1].H.z = zsplit;		// 0,2,4,6

    42 		childbb[(i<<1)+1].L.z = zsplit; // 1,3,5,7

    43 	}

    44 

    45 	// distribute shapes to children

    46 	ShapeList::iterator sh;

    47 	BBox shbb;

    48 	int child, both;

    49 	unsigned int shapenum = 0;

    50 	for (sh = shapes->begin(); sh != shapes->end(); sh++)

    51 	{

    52 		child = 0;

    53 		both = 0;

    54 		shbb = (*sh)->get_bbox();

    55 

    56 		if (shbb.L.x >= xsplit)

    57 			child |= 4; //right

    58 		else if (shbb.H.x > xsplit)

    59 			both |= 4; // both

    60 		// for left, do nothing

    61 

    62 		if (shbb.L.y >= ysplit)

    63 			child |= 2;

    64 		else if (shbb.H.y > ysplit)

    65 			both |= 2;

    66 

    67 		if (shbb.L.z >= zsplit)

    68 			child |= 1;

    69 		else if (shbb.H.z > zsplit)

    70 			both |= 1;

    71 

    72 		if (!both)

    73 		{

    74 			getChild(child)->addShape(*sh);

    75 			shapenum++;

    76 		}

    77 		else

    78 		{

    79 			// shape goes to more than one child

    80 			if (both == 7)

    81 			{

    82 				for (int i = 0; i < 8; i++)

    83 					getChild(i)->addShape(*sh);

    84 				shapenum += 8;

    85 			}

    86 			else if (both == 3 || both >= 5)

    87 			{

    88 				if (both == 3)

    89 				{

    90 					for (int i = 0; i < 4; i++)

    91 						getChild(child + i)->addShape(*sh);

    92 				}

    93 				else if (both == 5)

    94 				{

    95 					for (int i = 0; i < 4; i++)

    96 						getChild(child + i+(i>>1<<1))->addShape(*sh);

    97 				}

    98 				else if (both == 6)

    99 				{

   100 					for (int i = 0; i < 4; i++)

   101 						getChild(child + (i<<1))->addShape(*sh);

   102 				}

   103 				shapenum += 4;

   104 			}

   105 			else

   106 			{

   107 				getChild(child)->addShape(*sh);

   108 				getChild(child+both)->addShape(*sh);

   109 				shapenum += 2;

   110 			}

   111 		}

   112 	}

   113 

   114 	if (shapes->size() <= 8 && shapenum > 2*shapes->size())

   115 	{

   116 		// bad subdivision, revert

   117 		delete[] children;

   118 		return;

   119 	}

   120 

   121 	// remove shapes and set this node to non-leaf

   122 	delete shapes;

   123 	shapes = NULL;

   124 

   125 	// recursive subdivision

   126 	for (int i = 0; i < 8; i++)

   127 		children[i].subdivide(childbb[i], maxdepth-1);

   128 }

   129 

   130 void Octree::build()

   131 {

   132 	dbgmsg(1, "* building octree\n");

   133 	root = new OctreeNode();

   134 	ShapeList::iterator shape;

   135 	for (shape = shapes.begin(); shape != shapes.end(); shape++)

   136 		root->addShape(*shape);

   137 

   138 	root->subdivide(bbox, max_depth);

   139 	built = true;

   140 }

   141 

   142 static inline int first_node(const Float tx0, const Float ty0, const Float tz0,

   143 	const Float txm, const Float tym, const Float tzm)

   144 {

   145 	int res = 0;

   146 	if (tx0 > ty0)

   147 	{

   148 		if (tx0 > tz0)

   149 		{ // YZ

   150 			if (tym < tx0)

   151 				res |= 2;

   152 			if (tzm < tx0)

   153 				res |= 1;

   154 		}

   155 		else

   156 		{ // XY

   157 			if (txm < tz0)

   158 				res |= 4;

   159 			if (tym < tz0)

   160 				res |= 2;

   161 		}

   162 	}

   163 	else

   164 	{

   165 		if (ty0 > tz0)

   166 		{ // XZ

   167 			if (txm < ty0)

   168 				res |= 4;

   169 			if (tzm < ty0)

   170 				res |= 1;

   171 			return res;

   172 		}

   173 		else

   174 		{ // XY

   175 			if (txm < tz0)

   176 				res |= 4;

   177 			if (tym < tz0)

   178 				res |= 2;

   179 		}

   180 	}

   181 	return res;

   182 }

   183 

   184 static inline int next_node(const Float txm, const int xnode,

   185 	const Float tym, const int ynode, const Float tzm, const int znode)

   186 {

   187 	if (txm < tym)

   188 	{

   189 		if (txm < tzm)

   190 			return xnode;

   191 		else

   192 			return znode;

   193 	}

   194 	else

   195 	{

   196 		if (tym < tzm)

   197 			return ynode;

   198 		else

   199 			return znode;

   200 	}

   201 }

   202 

   203 static Shape *proc_subtree(const int a, const Float tx0, const Float ty0, const Float tz0,

   204 	const Float tx1, const Float ty1, const Float tz1, OctreeNode *node,

   205 	const Shape *origin_shape, const Ray &ray, Float &nearest_distance)

   206 {

   207 	Float txm, tym, tzm;

   208 	int curr_node;

   209 

   210 	// if ray does not intersect this node

   211 	if (tx1 < 0.0 || ty1 < 0.0 || tz1 < 0.0)

   212 		return NULL;

   213 

   214 	if (node->isLeaf())

   215 	{

   216 		Shape *nearest_shape = NULL;

   217 		ShapeList::iterator shape;

   218 		Float mindist = max(max(tx0,ty0),tz0);

   219 		Float dist = min(min(min(tx1,ty1),tz1),nearest_distance);

   220 		for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++)

   221 			if (*shape != origin_shape && (*shape)->intersect(ray, dist) && dist >= mindist)

   222 			{

   223 				nearest_shape = *shape;

   224 				nearest_distance = dist;

   225 			}

   226 		return nearest_shape;

   227 	}

   228 

   229 	txm = 0.5 * (tx0+tx1);

   230 	tym = 0.5 * (ty0+ty1);

   231 	tzm = 0.5 * (tz0+tz1);

   232 

   233 	curr_node = first_node(tx0,ty0,tz0,txm,tym,tzm);

   234 	Shape *shape = NULL;

   235 	while (curr_node < 8)

   236 	{

   237 		switch (curr_node)

   238 		{

   239 			case 0:

   240 				shape =proc_subtree (a,tx0,ty0,tz0,txm,tym,tzm,node->getChild(a), origin_shape, ray, nearest_distance);

   241 				curr_node = next_node(txm, 4, tym, 2, tzm, 1);

   242 				break;

   243 			case 1:

   244 				shape =proc_subtree (a,tx0,ty0,tzm,txm,tym,tz1,node->getChild(1^a), origin_shape, ray, nearest_distance);

   245 				curr_node = next_node(txm, 5, tym, 3, tz1, 8);

   246 				break;

   247 			case 2:

   248 				shape =proc_subtree (a,tx0,tym,tz0,txm,ty1,tzm,node->getChild(2^a), origin_shape, ray, nearest_distance);

   249 				curr_node = next_node(txm, 6, ty1, 8, tzm, 3);

   250 				break;

   251 			case 3:

   252 				shape =proc_subtree (a,tx0,tym,tzm,txm,ty1,tz1,node->getChild(3^a), origin_shape, ray, nearest_distance);

   253 				curr_node = next_node(txm, 7, ty1, 8, tz1, 8);

   254 				break;

   255 			case 4:

   256 				shape =proc_subtree (a,txm,ty0,tz0,tx1,tym,tzm,node->getChild(4^a), origin_shape, ray, nearest_distance);

   257 				curr_node = next_node(tx1, 8, tym, 6, tzm, 5);

   258 				break;

   259 			case 5:

   260 				shape =proc_subtree (a,txm,ty0,tzm,tx1,tym,tz1,node->getChild(5^a), origin_shape, ray, nearest_distance);

   261 				curr_node = next_node(tx1, 8, tym, 7, tz1, 8);

   262 				break;

   263 			case 6:

   264 				shape =proc_subtree (a,txm,tym,tz0,tx1,ty1,tzm,node->getChild(6^a), origin_shape, ray, nearest_distance);

   265 				curr_node = next_node(tx1, 8, ty1, 8, tzm, 7);

   266 				break;

   267 			case 7:

   268 				shape =proc_subtree (a,txm,tym,tzm,tx1,ty1,tz1,node->getChild(7^a), origin_shape, ray, nearest_distance);

   269 				curr_node = 8;

   270 				break;

   271 		}

   272 		if (shape != NULL)

   273 			return shape;

   274 	}

   275 	return NULL;

   276 }

   277 

   278 /*

   279 traversal algorithm paper as described in paper

   280 "An Efficient Parametric Algorithm for Octree Traversal"

   281 by J. Revelles, C. Urena and M. Lastra.

   282 */

   283 Shape * Octree::nearest_intersection(const Shape *origin_shape, const Ray &ray,

   284 		Float &nearest_distance)

   285 {

   286 	/* if we have no tree, fall back to naive test */

   287 	if (!built)

   288 		return Container::nearest_intersection(origin_shape, ray, nearest_distance);

   289 

   290 	int a = 0;

   291 	Vector3 ro = ray.o;

   292 	Vector3 rdir = ray.dir;

   293 

   294 	if (rdir.x < 0.0)

   295 	{

   296 		ro.x = (bbox.L.x+bbox.H.x) - ro.x;

   297 		rdir.x = -rdir.x;

   298 		a |= 4;

   299 	}

   300 	if (rdir.y < 0.0)

   301 	{

   302 		ro.y = (bbox.L.y+bbox.H.y) - ro.y;

   303 		rdir.y = -rdir.y;

   304 		a |= 2;

   305 	}

   306 	if (rdir.z < 0.0)

   307 	{

   308 		ro.z = (bbox.L.z+bbox.H.z) - ro.z;

   309 		rdir.z = -rdir.z;

   310 		a |= 1;

   311 	}

   312 	Float tx0 = (bbox.L.x - ro.x) / rdir.x;

   313 	Float tx1 = (bbox.H.x - ro.x) / rdir.x;

   314 	Float ty0 = (bbox.L.y - ro.y) / rdir.y;

   315 	Float ty1 = (bbox.H.y - ro.y) / rdir.y;

   316 	Float tz0 = (bbox.L.z - ro.z) / rdir.z;

   317 	Float tz1 = (bbox.H.z - ro.z) / rdir.z;

   318 

   319 	//Octree *node = root;

   320 	if (max(max(tx0,ty0),tz0) < min (min(tx1,ty1),tz1))

   321 		return proc_subtree(a,tx0,ty0,tz0,tx1,ty1,tz1,root,

   322 			origin_shape, ray, nearest_distance);

   323 	else

   324 		return NULL;

   325 }