/*
 * Pyrit Ray Tracer
 * file: octree.cc
 *
 * Radek Brich, 2006-2007
 */

#include "octree.h"

OctreeNode::~OctreeNode()
{
	if (shapes != NULL)
		delete shapes;
	else
		delete[] children;
}

void OctreeNode::subdivide(BBox bbox, int maxdepth)
{
	if (maxdepth <= 0 || shapes->size() <= 4)
		return;

	// make children
	children = new OctreeNode[8];

	// evaluate centres for axes
	const Float xsplit = (bbox.L.x + bbox.H.x)*0.5;
	const Float ysplit = (bbox.L.y + bbox.H.y)*0.5;
	const Float zsplit = (bbox.L.z + bbox.H.z)*0.5;

	// 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 = shapes->begin(); sh != shapes->end(); sh++)
	{
		for (int i = 0; i < 8; i++)
			if ((*sh)->intersect_bbox(childbb[i]))
			{
				getChild(i)->addShape(*sh);
				shapenum++;
			}
	}

	if (shapes->size() <= 8 && shapenum > 2*shapes->size())
	{
		// bad subdivision, revert
		delete[] children;
		return;
	}

	// remove shapes and set this node to non-leaf
	delete shapes;
	shapes = NULL;

	// recursive subdivision
	for (int i = 0; i < 8; i++)
		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;
}

static inline int first_node(const Float tx0, const Float ty0, const Float tz0,
	const Float txm, const Float tym, const Float tzm)
{
	int res = 0;
	if (tx0 > ty0)
	{
		if (tx0 > tz0)
		{ // YZ
			if (tym < tx0)
				res |= 2;
			if (tzm < tx0)
				res |= 1;
		}
		else
		{ // XY
			if (txm < tz0)
				res |= 4;
			if (tym < tz0)
				res |= 2;
		}
	}
	else
	{
		if (ty0 > tz0)
		{ // XZ
			if (txm < ty0)
				res |= 4;
			if (tzm < ty0)
				res |= 1;
			return res;
		}
		else
		{ // XY
			if (txm < tz0)
				res |= 4;
			if (tym < tz0)
				res |= 2;
		}
	}
	return res;
}

static inline 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;
	}
}

static Shape *proc_subtree(const int a, const Float tx0, const Float ty0, const Float tz0,
	const Float tx1, const Float ty1, const Float tz1, OctreeNode *node,
	const Shape *origin_shape, const Ray &ray, Float &nearest_distance)
{
	Float txm, tym, tzm;
	int curr_node;

	// if ray does not intersect this node
	if (tx1 < 0.0 || ty1 < 0.0 || tz1 < 0.0)
		return NULL;

	if (node->isLeaf())
	{
		Shape *nearest_shape = NULL;
		ShapeList::iterator shape;
		Float mindist = max(max(tx0,ty0),tz0);
		Float dist = min(min(min(tx1,ty1),tz1),nearest_distance);
		for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++)
			if (*shape != origin_shape && (*shape)->intersect(ray, dist) && dist >= mindist)
			{
				nearest_shape = *shape;
				nearest_distance = dist;
			}
		return nearest_shape;
	}

	txm = 0.5 * (tx0+tx1);
	tym = 0.5 * (ty0+ty1);
	tzm = 0.5 * (tz0+tz1);

	curr_node = first_node(tx0,ty0,tz0,txm,tym,tzm);
	Shape *shape = NULL;
	while (curr_node < 8)
	{
		switch (curr_node)
		{
			case 0:
				shape =proc_subtree (a,tx0,ty0,tz0,txm,tym,tzm,node->getChild(a), origin_shape, ray, nearest_distance);
				curr_node = next_node(txm, 4, tym, 2, tzm, 1);
				break;
			case 1:
				shape =proc_subtree (a,tx0,ty0,tzm,txm,tym,tz1,node->getChild(1^a), origin_shape, ray, nearest_distance);
				curr_node = next_node(txm, 5, tym, 3, tz1, 8);
				break;
			case 2:
				shape =proc_subtree (a,tx0,tym,tz0,txm,ty1,tzm,node->getChild(2^a), origin_shape, ray, nearest_distance);
				curr_node = next_node(txm, 6, ty1, 8, tzm, 3);
				break;
			case 3:
				shape =proc_subtree (a,tx0,tym,tzm,txm,ty1,tz1,node->getChild(3^a), origin_shape, ray, nearest_distance);
				curr_node = next_node(txm, 7, ty1, 8, tz1, 8);
				break;
			case 4:
				shape =proc_subtree (a,txm,ty0,tz0,tx1,tym,tzm,node->getChild(4^a), origin_shape, ray, nearest_distance);
				curr_node = next_node(tx1, 8, tym, 6, tzm, 5);
				break;
			case 5:
				shape =proc_subtree (a,txm,ty0,tzm,tx1,tym,tz1,node->getChild(5^a), origin_shape, ray, nearest_distance);
				curr_node = next_node(tx1, 8, tym, 7, tz1, 8);
				break;
			case 6:
				shape =proc_subtree (a,txm,tym,tz0,tx1,ty1,tzm,node->getChild(6^a), origin_shape, ray, nearest_distance);
				curr_node = next_node(tx1, 8, ty1, 8, tzm, 7);
				break;
			case 7:
				shape =proc_subtree (a,txm,tym,tzm,tx1,ty1,tz1,node->getChild(7^a), origin_shape, ray, nearest_distance);
				curr_node = 8;
				break;
		}
		if (shape != NULL)
			return shape;
	}
	return NULL;
}

/*
traversal algorithm paper as described in paper
"An Efficient Parametric Algorithm for Octree Traversal"
by J. Revelles, C. Urena and M. Lastra.
*/
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);

	int a = 0;
	Vector3 ro = ray.o;
	Vector3 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;
	}
	Float tx0 = (bbox.L.x - ro.x) / rdir.x;
	Float tx1 = (bbox.H.x - ro.x) / rdir.x;
	Float ty0 = (bbox.L.y - ro.y) / rdir.y;
	Float ty1 = (bbox.H.y - ro.y) / rdir.y;
	Float tz0 = (bbox.L.z - ro.z) / rdir.z;
	Float tz1 = (bbox.H.z - ro.z) / rdir.z;

	if (max(max(tx0,ty0),tz0) < min (min(tx1,ty1),tz1))
		return proc_subtree(a,tx0,ty0,tz0,tx1,ty1,tz1,root,
			origin_shape, ray, nearest_distance);
	else
		return NULL;
}