/*
 * 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;
	}
}

const 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
	OctreeTravState *st_cur = st;

#	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;

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