/*******************************************************************************
 * bsphere.cpp
 *
 * This module implements the bounding sphere calculations.
 *
 * ---------------------------------------------------------------------------
 * Persistence of Vision Ray Tracer ('POV-Ray') version 3.7.
 * Copyright 1991-2013 Persistence of Vision Raytracer Pty. Ltd.
 *
 * POV-Ray is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * POV-Ray is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 * ---------------------------------------------------------------------------
 * POV-Ray is based on the popular DKB raytracer version 2.12.
 * DKBTrace was originally written by David K. Buck.
 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
 * ---------------------------------------------------------------------------
 * $File: //depot/public/povray/3.x/source/backend/bounding/bsphere.cpp $
 * $Revision: #1 $
 * $Change: 6069 $
 * $DateTime: 2013/11/06 11:59:40 $
 * $Author: chrisc $
 *******************************************************************************/

// frame.h must always be the first POV file included (pulls in platform config)
#include "backend/frame.h"
#include "backend/bounding/bsphere.h"
#include "backend/math/vector.h"

#include <algorithm>

// this must be the last file included
#include "base/povdebug.h"

namespace pov
{

/*****************************************************************************
* Local preprocessor defines
******************************************************************************/

const int BRANCHING_FACTOR = 4;


/*****************************************************************************
* Static functions
******************************************************************************/

static void merge_spheres (VECTOR C, DBL *r, const VECTOR C1, DBL r1, const VECTOR C2, DBL r2);
static void recompute_bound (BSPHERE_TREE *Node);
static int sort_and_split(BSPHERE_TREE **Root, BSPHERE_TREE ***Elements, int *nElem, int first, int last, int& maxelements);

/*****************************************************************************
*
* FUNCTION
*
*   merge_spheres
*
* INPUT
*
*   C, r           - Center and radius^2 of new sphere
*   C1, r1, C2, r2 - Centers and radii^2 of spheres to merge
*
* OUTPUT
*
*   C, r
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Calculate a sphere that encloses two given spheres.
*
* CHANGES
*
*   Jul 1994 : Creation.
*
*   Oct 1994 : Added test for enclosed spheres. Calculate optimal sphere. [DB]
*
******************************************************************************/

static void merge_spheres(VECTOR C, DBL *r, const VECTOR C1, DBL r1, const VECTOR C2, DBL r2)
{
	DBL l, r1r, r2r, k1, k2;
	VECTOR D;

	VSub(D, C1, C2);

	VLength(l, D);

	/* Check if one sphere encloses the other. */

	r1r = sqrt(r1);
	r2r = sqrt(r2);

	if (l + r1r <= r2r)
	{
		Assign_Vector(C, C2);

		*r = r2;

		return;
	}

	if (l + r2r <= r1r)
	{
		Assign_Vector(C, C1);

		*r = r1;

		return;
	}

	k1 = (1.0 + (r1r - r2r) / l) / 2.0;
	k2 = (1.0 + (r2r - r1r) / l) / 2.0;

	VLinComb2(C, k1, C1, k2, C2);

	*r = Sqr((l + r1r + r2r) / 2.0);
}



/*****************************************************************************
*
* FUNCTION
*
*   recompute_bound
*
* INPUT
*
*   Node - Pointer to node
*
* OUTPUT
*
*   Node
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Recompute the bounding sphere of a given node in the bounding hierarchy,
*   i. e. find the bounding sphere that encloses the bounding spheres
*   of all nodes.
*
*   NOTE: The sphere found is probably not the tightest sphere possible!
*
* CHANGES
*
*   Jul 1994 : Creation.
*
*   Oct 1994 : Improved bounding sphere calculation. [DB]
*
******************************************************************************/

static void recompute_bound(BSPHERE_TREE *Node)
{
	short i;
	DBL r2;
	VECTOR C;

	Assign_Vector(C, Node->Node[0]->C);

	r2 = Node->Node[0]->r2;

	for (i = 1; i < Node->Entries; i++)
	{
		merge_spheres(C, &r2, C, r2, Node->Node[i]->C, Node->Node[i]->r2);
	}

	Assign_Vector(Node->C, C);

	Node->r2 = r2;
}



/*****************************************************************************
*
* FUNCTION
*
*   comp_elements
*
* INPUT
*
*   in_a, in_b - elements to compare
*
* OUTPUT
*
* RETURNS
*
*   int - result of comparison
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/

template<int Axis>
int CDECL comp_elements(const void *in_a, const void *in_b)
{
	DBL am, bm;
	typedef const BSPHERE_TREE *CONST_BSPHERE_TREE_PTR;

	am = (*reinterpret_cast<const CONST_BSPHERE_TREE_PTR *>(in_a))->C[Axis];
	bm = (*reinterpret_cast<const CONST_BSPHERE_TREE_PTR *>(in_b))->C[Axis];

	if (am < bm - EPSILON)
		return (-1);
	else
	{
		if (am > bm + EPSILON)
			return (1);
		else
			return (0);
	}
}



/*****************************************************************************
*
* FUNCTION
*
*   find_axis
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/

static int find_axis(BSPHERE_TREE **Elements, int first, int  last)
{
	int which = X;
	int i;
	DBL e, d = - BOUND_HUGE;
	VECTOR C, mins, maxs;

	Make_Vector(mins,  BOUND_HUGE,  BOUND_HUGE,  BOUND_HUGE);
	Make_Vector(maxs, -BOUND_HUGE, -BOUND_HUGE, -BOUND_HUGE);

	for (i = first; i < last; i++)
	{
		Assign_Vector(C, Elements[i]->C);

		mins[X] = min(mins[X], C[X]);
		maxs[X] = max(maxs[X], C[X]);

		mins[Y] = min(mins[Y], C[Y]);
		maxs[Y] = max(maxs[Y], C[Y]);

		mins[Z] = min(mins[Z], C[Z]);
		maxs[Z] = max(maxs[Z], C[Z]);
	}

	e = maxs[X] - mins[X];

	if (e > d)
	{
		d = e;  which = X;
	}

	e = maxs[Y] - mins[Y];

	if (e > d)
	{
		d = e;  which = Y;
	}

	e = maxs[Z] - mins[Z];

	if (e > d)
	{
		which = Z;
	}

	return (which);
}



/*****************************************************************************
*
* FUNCTION
*
*   build_area_table
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Generates a table of bounding sphere surface areas.
*
* CHANGES
*
*   Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/

static void build_area_table(BSPHERE_TREE **Elements, int a, int  b, DBL *areas)
{
	int i, imin, dir;
	DBL r2;
	VECTOR C;

	if (a < b)
	{
		imin = a;  dir = 1;
	}
	else
	{
		imin = b;  dir = -1;
	}

	Assign_Vector(C, Elements[a]->C);

	r2 = Elements[a]->r2;

	for (i = a; i != (b + dir); i += dir)
	{
		merge_spheres(C, &r2, C, r2, Elements[i]->C, Elements[i]->r2);

		areas[i-imin] = r2;
	}
}



/*****************************************************************************
*
* FUNCTION
*
*   sort_and_split
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/

static int sort_and_split(BSPHERE_TREE **Root, BSPHERE_TREE ***Elements, int *nElem, int first, int last, int& maxelements)
{
	int size, i, best_loc;
	DBL *area_left, *area_right;
	DBL best_index, new_index;
	BSPHERE_TREE *cd;

	int Axis = find_axis(*Elements, first, last);

	size = last - first;

	if (size <= 0)
	{
		return (1);
	}

	/*
	 * Actually, we could do this faster in several ways. We could use a
	 * logn algorithm to find the median along the given axis, and then a
	 * linear algorithm to partition along the axis. Oh well.
	 */

	switch(Axis)
	{
		case X:
			QSORT(reinterpret_cast<void *>(*Elements + first), size, sizeof(BSPHERE_TREE *), comp_elements<X>);
			break;
		case Y:
			QSORT(reinterpret_cast<void *>(*Elements + first), size, sizeof(BSPHERE_TREE *), comp_elements<Y>);
			break;
		case Z:
			QSORT(reinterpret_cast<void *>(*Elements + first), size, sizeof(BSPHERE_TREE *), comp_elements<Z>);
			break;
	}

	/*
	 * area_left[] and area_right[] hold the surface areas of the bounding
	 * boxes to the left and right of any given point. E.g. area_left[i] holds
	 * the surface area of the bounding box containing Elements 0 through i and
	 * area_right[i] holds the surface area of the box containing Elements
	 * i through size-1.
	 */

	area_left  = reinterpret_cast<DBL *>(POV_MALLOC(size * sizeof(DBL), "blob bounding hierarchy"));
	area_right = reinterpret_cast<DBL *>(POV_MALLOC(size * sizeof(DBL), "blob bounding hierarchy"));

	/* Precalculate the areas for speed. */

	build_area_table(*Elements, first, last - 1, area_left);
	build_area_table(*Elements, last - 1, first, area_right);

	best_index = area_right[0] * (size - 3.0);

	best_loc = - 1;

	/*
	 * Find the most effective point to split. The best location will be
	 * the one that minimizes the function N1*A1 + N2*A2 where N1 and N2
	 * are the number of objects in the two groups and A1 and A2 are the
	 * surface areas of the bounding boxes of the two groups.
	 */

	for (i = 0; i < size - 1; i++)
	{
		new_index = (i + 1) * area_left[i] + (size - 1 - i) * area_right[i + 1];

		if (new_index < best_index)
		{
			best_index = new_index;
			best_loc = i + first;
		}
	}

	POV_FREE(area_left);
	POV_FREE(area_right);

	/*
	 * Stop splitting if the BRANCHING_FACTOR is reached or
	 * if splitting stops being effective.
	 */

	if ((size <= BRANCHING_FACTOR) || (best_loc < 0))
	{
		cd = reinterpret_cast<BSPHERE_TREE *>(POV_MALLOC(sizeof(BSPHERE_TREE), "blob bounding hierarchy"));

		cd->Entries = (short)size;

		cd->Node = reinterpret_cast<BSPHERE_TREE **>(POV_MALLOC(size*sizeof(BSPHERE_TREE *), "blob bounding hierarchy"));

		for (i = 0; i < size; i++)
		{
			cd->Node[i] = (*Elements)[first+i];
		}

		recompute_bound(cd);

		*Root = cd;

		if (*nElem >= maxelements)
		{
			/* Prim array overrun, increase array by 50%. */

			maxelements = 1.5 * maxelements;

			/* For debugging only. */

// TODO FIXME			Debug_Info("Reallocing elements to %d\n", maxelements);

			*Elements = reinterpret_cast<BSPHERE_TREE **>(POV_REALLOC(*Elements, maxelements * sizeof(BSPHERE_TREE *), "bounding slabs"));
		}

		(*Elements)[*nElem] = cd;

		(*nElem)++;

		return (1);
	}
	else
	{
		sort_and_split(Root, Elements, nElem, first, best_loc + 1, maxelements);

		sort_and_split(Root, Elements, nElem, best_loc + 1, last, maxelements);

		return (0);
	}
}



/*****************************************************************************
*
* FUNCTION
*
*   Build_Bounding_Sphere_Hierarchy
*
* INPUT
*
*   nElem    - number of elements in Elements
*   Elements - array containing nElem elements
*
* OUTPUT
*
*   Root     - root node of the hierarchy
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Create the bounding sphere hierarchy for a given set of elements.
*   One element consists of an element number (index into the array
*   of elements; e.g. blob components, triangles) and a bounding
*   sphere enclosing this element (center and squared radius).
*
* CHANGES
*
*   Oct 1994 : Creation. (Derived from the bounding slab creation code)
*
******************************************************************************/

void Build_Bounding_Sphere_Hierarchy(BSPHERE_TREE **Root, int nElem, BSPHERE_TREE ***Elements)
{
	int low, high;

	// This is a resonable guess at the number of elements needed.
	// This array will be reallocated as needed if it isn't.
	int maxelements = 2 * nElem;

	// Do a sort on the elements, with the end result being
	// a tree of objects sorted along the x, y, and z axes.

	if (nElem > 0)
	{
		low  = 0;
		high = nElem;

		while (sort_and_split(Root, Elements, &nElem, low, high, maxelements) == 0)
		{
			low  = high;
			high = nElem;
		}
	}
}



/*****************************************************************************
*
* FUNCTION
*
*   Destroy_Bounding_Sphere_Hierarchy
*
* INPUT
*
*   Node - Pointer to current node
*
* OUTPUT
*
*   Node
*
* RETURNS
*
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Destroy bounding sphere hierarchy.
*
* CHANGES
*
*   Aug 1994 : Creation.
*
*   Dec 1994 : Fixed memory leakage. [DB]
*
******************************************************************************/

void Destroy_Bounding_Sphere_Hierarchy(BSPHERE_TREE *Node)
{
	short i;

	if (Node != NULL)
	{
		if (Node->Entries > 0)
		{
			/* This is a node. Free sub-nodes. */

			for (i = 0; i < Node->Entries; i++)
			{
				Destroy_Bounding_Sphere_Hierarchy(Node->Node[i]);
			}

			POV_FREE(Node->Node);
		}

		POV_FREE(Node);
	}
}

}