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

/* intersect.c
 *
 * Routines to perform intersection calculations between vectors and
 * solids.
 */

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

#include "basetype.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"
#include "list.h"
#include "hitlist.h"
#include "extent.h"
#include "hitdata.h"
#include "intersect.h"


/* Intersect
 *
 * Calculate intersction of ray with solid.  Discards any intersections
 * that are too close.
 *
 * Returns TRUE if intersection exists.
 * Returns information on intersection in hit_info.
 */

bool_t
Intersect(hit_data, ray, solid)
hit_data_t	*hit_data;
ray_t		*ray;
solid_t		*solid;
{
	hit_list_t	*hit_list;
	bool_t		have_hit;

	hit_list = Intersections(ray, solid);
	have_hit = (hit_list != HitListNull);

	if (have_hit && hit_data != HitDataNull)
		HitDataCalculate(hit_data, hit_list->elem);

	HitListFree(hit_list);
	return(have_hit);
}


/* Intersections
 *
 * Generates all the intersections (save those that are too near)
 * between a ray and a solid.
 */

hit_list_t *
Intersections(ray, solid)
ray_t	*ray;
solid_t	*solid;
{
        hit_list_t      *hit_list;

        hit_list = SolidIntersect(ray, solid);

        while (hit_list != HitListNull && hit_list->elem->t < REAL_TINY)
	{
		hit_list_t	*tmp;

		tmp = hit_list;
                hit_list = hit_list->ptr;
		tmp->ptr = HitListNull;
		HitListFree(tmp);
	}

	return(hit_list);
}


/* SolidIntersect
 *
 * Calculate all intersections of a ray with a solid.
 *
 * A join will have intersections from both of its sub-solids.
 * A meet has intersections from the left sub-solid, so long as
 * they are within the right, and similarly for the right sub-solid.
 * An atom has intersections from the quadric that describes it.
 */

hit_list_t *
SolidIntersect(ray, solid)
ray_t		*ray;
solid_t		*solid;
{
	solid_body_t	*body;
	hit_list_t	*list_left, *list_right;

	if (solid->bbox != BboxNull &&
	    !BboxIntersect(ray, solid->bbox))
			return(HitListNull);

	body = &solid->body;

	switch (solid->type)
	{
	case JOIN:
		list_left  = SolidIntersect(ray, body->join.left);
		list_right = SolidIntersect(ray, body->join.right);
		return(HitListMerge(list_left, list_right));

	case MEET:
		list_left  = SolidIntersect(ray, body->meet.left);
		list_right = SolidIntersect(ray, body->meet.right);
		list_left  = FilterInside(list_left, body->meet.right);
		list_right = FilterInside(list_right, body->meet.left);
		return(HitListMerge(list_left, list_right));

	case ATOM:
		return(QuadricIntersect(ray, body->atom.quadric, solid));

	default:
		FatalError(__FILE__, __LINE__, "SolidIntersect: Invalid solid");
	}
	/*NOTREACHED*/
}


/* AddHit
 *
 * Calculates and adds the intersection at parameter t into the
 * supplied hitlist.
 */

hit_list_t *
AddHit(t, ray, solid, hit_list)
real_t		t;
ray_t		*ray;
solid_t		*solid;
hit_list_t	*hit_list;
{
	point_t		*point;
	hit_info_t	*hit_info;

	/* Check to see if it's behind */

	if (t < REAL_ZERO)
		return(hit_list);

	point = RayInstantiate(t, ray);
	hit_info = HitInfo(solid, point, t);

	return(HitListInsert(hit_info, hit_list));
}


/* QuadricIntersect
 *
 * This function performs the intersection calculations between
 * a ray and a quadric surface.
 *
 * The solution is described elsewhere.  There are three main
 * cases.
 *
 * 1) The surface is absurd - there are no intersections.
 * 2) The surface is flat   - there may be one intersection.
 * 3) The surface is curved - there may be up to two intersections.
 */

hit_list_t *
QuadricIntersect(ray, quadric, solid)
ray_t		*ray;
quadric_t	*quadric;
solid_t		*solid;
{
	real_t		a, b, c, d, e, f, g, h, i, j;
	real_t		x0, y0, z0;
	real_t		vx, vy, vz;
	real_t		A, B, C, D;
	real_t		t1, t2;
	hit_list_t	*hit_list;

	QuadricCoefs(quadric, &a, &b, &c, &d, &e, &f, &g, &h, &i, &j);
	VectorGet(ray->vector, vx, vy, vz);
	PointGet(ray->point, x0, y0, z0);

	/* Calculate t*t coefficient */

	A = a*vx*vx + e*vy*vy + h*vz*vz +
	    REAL_TWO*b*vx*vy + REAL_TWO*c*vx*vz + REAL_TWO*f*vy*vz;

	/* Calculate t coefficient */

	B = d*vx + g*vy + i*vz;

	if (x0 != REAL_ZERO) B += x0 * (a*vx + b*vy + c*vz);
	if (y0 != REAL_ZERO) B += y0 * (e*vy + f*vz + b*vx);
	if (z0 != REAL_ZERO) B += z0 * (h*vz + c*vx + f*vy);

	B *= REAL_TWO;

	/* * Calculate constant coefficient */

	C = j;

	if (x0 != REAL_ZERO) C += x0 * (a*x0 + REAL_TWO * (b*y0 + c*z0 + d));
	if (y0 != REAL_ZERO) C += y0 * (e*y0 + REAL_TWO * (f*z0 + g));
	if (z0 != REAL_ZERO) C += z0 * (h*z0 + REAL_TWO * i);

	/* Find order of equation */

	if (A == REAL_ZERO)
	{
		/* Test for constant */

		if (B == REAL_ZERO)
			return(HitListNull);

		/* Linear */

		t1 = (real_t) (-C / B);

		return(AddHit(t1, ray, solid, HitListNull));
	}
	else
	{
		/* Quadratic */

		/* Calculate discriminant */

		D = B*B - REAL_FOUR*A*C;

		/* Check there is a solution */

		if (D < REAL_ZERO)
			return(HitListNull);

		/* Calculate it */

		D = (real_t) sqrt((double) D);
		A *= REAL_TWO;
		t1 = (-B + D) / A;
		t2 = (-B - D) / A;

		/* Calculate solutions */

		hit_list = AddHit(t1, ray, solid, HitListNull);
		hit_list = AddHit(t2, ray, solid, hit_list);

		return(hit_list);
	}
}