/*

 * octree.cc: Octree class

 *

 * This file is part of Pyrit Ray Tracer.

 *

 * Copyright 2007, 2008  Radek Brich

 *

 * Permission is hereby granted, free of charge, to any person obtaining a copy

 * of this software and associated documentation files (the "Software"), to deal

 * in the Software without restriction, including without limitation the rights

 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

 * copies of the Software, and to permit persons to whom the Software is

 * furnished to do so, subject to the following conditions:

 *

 * The above copyright notice and this permission notice shall be included in

 * all copies or substantial portions of the Software.

 *

 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN

 * THE SOFTWARE.

 */

#include "octree.h"

OctreeNode::~OctreeNode()

{

	if (isLeaf())

	{

		leaf = leaf^1; // zero leaf bit

		delete shapes;

	}

	else

		delete[] children;

}

void OctreeNode::subdivide(const BBox &bbox, int maxdepth)

{

	ShapeList *l_shapes = getShapes();

	// prepare children (this also sets this node as non-leaf)

	makeChildren();

	// evaluate centres for axes

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

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

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

	// set bounding boxes for children

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

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

	{

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

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

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

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

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

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

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

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

	}

	// distribute shapes to children

	ShapeList::iterator sh;

	unsigned int shapenum = 0;

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

	{

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

			if ((*sh)->intersect_bbox(childbb[i]))

			{

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

				shapenum++;

			}

	}

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

	|| shapenum >= 6*l_shapes->size())

	{

		// bad subdivision, revert

		delete[] children;

		setShapes(l_shapes);

		return;

	}

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

	delete l_shapes;

	// recursive subdivision

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

		if (maxdepth > 1 && getChild(i)->getShapes()->size() > 4)

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

}

void Octree::build()

{

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

	root = new OctreeNode();

	ShapeList::iterator shape;

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

		root->addShape(*shape);

	root->subdivide(bbox, max_depth);

	built = true;

}

/*******************************************************

octree traversal algorithm as described in paper

"An Efficient Parametric Algorithm for Octree Traversal"

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

see revision 37 for original recursive version

*******************************************************/

struct OctreeTravState

{

	Float tx0,ty0,tz0,tx1,ty1,tz1,txm,tym,tzm;

	OctreeNode *node;

	int next;

	OctreeTravState() {};

	OctreeTravState(

		const Float atx0, const Float aty0, const Float atz0,

		const Float atx1, const Float aty1, const Float atz1,

		const Float atxm, const Float atym, const Float atzm,

		OctreeNode *const anode, const int anext):

		tx0(atx0), ty0(aty0), tz0(atz0), tx1(atx1), ty1(aty1), tz1(atz1),

		txm(atxm), tym(atym), tzm(atzm), node(anode), next(anext) {};

};

inline const int &next_node(const Float &txm, const int &xnode,

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

{

	if (txm < tym)

	{

		if (txm < tzm)

			return xnode;

		else

			return znode;

	}

	else

	{

		if (tym < tzm)

			return ynode;

		else

			return znode;

	}

}

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

		Float &nearest_distance)

{

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

	if (!built)

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

#ifdef MSVC

	// MSVC wants constant expression here... hope it won't overflow :)

	OctreeTravState st[20];

#else

	OctreeTravState st[max_depth+1];

#endif

#	define node	st_cur->node

#	define tx0	st_cur->tx0

#	define ty0	st_cur->ty0

#	define tz0	st_cur->tz0

#	define tx1	st_cur->tx1

#	define ty1	st_cur->ty1

#	define tz1	st_cur->tz1

#	define txm	st_cur->txm

#	define tym	st_cur->tym

#	define tzm	st_cur->tzm

	int a = 0;

	Vector ro(ray.o);

	Vector rdir(ray.dir);

	if (rdir.x < 0.0)

	{

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

		rdir.x = -rdir.x;

		a |= 4;

	}

	if (rdir.y < 0.0)

	{

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

		rdir.y = -rdir.y;

		a |= 2;

	}

	if (rdir.z < 0.0)

	{

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

		rdir.z = -rdir.z;

		a |= 1;

	}

	if (rdir.x == 0.0f) rdir.x = Eps;

	if (rdir.y == 0.0f) rdir.y = Eps;

	if (rdir.z == 0.0f) rdir.z = Eps;

	rdir.x = 1.0f/rdir.x;

	rdir.y = 1.0f/rdir.y;

	rdir.z = 1.0f/rdir.z;

	tx0 = (bbox.L.x - ro.x) * rdir.x;

	tx1 = (bbox.H.x - ro.x) * rdir.x;

	ty0 = (bbox.L.y - ro.y) * rdir.y;

	ty1 = (bbox.H.y - ro.y) * rdir.y;

	tz0 = (bbox.L.z - ro.z) * rdir.z;

	tz1 = (bbox.H.z - ro.z) * rdir.z;

	if (max3(tx0,ty0,tz0) >= min3(tx1,ty1,tz1))

		return NULL;

	node = root;

	st_cur->next = -1;

	Shape *nearest_shape = NULL;

	for (;;)

	{

		if (st_cur->next == -1)

		{

			st_cur->next = 8;

			// if ray does intersect this node

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

			{

				if (node->isLeaf())

				{

					ShapeList::iterator shape;

					//register Float mindist = max3(tx0,ty0,tz0);

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

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

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

						{

							nearest_shape = *shape;

							nearest_distance = dist;

						}

					if (nearest_shape != NULL)

						return nearest_shape;

				}

				else

				{

					txm = 0.5f * (tx0+tx1);

					tym = 0.5f * (ty0+ty1);

					tzm = 0.5f * (tz0+tz1);

					// first node

					st_cur->next = 0;

					if (tx0 > ty0)

					{

						if (tx0 > tz0)

						{ // YZ

							if (tym < tx0)

								st_cur->next |= 2;

							if (tzm < tx0)

								st_cur->next |= 1;

						}

						else

						{ // XY

							if (txm < tz0)

								st_cur->next |= 4;

							if (tym < tz0)

								st_cur->next |= 2;

						}

					}

					else

					{

						if (ty0 > tz0)

						{ // XZ

							if (txm < ty0)

								st_cur->next |= 4;

							if (tzm < ty0)

								st_cur->next |= 1;

						}

						else

						{ // XY

							if (txm < tz0)

								st_cur->next |= 4;

							if (tym < tz0)

								st_cur->next |= 2;

						}

					}

				}

			}

		}

		while (st_cur->next == 8)

		{

			// pop state from stack

			if (st_cur == st)

				return NULL; // nothing to pop, finish

			--st_cur;

		}

		// push current state

		*(st_cur+1) = *st_cur;

		++st_cur;

		switch (st_cur->next)

		{

			case 0:

				tx1 = txm;

				ty1 = tym;

				tz1 = tzm;

				node = node->getChild(a);

				(st_cur-1)->next = next_node(txm, 4, tym, 2, tzm, 1);

				break;

			case 1:

				tz0 = tzm;

				tx1 = txm;

				ty1 = tym;

				node = node->getChild(1^a);

				(st_cur-1)->next = next_node(txm, 5, tym, 3, tz1, 8);

				break;

			case 2:

				ty0 = tym;

				tx1 = txm;

				tz1 = tzm;

				node = node->getChild(2^a);

				(st_cur-1)->next = next_node(txm, 6, ty1, 8, tzm, 3);

				break;

			case 3:

				ty0 = tym;

				tz0 = tzm;

				tx1 = txm;

				node = node->getChild(3^a);

				(st_cur-1)->next = next_node(txm, 7, ty1, 8, tz1, 8);

				break;

			case 4:

				tx0 = txm;

				ty1 = tym;

				tz1 = tzm;

				node = node->getChild(4^a);

				(st_cur-1)->next = next_node(tx1, 8, tym, 6, tzm, 5);

				break;

			case 5:

				tx0 = txm;

				tz0 = tzm;

				ty1 = tym;

				node = node->getChild(5^a);

				(st_cur-1)->next = next_node(tx1, 8, tym, 7, tz1, 8);

				break;

			case 6:

				tx0 = txm;

				ty0 = tym;

				tz1 = tzm;

				node = node->getChild(6^a);

				(st_cur-1)->next = next_node(tx1, 8, ty1, 8, tzm, 7);

				break;

			case 7:

				tx0 = txm;

				ty0 = tym;

				tz0 = tzm;

				node = node->getChild(7^a);

				(st_cur-1)->next = 8;

				break;

		}

		st_cur->next = -1;

	}

}