v21i013: A ray tracing program, Part06/08

[Rich Salz]

Submitted-by: Craig Kolb <craig at weedeater.math.yale.edu>
Posting-number: Volume 21, Issue 13
Archive-name: rayshade/part06

#! /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 <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of archive 6 (of 8)."
# Contents:  src/hf.c src/input_yacc.y src/raytrace.c
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'src/hf.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'src/hf.c'\"
else
echo shar: Extracting \"'src/hf.c'\" \(17271 characters\)
sed "s/^X//" >'src/hf.c' <<'END_OF_FILE'
X/*
X * hf.c
X *
X * Copyright (C) 1989, Craig E. Kolb
X *
X * This software may be freely copied, modified, and redistributed,
X * provided that this copyright notice is preserved on all copies.
X *
X * There is no warranty or other guarantee of fitness for this software,
X * it is provided solely .  Bug reports or fixes may be sent
X * to the author, who may or may not act on them as he desires.
X *
X * You may not include this software in a program or other software product
X * without supplying the source, or without informing the end-user that the
X * source is available for no extra charge.
X *
X * If you modify this software, you should include a notice giving the
X * name of the person performing the modification, the date of modification,
X * and the reason for such modification.
X *
X * $Id: hf.c,v 3.0 89/10/27 02:05:51 craig Exp $
X *
X * $Log:	hf.c,v $
X * Revision 3.0  89/10/27  02:05:51  craig
X * Baseline for first official release.
X * 
X */
X#include <stdio.h>
X#include <math.h>
X#include "typedefs.h"
X#include "funcdefs.h"
X#include "constants.h"
X
X/*
X * Any height values <= HF_UNSET is not considered to be part of the
X * height field. Any trianges containing such a vertex will not be
X * rendered.  This allows one to render non-square height fields.
X */
X#define HF_UNSET		(-1000.)
X/*
X * Number of datapoints in a single cell.  If you've got loads of memory,
X * decrease this number.  The 'optimal' number is quite system-dependent,
X * but something around 20 seems to work well.  For systems which don't
X * have much memory, this constant should be greater than or equal to
X * the largest height field which will be rendered, thus making the
X * algorithm non-hierarchical.
X */
X#define BESTSIZE		16
X/*
X * Size of triangle cache.
X */
X#define CACHESIZE		6
X/*
X * Used to differentiate between the two triangles used to represent a cell:
X *	a------d
X *      |\     |
X *      | \TRI1|	TRI1 == c-->d-->a-->c
X *      |  \   |
X *      |   \  |
X *	|    \ |
X *      |TRI2 \|	TRI2 == c-->a-->b-->d
X *      b------c
X */
X#define TRI1			1
X#define TRI2			2
X
Xdouble DDA2D(), CheckCell(), intHftri();
Xfloat minalt(), maxalt();
X
Xtypedef struct {
X	int stepX, stepY;
X	double tDX, tDY;
X	float minz, maxz;
X	int outX, outY;
X	Vector cp, pDX, pDY;
X} Trav2D;
X
XhfTri *CreateHfTriangle(), *GetQueuedTri();
X
XObject *
Xmakhf(surf, filename)
Xchar *surf, *filename;
X{
X	Hf *hf;
X	Object *newobj;
X	Primitive *prim;
X	FILE *fp;
X	float val, *maxptr, *minptr;
X	int i, j;
X
X	fp = fopen(filename, "r");
X	if (fp == (FILE *)NULL)
X		yyerror("Cannot open heightfield file");
X
X	prim = mallocprim();
X	newobj = new_object(NULL, HF, (char *)prim, (Trans *)NULL);
X	prim->type = HF;
X	prim->surf = find_surface(surf);
X	hf = (Hf *)Malloc(sizeof(Hf));
X	prim->objpnt.p_hf = hf;
X	/*
X	 * Make the following an option someday.
X	 */
X	hf->BestSize = BESTSIZE;
X	/*
X	 * Store the inverse for faster computation.
X	 */
X	hf->iBestSize = 1. / (float)hf->BestSize;
X	/*
X	 * Get HF size.
X	 */
X	if (fread(&hf->size, sizeof(int), 1, fp) == 0)
X		yyerror("Cannot read height field size.");
X
X	hf->data = (float **)share_malloc(hf->size * sizeof(float *));
X	for (i = 0; i < hf->size; i++) {
X		hf->data[i] = (float *)share_malloc(hf->size * sizeof(float));
X		/*
X		 * Read in row of HF data.
X		 */
X		if (fread(hf->data[i],sizeof(float),hf->size,fp) != hf->size)
X			yyerror("Not enough heightfield data.");
X		for (j = 0; j < hf->size; j++) {
X			val = hf->data[i][j];
X			if (val <= HF_UNSET) {
X				hf->data[i][j] = HF_UNSET;
X				/*
X				 * Don't include the point in min/max
X				 * calculations.
X				 */
X				continue;
X			}
X			if (val > hf->maxz)
X				hf->maxz = val;
X			if (val < hf->minz)
X				hf->minz = val;
X		}
X	}
X	fclose(fp);
X	/*
X	 * Allocate levels of grid.  hf->levels = log base BestSize of hf->size
X	 */
X	for (i = hf->size, hf->levels = 0; i > hf->BestSize; i /= hf->BestSize,
X				hf->levels++)
X			;
X	hf->levels++;
X	hf->qsize = CACHESIZE;
X	hf->q = (hfTri **)Calloc(hf->qsize, sizeof(hfTri *));
X	hf->qtail = 0;
X
X	hf->lsize = (int *)share_malloc(hf->levels * sizeof(int));
X	hf->spacing = (float *)share_malloc(hf->levels * sizeof(float));
X	hf->boundsmax = (float ***)share_malloc(hf->levels * sizeof(float **));
X	hf->boundsmin = (float ***)share_malloc(hf->levels * sizeof(float **));
X
X	hf->spacing[0] = hf->size -1;
X	hf->lsize[0] = (int)hf->spacing[0];
X	hf->boundsmax[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
X	hf->boundsmin[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
X	/*
X	 * Compute initial bounding boxes
X	 */
X	for (i = 0; i < hf->lsize[0]; i++) {
X		hf->boundsmax[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
X		hf->boundsmin[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
X		maxptr = hf->boundsmax[0][i];
X		minptr = hf->boundsmin[0][i];
X		for (j = 0; j < hf->lsize[0]; j++) {
X			*maxptr++ = maxalt(i, j, hf->data) + EPSILON;
X			*minptr++ = minalt(i, j, hf->data) - EPSILON;
X		}
X	}
X
X	for (i = 1; i < hf->levels; i++) {
X		hf->spacing[i] = hf->spacing[i-1] * hf->iBestSize;
X		hf->lsize[i] = (int)hf->spacing[i];
X		if (hf->lsize[i-1] % hf->BestSize != 0)
X			hf->lsize[i]++;
X		hf->boundsmax[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
X		hf->boundsmin[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
X		for (j = 0; j < hf->lsize[i]; j++) {
X			hf->boundsmax[i][j] = (float *)share_malloc(hf->lsize[i] *
X							sizeof(float));
X			hf->boundsmin[i][j] = (float *)share_malloc(hf->lsize[i] *
X							sizeof(float));
X		}
X		integrate_grid(hf, i);
X	}
X
X	hf->boundbox[LOW][X] = hf->boundbox[LOW][Y] = 0;
X	hf->boundbox[HIGH][X] = hf->boundbox[HIGH][Y] = 1;
X	hf->boundbox[LOW][Z] = hf->minz;
X	hf->boundbox[HIGH][Z] = hf->maxz;
X
X	return newobj;
X}
X
X/*
X * Intersect ray with height field.
X */
Xdouble
Xinthf(pos, ray, obj)
XVector *pos, *ray;
XPrimitive *obj;
X{
X	Hf *hf;
X	Ray tmpray;
X	Vector hitpos;
X	double offset, dist;
X	Trav2D trav;
X	extern unsigned long primtests[];
X
X	primtests[HF]++;
X	hf = obj->objpnt.p_hf;
X
X	/*
X	 * Find where we hit the hf cube.
X	 */
X	if (OutOfBounds(pos, hf->boundbox)) {
X		tmpray.pos = *pos;
X		tmpray.dir = *ray;
X		offset = IntBounds(&tmpray, hf->boundbox);
X		if (offset <= 0.)
X			/* Never hit hf. */
X			return 0.;
X		hitpos.x = pos->x + ray->x * offset;
X		hitpos.y = pos->y + ray->y * offset;
X		hitpos.z = pos->z + ray->z * offset;
X	} else
X		hitpos = *pos;
X	/*
X	 * Find out which cell "hitpoint" is in.  This could be
X	 * causing problems!
X	 */
X	if (equal(hitpos.x, 1.))
X		hitpos.x -= EPSILON;
X	if (equal(hitpos.y, 1.))
X		hitpos.y -= EPSILON;
X
X	if (ray->x < 0.) {
X		trav.stepX = trav.outX = -1;
X		trav.tDX = -1. / (ray->x * hf->spacing[hf->levels -1]);
X	} else if (ray->x > 0.) {
X		trav.stepX = 1;
X		trav.outX = hf->lsize[hf->levels -1];
X		/*
X		 * (1./size) / ray
X		 */
X		trav.tDX = 1. / (ray->x * hf->spacing[hf->levels -1]);
X	}
X
X	if (ray->y < 0.) {
X		trav.stepY = trav.outY = -1;
X		trav.tDY = -1. / (ray->y * hf->spacing[hf->levels -1]);
X	} else if (ray->y > 0.) {
X		trav.stepY = 1;
X		trav.outY = hf->lsize[hf->levels -1];
X		trav.tDY = 1. / (ray->y * hf->spacing[hf->levels -1]);
X	}
X
X	trav.pDX.x = ray->x * trav.tDX;
X	trav.pDX.y = ray->y * trav.tDX;
X	trav.pDX.z = ray->z * trav.tDX;
X	trav.pDY.x = ray->x * trav.tDY;
X	trav.pDY.y = ray->y * trav.tDY;
X	trav.pDY.z = ray->z * trav.tDY;
X
X	trav.cp = hitpos;
X	trav.minz = hf->minz;
X	trav.maxz = hf->maxz;
X	dist = DDA2D(hf, pos, ray, hf->levels -1, &trav);
X	return dist;
X}
X
X/*
X * Traverse the grid using a modified DDA algorithm.  If the extent of
X * the ray over a cell intersects the bounding volume defined by the
X * four corners of the cell, either recurse or perform ray/surface
X * intersection test.
X */
Xdouble
XDDA2D(hf, pos, ray, level, trav)
XHf *hf;
XVector *pos, *ray;
Xint level;
XTrav2D *trav;
X{
X	int x, y, size, posZ;
X	float **boundsmin, **boundsmax, spacing;
X	double dist, tX, tY;
X	Trav2D newtrav;
X	Vector nxp, nyp;
X
X	size = hf->lsize[level];
X	spacing = hf->spacing[level];
X
X	posZ = (ray->z > 0.);
X
X	x = trav->cp.x * hf->spacing[level];
X	y = trav->cp.y * hf->spacing[level];
X	boundsmax = hf->boundsmax[level];
X	boundsmin = hf->boundsmin[level];
X
X	if (trav->outX > size) trav->outX = size;
X	if (trav->outY > size) trav->outY = size;
X	if (trav->outX < 0) trav->outX = -1;
X	if (trav->outY < 0) trav->outY = -1;
X
X	if (ray->x < 0.) {
X		tX = (x /spacing - trav->cp.x) / ray->x;
X	} else if (ray->x > 0.)
X		tX = ((x+1)/spacing - trav->cp.x) / ray->x;
X	else
X		tX = FAR_AWAY;
X	if (ray->y < 0.) {
X		tY = (y /spacing - trav->cp.y) / ray->y;
X	} else if (ray->y > 0.)
X		tY = ((y+1)/spacing - trav->cp.y) / ray->y;
X	else
X		tY = FAR_AWAY;
X
X	nxp.x = trav->cp.x + tX * ray->x;
X	nxp.y = trav->cp.y + tX * ray->y;
X	nxp.z = trav->cp.z + tX * ray->z;
X
X	nyp.x = trav->cp.x + tY * ray->x;
X	nyp.y = trav->cp.y + tY * ray->y;
X	nyp.z = trav->cp.z + tY * ray->z;
X
X	do {
X		if (tX < tY) {
X			if ((posZ && trav->cp.z <= boundsmax[y][x] &&
X			     nxp.z >= boundsmin[y][x]) ||
X			    (!posZ && trav->cp.z >= boundsmin[y][x] &&
X			     nxp.z <= boundsmax[y][x])) {
X				if (level) {
X					/*
X					 * Recurve -- compute constants
X					 * needed for next level.
X					 * Nicely enough, this just
X					 * involves a few multiplications.
X					 */
X					newtrav = *trav;
X					newtrav.tDX *= hf->iBestSize;
X					newtrav.tDY *= hf->iBestSize;
X					newtrav.maxz = boundsmax[y][x];
X					newtrav.minz = boundsmin[y][x];
X					if (ray->x < 0.)
X						newtrav.outX=hf->BestSize*x-1;
X					else
X						newtrav.outX=hf->BestSize*(x+1);
X					if (ray->y < 0.)
X						newtrav.outY=hf->BestSize*y-1;
X					else
X						newtrav.outY=hf->BestSize*(y+1);
X					newtrav.pDX.x *= hf->iBestSize;
X					newtrav.pDX.y *= hf->iBestSize;
X					newtrav.pDX.z *= hf->iBestSize;
X					newtrav.pDY.x *= hf->iBestSize;
X					newtrav.pDY.y *= hf->iBestSize;
X					newtrav.pDY.z *= hf->iBestSize;
X					dist = DDA2D(hf,pos,ray,level-1,
X						&newtrav);
X				} else
X					dist = CheckCell(x,y,hf,ray,pos);
X				if (dist > EPSILON)
X					return dist;
X			}
X			x += trav->stepX;		/* Move in X */
X			if (x == trav->outX)		/* If outside, quit */
X				return 0.;
X			tX += trav->tDX;	/* Update position on ray */
X			trav->cp = nxp;		/* cur pos gets next pos */
X			nxp.x += trav->pDX.x;	/* Compute next pos */
X			nxp.y += trav->pDX.y;
X			nxp.z += trav->pDX.z;
X		} else {
X			if ((posZ && trav->cp.z <= boundsmax[y][x] &&
X			     nyp.z >= boundsmin[y][x]) ||
X			    (!posZ && trav->cp.z >= boundsmin[y][x] &&
X			     nyp.z <= boundsmax[y][x])) {
X				if (level) {
X					/* Recurse */
X					newtrav = *trav;
X					newtrav.tDX *= hf->iBestSize;
X					newtrav.tDY *= hf->iBestSize;
X					newtrav.maxz = boundsmax[y][x];
X					newtrav.minz = boundsmin[y][x];
X					if (ray->x < 0.)
X						newtrav.outX=hf->BestSize*x-1;
X					else
X						newtrav.outX=hf->BestSize*(x+1);
X					if (ray->y < 0.)
X						newtrav.outY=hf->BestSize*y-1;
X					else
X						newtrav.outY=hf->BestSize*(y+1);
X					newtrav.pDX.x *= hf->iBestSize;
X					newtrav.pDX.y *= hf->iBestSize;
X					newtrav.pDX.z *= hf->iBestSize;
X					newtrav.pDY.x *= hf->iBestSize;
X					newtrav.pDY.y *= hf->iBestSize;
X					newtrav.pDY.z *= hf->iBestSize;
X					dist = DDA2D(hf,pos,ray,level-1,
X							&newtrav);
X				} else
X					dist = CheckCell(x,y,hf,ray,pos);
X				if (dist > EPSILON)
X					return dist;
X			}
X			y += trav->stepY;
X			if (y == trav->outY)
X				return 0.;
X			tY += trav->tDY;
X			trav->cp = nyp;
X			nyp.x += trav->pDY.x;
X			nyp.y += trav->pDY.y;
X			nyp.z += trav->pDY.z;
X		}
X	} while ((trav->cp.x <= 1. && trav->cp.y <= 1.) &&
X		 ((posZ && trav->cp.z <= trav->maxz) ||
X		 (!posZ && trav->cp.z >= trav->minz)));
X
X		/*
X		 * while ((we're inside the horizontal bounding box)
X		 *		(usually caught by outX & outY, but
X		 *		 it's possible to go "too far" due to
X		 *		 the fact that our levels of grids do
X		 *		 not "nest" exactly if gridsize%BestSize != 0)
X		 *	  and
X		 *	  ((if ray->z is positive and we haven't gone through
X		 *	   the upper bounding plane) or
X		 *	  (if ray->z is negative and we haven't gone through
X		 *	   the lower bounding plane)));
X		 */
X
X	return 0.;
X}
X
X/*
X * Check for ray/cell intersection
X */
Xdouble
XCheckCell(x, y, hf, ray, pos)
Xint x, y;
XHf *hf;
XVector *ray, *pos;
X{
X	hfTri *tri1, *tri2;
X	double d1, d2;
X
X	d1 = d2 = FAR_AWAY;
X
X	if (tri1 = CreateHfTriangle(hf, x, y, x+1, y, x, y+1, TRI1))
X		d1 = intHftri(ray, pos, tri1);
X	if (tri2 = CreateHfTriangle(hf, x+1, y, x+1, y+1, x, y+1, TRI2))
X		d2 = intHftri(ray, pos, tri2);
X
X	if (d1 == FAR_AWAY && d2 == FAR_AWAY)
X		return 0;
X
X	if (d1 < d2) {
X		hf->hittri = *tri1;
X		return d1;
X	}
X
X	hf->hittri = *tri2;
X	return d2;
X}
X
XhfTri *
XCreateHfTriangle(hf, x1, y1, x2, y2, x3, y3, which)
XHf *hf;
Xint x1, y1, x2, y2, x3, y3, which;
X{
X	hfTri *tri;
X	double xid, yid;
X	Vector tmp1, tmp2;
X
X	/*
X	 * Don't use triangles with "unset" vertices.
X	 */
X	if (hf->data[y1][x1] == HF_UNSET ||
X	    hf->data[y2][x2] == HF_UNSET ||
X	    hf->data[y3][x3] == HF_UNSET)
X		return (hfTri *)0;
X
X	xid = (double)x1 / (double)(hf->size -1);
X	yid = (double)y1 / (double)(hf->size -1);
X
X	if ((tri = GetQueuedTri(hf, xid, yid, which)) != (hfTri *)0)
X		return tri;
X
X	tri = (hfTri *)Malloc(sizeof(hfTri));
X
X	tri->type = which;
X	tri->v1.x = xid;
X	tri->v1.y = yid;
X	tri->v1.z = hf->data[y1][x1];
X	tri->v2.x = (double)x2 / (double)(hf->size-1);
X	tri->v2.y = (double)y2 / (double)(hf->size-1);
X	tri->v2.z = hf->data[y2][x2];
X	tri->v3.x = (double)x3 / (double)(hf->size-1);
X	tri->v3.y = (double)y3 / (double)(hf->size-1);
X	tri->v3.z = hf->data[y3][x3];
X
X	tmp1.x = tri->v2.x - tri->v1.x;
X	tmp1.y = tri->v2.y - tri->v1.y;
X	tmp1.z = tri->v2.z - tri->v1.z;
X	tmp2.x = tri->v3.x - tri->v1.x;
X	tmp2.y = tri->v3.y - tri->v1.y;
X	tmp2.z = tri->v3.z - tri->v1.z;
X
X	(void)crossp(&tri->norm, &tmp1, &tmp2);
X
X	tri->d = -dotp(&tri->v1, &tri->norm);
X
X	QueueTri(hf, tri);
X	return tri;
X}
X
X/*
X * Intersect ray with right isoscoles triangle, the hypotenuse of which
X * has slope of -1.
X */
Xdouble
XintHftri(ray, pos, tri)
XhfTri *tri;
XVector *pos, *ray;
X{
X	double u, v, dist, xpos, ypos;
X
X	u = dotp(&tri->norm, pos) + tri->d;
X	v = dotp(&tri->norm, ray);
X
X	if ((u <= 0. || v > -EPSILON) && (u >= 0. && v < EPSILON))
X		return FAR_AWAY;
X
X	dist = -u / v;
X
X	if (dist < EPSILON)
X		return FAR_AWAY;
X
X	xpos = pos->x + dist*ray->x;
X	ypos = pos->y + dist*ray->y;
X
X	if (tri->type == TRI1 && xpos >= tri->v1.x && ypos >= tri->v1.y &&
X		xpos + ypos <= tri->v2.x + tri->v2.y)
X		return dist;
X	if (tri->type == TRI2 && xpos <= tri->v2.x && ypos <= tri->v2.y &&
X		xpos + ypos >= tri->v1.x + tri->v1.y)
X		return dist;
X	return FAR_AWAY;
X}
X
X/*
X * Compute normal to height field.
X */
Xnrmhf(pos, prim, nrm)
XVector *pos, *nrm;
XPrimitive *prim;
X{
X	*nrm = prim->objpnt.p_hf->hittri.norm;
X}
X
X/*
X * Compute heightfield bounding box.
X */
Xhfextent(obj, bounds)
XPrimitive *obj;
Xdouble bounds[2][3];
X{
X	Hf *hf = obj->objpnt.p_hf;
X
X	/*
X	 * By default, height fields are centered at (0.5, 0.5, 0.)
X	 */
X	bounds[LOW][X] = bounds[LOW][Y] = 0;
X	bounds[HIGH][X] = bounds[HIGH][Y] = 1;
X	bounds[LOW][Z] = hf->minz;
X	bounds[HIGH][Z] = hf->maxz;
X}
X
X/*
X * Build min/max altitude value arrays for the given grid level.
X */
Xintegrate_grid(hf, level)
XHf *hf;
Xint level;
X{
X	int i, j, k, l, ii, ji;
X	float max_alt, min_alt;
X	float **maxinto, **mininto, **frommax, **frommin, *minptr, *maxptr;
X	int insize, fromsize, fact;
X
X	maxinto = hf->boundsmax[level];
X	mininto = hf->boundsmin[level];
X	insize = hf->lsize[level];
X	frommax = hf->boundsmax[level-1];
X	frommin = hf->boundsmin[level-1];
X	fact = hf->BestSize;
X	fromsize = hf->lsize[level-1];
X
X	ii = 0;
X
X	for (i = 0; i < insize; i++) {
X		ji = 0;
X		for (j = 0; j < insize; j++) {
X			max_alt = HF_UNSET;
X			min_alt = -HF_UNSET;
X			for (k = 0; k <= fact; k++) {
X				if (ii+k >= fromsize)
X					continue;
X				maxptr = &frommax[ii+k][ji];
X				minptr = &frommin[ii+k][ji];
X				for (l = 0; l <= fact; l++,maxptr++,minptr++) {
X					if (ji+l >= fromsize)
X						continue;
X					if (*maxptr > max_alt)
X						max_alt = *maxptr;
X					if (*minptr < min_alt)
X						min_alt = *minptr;
X				}
X			}
X			maxinto[i][j] = max_alt + EPSILON;
X			mininto[i][j] = min_alt - EPSILON;
X			ji += fact;
X		}
X		ii += fact;
X	}
X}
X
X/*
X * Place the given triangle in the triangle cache.
X */
XQueueTri(hf, tri)
XHf *hf;
XhfTri *tri;
X{
X	if (hf->q[hf->qtail])		/* Free old triangle data */
X		free((char *)hf->q[hf->qtail]);
X	hf->q[hf->qtail] = tri;		/* Put on tail */
X	hf->qtail = (hf->qtail + 1) % hf->qsize;	/* Increment tail */
X}
X
X/*
X * Search list of cached trianges to see if this triangle has been
X * cached.  If so, return a pointer to it.  If not, return null pointer.
X */
XhfTri *
XGetQueuedTri(hf, x, y, flag)
XHf *hf;
Xdouble x, y;
Xint flag;
X{
X	register int i;
X	register hfTri **tmp;
X
X	for (i = 0, tmp = hf->q; i < hf->qsize; i++, tmp++) {
X		if (*tmp && (*tmp)->v1.x == x && (*tmp)->v1.y == y &&
X		    (*tmp)->type == flag)
X			return *tmp;	/* vertices & flag match, return it */
X	}
X
X	return (hfTri *)0;
X}
X
X/*
X * Return maximum height of cell indexed by y,x.  This could be done
X * as a macro, but many C compliers will choke on it.
X */
Xfloat
Xminalt(y,x,data)
Xint x, y;
Xfloat **data;
X{
X	float  min_alt;
X
X	min_alt = min(data[y][x], data[y+1][x]);
X	min_alt = min(min_alt, data[y][x+1]);
X	min_alt = min(min_alt, data[y+1][x+1]);
X	return min_alt;
X}
X
X/*
X * Return maximum cell height, as above.
X */
Xfloat
Xmaxalt(y,x,data)
Xint x, y;
Xfloat **data;
X{
X	float  max_alt;
X
X	max_alt = max(data[y][x], data[y+1][x]);
X	max_alt = max(max_alt, data[y][x+1]);
X	max_alt = max(max_alt, data[y+1][x+1]);
X	return max_alt;
X}
END_OF_FILE
if test 17271 -ne `wc -c <'src/hf.c'`; then
    echo shar: \"'src/hf.c'\" unpacked with wrong size!
fi
# end of 'src/hf.c'
fi
if test -f 'src/input_yacc.y' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'src/input_yacc.y'\"
else
echo shar: Extracting \"'src/input_yacc.y'\" \(14190 characters\)
sed "s/^X//" >'src/input_yacc.y' <<'END_OF_FILE'
X/* input_yacc.y								    */
X/*									    */
X/* Copyright (C) 1989, Craig E. Kolb					    */
X/*									    */
X/* This software may be freely copied, modified, and redistributed,	    */
X/* provided that this copyright notice is preserved on all copies.	    */
X/* 									    */
X/* There is no warranty or other guarantee of fitness for this software,    */
X/* it is provided solely "as is".  Bug reports or fixes may be sent	    */
X/* to the author, who may or may not act on them as he desires.		    */
X/*									    */
X/* You may not include this software in a program or other software product */
X/* without supplying the source, or without informing the end-user that the */
X/* source is available for no extra charge.				    */
X/*									    */
X/* If you modify this software, you should include a notice giving the	    */
X/* name of the person performing the modification, the date of modification,*/
X/* and the reason for such modification.				    */
X/*									    */
X/* $Id: input_yacc.y,v 3.0 89/10/27 02:05:53 craig Exp $ */
X%{
X#include <stdio.h>
X#include "constants.h"
X#include "typedefs.h"
X#include "funcdefs.h"
X#include "texture.h"
X
Xint Npoints=0, CurXSize, CurYSize, CurZSize;
XObject *LastObj = (Object *)0;
XObjList *CurObj, *ListTmp;
XSurface *stmp;
XTexture *CurText;
XTransInfo *CurTrans = (TransInfo *)0, CurITrans;
XPointList *Polypoints, *Point;
Xextern FILE *yyin;
Xextern Object *World;
Xextern int WorldXSize, WorldYSize, WorldZSize, nlight, Xres, Yres, maxlevel;
Xextern int yylineno, Jittered, JitSamples, pixel_div;
Xextern double hfov, vfov, RedContrast, GreenContrast, BlueContrast;
Xextern double TreeCutoff;
Xextern Vector eyep, lookp, up;
Xextern char outfilename[];
Xextern Color background;
Xextern SurfaceList *Surfaces;
Xextern Light light[];
Xextern Fog *GlobalFog;
Xextern Mist *GlobalMist;
X%}
X%union {
X	char *c;
X	int i;
X	double d;
X	Vector v;
X	Color col;
X	struct Texture *text;
X}
X%token <i> tINT
X%token <d> tFLOAT
X%token <c> tSTRING
X%token tADAPTIVE tBACKGROUND tBLOTCH tBOX tBUMP tCONE tCYL tDIRECTIONAL
X%token tENDDEF tEXTENDED tEYEP tFBM tFBMBUMP tFOG tFOV tGRID
X%token tHEIGHTFIELD tJITTERED tLIGHT tLIST tLOOKP tMARBLE tMAXDEPTH tMIST
X%token tOBJECT tOUTFILE
X%token tPLANE tPOINT tPOLY tROTATE tSAMPLES
X%token tSCALE tSCREEN tSPHERE tSTARTDEF tSUPERQ tSURFACE tRESOLUTION
X%token tTHRESH tTRANSLATE tTRANSFORM tTRIANGLE tUP tENDFILE
X%token tTEXTURE tCHECKER tWOOD tCONTRAST tCUTOFF
X%type <d> Fnumber
X%type <c> String
X%type <v> Vector
X%type <col> Color
X%type <text> Texturetype
X%%
XItems		: /* empty */
X		| Items Item
X		;
XItem		: Eyep
X		| Lookp
X		| Up
X		| Fov
X		| Screen
X		| Maxdepth
X		| Samples
X		| Jittered
X		| Adaptive
X		| Contrast
X		| Cutoff
X		| Background
X		| Light
X		| Primitive
X		| Child
X		| Surface
X		| Outfile
X		| List
X		| Grid
X		| Object
X		| Fog
X		| Mist
X		| tENDFILE		/* For backward compatibility */
X		;
XList		: tLIST
X		{
X			if (CurObj->data)
X				CurObj->data->type = LIST;
X			else
X				World->type = LIST;
X		}
X		;
XGrid		: tGRID tINT tINT tINT
X		{
X			if (CurObj->data) {
X				CurObj->data->type = GRID;
X				CurXSize = $2;
X				CurYSize = $3;
X				CurZSize = $4;
X			} else {
X				World->type = GRID;
X				WorldXSize = $2;
X				WorldYSize = $3;
X				WorldZSize = $4;
X			}
X		}
X		;
XPrimitive	: Prim Textures
X		{
X			if (LastObj)
X				/* User may have botched prim. def. */
X				LastObj->texture = CurText;
X			CurText = (Texture *)0;
X			LastObj = (Object *)0;
X		}
X		;
XPrim		: Primtype Transforms
X		{
X			if (LastObj != (Object *)0) {
X				/*
X				 * Compute boundings box of primitive.
X				 * if box's low X is > its high X,
X				 * then the prim is unbounded.
X				 */
X				set_prim_bounds(LastObj);
X				/*
X				 * Add primitive to current object
X				 * list.  When the obj. def. is complete
X				 * make_list() will calculate the
X				 * bounding box of the entire object.
X				 */
X				add_prim(LastObj, CurObj->data);
X				if (CurTrans) {
X					/*
X					 * Compute the bounding box of the
X					 * transformed object.
X					 */
X					transform_bounds(CurTrans,
X							LastObj->bounds);
X					invert_trans(&CurITrans, CurTrans);
X					LastObj->trans = new_trans(CurTrans,
X								&CurITrans);
X					free((char *)CurTrans);
X					CurTrans = (TransInfo *)0;
X				}
X			} else {
X				/*
X				 * There was something wrong with the def.
X				 */
X				if (CurTrans) {
X					free((char *)CurTrans);
X					CurTrans = (TransInfo *)0;
X				}
X			}
X		}
X		;
XPrimtype	: Plane
X		| Sphere
X		| Box
X		| Triangle
X		| Cylinder
X		| Cone
X		| Superq
X		| Poly
X		| HeightField
X		;
XObject		: Objectdef Textures
X		{
X			CurObj->data->texture = CurText;
X			CurText = (Texture *)0;
X			free((char *)CurObj);
X			CurObj = CurObj->next;
X		}
X		;
XObjectdef	: Startdef Objdefs tENDDEF
X		{
X			/*
X			 * Object definition.
X			 */
X			LastObj = (Object *)NULL;
X			if (CurObj->data->data == (char *)0) {
X				fprintf(stderr,"Warning: null object defined");
X				fprintf(stderr," (line %d)\n",yylineno);
X			} else {
X				if (CurObj->data->type == GRID) {
X					list2grid(CurObj->data, CurXSize,
X						CurYSize, CurZSize);
X				} else {
X					/*
X				 	 * Object is a list -- transform the
X				 	 * linked list (ObjList) into a List.
X				 	 */
X					make_list(CurObj->data);
X				}
X				/*
X			 	 * Add this new object to the list of
X			 	 * defined objects.
X			 	 */
X				add_to_objects(CurObj->data);
X			}
X		}
X		;
XStartdef	: tSTARTDEF String
X		/*
X		 * define <name>
X		 */
X		{
X			/*
X			 * Once we know the bounding box of this object
X			 * (and the user hasn't specified that this object
X			 * is stored in a List), then we enlist everything
X			 * it contains.
X			 * The new object's DATA field points to an ObjList
X			 * until the definition complete, when the ObjList
X			 * is then turned into either a Grid or a List.
X			 */
X			ListTmp = (ObjList *)Malloc(sizeof(ObjList));
X			ListTmp->data = new_object($2, LIST, (char *)NULL,
X						(Trans *)NULL);
X			ListTmp->next = CurObj;
X			CurObj = ListTmp;
X		}
X		;
XObjdefs		: Objdefs Objdef
X		|
X		;
XObjdef		: Primitive
X		| Surface
X		| Child
X		| List
X		| Grid
X		| Object
X		;
XTextures	: Textures Texture
X		|
X		;
XTexture		: tTEXTURE Texturetype Transforms
X		{
X			/*
X			 * Set transformation information.
X			 */
X			if (CurTrans) {
X				invert_trans(&CurITrans, CurTrans);
X				$2->trans = new_trans(CurTrans,&CurITrans);
X				free((char *)CurTrans);
X				CurTrans = (TransInfo *)NULL;
X			}
X			/*
X			 * Walk to the end of list of textures and
X			 * append new texture.  This is done so that
X			 * textures are implied in the expected order.
X			 */
X			{
X				Texture *tp;
X
X				$2->next = (Texture *)0;
X
X				if (CurText) {
X					for (tp=CurText;tp->next;tp=tp->next)
X							;
X					tp->next = $2;
X
X				} else {
X					CurText = $2;
X				}
X			}
X		}
X		;
XTexturetype	: tCHECKER String
X		{
X			$$ = NewCheckText($2);
X		}
X		| tBLOTCH Fnumber String
X		{
X			$$ = NewBlotchText($2, $3);
X		}
X		| tBUMP Fnumber
X		{
X			$$ = NewBumpText($2);
X		}
X		| tMARBLE
X		{
X			$$ = NewMarbleText((char *)NULL);
X		}
X		| tMARBLE String
X		{
X			$$ = NewMarbleText($2);
X		}
X		| tFBM Fnumber Fnumber Fnumber Fnumber tINT Fnumber
X		{
X			$$ = NewfBmText($2, $3, $4, $5, $6, $7, (char *)0);
X		}
X		| tFBM Fnumber Fnumber Fnumber Fnumber tINT Fnumber tSTRING
X		{
X			$$ = NewfBmText($2, $3, $4, $5, $6, $7, $8);
X		}
X		| tFBMBUMP Fnumber Fnumber Fnumber Fnumber tINT
X		{
X			$$ = NewfBmBumpText($2, $3, $4, $5, $6);
X		}
X		| tWOOD
X		{
X			$$ = NewWoodText();
X		}
X		;
XChild		: Childdef Textures
X		{
X			LastObj->texture = CurText;
X			CurText = (Texture *)0;
X			LastObj = (Object *)0;
X		}
X		;
XChilddef	: tOBJECT String Transforms
X		{
X			LastObj = add_child_named($2, CurObj->data);
X			if (CurTrans) {
X				transform_bounds(CurTrans, LastObj->bounds);
X				invert_trans(&CurITrans, CurTrans);
X				if (LastObj->trans) {
X					mmult(&LastObj->trans->obj2world,
X						CurTrans,
X						&LastObj->trans->obj2world);
X					mmult(&LastObj->trans->world2obj,
X						&CurITrans,
X						&LastObj->trans->world2obj);
X				} else
X					LastObj->trans = new_trans(CurTrans,
X								&CurITrans);
X				free((char *)CurTrans);
X				CurTrans = (TransInfo *)NULL;
X			}
X		}
X		;
XTransforms	: Transforms Transform
X		| /* empty */
X		;
XTransform	: tTRANSLATE Vector
X		{
X			if (CurTrans == (TransInfo *)0)
X				CurTrans = new_transinfo();
X			translate(CurTrans, &($2));
X		}
X		| tROTATE Vector Fnumber
X		{
X			if (CurTrans == (TransInfo *)0)
X				CurTrans = new_transinfo();
X
X			rotate(CurTrans, &($2), deg2rad($3));
X		}
X		| tSCALE Fnumber Fnumber Fnumber
X		{
X			if (CurTrans == (TransInfo *)0)
X				CurTrans = new_transinfo();
X			scale(CurTrans, $2, $3, $4);
X		}
X		| tTRANSFORM Fnumber Fnumber Fnumber
X				Fnumber Fnumber Fnumber
X				Fnumber Fnumber Fnumber
X		{
X			if (CurTrans == (TransInfo *)0)
X				CurTrans = new_transinfo();
X			explicit_trans(CurTrans,
X				$2, $3, $4, $5, $6, $7, $8, $9, $10,
X				0., 0., 0., CurTrans);
X		}
X		| tTRANSFORM Fnumber Fnumber Fnumber
X				Fnumber Fnumber Fnumber
X				Fnumber Fnumber Fnumber
X				Fnumber Fnumber Fnumber
X		{
X			if (CurTrans == (TransInfo *)0)
X				CurTrans = new_transinfo();
X			explicit_trans(CurTrans,
X				$2, $3, $4, $5, $6, $7, $8, $9, $10,
X				$11, $12, $13,CurTrans);
X		}
X		;
XEyep		: tEYEP Vector Transforms
X		{
X			eyep = $2;
X			if (CurTrans) {
X				transform_point(&eyep, CurTrans);
X				free((char *)CurTrans);
X				CurTrans = (TransInfo *)0;
X			}
X		}
X		;
XLookp		: tLOOKP Vector
X		{
X			lookp = $2;
X		}
X		;
XUp		: tUP Vector
X		{
X			up = $2;
X		}
X		;
XFov		: tFOV Fnumber Fnumber
X		{
X			hfov = $2; vfov = $3;
X		}
X		| tFOV Fnumber
X		{
X			hfov = $2;
X		}
X		;
XSamples		: tSAMPLES tINT
X		{
X			if (JitSamples == UNSET)
X				JitSamples = $2;
X		}
X		;
XAdaptive	: tADAPTIVE tINT
X		{
X			if (pixel_div == UNSET && !Jittered)
X				pixel_div = $2;
X		}
X		;
XContrast	: tCONTRAST Fnumber Fnumber Fnumber
X		{
X			if (RedContrast == UNSET ||
X			    GreenContrast == UNSET ||
X			    BlueContrast == UNSET) {
X				RedContrast = $2;
X				GreenContrast = $3;
X				BlueContrast = $4;
X			}
X		}
X		;
XCutoff		: tCUTOFF Fnumber
X		{
X			if (TreeCutoff == UNSET)
X				TreeCutoff = $2;
X		}
X		;
XJittered	: tJITTERED
X		{
X			if (pixel_div == UNSET)
X				Jittered = TRUE;
X		}
X		;
XScreen		: tSCREEN tINT tINT
X		{
X			if (Xres == UNSET || Yres == UNSET) {
X				Xres = $2;
X				Yres = $3;
X			}
X		}
X		| tRESOLUTION tINT tINT
X		{
X			if (Xres == UNSET || Yres == UNSET) {
X				Xres = $2;
X				Yres = $3;
X			}
X		}
X		;
XMaxdepth	: tMAXDEPTH tINT
X		{
X			maxlevel = $2;
X		}
X		;
XBackground	: tBACKGROUND Color
X		{
X			background = $2;
X		}
X		;
XLight		: Lightdef tPOINT Vector
X		{
X			light[nlight].pos = $3;
X			light[nlight].type = LOCAL;
X			nlight++;
X		}
X		| Lightdef tDIRECTIONAL Vector
X		{
X			(void)normalize(&($3));
X			light[nlight].pos = $3;
X			light[nlight].type = DIRECTIONAL;
X			nlight++;
X		}
X		| Lightdef tEXTENDED Vector Fnumber
X		{
X			light[nlight].pos = $3;
X			light[nlight].radius = $4;
X			light[nlight].type = EXTENDED;
X			nlight++;
X		}
X		;
XLightdef	: tLIGHT Fnumber
X		{
X			if (nlight == LIGHTS)
X				yyerror("Too many lights.");
X			light[nlight].color.r = $2;
X			light[nlight].color.g = $2;
X			light[nlight].color.b = $2;
X		}
X		| tLIGHT Color
X		{
X			if (nlight == LIGHTS)
X				yyerror("Too many lights.\n");
X			light[nlight].color = $2;
X		}
X		;
XSurface		: tSURFACE String
X			Color Color Color
X			Fnumber Fnumber Fnumber Fnumber
X		{
X			/*
X			 * surface name
X			 * 	amb
X			 * 	diff
X			 * 	spec coef reflect refract krefract
X			 */
X			stmp = make_surface($2, &($3), &($4), &($5), $6, $7,
X					$8, $9, 0., 0.);
X			Surfaces = add_surface(stmp, Surfaces);
X
X		}
X		| tSURFACE String Color Color Color
X			Fnumber Fnumber Fnumber Fnumber Fnumber Fnumber
X		{
X			/*
X			 * surface name
X			 * 	amb
X			 * 	diff
X			 * 	spec coef reflect refract krefract
X			 */
X			stmp = make_surface($2, &($3), &($4), &($5), $6, $7,
X					$8, $9, $10, $11);
X			Surfaces = add_surface(stmp, Surfaces);
X
X		}
X		;
XHeightField	: tHEIGHTFIELD tSTRING tSTRING
X		{
X			LastObj = makhf($2, $3);
X		}
X		;
XPoly		: tPOLY tSTRING Polypoints
X		{
X			LastObj = makpoly($2, Polypoints, Npoints);
X			Polypoints = (PointList *)0;
X			Npoints = 0;
X		}
X		;
XPolypoints	: /* empty */
X		| Polypoints Polypoint
X		;
XPolypoint	: Vector
X		{
X			Point = (PointList *)Malloc(sizeof(PointList));
X			Point->vec = $1;
X			Point->next = Polypoints;
X			Polypoints = Point;
X			Npoints++;
X		}
X		;
XCone		: tCONE tSTRING Vector Vector Fnumber Fnumber
X		{
X			CurTrans = new_transinfo();
X			LastObj = makcone($2, &($3), &($4), $5, $6, CurTrans);
X		}
X		;
XCylinder	: tCYL tSTRING Vector Vector Fnumber
X		{
X			/*
X			 * Cylinders automagically define a
X			 * transformation matrix.
X			 */
X			CurTrans = new_transinfo();
X			LastObj = makcyl($2, &($3), &($4), $5, CurTrans);
X		}
X		;
XSphere		: tSPHERE tSTRING Fnumber Vector
X		{
X			LastObj = maksph($2, $3, &($4));
X		}
X		;
XBox		: tBOX tSTRING
X			Fnumber Fnumber Fnumber
X			Fnumber Fnumber Fnumber
X		{
X			LastObj = makbox($2, $3, $4, $5, $6, $7, $8);
X		}
X		;
XTriangle	: tTRIANGLE tSTRING Vector Vector Vector
X		{
X			LastObj = maktri(TRIANGLE, $2, &($3), &($4), &($5),
X					(Vector *)0, (Vector *)0, (Vector *)0);
X		}
X		| tTRIANGLE tSTRING Vector Vector Vector Vector Vector Vector
X		{
X			LastObj = maktri(PHONGTRI, $2, &($3), &($5),
X					&($7), &($4), &($6), &($8));
X		}
X		;
XSuperq		: tSUPERQ tSTRING
X			Fnumber Fnumber Fnumber
X			Fnumber Fnumber Fnumber
X			Fnumber
X		{
X			LastObj = maksup($2, $3, $4, $5, $6, $7, $8, $9);
X		}
X		;
XPlane		: tPLANE tSTRING Vector Vector
X		{
X			LastObj = makplane($2, &($3), &($4));
X		}
X		;
XOutfile		: tOUTFILE String
X		{
X			if (*outfilename != (char)NULL)
X				fprintf(stderr,"Ignoring output name \"%s\"\n",
X					$2);
X			else
X				strcpy(outfilename, $2);
X		}
X		;
XMist		: tMIST Color Color Fnumber Fnumber
X		{
X			GlobalMist = (Mist *)Malloc(sizeof(Mist));
X			GlobalMist->color = $2;
X			GlobalMist->trans = $3;
X			GlobalMist->zero = $4;
X			GlobalMist->scale = 1. / $5;
X		}
X		;
XFog		: tFOG Fnumber Color
X		{
X			GlobalFog = (Fog *)Malloc(sizeof(Fog));
X			GlobalFog->trans = 1./$2;
X			GlobalFog->color = $3;
X		}
X		;
XColor		: Fnumber Fnumber Fnumber
X		{
X			$$.r = $1;
X			$$.g = $2;
X			$$.b = $3;
X		}
X		;
XVector		: Fnumber Fnumber Fnumber
X		{
X			$$.x = $1;
X			$$.y = $2;
X			$$.z = $3;
X		}
X		;
XFnumber		: tFLOAT
X		{ $$ = $1;}
X		| tINT
X		{ $$ = $1;}
X		;
XString		: tSTRING
X		{ $$ = $1;}
X%%
Xyyerror(s)
Xchar *s;
X{
X	extern char tmpname[];
X
X	fprintf(stderr,"rayshade: line %d: %s\n", yylineno, s);
X	unlink(tmpname);
X	exit(2);
X}
X
END_OF_FILE
if test 14190 -ne `wc -c <'src/input_yacc.y'`; then
    echo shar: \"'src/input_yacc.y'\" unpacked with wrong size!
fi
# end of 'src/input_yacc.y'
fi
if test -f 'src/raytrace.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'src/raytrace.c'\"
else
echo shar: Extracting \"'src/raytrace.c'\" \(15323 characters\)
sed "s/^X//" >'src/raytrace.c' <<'END_OF_FILE'
X/*
X * raytrace.c
X *
X * Copyright (C) 1989, Craig E. Kolb
X *
X * This software may be freely copied, modified, and redistributed,
X * provided that this copyright notice is preserved on all copies.
X *
X * There is no warranty or other guarantee of fitness for this software,
X * it is provided solely .  Bug reports or fixes may be sent
X * to the author, who may or may not act on them as he desires.
X *
X * You may not include this software in a program or other software product
X * without supplying the source, or without informing the end-user that the
X * source is available for no extra charge.
X *
X * If you modify this software, you should include a notice giving the
X * name of the person performing the modification, the date of modification,
X * and the reason for such modification.
X *
X * $Id: raytrace.c,v 3.0 89/10/27 02:06:02 craig Exp $
X *
X * $Log:	raytrace.c,v $
X * Revision 3.0  89/10/27  02:06:02  craig
X * Baseline for first official release.
X * 
X */
X
X/*
X * This module could use some work.  In particular, a better antialiasing
X * scheme would be nice.
X */
X
X#ifdef LINDA
X/*
X * If you're using linda, be sure to *move* this file to "raytrace.cl"
X * before compiling.
X */
X#include "linda.h"
X#endif
X#include <math.h>
X#include <stdio.h>
X#include "typedefs.h"
X#include "constants.h"
X#include "funcdefs.h"
X#include "raytrace.h"
X
XColor	*pixel_buf[2], background;	/* Point buffer, background color */
XColor	*out_buf;			/* Output pixel buffer */
Xint	pixel_div = UNSET;		/* max times to subdivide pixel */
Xint	JitSamples = UNSET;		/* sqrt of # jittered samples/pixel */
Xdouble	JitterWeight;			/* 1. / (Total Samples Taken) */
Xint	Jittered;			/* Use distributed raytracing */
Xint	*SampleNumbers;
Xint	SampleNumber;
Xint	StartLine = UNSET;
Xdouble	RedContrast = UNSET, GreenContrast = UNSET, BlueContrast = UNSET;
Xdouble	SampleSpacing;
X
X#ifdef LINDA
Xextern	int Workers;
Xint adapt_worker(), dist_worker();
X#endif
X
Xextern	int Xres, Yres;
Xextern	Vector eyep, firstray;
XRay	TopRay;				/* Top-level ray. */
X
X#ifdef LINDA
X
Xreal_main(argc, argv)
Xint argc;
Xchar **argv;
X{
X	/*
X	 * Linda wants all routines, including lmain(), to be in a single
X	 * file.  So, we have to perform this bit of silliness.
X	 */
X	return rayshade_main(argc, argv);
X}
X
X#endif
X
Xraytrace()
X{
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 = eyep;
X	TopRay.media = (SurfaceList *)0;
X	TopRay.shadow = FALSE;
X
X	out_buf = (Color *)Malloc(Xres * sizeof(Color));
X
X	if (Jittered)
X		distributed_trace();
X	else
X		adaptive_trace();
X}
X
X/*
X * Raytrace an image using "distributed" raytracing.
X */
Xdistributed_trace()
X{
X	register int y;
X	extern FILE *fstats;
X	extern unsigned long EyeRays;
X	extern int Verbose;
X
X	switch (JitSamples) {
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		default:
X			fprintf(stderr,"Sorry, %d rays/pixel not supported.\n",
X				JitSamples*JitSamples);
X			exit(2);
X	}
X
X	JitterWeight= 1. / (JitSamples * JitSamples);
X	SampleSpacing = 1. / JitSamples;
X
X#ifdef LINDA
X	/*
X	 * Linda strategy:
X	 * There is one tuple, named "scaninfo", which holds the number of
X	 * the next line to be worked on.  A worker ins the scaninfo tuple,
X	 * notes its value, and increments it before outing it again.
X	 * The supervisor ins finished scanlines in order and writes them to
X	 * the output file.
X	 */
X	fprintf(fstats,"Using %d workers.\n",Workers);
X	out("scaninfo", StartLine);
X	for (y = 0; y < Workers; y++)
X		eval("worker", dist_worker());
X	for (y = StartLine; y >= 0 ; y--) {
X		in("result", y, ? out_buf);
X		if (Verbose)
X			fprintf(stderr,"Supervisor: inned %d\n",y);
X		if (y % 10 == 0)
X			fprintf(fstats, "Finished line %d.\n",y);
X		outline(out_buf);
X	}
X	for (y = 0; y < Workers; y++)
X		in("worker", ? int);
X#else
X	/*
X	 * Trace each scanline, writing results to output file.
X	 */
X	for (y = StartLine; y >= 0; y--) {
X		trace_jit_line(y, out_buf);
X		outline(out_buf);
X		if (y % 10 == 0) {
X			fprintf(fstats,"Finished line %d (%ld rays)\n",y,
X							EyeRays);
X			fflush(fstats);
X		}
X	}
X#endif
X}
X
Xtrace_jit_line(line, buf)
Xint line;
XColor *buf;
X{
X	register int x;
X	for (x = 0; x < Xres; x++)
X		trace_jit_pixel(x, line, buf++);
X}
X
Xtrace_jit_pixel(xp, yp, color)
Xint xp, yp;
XColor *color;
X{
X	Color tmp;
X	double x, y, xpos, ypos;
X	int i, j, index;
X
X	ypos = (double)yp - 0.5;
X	color->r = color->g = color->b = 0.;
X	index = 0;
X	for (i = 0; i < JitSamples; i++, ypos += SampleSpacing) {
X		xpos = (double)xp - 0.5;
X		for (j = 0; j < JitSamples; j++, xpos += SampleSpacing) {
X			x = xpos + nrand() * SampleSpacing;
X			y = ypos + nrand() * SampleSpacing;
X			SampleNumber = SampleNumbers[index++];
X			trace_point(x, y, &tmp);
X			color->r += tmp.r;
X			color->g += tmp.g;
X			color->b += tmp.b;
X		}
X	}
X	ScaleColor(JitterWeight, *color, color);
X}
X
X/*
X * Raytrace an image using adaptive supersampling to perform antialising.
X */
Xadaptive_trace()
X{
X	register int line;
X	extern unsigned long EyeRays;
X	extern int maxlevel, Verbose;
X	extern FILE *fstats;
X
X	/*
X	 * In the adaptive supersampling case, Jitter, JitterWeight,
X	 * and SampleSpacing refer are used in sampling extended
X	 * light sources.  JitterWeight has a -4 in the denominator
X	 * due to the fact that we don't sample the corners of extended
X	 * sources when performing adaptive supersampling.
X	 */
X	JitterWeight = 1. / (JitSamples * JitSamples - 4);
X	SampleSpacing = 1. / JitSamples;
X
X	pixel_buf[0] = (Color *)Malloc(sizeof(Color) *
X				(unsigned)(Xres + 1));
X	pixel_buf[1] = (Color *)Malloc(sizeof(Color) *
X				(unsigned)(Xres + 1));
X	/*
X	 * Minimum pixel square size.
X	 */
X	Minsquare = 1./pow(2.,(double)pixel_div);
X	/*
X	 * At any time, there can be a maximum of 3 * pixel_div + 1
X	 * pixel squares on the stack.
X	 */
X	SquareStack = (pixel_square *)Malloc((unsigned)(3 * pixel_div + 1) *
X				sizeof(pixel_square));
X	/*
X	 * A pixel is treated as a square through whose corners rays
X	 * are traced.  If the color values at the corners are
X	 * "too different" (as measured by pixel_ok()), the square is
X	 * divided into four squares (tracing 5 additional rays)
X	 * and the process is repeated on the four new squares.
X	 * The color assigned to a square is the average of the 4 corners.
X	 * Note that this means that adjacent super-sampled pixels
X	 * cause the same ray to be traced more than once.
X	 * This scheme is adequate but far from wonderful.
X	 */
X#ifdef LINDA
X	/*
X	 * This is a bit more complicated than 'jittered' sampling,
X	 * as workers must know about other scanlines.
X	 *
X	 * The basic idea is to have a "scanline" tuple which
X	 * holds the last pixline handed out and the last scanline
X	 * completed.  This should be improved by keeping a tuple
X	 * containing the last completed scanline/pixline, and not
X	 * just the last ones assigned.  (The problem being that a
X	 * worker can try to in a pixline while another worker
X	 * is still working on it.)
X	 */
X	fprintf(fstats,"Using %d workers.\n",Workers);
X	out("scaninfo", StartLine+1, StartLine+2);
X	for (line = 0; line < Workers; line++)
X		eval("worker", adapt_worker());
X
X	/*
X	 * in() the scanlines in order, writing to output file as we go.
X	 * This loop can easily be modified to make the supervisor
X	 * do some work, too.
X	 */
X	for (line = StartLine; line >= 0;) {
X		if (!adapt_job(TRUE))
X			sleep(5);
X		while (inp("result", line, ? out_buf)) {
X			if (Verbose)
X				fprintf(stderr,"Supervisor: inned %d\n",line);
X			if (line % 10 == 0)
X				fprintf(fstats, "Finished line %d.\n",line);
X			outline(out_buf);
X			if (--line < 0)
X				break;
X		}
X	}
X	if (Verbose)
X		fprintf(stderr,"Finished -- inning workers.\n");
X	for (line = 0; line < Workers; line++)
X		in("worker", ? int);
X#else
X	line = StartLine + 1;
X	trace_line(line, &pixel_buf[line & 01][0]);
X	/*
X	 * Work bottom-to-top, as that's the way Utah-raster wants to
X	 * operate.
X	 */
X	for(line--;line >= 0;line--) {
X		trace_line(line, &pixel_buf[line & 01][0]);
X		subdivide_line(line, pixel_buf[line & 01],
X				     pixel_buf[(line+1) & 01],
X				     out_buf);
X		outline(out_buf);
X		if(line % 10 == 0) {
X			fprintf(fstats,"Finished line %d (%ld rays)\n",line,
X								EyeRays);
X			fflush(fstats);
X		}
X	}
X#endif
X}
X
X/*
X * Trace a line of sample points along "line".
X */
Xtrace_line(line, buf)
Xint line;
XColor *buf;
X{
X	register int x;
X	/*
X	 * We need to trace Xres + 1 rays.
X	 */
X	for(x = 0; x <= Xres;x++)
X		trace_point((double)x, (double)line, buf + x);
X}
X
X/*
X * Given the two arrays of sample points which define the upper and
X * lower edges of all the pixel squares in line "y," push each
X * square in turn on the pixelsquare stack and determine a color value
X * for the pixel by calling subdivide_square().
X */
Xsubdivide_line(y, upper, lower, buf)
Xint y;
XColor *upper, *lower, *buf;
X{
X	register int x;
X
X	/*
X	 * Upper is the array of
X	 * sample values which define the "upper" part (corners) of this
X	 * row, while lower are the "lower" corners.  For the
X	 * next (lower) row, the current "upper" becomes "lower".
X	 */
X	for(x = 0; x < Xres;x++) {
X		SquareStack[0].x = (float)x;
X		SquareStack[0].y = (float)y;
X		SquareStack[0].size = 1.0;
X		SquareStack[0].ul = upper[x];
X		SquareStack[0].ur = upper[x+1];
X		SquareStack[0].ll = lower[x];
X		SquareStack[0].lr = lower[x+1];
X		subdivide_square(&buf[x]);
X	}
X}
X
X/*
X * Bleck, but it saves us a function call and keeps the code much cleaner.
X */
X#define push_square(u,v,s,a,b,c,d) {\
X			Stackp->x = u; \
X			Stackp->y = v; \
X			Stackp->size = s; \
X			Stackp->ul = a; \
X			Stackp->ur = b; \
X			Stackp->ll = c; \
X			Stackp->lr = d; \
X			Stackp++;}
X
X/*
X * Subdivide a pixel square.
X */
Xsubdivide_square(color)
XColor *color;
X{
X	register pixel_square *Stackp;
X	register float x, y;
X	float size, halfsize;
X	double avfact, xdelta, ydelta;
X	Color ul, ur, ll, lr, u, d, l, r, c;
X	extern unsigned long SuperSampled;
X
X	color->r = color->g = color->b = 0.;
X	Stackp = SquareStack + 1;
X
X	do {
X		Stackp--;
X		size = Stackp->size;
X		ul = Stackp->ul;
X		ur = Stackp->ur;
X		ll = Stackp->ll;
X		lr = Stackp->lr;
X
X		if(size <= Minsquare || pixel_ok(&ul,&ur,&ll,&lr)) {
X			/*
X			 * The square is either the smallest allowed, or
X			 * the four corners of the square are similar.
X			 * Average the four corners (weighted by the
X			 * size of the square) to get this square's
X			 * contribution to the whole pixel's color.
X			 */
X			avfact = (size * size) * 0.25;
X			color->r += (ul.r + ur.r + ll.r + lr.r) * avfact;
X			color->g += (ul.g + ur.g + ll.g + lr.g) * avfact;
X			color->b += (ul.b + ur.b + ll.b + lr.b) * avfact;
X			continue;
X		}
X		/*
X		 * Subdivide into four squares -- trace 5 additional
X		 * rays and push the appropriate squares on the pixelsquare
X		 * stack.
X		 */
X		x = Stackp->x;
X		y = Stackp->y;
X		halfsize = size * 0.5;
X		xdelta = (double)(x + halfsize);
X		ydelta = (double)(y + halfsize);
X		trace_point(xdelta, (double)y, &u);
X		trace_point((double)x, ydelta, &l);
X		trace_point(xdelta, ydelta, &c);
X		trace_point((double)(x + size),ydelta, &r);
X		trace_point(xdelta, (double)(y + size), &d);
X		if(size == 1.)
X			SuperSampled++;
X		push_square(x, y, halfsize, ul, u, l, c);
X		push_square((float)xdelta, y, halfsize, u, ur, c, r);
X		push_square(x, (float)ydelta, halfsize, l, c, ll, d);
X		push_square((float)xdelta, (float)ydelta, halfsize,
X					c, r, d, lr);
X	} while (Stackp != SquareStack);
X}
X
X/*
X * Trace a ray through x, y on the screen, placing the result in "color."
X */
Xtrace_point(x, y, color)
Xdouble x, y;
XColor *color;
X{
X	double dist;
X	HitInfo hitinfo;
X	extern Vector scrnx, scrny;
X	extern unsigned long EyeRays;
X	extern double TraceRay();
X
X	/*
X	 * Calculate ray direction.
X	 */
X	EyeRays++;
X	TopRay.dir.x = firstray.x + x*scrnx.x - y*scrny.x;
X	TopRay.dir.y = firstray.y + x*scrnx.y - y*scrny.y;
X	TopRay.dir.z = firstray.z + x*scrnx.z - y*scrny.z;
X
X	(void)normalize(&TopRay.dir);
X
X	/*
X	 * Do the actual ray trace.
X	 */
X	dist = TraceRay((Primitive *)NULL, &TopRay, &hitinfo);
X	if (dist > 0.)
X		/*
X		 * There was a valid intersection.
X		 */
X		ShadeRay(&hitinfo, &TopRay, dist, &background, color, 1.0);
X	else
X		/* Use background color */
X		*color = background;
X}
X
X/*
X * Return TRUE if this pixel is okay and doesn't need to be supersampled,
X * FALSE otherwise.
X */
Xpixel_ok(w,x,y,z)
XColor *w, *x, *y, *z;
X{
X	double rmax, rmin, gmax, gmin, bmax, bmin;
X	double rsum, gsum, bsum;
X
X	/*
X	 * Find min & max R, G, & B.
X	 */
X	rmax = max(w->r, max(x->r, max(y->r, z->r)));
X	rmin = min(w->r, min(x->r, min(y->r, z->r)));
X	gmax = max(w->g, max(x->g, max(y->g, z->g)));
X	gmin = min(w->g, min(x->g, min(y->g, z->g)));
X	bmax = max(w->b, max(x->b, max(y->b, z->b)));
X	bmin = min(w->b, min(x->b, min(y->b, z->b)));
X
X	/*
X	 * Contrast is defined as (Max - Min) / (Max + Min) for each
X	 * of RG&B.  If any of these values is greater than the maximum
X	 * allowed, we have to supersample.
X	 */
X	rsum = rmax + rmin;
X	gsum = gmax + gmin;
X	bsum = bmax + bmin;
X	if ((rsum == 0. || (rmax - rmin) / rsum < RedContrast) &&
X	    (gsum == 0. || (bmax - bmin) / gsum < BlueContrast) &&
X	    (bsum == 0. || (gmax - gmin) / bsum < GreenContrast))
X		return TRUE;
X
X	return FALSE;
X}
X
X#ifdef LINDA
Xdist_worker()
X{
X	while (dist_job())
X		;
X	return;
X}
X
X/*
X * Worker trained to perform distributed ray-tracing.
X */
Xdist_job()
X{
X	int y;
X	extern int Verbose;
X
X	in("scaninfo", ? y);
X	if (y < 0) {
X		out("scaninfo", y);
X		return 0;
X	}
X	if (Verbose)
X		fprintf(stderr,"Worker: inned %d\n",y);
X	out("scaninfo", y-1);
X	trace_jit_line(y, out_buf);
X	if (Verbose)
X		fprintf(stderr,"Worker: outing %d\n",y);
X	out("result", y, out_buf : Xres);
X	return 1;
X}
X
Xadapt_worker()
X{
X	while (adapt_job(FALSE))
X		;
X	return;
X}
X
Xadapt_job(supervisor)
Xint supervisor;
X{
X	int lastpix, lastscan;
X	extern int Verbose;
X
X	in("scaninfo", ? lastpix, ? lastscan);
X	if (lastpix <= 0) {
X		out("scaninfo", lastpix, lastscan);
X		if (Verbose)
X			fprintf(stderr,"Worker:  all finished!\n");
X		return FALSE;
X	}
X
X	if (rdp("scanline", lastpix -1, ? pixel_buf[0]) &&
X	    inp("scanline", lastpix, ? pixel_buf[1])) {
X		lastpix--;
X		out("scaninfo", lastpix, lastscan);
X		if (Verbose)
X			fprintf(stderr,"%s: doing pixline %d\n",
X				supervisor ? "Supervisor" : "Worker",
X					lastpix);
X		subdivide_line(lastpix, pixel_buf[0], pixel_buf[1],
X					out_buf);
X		out("result", lastpix, out_buf : Xres);
X	} else if (supervisor) {
X		/*
X		 * Don't let the supervisor get caught up in
X		 * ray-tracing a whole scanline.  It might take
X		 * a long, long time, causing tuple-space to get
X		 * jammed with finished pixlines, and...
X		 */
X		if (Verbose)
X			fprintf(stderr,"Supervisor: nothing to do...\n");
X		out ("scaninfo", lastpix, lastscan);
X		return FALSE;
X	} else if (lastscan > 0) {
X		lastscan--;
X		out("scaninfo", lastpix, lastscan);
X		if (Verbose)
X			fprintf(stderr,"Worker: doing scan %d\n",
X					lastscan);
X		trace_line(lastscan, pixel_buf[0]);
X		out("scanline", lastscan, pixel_buf[0] : Xres + 1);
X	} else {
X		/*
X		 * Nothing to do until somebody finishes a scanline.
X		 */
X		if (Verbose) {
X			fprintf(stderr,"Worker idle... ");
X			fprintf(stderr,"pix = %d, scan = %d\n",
X					lastpix, lastscan);
X		}
X		out("scaninfo", lastpix, lastscan);
X		sleep(2);
X	}
X	return 1;
X}
X#endif
END_OF_FILE
if test 15323 -ne `wc -c <'src/raytrace.c'`; then
    echo shar: \"'src/raytrace.c'\" unpacked with wrong size!
fi
# end of 'src/raytrace.c'
fi
echo shar: End of archive 6 \(of 8\).
cp /dev/null ark6isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 8 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0


-- 
Please send comp.sources.unix-related mail to rsalz at uunet.uu.net.
Use a domain-based address or give alternate paths, or you may lose out.