#ifndef lint
static char *sccsid = "@(#)bbox.c	1.1 (Steve Hill) 3/9/90";
#endif

/* bbox.c
 *
 * Routines to calculate and combine bounding boxes.  A bounding box
 * consists of six planes - two perpendicular to each axis.  The
 * planes define the area within which a particular solid lies.
 *
 * Bounding boxes are used to speed up the ray tracing process.
 * They are simpler to perform intersection calculations upon than
 * general quadric equations.
 *
 * A bounding box is only calculated on demand.  Use them where there
 * are a large number of shapes in a small space.
 *
 * Restrictions:  Only ellipsoids have bounding boxes calculated
 * for them.  All other forms are considered to be unbounded.
 * Although this may not be so in practice, the general case is
 * too hard.  If you want a bounded cylinder, say, then intersect
 * it with a sphere or ellipsoid of the appropriate size ie. big
 * enough to encompass the whole shape.
 */

#include <stdio.h>
#include <math.h>


#include "basetype.h"
#include "malloc.h"
#include "cartesian.h"
#include "error.h"
#include "matrix.h"
#include "point.h"
#include "vector.h"
#include "indexlist.h"
#include "ray.h"
#include "quadric.h"
#include "colour.h"
#include "surface.h"
#include "bitmap.h"
#include "pattern.h"
#include "solid.h"
#include "bbox.h"


/* Bbox
 *
 * Allocate and initialise a bounding box structure.
 */

bbox_t *
Bbox(x_min, x_max, y_min, y_max, z_min, z_max)
real_t	x_min, x_max, y_min, y_max, z_min, z_max;
{
	bbox_t	*bbox;

	bbox = (bbox_t *) malloc(sizeof(bbox_t));
	if (bbox == BboxNull)
		FatalError(__FILE__, __LINE__, "Bbox: out of memory");

	bbox->x_min = x_min;
	bbox->y_min = y_min;
	bbox->z_min = z_min;
	bbox->x_max = x_max;
	bbox->y_max = y_max;
	bbox->z_max = z_max;

	return(bbox);
}


/* BboxCopy
 *
 * Create a duplicate of a bounding box.
 */

bbox_t *
BboxCopy(bbox)
bbox_t	*bbox;
{
	if (bbox == BboxNull)
		return(BboxNull);

	return(Bbox(bbox->x_min, bbox->x_max,
		    bbox->y_min, bbox->y_max,
		    bbox->z_min, bbox->z_max));
}


/* BboxLet
 *
 * Copy bbox2 into bbox1.
 */

bbox_t *
BboxLet(bbox1, bbox2)
bbox_t	*bbox1, *bbox2;
{
	bbox1->x_min = bbox2->x_min;
	bbox1->y_min = bbox2->y_min;
	bbox1->z_min = bbox2->z_min;
	bbox1->x_max = bbox2->x_max;
	bbox1->y_max = bbox2->y_max;
	bbox1->z_max = bbox2->z_max;

	return(bbox1);
}


/* BboxFree
 *
 * Free a bounding box.
 */

void
BboxFree(bbox)
bbox_t	*bbox;
{
	if (bbox == BboxNull)
		return;

	free((char *) bbox);
}


/* BboxPrint
 *
 * Print bounding box to specified file.
 */

void
BboxPrint(file, bbox)
FILE	*file;
bbox_t	*bbox;
{
	if (bbox == BboxNull)
	{
		fprintf(file, "Null Bbox\n");
		return;
	}
	fprintf(file, "X: %f <-> %f\n", bbox->x_min, bbox->x_max);
	fprintf(file, "Y: %f <-> %f\n", bbox->y_min, bbox->y_max);
	fprintf(file, "Z: %f <-> %f\n", bbox->z_min, bbox->z_max);
}


/* BboxJoin
 *
 * Join two bounding boxes.  The join of two boxes is the box
 * that is large enough to hold both of them.
 */

void
BboxJoin(bbox, bbox1, bbox2)
bbox_t	*bbox, *bbox1, *bbox2;
{
	bbox->x_min = RealMin(bbox1->x_min, bbox2->x_min);
	bbox->y_min = RealMin(bbox1->y_min, bbox2->y_min);
	bbox->z_min = RealMin(bbox1->z_min, bbox2->z_min);
	bbox->x_max = RealMax(bbox1->x_max, bbox2->x_max);
	bbox->y_max = RealMax(bbox1->y_max, bbox2->y_max);
	bbox->z_max = RealMax(bbox1->z_max, bbox2->z_max);
}


/* BboxMeet
 *
 * Meet two bounding boxes.  The meet of two boxes is the box
 * that encompasses their intersection.
 */

void
BboxMeet(bbox, bbox1, bbox2)
bbox_t	*bbox, *bbox1, *bbox2;
{
	bbox->x_min = RealMax(bbox1->x_min, bbox2->x_min);
	bbox->y_min = RealMax(bbox1->y_min, bbox2->y_min);
	bbox->z_min = RealMax(bbox1->z_min, bbox2->z_min);
	bbox->x_max = RealMin(bbox1->x_max, bbox2->x_max);
	bbox->y_max = RealMin(bbox1->y_max, bbox2->y_max);
	bbox->z_max = RealMin(bbox1->z_max, bbox2->z_max);
}


/* BboxAtom
 *
 * Find the bounding box of a quadric surface.
 * The routine does rather more work than is necessary, but may be
 * generalised in the future.
 */

void
BboxAtom(bbox, quadric, sort)
bbox_t		*bbox;
quadric_t	*quadric;
quadric_sort_t	sort;
{
	real_t	a, b, c, d, e, f, g, h, i, j;
	real_t	A, B, C, D, E, F, G, H, I, J;

	if (sort != ELLIPSOID)
	{
		bbox->x_min = -REAL_BIG;
		bbox->y_min = -REAL_BIG;
		bbox->z_min = -REAL_BIG;
		bbox->x_max = REAL_BIG;
		bbox->y_max = REAL_BIG;
		bbox->z_max = REAL_BIG;
		return;
	}

	QuadricCoefs(quadric, &a, &b, &c, &d, &e, &f, &g, &h, &i, &j);

	/* Calculate co-factors
	 */

	A = e*h*j - e*i*i - f*f*j + 2*f*g*i - g*g*h;
	B = -b*h*j + b*i*i + c*f*j - c*g*i - d*f*i + d*g*h;
	C = b*f*j - b*g*i - c*e*j + c*g*g + d*e*i - d*g*f;
	D = -b*f*i + b*g*h + c*e*i - c*g*f - d*e*h + d*f*f;

	E = a*h*j - a*i*i - c*c*j + 2*c*d*i - d*d*h;
	F = -a*f*j + a*g*i + c*b*j - c*d*g - d*b*i + d*d*f;
	G = a*f*i - a*g*h - c*b*i + c*c*g + d*b*h - d*c*f;
	H = a*e*j - a*g*g - b*b*j + 2*b*d*g - d*d*e;

	I = -a*e*i + a*g*f + b*b*i - b*c*g - d*b*f + d*c*e;
	J = a*e*h - a*f*f - b*b*h + 2*b*c*f - c*c*e;

	BboxQuadratic(&bbox->x_min, &bbox->x_max, J, -REAL_TWO*D, A);
	BboxQuadratic(&bbox->y_min, &bbox->y_max, J, -REAL_TWO*G, E);
	BboxQuadratic(&bbox->z_min, &bbox->z_max, J, -REAL_TWO*I, H);
}


/* BboxQuadratic
 *
 * Solve the quadratic equation:
 *
 *   2
 * ax  + bx + c = 0
 *
 * Return sorted results in min and max.  If the equation has imaginary
 * roots, then return -REAL_BIG, +REAL_BIG.
 */

void
BboxQuadratic(min, max, a, b, c)
real_t	*min, *max, a, b, c;
{
	real_t	d;

	if (a == REAL_ZERO)
	{
		if (b == REAL_ZERO)
		{
			*min = -REAL_BIG;
			*max = REAL_BIG;
		}
		else
			*min = *max = -c / b;

		return;
	}
	
	d = b*b - REAL_FOUR*a*c;

	if (d < REAL_ZERO)
	{
		*min = -REAL_BIG;
		*max = REAL_BIG;
	}
	else
	{
		real_t	sol1, sol2;

		sol1 = (-b + sqrt(d)) / (REAL_TWO * a);
		sol2 = (-b - sqrt(d)) / (REAL_TWO * a);

		RealSort(sol1, sol2);

		*min = sol1;
		*max = sol2;
	}
}


/* BboxCalculate
 *
 * Calculate bounding box for solid.
 */

void
BboxCalculate(bbox, solid)
bbox_t	*bbox;
solid_t	*solid;
{
	bbox_t		bbox1, bbox2;
	solid_body_t	*body;

	body = &solid->body;

	switch (solid->type)
	{
	case ATOM:
		BboxAtom(bbox, body->atom.quadric, body->atom.sort);
		break;

	case JOIN:
		BboxCalculate(&bbox1, body->join.left);
		BboxCalculate(&bbox2, body->join.right);
		BboxJoin(bbox, &bbox1, &bbox2);
		break;

	case MEET:
		BboxCalculate(&bbox1, body->meet.left);
		BboxCalculate(&bbox2, body->meet.right);
		BboxMeet(bbox, &bbox1, &bbox2);
		break;

	default:
		FatalError(__FILE__, __LINE__, "BboxCalculate: bad solid type");
	}

	if (solid->bounded == TRUE)
		solid->bbox = BboxCopy(bbox);
}


/* BboxIntersect
 *
 * Interesct a ray with a bounding box.  Returns TRUE if an
 * Intersection occurs, FALSE otherwise.
 */

bool_t
BboxIntersect(ray, bbox)
ray_t	*ray;
bbox_t	*bbox;
{
	vector_t	*vector;
	point_t		*point;
	real_t		t1, t2,
			tmax_x, tmin_x,
			tmax_y, tmin_y,
			tmax_z, tmin_z,
			tmax, tmin;

	vector = ray->vector;
	point  = ray->point;

	if (RealAbs(vector->x) < REAL_TINY)
	{
		if (point->x < bbox->x_min || point->x > bbox->x_max)
			return(FALSE);
		else
		{
			tmax_x = REAL_BIG;
			tmin_x = -REAL_BIG;
		}
	}
	else
	{
		t1 = (bbox->x_min - point->x) / vector->x;
		t2 = (bbox->x_max - point->x) / vector->x;

		tmax_x = RealMax(t1, t2);
		tmin_x = RealMin(t1, t2);
		if (tmax_x < REAL_ZERO)
			return(FALSE);
	}

	if (RealAbs(vector->y) < REAL_TINY)
	{
		if (point->y < bbox->y_min || point->y > bbox->y_max)
			return(FALSE);
		else
		{
			tmax_y = REAL_BIG;
			tmin_y = -REAL_BIG;
		}
	}
	else
	{
		t1 = (bbox->y_min - point->y) / vector->y;
		t2 = (bbox->y_max - point->y) / vector->y;

		tmax_y = RealMax(t1, t2);
		tmin_y = RealMin(t1, t2);
		if (tmax_y < REAL_ZERO)
			return(FALSE);
	}

	if (RealAbs(vector->z) < REAL_TINY)
	{
		if (point->z < bbox->z_min || point->z > bbox->z_max)
			return(FALSE);
		else
		{
			tmax_z = REAL_BIG;
			tmin_z = -REAL_BIG;
		}
	}
	else
	{
		t1 = (bbox->z_min - point->z) / vector->z;
		t2 = (bbox->z_max - point->z) / vector->z;

		tmax_z = RealMax(t1, t2);
		tmin_z = RealMin(t1, t2);
		if (tmax_z < REAL_ZERO)
			return(FALSE);
	}

	tmax = RealMin(RealMin(tmax_x, tmax_y), tmax_z);
	tmin = RealMax(RealMax(tmin_x, tmin_y), tmin_z);

	if (tmax < REAL_ZERO)
		return(FALSE);
	
	if (tmax < tmin)
		return(FALSE);

	return(TRUE);
}