2 #include <stack>

     5 

     6 class SortableShape

     7 {

     8 public:

     9 	Shape *shape;

    10 	BBox bbox;

    11 	short axis;

    12 	short mark;

    13 

    14 	SortableShape(Shape *aShape, short &aAxis): shape(aShape), axis(aAxis), mark(0)

    15 		{ bbox = shape->get_bbox(); };

    16 	friend bool operator<(const SortableShape& a, const SortableShape& b)

    17 		{ return a.bbox.L.cell[a.axis] < b.bbox.L.cell[b.axis]; };

    18 	void setAxis(short aAxis) { axis = aAxis; };

    19 	void setMark() { mark = 1; };

    20 	short hasMark() { return mark; };

    21 };

    22 

    23 class SortableShapeList: public vector<SortableShape>

    24 {

    25 public:

    26 	SortableShapeList(ShapeList &shapes, short axis)

    27 	{

    28 		ShapeList::iterator shape;

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

    30 			push_back(SortableShape(*shape, axis));

    31 	};

    32 };

    33 

    34 class SplitPos

    35 {

    36 public:

    37 	float pos;

    38 	int lnum, rnum;

    39 	SplitPos(): pos(0.0), lnum(0), rnum(0) {};

    40 	SplitPos(float &aPos): pos(aPos), lnum(0), rnum(0) {};

    41 	friend bool operator<(const SplitPos& a, const SplitPos& b)

    42 		{ return a.pos < b.pos; };

    43 	friend bool operator==(const SplitPos& a, const SplitPos& b)

    44 		{ return a.pos == b.pos; };

    45 };

    46 

    47 class SplitList: public vector<SplitPos>

    48 {

    49 };

    50 

    51 // stack element for kd-tree traversal

    52 struct StackElem

    53 {

    54 	KdNode* node; /* pointer to far child */

    55 	float t; /* the entry/exit signed distance */

    56 	Vector3 pb; /* the coordinates of entry/exit point */

    57 };

    58 

    59 // ----------------------------------------

   297 /* algorithm by Vlastimil Havran, Heuristic Ray Shooting Algorithms, appendix C */

   298 Shape *KdTree::nearest_intersection(const Shape *origin_shape, const Ray &ray,

   299 	float &nearest_distance)

   300 {

   301 	float a, b; /* entry/exit signed distance */

   302 	float t;    /* signed distance to the splitting plane */

   303 

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

   305 	if (!built)

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

   307 

   308 	if (!bbox.intersect(ray, a, b))

   309 		return NULL;

   310 

   311 	stack<StackElem*> st;

   312 

   313 	/* pointers to the far child node and current node */

   314 	KdNode *farchild, *node;

   315 	node = root; /* start from the kd-tree root node */

   316 

   317 	StackElem *enPt = new StackElem();

   318 	enPt->t = a; /* set the signed distance */

   319 	enPt->node = NULL;

   320 

   321 	/* distinguish between internal and external origin */

   322 	if (a >= 0.0) /* a ray with external origin */

   323 		enPt->pb = ray.a + ray.dir * a;

   324 	else /* a ray with internal origin */

   325 		enPt->pb = ray.a;

   326 

   327 	/* setup initial exit point in the stack */

   328 	StackElem *exPt = new StackElem();

   329 	exPt->t = b;

   330 	exPt->pb = ray.a + ray.dir * b;

   331 	exPt->node = NULL; /* set termination flag */

   332 

   333 	st.push(enPt);

   334 

   335 	/* loop, traverse through the whole kd-tree, until an object is intersected or ray leaves the scene */

   336 	while (node)

   337 	{

   338 		/* loop until a leaf is found */

   339 		while (!node->isLeaf())

   340 		{

   341 			/* retrieve position of splitting plane */

   342 			float splitVal = node->getSplit();

   343 			short axis = node->getAxis();

   344 

   345 			if (enPt->pb[axis] <= splitVal)

   346 			{

   347 				if (exPt->pb[axis] <= splitVal)

   348 				{ /* case N1, N2, N3, P5, Z2, and Z3 */

   349 					node = node->getLeftChild();

   350 					continue;

   351 				}

   352 				if (exPt->pb[axis] == splitVal)

   353 				{ /* case Z1 */

   354 					node = node->getRightChild();

   355 					continue;

   356 				}

   357 				/* case N4 */

   358 				farchild = node->getRightChild();

   359 				node = node->getLeftChild();

   360 			}

   361 			else

   362 			{ /* (enPt->pb[axis] > splitVal) */

   363 				if (splitVal < exPt->pb[axis])

   364 				{

   365 					/* case P1, P2, P3, and N5 */

   366 					node = node->getRightChild();

   367 					continue;

   368 				}

   369 				/* case P4*/

   370 				farchild = node->getLeftChild();

   371 				node = node->getRightChild();

   372 			}

   373 

   374 			/* case P4 or N4 . . . traverse both children */

   375 

   376 			/* signed distance to the splitting plane */

   377 			t = (splitVal - ray.a.cell[axis]) / ray.dir.cell[axis];

   378 

   379 			/* setup the new exit point */

   380 			st.push(exPt);

   381 			exPt = new StackElem();

   382 

   383 			/* push values onto the stack */

   384 			exPt->t = t;

   385 			exPt->node = farchild;

   386 			exPt->pb[axis] = splitVal;

   387 			exPt->pb[(axis+1)%3] = ray.a.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];

   388 			exPt->pb[(axis+2)%3] = ray.a.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];

   389 		} /* while */

   390 

   391 		/* current node is the leaf . . . empty or full */

   392 		/* "intersect ray with each object in the object list, discarding "

   393 		"those lying before stack[enPt].t or farther than stack[exPt].t" */

   394 		Shape *nearest_shape = NULL;

   395 		float dist = exPt->t;

   396 		ShapeList::iterator shape;

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

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

   399 			&& dist >= enPt->t)

   400 				nearest_shape = *shape;

   401 

   402 		if (nearest_shape)

   403 		{

   404 			nearest_distance = dist;

   405 			return nearest_shape;

   406 		}

   407 

   408 		/* pop from the stack */

   409 		enPt = exPt; /* the signed distance intervals are adjacent */

   410 

   411 		/* retrieve the pointer to the next node, it is possible that ray traversal terminates */

   412 		node = enPt->node;

   413 

   414 		// pop

   415 		exPt = st.top();

   416 		st.pop();

   417 	} /* while */

   418 

   419 	/* ray leaves the scene */

   420 	return NULL;

   421 }

   422