#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh 'libray/libcommon/transform.c' <<'END_OF_FILE'
X/*
X * transform.c
X *
X * Copyright (C) 1989, 1991, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: transform.c,v 4.0.1.1 91/09/29 15:37:06 cek Exp Locker: cek $
X *
X * $Log:	transform.c,v $
X * Revision 4.0.1.1  91/09/29  15:37:06  cek
X * patch1: CoordSysTransform did not detect up-Z == -1.
X * 
X * Revision 4.0  91/07/17  14:32:25  kolb
X * Initial version.
X * 
X */
X#include "common.h"
X
X/*
X * Matrices are indexed row-first; that is:
X * matrix[ROW][COLUMN]
X */
X/*
X * Allocate new structure that holds both object-to-world and
X * world-to-object space transformation structures.  It probably
X * should hold pointers to these structures.
X */
XTrans *
XTransCreate(tr, meth)
XTransRef tr;
XTransMethods *meth;
X{
X	Trans *res;
X
X	res = (Trans *)share_malloc(sizeof(Trans));
X	res->tr = tr;
X	res->methods = meth;
X	res->animated = FALSE;
X	res->assoc = (ExprAssoc *)NULL;
X	res->prev = res->next = (Trans *)NULL;
X	MatrixInit(&res->trans);
X	MatrixInit(&res->itrans);
X	return res;
X}
X
Xvoid
XTransFree(trans)
XTrans *trans;
X{
X	if (trans->tr)
X		free((voidstar)trans->tr);
X	free((voidstar)trans);
X}
X
Xvoid
XTransAssoc(trans, ptr, expr)
XTrans *trans;
XFloat *ptr;
XExpr *expr;
X{
X	ExprAssoc *assoc;
X
X	if (expr->timevary) {
X		/*
X		 * Gotta store the sucker.
X		 */
X		trans->assoc = AssocCreate(ptr, expr, trans->assoc);
X		trans->animated = TRUE;
X	} else {
X		*ptr = expr->value;
X	}
X	fflush(stderr);
X}
X
X/*
X * Allocate new transformation 'matrix'.
X */
XRSMatrix *
XMatrixCreate()
X{
X	RSMatrix *res;
X
X	res = (RSMatrix *)share_malloc(sizeof(RSMatrix));
X	MatrixInit(res);
X	return res;
X}
X
X/*
X * Multiply m1 and m2, copy result into "res".
X */
Xvoid
XMatrixMult(t1, t2, res)
XRSMatrix *t1, *t2, *res;
X{
X	register int i;
X	RSMatrix tmp;
X
X	for (i = 0; i < 3; i++) {
X		tmp.matrix[i][0] = t1->matrix[i][0] * t2->matrix[0][0] +
X			  	   t1->matrix[i][1] * t2->matrix[1][0] +
X			  	   t1->matrix[i][2] * t2->matrix[2][0];
X		tmp.matrix[i][1] = t1->matrix[i][0] * t2->matrix[0][1] +
X			  	   t1->matrix[i][1] * t2->matrix[1][1] +
X			  	   t1->matrix[i][2] * t2->matrix[2][1];
X		tmp.matrix[i][2] = t1->matrix[i][0] * t2->matrix[0][2] +
X			  	   t1->matrix[i][1] * t2->matrix[1][2] +
X		  		   t1->matrix[i][2] * t2->matrix[2][2];
X	}
X
X	tmp.translate.x = t1->translate.x * t2->matrix[0][0] +
X			  t1->translate.y * t2->matrix[1][0] +
X			  t1->translate.z * t2->matrix[2][0] + t2->translate.x;
X	tmp.translate.y = t1->translate.x * t2->matrix[0][1] +
X			  t1->translate.y * t2->matrix[1][1] +
X			  t1->translate.z * t2->matrix[2][1] + t2->translate.y;
X	tmp.translate.z = t1->translate.x * t2->matrix[0][2] +
X			  t1->translate.y * t2->matrix[1][2] +
X			  t1->translate.z * t2->matrix[2][2] + t2->translate.z;
X	MatrixCopy(&tmp, res);
X}
X
X/*
X * Return transformation information to map the "coordinate system"
X * with the given origin, "up" vector, radius, and up axis lengths to
X * one in which the "up" vector is the Z axis and the x/y/up axes
X * have unit length.  This is useful for transforming a general
X * form of a primitive into a canonical, Z-axis aligned, unit size
X * primitive, facilitating intersection testing.
X * Assumes that "up" is normalized.
X */
Xvoid
XCoordSysTransform(origin, up, r, len, trans)
XVector *origin, *up;
XFloat r, len;
XTrans *trans;
X{
X	RSMatrix tmp;
X	Vector atmp;
X
X	ScaleMatrix(r, r, len, &trans->trans);
X	if (1. - fabs(up->z) < EPSILON) {
X		atmp.x = 1.;
X		atmp.y = atmp.z = 0.;
X	} else {
X		atmp.x = up->y;
X		atmp.y = -up->x;
X		atmp.z= 0.;
X	}
X	/*
X	 * Might want to make sure that |up->z| is < 1.
X	 */
X	RotationMatrix(atmp.x, atmp.y, atmp.z, -acos(up->z), &tmp);
X	MatrixMult(&trans->trans, &tmp, &trans->trans);
X	TranslationMatrix(origin->x, origin->y, origin->z, &tmp);
X	MatrixMult(&trans->trans, &tmp, &trans->trans);
X	MatrixInvert(&trans->trans, &trans->itrans);
X}
X
Xvoid
XTransCopy(from, into)
XTrans *into, *from;
X{
X	MatrixCopy(&from->trans, &into->trans);
X	MatrixCopy(&from->itrans, &into->itrans);
X}
X
Xvoid
XTransInvert(from, into)
XTrans *into, *from;
X{
X	RSMatrix ttmp;
X	/*
X	 * In case into == from...
X	 */
X	if (from == into) {
X		ttmp = from->trans;
X		into->trans = from->itrans;
X		into->itrans = ttmp;
X	} else {
X		into->trans = from->itrans;
X		into->itrans = from->trans;
X	}
X}
X
X/*
X * Copy a given transformation structure.
X */
Xvoid
XMatrixCopy(from, into)
XRSMatrix *into, *from;
X{
X	into->matrix[0][0] = from->matrix[0][0];
X	into->matrix[0][1] = from->matrix[0][1];
X	into->matrix[0][2] = from->matrix[0][2];
X	into->matrix[1][0] = from->matrix[1][0];
X	into->matrix[1][1] = from->matrix[1][1];
X	into->matrix[1][2] = from->matrix[1][2];
X	into->matrix[2][0] = from->matrix[2][0];
X	into->matrix[2][1] = from->matrix[2][1];
X	into->matrix[2][2] = from->matrix[2][2];
X	into->translate = from->translate;
X}
X
Xvoid
XTransInit(trans)
XTrans *trans;
X{
X	MatrixInit(&trans->trans);
X	MatrixInit(&trans->itrans);
X}
X
Xvoid
XTransCompose(t1, t2, res)
XTrans *t1, *t2, *res;
X{
X	MatrixMult(&t1->trans, &t2->trans, &res->trans);
X	MatrixMult(&t2->itrans, &t1->itrans, &res->itrans);
X}
X
X/*
X * Initialize transformation structure.
X */
Xvoid
XMatrixInit(trans)
XRSMatrix *trans;
X{
X	trans->matrix[0][0] = trans->matrix[1][1] = trans->matrix[2][2] = 1.;
X	trans->matrix[0][1] = trans->matrix[0][2] = trans->matrix[1][0] =
X	trans->matrix[1][2] = trans->matrix[2][0] = trans->matrix[2][1] = 0.;
X	trans->translate.x = trans->translate.y = trans->translate.z = 0.;
X}
X
X/*
X * Calculate inverse of the given transformation structure.
X */
Xvoid
XMatrixInvert(trans, inverse)
XRSMatrix *inverse, *trans;
X{
X	RSMatrix ttmp;
X	int i;
X	Float d;
X	extern int yylineno;
X
X	ttmp.matrix[0][0] = trans->matrix[1][1]*trans->matrix[2][2] -
X			    trans->matrix[1][2]*trans->matrix[2][1];
X	ttmp.matrix[1][0] = trans->matrix[1][0]*trans->matrix[2][2] -
X			    trans->matrix[1][2]*trans->matrix[2][0];
X	ttmp.matrix[2][0] = trans->matrix[1][0]*trans->matrix[2][1] -
X			    trans->matrix[1][1]*trans->matrix[2][0];
X
X	ttmp.matrix[0][1] = trans->matrix[0][1]*trans->matrix[2][2] -
X			    trans->matrix[0][2]*trans->matrix[2][1];
X	ttmp.matrix[1][1] = trans->matrix[0][0]*trans->matrix[2][2] -
X			    trans->matrix[0][2]*trans->matrix[2][0];
X	ttmp.matrix[2][1] = trans->matrix[0][0]*trans->matrix[2][1] -
X			    trans->matrix[0][1]*trans->matrix[2][0];
X
X	ttmp.matrix[0][2] = trans->matrix[0][1]*trans->matrix[1][2] -
X			    trans->matrix[0][2]*trans->matrix[1][1];
X	ttmp.matrix[1][2] = trans->matrix[0][0]*trans->matrix[1][2] -
X			    trans->matrix[0][2]*trans->matrix[1][0];
X	ttmp.matrix[2][2] = trans->matrix[0][0]*trans->matrix[1][1] -
X			    trans->matrix[0][1]*trans->matrix[1][0];
X
X	d = trans->matrix[0][0]*ttmp.matrix[0][0] -
X	    trans->matrix[0][1]*ttmp.matrix[1][0] +
X	    trans->matrix[0][2]*ttmp.matrix[2][0];
X
X	if (fabs(d) < EPSILON*EPSILON)
X		RLerror(RL_PANIC, "Singular matrix.\n",yylineno);
X
X	ttmp.matrix[0][0] /= d;
X	ttmp.matrix[0][2] /= d;
X	ttmp.matrix[1][1] /= d;
X	ttmp.matrix[2][0] /= d;
X	ttmp.matrix[2][2] /= d;
X
X	d = -d;
X
X	ttmp.matrix[0][1] /= d;
X	ttmp.matrix[1][0] /= d;
X	ttmp.matrix[1][2] /= d;
X	ttmp.matrix[2][1] /= d;
X
X	ttmp.translate.x = -(ttmp.matrix[0][0]*trans->translate.x +
X			     ttmp.matrix[1][0]*trans->translate.y +
X			     ttmp.matrix[2][0]*trans->translate.z);
X	ttmp.translate.y = -(ttmp.matrix[0][1]*trans->translate.x +
X			     ttmp.matrix[1][1]*trans->translate.y +
X			     ttmp.matrix[2][1]*trans->translate.z);
X	ttmp.translate.z = -(ttmp.matrix[0][2]*trans->translate.x +
X			     ttmp.matrix[1][2]*trans->translate.y +
X			     ttmp.matrix[2][2]*trans->translate.z);
X
X	MatrixCopy(&ttmp, inverse);
X}
X
X/*
X * Apply a transformation to a point (translation affects the point).
X */
Xvoid
XPointTransform(vec, trans)
XVector *vec;
XRSMatrix *trans;
X{
X	Vector tmp;
X
X	tmp.x = vec->x * trans->matrix[0][0] + vec->y * trans->matrix[1][0] +
X			vec->z * trans->matrix[2][0] + trans->translate.x;
X	tmp.y = vec->x * trans->matrix[0][1] + vec->y * trans->matrix[1][1] +
X			vec->z * trans->matrix[2][1] + trans->translate.y;
X	tmp.z = vec->x * trans->matrix[0][2] + vec->y * trans->matrix[1][2] +
X			vec->z * trans->matrix[2][2] + trans->translate.z;
X	*vec = tmp;
X}
X
X/*
X * 'c1x' is the X (0th) component of the first column, and so on.
X */
Xvoid
XArbitraryMatrix(c1x, c2x, c3x, c1y, c2y, c3y, c1z, c2z, c3z, tx, ty, tz, trans)
XFloat c1x, c1y, c1z, c2x, c2y, c2z, c3x, c3y, c3z, tx, ty, tz;
XRSMatrix *trans;
X{
X	trans->matrix[0][0] = c1x;
X	trans->matrix[1][0] = c1y;
X	trans->matrix[2][0] = c1z;
X
X	trans->matrix[0][1] = c2x;
X	trans->matrix[1][1] = c2y;
X	trans->matrix[2][1] = c2z;
X
X	trans->matrix[0][2] = c3x;
X	trans->matrix[1][2] = c3y;
X	trans->matrix[2][2] = c3z;
X
X	trans->translate.x = tx;
X	trans->translate.y = ty;
X	trans->translate.z = tz;
X}
X
X/*
X * Apply transformation to a vector (translations have no effect).
X */
Xvoid
XVecTransform(vec, trans)
XVector *vec;
XRSMatrix *trans;
X{
X	Vector tmp;
X
X	tmp.x = vec->x*trans->matrix[0][0] +
X		vec->y*trans->matrix[1][0] + vec->z*trans->matrix[2][0];
X	tmp.y = vec->x*trans->matrix[0][1] +
X		vec->y*trans->matrix[1][1] + vec->z*trans->matrix[2][1];
X	tmp.z = vec->x*trans->matrix[0][2] +
X		vec->y*trans->matrix[1][2] + vec->z*trans->matrix[2][2];
X
X	*vec = tmp;
X}
X
X/*
X * Transform normal -- multiply by the transpose of the given
X * matrix (which is the inverse of the 'desired' transformation).
X */
Xvoid
XNormalTransform(norm, it)
XVector *norm;
XRSMatrix *it;
X{
X	Vector onorm;
X
X	onorm = *norm;
X
X	norm->x = onorm.x*it->matrix[0][0] + onorm.y*it->matrix[0][1] +
X				onorm.z*it->matrix[0][2];
X	norm->y = onorm.x*it->matrix[1][0] + onorm.y*it->matrix[1][1] +
X				onorm.z*it->matrix[1][2];
X	norm->z = onorm.x*it->matrix[2][0] + onorm.y*it->matrix[2][1] +
X				onorm.z*it->matrix[2][2];
X	(void)VecNormalize(norm);
X}
X
X/*
X * Transform "ray" by transforming the origin point and direction vector.
X */
XFloat
XRayTransform(ray, trans)
XRay *ray;
XRSMatrix *trans;
X{
X	PointTransform(&ray->pos, trans);
X	VecTransform(&ray->dir, trans);
X	return VecNormalize(&ray->dir);
X}
X
Xvoid
XTransPropagate(trans)
XTrans *trans;
X{
X	(*trans->methods->propagate)(trans->tr, &trans->trans, &trans->itrans);
X}
X
Xvoid
XTransResolveAssoc(trans)
XTrans *trans;
X{
X	Trans *curtrans;
X	ExprAssoc *curassoc;
X
X	for (curtrans = trans; curtrans; curtrans = curtrans->next) {
X		for (curassoc = curtrans->assoc; curassoc; curassoc = curassoc->next) {
X			*curassoc->lhs = ExprEval(curassoc->expr);
X		}
X		if (curtrans->assoc)
X			TransPropagate(curtrans);
X	}
X}
X
Xvoid
XTransComposeList(list, result)
XTrans *list, *result;
X{
X	TransCopy(list, result);
X	for (list = list->next; list; list = list->next)
X		TransCompose(list, result, result);
X}
END_OF_FILE
if test 10843 -ne `wc -c <'libray/libcommon/transform.c'`; then
    echo shar: \"'libray/libcommon/transform.c'\" unpacked with wrong size!
fi
# end of 'libray/libcommon/transform.c'
fi
if test -f 'libray/libobj/csg.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'libray/libobj/csg.c'\"
else
echo shar: Extracting \"'libray/libobj/csg.c'\" \(10585 characters\)
sed "s/^X//" >'libray/libobj/csg.c' <<'END_OF_FILE'
X/*
X * csg.c
X *
X * Copyright (C) 1991, Rod G. Bogart, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: csg.c,v 4.0 91/07/17 14:37:00 kolb Exp Locker: kolb $
X *
X * $Log:	csg.c,v $
X * Revision 4.0  91/07/17  14:37:00  kolb
X * Initial version.
X * 
X */
X#include "geom.h"
X#include "csg.h"
X
X#define csg_set_enter(l, f) ((l)->data[0].enter = (f) + 1)
X
Xstatic Methods *iCsgMethods = NULL;
Xstatic char csgName[] = "csg";
X
Xstatic int	CsgUnionInt(), CsgDifferenceInt(),
X		CsgIntersectInt();
Xstatic void	CsgHitlistCopy(), CsgSetBounds();
X
XCsg *
XCsgCreate(op)
Xint op;
X{
X	Csg *csg;
X
X	csg = (Csg *)share_malloc(sizeof(Csg));
X	csg->operator = op;
X	csg->obj1 = csg->obj2 = (Geom *)NULL;
X
X
X	switch(op) {
X		case CSG_UNION:
X			csg->intmeth = CsgUnionInt;
X			break;
X		case CSG_INTERSECT:
X			csg->intmeth = CsgIntersectInt;
X			break;
X		case CSG_DIFFERENCE:
X			csg->intmeth = CsgDifferenceInt;
X			break;
X		default:
X			RLerror(RL_ABORT, "Unknown csg op type %d?\n",op);
X	}
X	return csg;
X}
X
XMethods *
XCsgMethods()
X{
X	if (iCsgMethods== (Methods *)NULL) {
X		iCsgMethods = MethodsCreate();
X		iCsgMethods->create = (GeomCreateFunc *)CsgCreate;
X		iCsgMethods->convert = CsgConvert;
X		iCsgMethods->methods = CsgMethods;
X		iCsgMethods->name = CsgName;
X		iCsgMethods->intersect = CsgIntersect;
X		iCsgMethods->bounds = CsgBounds;
X		iCsgMethods->checkbounds = FALSE;
X		iCsgMethods->closed = TRUE;
X	}
X	return iCsgMethods;
X}
X
Xchar *
XCsgName()
X{
X	return csgName;
X}
X
Xcsg_intersect_objs(csg, ray, hit1, hit2, mindist, dist1, dist2)
XCsg *csg;
XRay *ray;
XHitList *hit1, *hit2;
XFloat mindist, *dist1, *dist2;
X{
X	int operator;
X
X	hit1->nodes = 0;
X	hit2->nodes = 0;
X	*dist1 = FAR_AWAY;
X	*dist2 = FAR_AWAY;
X	operator = csg->operator;
X
X	if (!intersect(csg->obj1, ray, hit1, mindist, dist1) &&
X	    ((operator == CSG_INTERSECT) || (operator == CSG_DIFFERENCE))) {
X		/*
X		 * Intersection and Difference cases: if you miss the first
X		 * object, you missed the whole thing.
X		 */
X		return FALSE;
X	}
X
X	if (!intersect(csg->obj2, ray, hit2, mindist, dist2) &&
X	    ((operator == CSG_INTERSECT) ||
X	     (hit1->nodes == 0) && (operator == CSG_UNION))) {
X		/*
X		 * Intersect case:  if you miss either object, you miss whole
X		 * Union case: if you miss both object, you miss whole
X		 */
X		return FALSE;
X	}
X
X	return TRUE;
X}
X
Xint
Xcsg_enter_obj(hitp)
XHitList *hitp;
X{
X	if (hitp->data[0].enter)
X		return hitp->data[0].enter - 1;
X
X	return PrimEnter(hitp->data[0].obj, &hitp->data[0].ray,
X			hitp->data[0].mindist, hitp->data[0].dist);
X}
X
Xstatic int
XCsgUnionInt(ray, hit1p, hit2p, dist1, dist2, hitclose, distclose)
XRay *ray;
XHitList *hit1p, *hit2p, **hitclose;
XFloat dist1, dist2, *distclose;
X{
X	Float distnext;
X	HitList hitnext, *hittmp;
X
X	while (TRUE) {
X		if (hit2p->nodes == 0 ||
X		    csg_enter_obj(hit2p)) {
X			/* return hit1 */
X			*hitclose = hit1p;
X			*distclose = dist1;
X			csg_set_enter(hit1p, csg_enter_obj(hit1p));
X			return TRUE;
X		} else {
X			distnext = FAR_AWAY;
X			hitnext.nodes = 0;
X			if (!intersect(hit1p->data[hit1p->nodes-1].obj,
X			    ray, &hitnext, dist2+EPSILON, &distnext)) {
X				/*
X				 * None of obj1 beyond, return hit2 (leaving)
X				 */
X				*hitclose = hit2p;
X				*distclose = dist2;
X				csg_set_enter(hit2p, FALSE);
X				return TRUE;
X			} else {
X				/*
X				 * Since hit1 is supposed to be the close one,
X				 * swap them and copy hitnext into hit2.
X	     			 */
X				hittmp = hit1p;
X				hit1p = hit2p;
X				hit2p = hittmp;
X				dist1 = dist2;
X				CsgHitlistCopy(&hitnext, hit2p);
X				dist2 = distnext;
X				/* and continue */
X			}
X		}
X	}
X}
X
Xstatic int
XCsgIntersectInt(ray, hit1p, hit2p, dist1, dist2, hitclose, distclose)
XRay *ray;
XHitList *hit1p, *hit2p, **hitclose;
XFloat dist1, dist2, *distclose;
X{
X	HitList *hittmp, hitnext;
X	Float distnext;
X
X	while (TRUE) {
X		if (!csg_enter_obj(hit2p)) {
X			/* Ray is leaving obj2 */
X			/* Return hit1 info */
X			*hitclose = hit1p;
X			*distclose = dist1;
X			csg_set_enter(hit1p, csg_enter_obj(hit1p));
X			return TRUE;
X		} else {
X			distnext = FAR_AWAY;
X			hitnext.nodes = 0;
X			if (!intersect(hit1p->data[hit1p->nodes-1].obj,
X			    ray, &hitnext, dist2+EPSILON, &distnext)) {
X				/*
X				 * None of obj1 beyond, so return miss
X				 */
X				return FALSE;
X			} else {
X				/*
X				 * Since hit1 is supposed to be the
X				 * close one, swap them and copy
X				 * hitnext into hit2.
X				 */
X				hittmp = hit1p;
X				hit1p = hit2p;
X				hit2p = hittmp;
X				dist1 = dist2;
X				CsgHitlistCopy(&hitnext, hit2p);
X				dist2 = distnext;
X				/* and continue */
X			}
X		}
X	}
X}
X
Xstatic int
XCsgDifferenceInt(ray, hit1p, hit2p, dist1, dist2, hitclose, distclose)
XRay *ray;
XHitList *hit1p, *hit2p, **hitclose;
XFloat dist1, dist2, *distclose;
X{
X	Float distnext;
X	HitList hitnext;
X
X	while (TRUE) {
X		if (dist1 < dist2) {
X			if (hit2p->nodes == 0 ||
X			    csg_enter_obj(hit2p)) {
X				/* return hit1 */
X				*hitclose = hit1p;
X				*distclose = dist1;
X				csg_set_enter(hit1p, csg_enter_obj(hit1p));
X				return TRUE;
X			} else {
X				distnext = FAR_AWAY;
X				hitnext.nodes = 0;
X				if (!intersect(hit1p->data[hit1p->nodes-1].obj,
X				    ray, &hitnext, dist2+EPSILON, &distnext)) {
X					/*
X					 * None of obj1 beyond, so
X					 * return miss
X					 */
X					return FALSE;
X				} else {
X					dist1 = distnext;
X					CsgHitlistCopy(&hitnext, hit1p);
X					/* and continue */
X				}
X			}
X		} else { /* dist1 <= dist2 */
X			if (hit1p->nodes == 0) {
X				/* return a miss */
X				return FALSE;
X			}
X			if (!csg_enter_obj(hit1p)) {
X				/*
X				 * return hit2, but invert hit2
X				 * Enter/Leave flag
X				 */
X				*hitclose = hit2p;
X				*distclose = dist2;
X				csg_set_enter(hit2p, !csg_enter_obj(hit2p));
X				return TRUE;
X			} else {
X				distnext = FAR_AWAY;
X				hitnext.nodes = 0;
X				if (!intersect(hit2p->data[hit2p->nodes-1].obj,
X				    ray, &hitnext, dist1+EPSILON, &distnext)) {
X					/*
X					 * None of obj2 beyond, so
X					 * return hit1
X					 */
X					*hitclose = hit1p;
X					*distclose = dist1;
X					/* we know we're entering obj1 */
X					csg_set_enter(hit1p, TRUE);
X					return TRUE;
X				} else {
X					dist2 = distnext;
X					CsgHitlistCopy(&hitnext, hit2p);
X					/* and continue */
X				}
X			}
X		}
X	}
X}
X
Xint
XCsgIntersect(csg, ray, hitlist, mindist, maxdist)
XCsg *csg;
XRay *ray;
XHitList *hitlist;
XFloat mindist, *maxdist;
X{
X	Float dist1, dist2, disttmp, distclose;
X	HitList hit1, hit2, *hit1p, *hit2p, *hitclose;
X
X	hit1p = &hit1;
X	hit2p = &hit2;
X	if (!csg_intersect_objs(csg, ray, hit1p, hit2p, mindist,
X	    &dist1, &dist2)) {
X		/* missed the csg object */
X		return FALSE;
X	}
X
X	if ((dist1 > dist2) &&
X	    (csg->operator == CSG_UNION || csg->operator == CSG_INTERSECT)) {
X		/* swap so 1 is closest (except in difference case) */
X		disttmp = dist2;  
X		dist2 = dist1;  
X		dist1 = disttmp;
X		hit1p = &hit2;  
X		hit2p = &hit1;
X	}
X
X	/*
X	 * Call appropriate intersection method.  If FALSE is return,
X	 * no hit of any kind was found.
X	 */
X	if (!(*csg->intmeth)(ray, hit1p, hit2p, dist1, dist2,
X	    &hitclose, &distclose))
X		return FALSE;
X
X	/*
X	 * At this time, the closest hit is in hitclose and
X	 * distclose.
X	 */
X	if (distclose < mindist || distclose > *maxdist)
X		return FALSE;
X
X	CsgHitlistCopy(hitclose, hitlist);
X	*maxdist = distclose;
X	return TRUE;
X}
X
Xstatic void
XCsgHitlistCopy(from, to)
XHitList *from, *to;
X{
X	int i;
X
X	to->nodes = from->nodes;
X	for (i = 0; i < from->nodes; i++)
X		to->data[i] = from->data[i];
X}
X
Xstatic void
XCsgSetBounds(csg, bounds)
XCsg *csg;
XFloat bounds[2][3];
X{
X	GeomComputeBounds(csg->obj1);
X	GeomComputeBounds(csg->obj2);
X
X	switch (csg->operator) {
X	case CSG_UNION:
X		bounds[LOW][X] = min(csg->obj1->bounds[LOW][X], csg->obj2->bounds[LOW][X]);
X		bounds[HIGH][X] = max(csg->obj1->bounds[HIGH][X], csg->obj2->bounds[HIGH][X]);
X		bounds[LOW][Y] = min(csg->obj1->bounds[LOW][Y], csg->obj2->bounds[LOW][Y]);
X		bounds[HIGH][Y] = max(csg->obj1->bounds[HIGH][Y], csg->obj2->bounds[HIGH][Y]);
X		bounds[LOW][Z] = min(csg->obj1->bounds[LOW][Z], csg->obj2->bounds[LOW][Z]);
X		bounds[HIGH][Z] = max(csg->obj1->bounds[HIGH][Z], csg->obj2->bounds[HIGH][Z]);
X		break;
X	case CSG_INTERSECT:
X		bounds[LOW][X] = max(csg->obj1->bounds[LOW][X], csg->obj2->bounds[LOW][X]);
X		bounds[HIGH][X] = min(csg->obj1->bounds[HIGH][X], csg->obj2->bounds[HIGH][X]);
X		bounds[LOW][Y] = max(csg->obj1->bounds[LOW][Y], csg->obj2->bounds[LOW][Y]);
X		bounds[HIGH][Y] = min(csg->obj1->bounds[HIGH][Y], csg->obj2->bounds[HIGH][Y]);
X		bounds[LOW][Z] = max(csg->obj1->bounds[LOW][Z], csg->obj2->bounds[LOW][Z]);
X		bounds[HIGH][Z] = min(csg->obj1->bounds[HIGH][Z], csg->obj2->bounds[HIGH][Z]);
X		break;
X	case CSG_DIFFERENCE:
X		bounds[LOW][X] = csg->obj1->bounds[LOW][X];
X		bounds[HIGH][X] = csg->obj1->bounds[HIGH][X];
X		bounds[LOW][Y] = csg->obj1->bounds[LOW][Y];
X		bounds[HIGH][Y] = csg->obj1->bounds[HIGH][Y];
X		bounds[LOW][Z] = csg->obj1->bounds[LOW][Z];
X		bounds[HIGH][Z] = csg->obj1->bounds[HIGH][Z];
X		break;
X	default:
X		RLerror(RL_ABORT, "Unknown csg operator type %d?\n",
X			csg->operator);
X	}
X}
X
X/*
X * Return index of hitlist node closest to the leaf and not below any
X * CSG object.
X */
Xint
XFirstCSGGeom(hitlist)
XHitList *hitlist;
X{
X	int i;
X
X	/*
X	 * UUUUGLY -- detect if obj is a CsgGeom by comparing
X	 * methods with iCsgMethods.
X	 */
X	for (i = hitlist->nodes -1; i; i--)
X		if (hitlist->data[i].obj->methods == iCsgMethods)
X			return i;
X	return 0;
X}
X
Xvoid
XCsgBounds(csg, bounds)
XCsg *csg;
XFloat bounds[2][3];
X{
X	CsgSetBounds(csg, csg->bounds);
X	BoundsCopy(csg->bounds, bounds);
X}
X
Xint
XCsgConvert(csg, list)
XCsg *csg;
XGeom *list;
X{
X	static int OpenAdvised = FALSE;
X	/*
X	 * Currently, this only handles two objects.
X	 * Will be fixed in the future.
X	 * No really we promise.
X	 */
X	if (!list || !list->next) {
X		RLerror(RL_WARN, "CSG needs at least two objects.\n");
X		return 0;
X	}
X	if (list->next->next) {
X		RLerror(RL_WARN, "Currently, CSG only handles two objects.\n");
X		return 0;
X	}
X	/*
X	 * Things are put into lists backwards....
X	 */
X	csg->obj2 = list;
X	csg->obj1 = list->next;
X	if ((!csg->obj1->methods->closed || !csg->obj2->methods->closed) &&
X	    !OpenAdvised) {
X		RLerror(RL_ADVISE,
X			"Performing CSG with non-closed object(s).\n");
X		OpenAdvised = TRUE;
X	}
X	return csg->obj1->prims + csg->obj2->prims;
X}
X
Xvoid
XCsgMethodRegister(meth)
XUserMethodType meth;
X{
X	if (iCsgMethods)
X		iCsgMethods->user = meth;
X}
END_OF_FILE
if test 10585 -ne `wc -c <'libray/libobj/csg.c'`; then
    echo shar: \"'libray/libobj/csg.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/csg.c'
fi
if test -f 'libray/libobj/grid.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'libray/libobj/grid.c'\"
else
echo shar: Extracting \"'libray/libobj/grid.c'\" \(11807 characters\)
sed "s/^X//" >'libray/libobj/grid.c' <<'END_OF_FILE'
X/*
X * grid.c
X *
X * Copyright (C) 1989, 1991, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: grid.c,v 4.0.1.1 91/10/04 15:55:37 cek Exp Locker: cek $
X *
X * $Log:	grid.c,v $
X * Revision 4.0.1.1  91/10/04  15:55:37  cek
X * patch1: Removed straight floating point comparisons.
X * 
X * Revision 4.0  91/07/17  14:38:02  kolb
X * Initial version.
X * 
X */
X#include "geom.h"
X#include "grid.h"
X
Xstatic Methods *iGridMethods = NULL;
Xstatic char gridName[] = "grid";
X
Xstatic unsigned long raynumber = 1;		/* Current "ray number". */
X						/* (should be "grid number") */
Xstatic void engrid(), GridFreeVoxels();
Xstatic int pos2grid(), CheckVoxel();
X
XGrid *
XGridCreate(x, y, z)
Xint x, y, z;
X{
X	Grid *grid;
X
X	if (x < 1 || y < 1 || z < 1) {
X		RLerror(RL_WARN, "Invalid grid specification.\n");
X		return (Grid *)NULL;
X	}
X	grid = (Grid *)share_calloc(1, sizeof(Grid));
X	grid->xsize = x;
X	grid->ysize = y;
X	grid->zsize = z;
X	return grid;	
X}
X
Xchar *
XGridName()
X{
X	return gridName;
X}
X
X/*
X * Intersect ray with grid, returning distance from "pos" to
X * nearest intersection with an object in the grid.  Returns 0.
X * if no intersection.
X */
Xint
XGridIntersect(grid, ray, hitlist, mindist, maxdist)
XGrid *grid;
XRay *ray;
XHitList *hitlist;
XFloat mindist, *maxdist;
X{
X	GeomList *list;
X	Geom *obj;
X	int hit;
X	Float offset, tMaxX, tMaxY, tMaxZ;
X	Float tDeltaX, tDeltaY, tDeltaZ, *raybounds[2][3];
X	int stepX, stepY, stepZ, outX, outY, outZ, x, y, z;
X	Vector curpos, nXp, nYp, nZp, np, pDeltaX, pDeltaY, pDeltaZ;
X	unsigned long counter;
X
X	hit = FALSE;
X	/*
X	 * Check unbounded objects.
X	 */
X	for (obj = grid->unbounded ; obj; obj = obj->next) {
X		if (intersect(obj, ray, hitlist, mindist, maxdist))
X			hit = TRUE;
X	}
X
X	/*
X	 * If outside of the bounding box, check to see if we
X	 * hit it.
X	 */
X	VecAddScaled(ray->pos, mindist, ray->dir, &curpos);
X	if (OutOfBounds(&curpos, grid->bounds)) {
X		offset = *maxdist;
X		if (!BoundsIntersect(ray, grid->bounds, mindist, &offset))
X			/*
X			 * Ray never hit grid space.
X			 */
X			return hit;
X		/*
X		 * else
X		 *	The ray enters voxel space before it hits
X		 * 	an unbounded object.
X		 */
X		VecAddScaled(ray->pos, offset, ray->dir, &curpos);
X	} else
X		offset = mindist;
X
X	counter = raynumber++;
X
X	/*
X	 * tMaxX is the absolute distance from the ray origin we must move
X	 *		to get to the next voxel in the X
X	 *		direction.  It is incrementally updated
X	 * 		by DDA as we move from voxel-to-voxel.
X	 * tDeltaX is the relative amount along the ray we move to
X	 *		get to the next voxel in the X direction. Thus,
X	 *		when we decide to move in the X direction,
X	 * 		we increment tMaxX by tDeltaX.
X	 */
X	x = x2voxel(grid, curpos.x);
X	if (x == grid->xsize)
X		x--;
X	if (fabs(ray->dir.x) < EPSILON) {
X		tMaxX = FAR_AWAY;
X		raybounds[LOW][X] = &curpos.x;
X		raybounds[HIGH][X] = &np.x;
X		tDeltaX = 0.;
X	} else if (ray->dir.x < 0.) {
X		tMaxX = offset + (voxel2x(grid, x) - curpos.x) / ray->dir.x;
X		tDeltaX = grid->voxsize[X] / - ray->dir.x;
X		stepX = outX = -1;
X		raybounds[LOW][X] = &np.x;
X		raybounds[HIGH][X] = &curpos.x;
X	} else {
X		tMaxX = offset + (voxel2x(grid, x+1) - curpos.x) / ray->dir.x;
X		tDeltaX = grid->voxsize[X] / ray->dir.x;
X		stepX = 1;
X		outX = grid->xsize;
X		raybounds[LOW][X] = &curpos.x;
X		raybounds[HIGH][X] = &np.x;
X	}
X
X	y = y2voxel(grid, curpos.y);
X	if (y == grid->ysize)
X		y--;
X
X	if (fabs(ray->dir.y) < EPSILON) {
X		tMaxY = FAR_AWAY;
X		raybounds[LOW][Y] = &curpos.y;
X		raybounds[HIGH][Y] = &np.y;
X		tDeltaY = 0.;
X	} else if (ray->dir.y < 0.) {
X		tMaxY = offset + (voxel2y(grid, y) - curpos.y) / ray->dir.y;
X		tDeltaY = grid->voxsize[Y] / - ray->dir.y;
X		stepY = outY = -1;
X		raybounds[LOW][Y] = &np.y;
X		raybounds[HIGH][Y] = &curpos.y;
X	} else {
X		tMaxY = offset + (voxel2y(grid, y+1) - curpos.y) / ray->dir.y;
X		tDeltaY = grid->voxsize[Y] / ray->dir.y;
X		stepY = 1;
X		outY = grid->ysize;
X		raybounds[LOW][Y] = &curpos.y;
X		raybounds[HIGH][Y] = &np.y;
X	}
X
X	z = z2voxel(grid, curpos.z);
X	if (z == grid->zsize)
X		z--;
X	if (fabs(ray->dir.z) < EPSILON) {
X		tMaxZ = FAR_AWAY;
X		raybounds[LOW][Z] = &curpos.z;
X		raybounds[HIGH][Z] = &np.z;
X		tDeltaZ = 0.;
X	} else if (ray->dir.z < 0.) {
X		tMaxZ = offset + (voxel2z(grid, z) - curpos.z) / ray->dir.z;
X		tDeltaZ = grid->voxsize[Z] / - ray->dir.z;
X		stepZ = outZ = -1;
X		raybounds[LOW][Z] = &np.z;
X		raybounds[HIGH][Z] = &curpos.z;
X	} else {
X		tMaxZ = offset + (voxel2z(grid, z+1) - curpos.z) / ray->dir.z;
X		tDeltaZ = grid->voxsize[Z] / ray->dir.z;
X		stepZ = 1;
X		outZ = grid->zsize;
X		raybounds[LOW][Z] = &curpos.z;
X		raybounds[HIGH][Z] = &np.z;
X	}
X
X	VecScale(tDeltaX, ray->dir, &pDeltaX);
X	VecScale(tDeltaY, ray->dir, &pDeltaY);
X	VecScale(tDeltaZ, ray->dir, &pDeltaZ);
X
X	VecAddScaled(ray->pos, tMaxX, ray->dir, &nXp);
X	VecAddScaled(ray->pos, tMaxY, ray->dir, &nYp);
X	VecAddScaled(ray->pos, tMaxZ, ray->dir, &nZp);
X
X	while (TRUE) {
X		list = grid->cells[x][y][z];
X		if (tMaxX < tMaxY && tMaxX < tMaxZ) {
X			if (list) {
X				np = nXp;
X			    	if (CheckVoxel(list,ray,raybounds,
X			    	    hitlist,counter,offset,maxdist))
X					hit = TRUE;
X			}
X			x += stepX;
X			if (*maxdist < tMaxX || x == outX)
X				break;
X			tMaxX += tDeltaX;
X			curpos = nXp;
X			nXp.x += pDeltaX.x;
X			nXp.y += pDeltaX.y;
X			nXp.z += pDeltaX.z;
X		} else if (tMaxZ < tMaxY) {
X			if (list) {
X				np = nZp;
X			    	if (CheckVoxel(list,ray, raybounds,
X			    	    hitlist,counter,offset,maxdist))
X					hit = TRUE;
X			}
X			z += stepZ;
X			if (*maxdist < tMaxZ || z == outZ)
X				break;
X			tMaxZ += tDeltaZ;
X			curpos = nZp;
X			nZp.x += pDeltaZ.x;
X			nZp.y += pDeltaZ.y;
X			nZp.z += pDeltaZ.z;
X		} else {
X			if (list) {
X				np = nYp;
X			    	if (CheckVoxel(list,ray,raybounds,
X			    	    hitlist,counter,offset,maxdist))
X					hit = TRUE;
X			}
X			y += stepY;
X			if (*maxdist < tMaxY || y == outY)
X				break;
X			tMaxY += tDeltaY;
X			curpos = nYp;
X			nYp.x += pDeltaY.x;
X			nYp.y += pDeltaY.y;
X			nYp.z += pDeltaY.z;
X		}
X	}
X	return hit;
X}
X
X/*
X * Intersect ray with objects in grid cell.  Note that there are a many ways
X * to speed up this routine, all of which uglify the code to a large extent.
X */
Xstatic int
XCheckVoxel(list,ray,raybounds,hitlist,counter,mindist,maxdist)
XGeomList *list;
XRay *ray;
XFloat *raybounds[2][3];
XHitList *hitlist;
Xunsigned long counter;
XFloat mindist, *maxdist;
X{
X	Geom *obj;
X	int hit;
X	Float lx, hx, ly, hy, lz, hz;
X
X	lx = *raybounds[LOW][X];
X	hx = *raybounds[HIGH][X];
X	ly = *raybounds[LOW][Y];
X	hy = *raybounds[HIGH][Y];
X	lz = *raybounds[LOW][Z];
X	hz = *raybounds[HIGH][Z];
X
X	hit = FALSE;
X
X	do {
X		obj = list->obj;
X		/*
X		 * If object's counter is greater than or equal to the
X		 * number associated with the current grid,
X		 * don't bother checking again.  In addition, if the
X		 * bounding box of the ray's extent in the voxel does
X		 * not intersect the bounding box of the object, don't bother.
X		 */
X#ifdef SHAREDMEM
X		if (*obj->counter < counter &&
X#else
X		if (obj->counter < counter &&
X#endif
X		    obj->bounds[LOW][X] <= hx  &&
X		    obj->bounds[HIGH][X] >= lx &&
X		    obj->bounds[LOW][Y] <= hy  &&
X		    obj->bounds[HIGH][Y] >= ly &&
X		    obj->bounds[LOW][Z] <= hz  &&
X		    obj->bounds[HIGH][Z] >= lz) {
X#ifdef SHAREDMEM
X			*obj->counter = counter;
X#else
X			obj->counter = counter;
X#endif
X			if (intersect(obj, ray, hitlist, mindist, maxdist))
X				hit = TRUE;
X		}
X	} while ((list = list->next) != (GeomList *)0);
X
X	return hit;
X}
X
Xint
XGridConvert(grid, objlist)
XGrid *grid;
XGeom *objlist;
X{
X	int num;
X
X	/*
X	 * Keep linked list of all bounded objects in grid; it may come
X	 * in handy.
X	 */
X	grid->objects = objlist;
X	for (num = 0; objlist; objlist = objlist->next)
X		num += objlist->prims;
X
X	return num;
X}
X
Xvoid
XGridBounds(grid, bounds)
XGrid *grid;
XFloat bounds[2][3];
X{
X	Geom *obj, *ltmp;
X	int x, y, i;
X
X	BoundsInit(bounds);
X	/*
X	 * For each object on the list,
X	 * compute its bounds...
X	 */
X	/*
X	 * Find bounding box of bounded objects and get list of
X	 * unbounded objects.
X	 */
X	grid->unbounded = GeomComputeAggregateBounds(&grid->objects,
X				grid->unbounded, grid->bounds);
X
X	BoundsCopy(grid->bounds, bounds);
X
X	grid->voxsize[X] = (grid->bounds[HIGH][X]-grid->bounds[LOW][X])/
X				grid->xsize;
X	grid->voxsize[Y] = (grid->bounds[HIGH][Y]-grid->bounds[LOW][Y])/
X				grid->ysize;
X	grid->voxsize[Z] = (grid->bounds[HIGH][Z]-grid->bounds[LOW][Z])/
X				grid->zsize;
X
X	if (grid->cells == (GeomList ****)NULL) {
X		/*
X	 	 * Allocate voxels.
X	 	 */
X		grid->cells = (GeomList ****)share_malloc(grid->xsize *
X					sizeof(GeomList ***));
X		for (x = 0; x < grid->xsize; x++) {
X			grid->cells[x] = (GeomList ***)share_malloc(grid->ysize *
X				sizeof(GeomList **));
X			for (y = 0; y < grid->ysize; y++)
X				grid->cells[x][y] = (GeomList **)share_calloc(
X					(unsigned)grid->zsize,sizeof(GeomList *));
X		}
X	} else {
X		/*
X		 * New frame...
X		 * Free up the objlists in each voxel.
X		 */
X		GridFreeVoxels(grid);
X	}
X
X	/*
X	 * objlist now holds a linked list of bounded objects.
X	 */
X	for (ltmp = grid->objects; ltmp != (Geom *)0; ltmp = ltmp->next)
X		engrid(ltmp, grid);
X}
X
Xstatic void
XGridFreeVoxels(grid)
XGrid *grid;
X{
X	int x, y, z;
X	GeomList *cell, *next;
X
X	for (x = 0; x < grid->xsize; x++) {
X		for (y = 0; y < grid->ysize; y++) {
X			for (z = 0; z < grid->zsize; z++) {
X				for (cell = grid->cells[x][y][z]; cell; cell = next) {
X					next = cell->next;
X					free((voidstar)cell);
X				}
X				grid->cells[x][y][z] = (GeomList *)NULL;
X			}
X		}
X	}
X}
X
XMethods *
XGridMethods()
X{
X	if (iGridMethods == (Methods *)NULL) {
X		iGridMethods = MethodsCreate();
X		iGridMethods->methods = GridMethods;
X		iGridMethods->create = (GeomCreateFunc *)GridCreate;
X		iGridMethods->intersect = GridIntersect;
X		iGridMethods->name = GridName;
X		iGridMethods->convert = GridConvert;
X		iGridMethods->bounds = GridBounds;
X		iGridMethods->checkbounds = FALSE;
X		iGridMethods->closed = TRUE;
X	}
X	return iGridMethods;
X}
X
X/*
X * Place an object in a grid.
X */
Xstatic void
Xengrid(obj, grid)
XGeom *obj;
XGrid *grid;
X{
X	int x, y, z, low[3], high[3];
X	GeomList *ltmp;
X
X	/*
X	 * This routine should *never* be passed an unbounded object, but...
X	 */
X	if (!pos2grid(grid, obj->bounds[LOW], low) ||
X	    !pos2grid(grid, obj->bounds[HIGH], high) ||
X	    obj->bounds[LOW][X] > obj->bounds[HIGH][X]) {
X		/*
X		 * Geom is partially on wholly outside of
X		 * grid -- this should never happen, but just
X		 * in case...
X		 */
X		RLerror(RL_ABORT, "Engrid got an unbounded object?!\n");
X		return;
X	    }
X
X	/*
X	 * For each voxel that intersects the object's bounding
X	 * box, add pointer to this object to voxel's linked list.
X 	 */
X	for (x = low[X]; x <= high[X]; x++) {
X		for (y = low[Y]; y <= high[Y]; y++) {
X			for (z = low[Z]; z <= high[Z]; z++) {
X				ltmp = (GeomList *)share_malloc(sizeof(GeomList));
X				ltmp->obj = obj;
X				ltmp->next = grid->cells[x][y][z];
X				grid->cells[x][y][z] = ltmp;
X			}
X		}
X	}
X}
X
X/*
X * Convert 3D point to index into grid's voxels.
X */
Xstatic int
Xpos2grid(grid, pos, index)
XGrid *grid;
XFloat pos[3];
Xint index[3];
X{
X	index[X] = (int)(x2voxel(grid, pos[0]));
X	index[Y] = (int)(y2voxel(grid, pos[1]));
X	index[Z] = (int)(z2voxel(grid, pos[2]));
X
X	if (index[X] == grid->xsize)
X		index[X]--;
X	if (index[Y] == grid->ysize)
X		index[Y]--;
X	if (index[Z] == grid->zsize)
X		index[Z]--;
X
X	if (index[X] < 0 || index[X] >= grid->xsize ||
X	    index[Y] < 0 || index[Y] >= grid->ysize ||
X	    index[Z] < 0 || index[Z] >= grid->zsize)
X		return FALSE;
X	return TRUE;
X}
X
Xvoid
XGridMethodRegister(meth)
XUserMethodType meth;
X{
X	if (iGridMethods)
X		iGridMethods->user = meth;
X}
END_OF_FILE
if test 11807 -ne `wc -c <'libray/libobj/grid.c'`; then
    echo shar: \"'libray/libobj/grid.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/grid.c'
fi
if test -f 'rayshade/raytrace.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'rayshade/raytrace.c'\"
else
echo shar: Extracting \"'rayshade/raytrace.c'\" \(11927 characters\)
sed "s/^X//" >'rayshade/raytrace.c' <<'END_OF_FILE'
X/*
X * raytrace.c
X *
X * Copyright (C) 1989, 1991, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: raytrace.c,v 4.0.1.1 92/01/10 17:13:02 cek Exp Locker: cek $
X *
X * $Log:	raytrace.c,v $
X * Revision 4.0.1.1  92/01/10  17:13:02  cek
X * patch3: Made status report print actual scanline number.
X * 
X * Revision 4.0  91/07/17  14:50:49  kolb
X * Initial version.
X * 
X */
X
X#include "rayshade.h"
X#include "libsurf/atmosphere.h"
X#include "libsurf/surface.h"
X#include "libcommon/sampling.h"
X#include "options.h"
X#include "stats.h"
X#include "raytrace.h"
X#include "viewing.h"
X
X#define UNSAMPLED	-1
X#define SUPERSAMPLED	-2
X
Xtypedef struct {
X	Pixel	*pix;	/* Pixel values */
X	int	*samp;	/* Sample number */
X} Scanline;
X
Xstatic int		*SampleNumbers;
Xstatic void	RaytraceInit();
X
Xstatic Ray	TopRay;				/* Top-level ray. */
XFloat		SampleTime();
X
XPixel		WhitePix = {1., 1., 1., 1.},
X		BlackPix = {0., 0., 0., 0.};
X
X/*
X * "Dither matrices" used to encode the 'number' of a ray that samples a
X * particular portion of a pixel.  Hand-coding is ugly, but...
X */
Xstatic int OneSample[1] = 	{0};
Xstatic int TwoSamples[4] =	{0, 2,
X				 3, 1};
Xstatic int ThreeSamples[9] =	{0, 2, 7,
X				 6, 5, 1,
X				 3, 8, 4};
Xstatic int FourSamples[16] =	{ 0,  8,  2, 10,
X				 12,  4, 14,  6,
X				  3, 11,  1,  9,
X				 15,  7, 13,  5};
Xstatic int FiveSamples[25] =	{ 0,  8, 23, 17,  2,
X				 19, 12,  4, 20, 15,
X				  3, 21, 16,  9,  6,
X				 14, 10, 24,  1, 13,
X				 22,  7, 18, 11,  5};
Xstatic int SixSamples[36] =	{ 6, 32,  3, 34, 35,  1,
X				  7, 11, 27, 28,  8, 30,
X				 24, 14, 16, 15, 23, 19,
X				 13, 20, 22, 21, 17, 18,
X				 25, 29, 10,  9, 26, 12,
X				 36,  5, 33,  4,  2, 31};
Xstatic int SevenSamples[49] =	{22, 47, 16, 41, 10, 35,  4,
X				  5, 23, 48, 17, 42, 11, 29,
X				 30,  6, 24, 49, 18, 36, 12,
X				 13, 31,  7, 25, 43, 19, 37,
X				 38, 14, 32,  1, 26, 44, 20,
X				 21, 39,  8, 33,  2, 27, 45,
X				 46, 15, 40,  9, 34,  3, 28};
Xstatic int EightSamples[64] =	{ 8, 58, 59,  5,  4, 62, 63,  1,
X				 49, 15, 14, 52, 53, 11, 10, 56,
X				 41, 23, 22, 44, 45, 19, 18, 48,
X				 32, 34, 35, 29, 28, 38, 39, 25,
X				 40, 26, 27, 37, 36, 30, 31, 33,
X				 17, 47, 46, 20, 21, 43, 42, 24,
X				  9, 55, 54, 12, 13, 51, 50, 16,
X				 64,  2,  3, 61, 60,  6,  7, 57};
X
Xvoid	AdaptiveRefineScanline(), FullySamplePixel(), FullySampleScanline(),
X	SingleSampleScanline();
Xstatic int	ExcessiveContrast();
Xstatic Scanline scan0, scan1, scan2;
X
X
Xvoid
Xraytrace(argc, argv)
Xint argc;
Xchar **argv;
X{
X	int y, *tmpsamp;
X	Pixel *tmppix;
X	Float usertime, systime, lasttime;
X
X	/*
X	 * If this is the first frame,
X	 * allocate scanlines, etc.
X	 */
X	if (Options.framenum == Options.startframe)
X		RaytraceInit();
X	/*
X	 * The top-level ray TopRay always has as its origin the
X	 * eye position and as its medium NULL, indicating that it
X	 * is passing through a medium with index of refraction
X	 * equal to DefIndex.
X	 */
X	TopRay.pos = Camera.pos;
X	TopRay.media = (Medium *)0;
X	TopRay.depth = 0;
X
X	/*
X	 * Always fully sample the bottom and top rows and the left
X	 * and right column of pixels.  This minimizes artifacts that
X	 * may arise when piecing together images.
X	 */
X	FullySampleScanline(0, &scan0);
X
X	SingleSampleScanline(1, &scan1);
X	FullySamplePixel(0, 1, &scan1.pix[0], &scan1.samp[0]);
X	FullySamplePixel(Screen.xsize -1, 1, &scan1.pix[Screen.xsize -1],
X		&scan1.samp[Screen.xsize -1]);
X
X	lasttime = 0;
X	for (y = 1; y < Screen.ysize; y++) {
X		SingleSampleScanline(y+1, &scan2);
X		FullySamplePixel(0, y+1, &scan2.pix[0], &scan2.samp[0]);
X		FullySamplePixel(Screen.xsize -1, y+1,
X			&scan2.pix[Screen.xsize -1],
X			&scan2.samp[Screen.xsize -1]);
X
X		if (Sampling.sidesamples > 1)
X			AdaptiveRefineScanline(y,&scan0,&scan1,&scan2);
X
X		PictureWriteLine(scan0.pix);
X
X		tmppix = scan0.pix;
X		tmpsamp = scan0.samp;
X		scan0.pix = scan1.pix;
X		scan0.samp = scan1.samp;
X		scan1.pix = scan2.pix;
X		scan1.samp = scan2.samp;
X		scan2.pix = tmppix;
X		scan2.samp = tmpsamp;
X
X		if ((y+Screen.miny-1) % Options.report_freq == 0) {
X			fprintf(Stats.fstats,"Finished line %d (%lu rays",
X						y+Screen.miny-1,
X						Stats.EyeRays);
X			if (Options.verbose) {
X				/*
X				 * Report total CPU and split times.
X				 */
X				RSGetCpuTime(&usertime, &systime);
X				fprintf(Stats.fstats,", %2.2f sec,",
X						usertime+systime);
X				fprintf(Stats.fstats," %2.2f split",
X						usertime+systime-lasttime);
X				lasttime = usertime+systime;
X			}
X			fprintf(Stats.fstats,")\n");
X			(void)fflush(Stats.fstats);
X		}
X
X	}
X	/*
X	 * Supersample last scanline.
X	 */
X	for (y = 1; y < Screen.xsize -1; y++) {
X		if (scan0.samp[y] != SUPERSAMPLED)
X			FullySamplePixel(y, Screen.ysize -1,
X				&scan0.pix[y],
X				&scan0.samp[y]);
X	}
X	PictureWriteLine(scan0.pix);
X}
X
Xvoid
XSingleSampleScanline(line, data)
Xint line;
XScanline *data;
X{
X	Float upos, vpos, yp;
X	int x, usamp, vsamp;
X	Pixel tmp;
X
X	yp = line + Screen.miny - 0.5*Sampling.filterwidth;
X	for (x = 0; x < Screen.xsize; x++) {
X		/*
X		 * Pick a sample number...
X		 */
X		data->samp[x] = nrand() * Sampling.totsamples;
X		/*
X		 * Take sample corresponding to sample #.
X		 */
X		usamp = data->samp[x] % Sampling.sidesamples;
X		vsamp = data->samp[x] / Sampling.sidesamples;
X
X		vpos = yp + vsamp * Sampling.filterdelta;
X		upos = x + Screen.minx - 0.5*Sampling.filterwidth +
X				usamp*Sampling.filterdelta;
X		if (Options.jitter) {
X			vpos += nrand()*Sampling.filterdelta;
X			upos += nrand()*Sampling.filterdelta;
X		}
X		TopRay.time = SampleTime(SampleNumbers[data->samp[x]]);
X		SampleScreen(upos, vpos, &TopRay,
X			&data->pix[x], SampleNumbers[data->samp[x]]);
X		if (Options.samplemap)
X			data->pix[x].alpha = 0;
X	}
X}
X
Xvoid
XFullySampleScanline(line, data)
Xint line;
XScanline *data;
X{
X	int x;
X
X	for (x = 0; x < Screen.xsize; x++) {
X		data->samp[x] = UNSAMPLED;
X		FullySamplePixel(x, line, &data->pix[x], &data->samp[x]);
X	}
X}
X
Xvoid
XFullySamplePixel(xp, yp, pix, prevsamp)
Xint xp, yp;
XPixel *pix;
Xint *prevsamp;
X{
X	Float upos, vpos, u, v;
X	int x, y, sampnum;
X	Pixel ctmp;
X
X	if (*prevsamp == SUPERSAMPLED)
X		return;	/* already done */
X
X	Stats.SuperSampled++;
X	if (*prevsamp == UNSAMPLED) {
X		/*
X		 * No previous sample; initialize to black.
X		 */
X		pix->r = pix->g = pix->b = pix->alpha = 0.;
X	} else {
X		if (Sampling.sidesamples == 1) {
X			*prevsamp = SUPERSAMPLED;
X			return;
X		}
X		x = *prevsamp % Sampling.sidesamples;
X		y = *prevsamp / Sampling.sidesamples;
X		pix->r *= Sampling.filter[x][y];
X		pix->g *= Sampling.filter[x][y];
X		pix->b *= Sampling.filter[x][y];
X		pix->alpha *= Sampling.filter[x][y];
X	}
X
X	sampnum = 0;
X	xp += Screen.minx;
X	vpos = Screen.miny + yp - 0.5*Sampling.filterwidth;
X	for (y = 0; y < Sampling.sidesamples; y++,
X	     vpos += Sampling.filterdelta) {
X		upos = xp - 0.5*Sampling.filterwidth;
X		for (x = 0; x < Sampling.sidesamples; x++,
X		     upos += Sampling.filterdelta) {
X			if (sampnum != *prevsamp) {
X				if (Options.jitter) {
X					u = upos + nrand()*Sampling.filterdelta;
X					v = vpos + nrand()*Sampling.filterdelta;
X				} else {
X					u = upos;
X					v = vpos;
X				}
X				TopRay.time = SampleTime(SampleNumbers[sampnum]);
X				SampleScreen(u, v, &TopRay, &ctmp,
X					SampleNumbers[sampnum]);
X				pix->r += ctmp.r*Sampling.filter[x][y];
X				pix->g += ctmp.g*Sampling.filter[x][y];
X				pix->b += ctmp.b*Sampling.filter[x][y];
X				pix->alpha += ctmp.alpha*Sampling.filter[x][y];
X			}
X			if (++sampnum == Sampling.totsamples)
X				sampnum = 0;
X		}
X	}
X
X	if (Options.samplemap)
X		pix->alpha = 255;
X
X	*prevsamp = SUPERSAMPLED;
X}
X
Xvoid
XAdaptiveRefineScanline(y, scan0, scan1, scan2)
Xint y;
XScanline *scan0, *scan1, *scan2;
X{
X	int x, done;
X
X	/*
X	 * Walk down scan1, looking at 4-neighbors for excessive contrast.
X	 * If found, supersample *all* neighbors not already supersampled.
X	 * The process is repeated until either there are no
X	 * high-contrast regions or all such regions are already supersampled.
X	 */
X
X	do {
X		done = TRUE;
X		for (x = 1; x < Screen.xsize -1; x++) {
X			/*
X		 	 * Find min and max RGB for area we care about
X			 */
X			if (ExcessiveContrast(x, scan0->pix, scan1->pix,
X			    scan2->pix)) {
X				if (scan1->samp[x-1] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x-1, y,
X						&scan1->pix[x-1],
X						&scan1->samp[x-1]);
X				}
X				if (scan0->samp[x] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x, y-1,
X						&scan0->pix[x],
X						&scan0->samp[x]);
X				}
X				if (scan1->samp[x+1] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x+1, y,
X						&scan1->pix[x+1],
X						&scan1->samp[x+1]);
X				}
X				if (scan2->samp[x] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x, y+1,
X						&scan2->pix[x],
X						&scan2->samp[x]);
X				}
X				if (scan1->samp[x] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x, y,
X						&scan1->pix[x],
X						&scan1->samp[x]);
X				}
X			}
X		}
X	} while (!done);
X}
X
Xstatic int
XExcessiveContrast(x, pix0, pix1, pix2)
Xint x;
XPixel *pix0, *pix1, *pix2;
X{
X	Float mini, maxi, sum, diff;
X
X	maxi = max(pix0[x].r, pix1[x-1].r);
X	if (pix1[x].r > maxi) maxi = pix1[x].r;
X	if (pix1[x+1].r > maxi) maxi = pix1[x+1].r;
X	if (pix2[x].r > maxi) maxi = pix2[x].r;
X
X	mini = min(pix0[x].r, pix1[x-1].r);
X	if (pix1[x].r < mini) mini = pix1[x].r;
X	if (pix1[x+1].r < mini) mini = pix1[x+1].r;
X	if (pix2[x].r < mini) mini = pix2[x].r;
X
X	diff = maxi - mini;
X	sum = maxi + mini;
X	if (sum > EPSILON && diff/sum > Options.contrast.r)
X		return TRUE;
X
X	maxi = max(pix0[x].g, pix1[x-1].g);
X	if (pix1[x].g > maxi) maxi = pix1[x].g;
X	if (pix1[x+1].g > maxi) maxi = pix1[x+1].g;
X	if (pix2[x].g > maxi) maxi = pix2[x].g;
X
X	mini = min(pix0[x].g, pix1[x-1].g);
X	if (pix1[x].g < mini) mini = pix1[x].g;
X	if (pix1[x+1].g < mini) mini = pix1[x+1].g;
X	if (pix2[x].g < mini) mini = pix2[x].g;
X
X	diff = maxi - mini;
X	sum = maxi + mini;
X
X	if (sum > EPSILON && diff/sum > Options.contrast.g)
X		return TRUE;
X
X	maxi = max(pix0[x].b, pix1[x-1].b);
X	if (pix1[x].b > maxi) maxi = pix1[x].b;
X	if (pix1[x+1].b > maxi) maxi = pix1[x+1].b;
X	if (pix2[x].b > maxi) maxi = pix2[x].b;
X
X	mini = min(pix0[x].b, pix1[x-1].b);
X	if (pix1[x].b < mini) mini = pix1[x].b;
X	if (pix1[x+1].b < mini) mini = pix1[x+1].b;
X	if (pix2[x].b < mini) mini = pix2[x].b;
X
X	diff = maxi - mini;
X	sum = maxi + mini;
X	if (sum > EPSILON && diff/sum > Options.contrast.b)
X		return TRUE;
X
X	return FALSE;
X}
X
XFloat
XSampleTime(sampnum)
Xint sampnum;
X{
X	Float window, jitter = 0.0, res;
X
X	if (Options.shutterspeed <= 0.)
X		return Options.framestart;
X	if (Options.jitter)
X		jitter = nrand();
X	window = Options.shutterspeed / Sampling.totsamples;
X	res = Options.framestart + window * (sampnum + jitter);
X	TimeSet(res);
X	return res;
X}
X
Xstatic void
XRaytraceInit()
X{
X
X	switch (Sampling.sidesamples) {
X		case 1:
X			SampleNumbers = OneSample;
X			break;
X		case 2:
X			SampleNumbers = TwoSamples;
X			break;
X		case 3:
X			SampleNumbers = ThreeSamples;
X			break;
X		case 4:
X			SampleNumbers = FourSamples;
X			break;
X		case 5:
X			SampleNumbers = FiveSamples;
X			break;
X		case 6:
X			SampleNumbers = SixSamples;
X			break;
X		case 7:
X			SampleNumbers = SevenSamples;
X			break;
X		case 8:
X			SampleNumbers = EightSamples;
X			break;
X		default:
X			RLerror(RL_PANIC,
X				"Sorry, %d rays/pixel not supported.\n",
X					Sampling.totsamples);
X	}
X
X	/*
X 	 * Allocate pixel arrays and arrays to store sampling info.
X 	 */
X	scan0.pix = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
X	scan1.pix = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
X	scan2.pix = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
X
X	scan0.samp = (int *)Malloc(Screen.xsize * sizeof(int));
X	scan1.samp = (int *)Malloc(Screen.xsize * sizeof(int));
X	scan2.samp = (int *)Malloc(Screen.xsize * sizeof(int));
X}
END_OF_FILE
if test 11927 -ne `wc -c <'rayshade/raytrace.c'`; then
    echo shar: \"'rayshade/raytrace.c'\" unpacked with wrong size!
fi
# end of 'rayshade/raytrace.c'
fi
echo shar: End of archive 14 \(of 19\).
cp /dev/null ark14isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 19 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0