#! /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 16 (of 19)."
# Contents:  libray/libobj/hf.c raypaint/render.c
# Wrapped by kolb@woody on Wed Jul 17 17:56:55 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'libray/libobj/hf.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'libray/libobj/hf.c'\"
else
echo shar: Extracting \"'libray/libobj/hf.c'\" \(17333 characters\)
sed "s/^X//" >'libray/libobj/hf.c' <<'END_OF_FILE'
X/*
X * hf.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: hf.c,v 4.0 91/07/17 14:38:15 kolb Exp Locker: kolb $
X *
X * $Log:	hf.c,v $
X * Revision 4.0  91/07/17  14:38:15  kolb
X * Initial version.
X * 
X */
X#include "geom.h"
X#include "hf.h"
X
Xstatic Methods *iHfMethods = NULL;
Xstatic char hfName[] = "heighfield";
X
Xstatic void integrate_grid(), QueueTri();
Xstatic int DDA2D(), CheckCell();
Xstatic Float intHftri();
Xstatic float minalt(), maxalt();
X
Xtypedef struct {
X	int stepX, stepY;
X	Float tDX, tDY;
X	float minz, maxz;
X	int outX, outY;
X	Vector cp, pDX, pDY;
X} Trav2D;
X
XhfTri *CreateHfTriangle(), *GetQueuedTri();
X
Xunsigned long HFTests, HFHits;
X
XHf *
XHfCreate(filename)
Xchar *filename;
X{
X	Hf *hf;
X	FILE *fp;
X	float val, *maxptr, *minptr;
X	int i, j;
X
X	fp = fopen(filename, "r");
X	if (fp == (FILE *)NULL) {
X		RLerror(RL_ABORT, "Cannot open heightfield file \"%s\".",
X				filename);
X		return (Hf *)NULL;
X	}
X
X	hf = (Hf *)Malloc(sizeof(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((char *)&hf->size, sizeof(int), 1, fp) == 0) {
X		RLerror(RL_ABORT, "Cannot read height field size.");
X		return (Hf *)NULL;
X	}
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((char *)hf->data[i],sizeof(float),hf->size,fp)
X		    != hf->size) {
X			RLerror(RL_ABORT, "Not enough heightfield data.");
X			return (Hf *)NULL;
X		}
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	(void)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((unsigned)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 ((Float)hf->lsize[i] != hf->spacing[i])
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 hf;
X}
X
XMethods *
XHfMethods()
X{
X	if (iHfMethods == (Methods *)NULL) {
X		iHfMethods = MethodsCreate();
X		iHfMethods->create = (GeomCreateFunc *)HfCreate;
X		iHfMethods->methods = HfMethods;
X		iHfMethods->name = HfName;
X		iHfMethods->intersect = HfIntersect;
X		iHfMethods->normal = HfNormal;
X		iHfMethods->uv = HfUV;
X		iHfMethods->bounds = HfBounds;
X		iHfMethods->stats = HfStats;
X		iHfMethods->checkbounds = TRUE;
X		iHfMethods->closed = FALSE;
X	}
X	return iHfMethods;
X}
X
X/*
X * Intersect ray with height field.
X */
Xint
XHfIntersect(hf, ray, mindist, maxdist)
XHf *hf;
XRay *ray;
XFloat mindist, *maxdist;
X{
X	Vector hitpos;
X	Float offset;
X	Trav2D trav;
X
X	HFTests++;
X
X	/*
X	 * Find where we hit the hf cube.
X	 */
X	VecAddScaled(ray->pos, mindist, ray->dir, &hitpos);
X	if (OutOfBounds(&hitpos, hf->boundbox)) {
X		offset = *maxdist;
X		if (!BoundsIntersect(ray, hf->boundbox, mindist, &offset))
X			return FALSE;
X		hitpos.x = ray->pos.x + ray->dir.x * offset;
X		hitpos.y = ray->pos.y + ray->dir.y * offset;
X		hitpos.z = ray->pos.z + ray->dir.z * offset;
X	} else
X		hitpos = ray->pos;
X	/*
X	 * Find out in which cell "hitpoint" is.
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->dir.x < 0.) {
X		trav.stepX = trav.outX = -1;
X		trav.tDX = -1. / (ray->dir.x * hf->spacing[hf->levels -1]);
X	} else if (ray->dir.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->dir.x * hf->spacing[hf->levels -1]);
X	}
X
X	if (ray->dir.y < 0.) {
X		trav.stepY = trav.outY = -1;
X		trav.tDY = -1. / (ray->dir.y * hf->spacing[hf->levels -1]);
X	} else if (ray->dir.y > 0.) {
X		trav.stepY = 1;
X		trav.outY = hf->lsize[hf->levels -1];
X		trav.tDY = 1. / (ray->dir.y * hf->spacing[hf->levels -1]);
X	}
X
X	trav.pDX.x = ray->dir.x * trav.tDX;
X	trav.pDX.y = ray->dir.y * trav.tDX;
X	trav.pDX.z = ray->dir.z * trav.tDX;
X	trav.pDY.x = ray->dir.x * trav.tDY;
X	trav.pDY.y = ray->dir.y * trav.tDY;
X	trav.pDY.z = ray->dir.z * trav.tDY;
X
X	trav.cp = hitpos;
X	trav.minz = hf->minz;
X	trav.maxz = hf->maxz;
X	if (DDA2D(hf, &ray->pos, &ray->dir, hf->levels -1, &trav, maxdist)) {
X		HFHits++;
X		return TRUE;
X	}
X	return FALSE;
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 */
Xstatic int
XDDA2D(hf, pos, ray, level, trav, maxdist)
XHf *hf;
XVector *pos, *ray;
Xint level;
XTrav2D *trav;
XFloat *maxdist;
X{
X	int x, y, size, posZ;
X	float **boundsmin, **boundsmax, spacing;
X	Float 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	if (x == size)
X		x--;
X	y = trav->cp.y * hf->spacing[level];
X	if (y == size)
X		y--;
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					 * Recurse -- 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					if (DDA2D(hf,pos,ray,level-1,
X						&newtrav, maxdist))
X						return TRUE;
X				} else if (CheckCell(x,y,hf,ray,pos,maxdist))
X					return TRUE;
X			}
X			x += trav->stepX;		/* Move in X */
X			if (*maxdist < tX || x == trav->outX)
X				/* If outside, quit */
X				return FALSE;
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					if (DDA2D(hf,pos,ray,level-1,
X							&newtrav, maxdist))
X						return TRUE;
X				} else if (CheckCell(x,y,hf,ray,pos,maxdist))
X					return TRUE;
X			}
X			y += trav->stepY;
X			if (*maxdist < tY || y == trav->outY)
X				return FALSE;
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 FALSE;
X}
X
X/*
X * Check for ray/cell intersection
X */
Xstatic int
XCheckCell(x, y, hf, ray, pos, maxdist)
Xint x, y;
XHf *hf;
XVector *ray, *pos;
XFloat *maxdist;
X{
X	hfTri *tri1, *tri2;
X	Float 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 FALSE;
X
X	if (d1 < d2) {
X		if (d1 < *maxdist) {
X			hf->hittri = *tri1;
X			*maxdist = d1;
X			return TRUE;
X		}
X		return FALSE;
X	}
X
X	if (d2 < *maxdist) {
X		hf->hittri = *tri2;
X		*maxdist = d2;
X		return TRUE;
X	}
X	return FALSE;
X}
X
Xstatic hfTri *
XCreateHfTriangle(hf, x1, y1, x2, y2, x3, y3, which)
XHf *hf;
Xint x1, y1, x2, y2, x3, y3, which;
X{
X	hfTri *tri;
X	Float 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 = (Float)x1 / (Float)(hf->size -1);
X	yid = (Float)y1 / (Float)(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 = (Float)x2 / (Float)(hf->size-1);
X	tri->v2.y = (Float)y2 / (Float)(hf->size-1);
X	tri->v2.z = hf->data[y2][x2];
X	tri->v3.x = (Float)x3 / (Float)(hf->size-1);
X	tri->v3.y = (Float)y3 / (Float)(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)VecNormCross(&tmp1, &tmp2, &tri->norm);
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 */
Xstatic Float
XintHftri(ray, pos, tri)
XhfTri *tri;
XVector *pos, *ray;
X{
X	Float 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 */
X/*ARGSUSED*/
Xint
XHfNormal(hf, pos, nrm, gnrm)
XHf *hf;
XVector *pos, *nrm, *gnrm;
X{
X	*gnrm = *nrm = hf->hittri.norm;
X	return FALSE;
X}
X
X/*ARGSUSED*/
Xvoid
XHfUV(hf, pos, norm, uv, dpdu, dpdv)
XHf *hf;
XVector *pos, *norm, *dpdu, *dpdv;
XVec2d *uv;
X{
X	uv->u = pos->x;
X	uv->v = pos->y;
X	if (dpdu) {
X		dpdu->x = 1.;
X		dpdu->y = dpdv->z = 0.;
X		dpdv->x = dpdv->z = 0.;
X		dpdv->y = 1.;
X	}
X}
X
X/*
X * Compute heightfield bounding box.
X */
Xvoid
XHfBounds(hf, bounds)
XHf *hf;
XFloat bounds[2][3];
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 */
Xstatic void
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 */
Xstatic void
XQueueTri(hf, tri)
XHf *hf;
XhfTri *tri;
X{
X	if (hf->q[hf->qtail])		/* Free old triangle data */
X		free((voidstar)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 */
Xstatic hfTri *
XGetQueuedTri(hf, x, y, flag)
XHf *hf;
XFloat 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 */
Xstatic float
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 */
Xstatic float
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}
X
Xchar *
XHfName()
X{
X	return hfName;
X}
X
Xvoid
XHfStats(tests, hits)
Xunsigned long *tests, *hits;
X{
X	*tests = HFTests;
X	*hits = HFHits;
X}
X
Xvoid
XHfMethodRegister(meth)
XUserMethodType meth;
X{
X	if (iHfMethods)
X		iHfMethods->user = meth;
X}
END_OF_FILE
if test 17333 -ne `wc -c <'libray/libobj/hf.c'`; then
    echo shar: \"'libray/libobj/hf.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/hf.c'
fi
if test -f 'raypaint/render.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'raypaint/render.c'\"
else
echo shar: Extracting \"'raypaint/render.c'\" \(17440 characters\)
sed "s/^X//" >'raypaint/render.c' <<'END_OF_FILE'
X/*
X * render.c
X *
X * Copyright (C) 1989, 1991 Craig E. Kolb, Rod G. Bogart
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: render.c,v 4.0 91/07/17 17:37:09 kolb Exp Locker: kolb $
X *
X * $Log:	render.c,v $
X * Revision 4.0  91/07/17  17:37:09  kolb
X * Initial version.
X * 
X */
X
X#include "rayshade.h"
X#include "libcommon/sampling.h"
X#include "libsurf/atmosphere.h"
X#include "viewing.h"
X#include "options.h"
X#include "stats.h"
X#include "picture.h"
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		*SampleNumbers;
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#define RFAC	0.299
X#define GFAC	0.587
X#define BFAC	0.114
X
X#define NOT_CLOSED 0
X#define CLOSED_PARTIAL   1
X#define CLOSED_SUPER 2
X/*
X * If a region has area < MINAREA, it is considered to be "closed",
X * (a permanent leaf).  Regions that meet this criterion
X * are drawn pixel-by-pixel rather.
X */
X#define MINAREA		9
X
X#define SQ_AREA(s)	(((s)->xsize+1)*((s)->ysize+1))
X
X#define PRIORITY(s)	((s)->var * SQ_AREA(s))
X
X#define INTENSITY(p)	((RFAC*(p)[0] + GFAC*(p)[1] + BFAC*(p)[2])/255.)
X
X#define OVERLAPS_RECT(s) (!Rectmode || \
X				((s)->xpos <= Rectx1 && \
X				 (s)->ypos <= Recty1 && \
X				 (s)->xpos+(s)->xsize >= Rectx0 && \
X				 (s)->ypos+(s)->ysize >= Recty0))
X
Xtypedef unsigned char RGB[3];
X
Xstatic RGB **Image;
Xstatic char **SampleMap;
X
X/*
X * Sample square
X */
Xtypedef struct SSquare {
X	short xpos, ypos, xsize, ysize;
X	short depth;
X	short leaf, closed;
X	float mean, var, prio;
X	struct SSquare *child[4], *parent;
X} SSquare;
X
XSSquare *SSquares, *SSquareCreate(), *SSquareInstall(), *SSquareSelect(),
X	*SSquareFetchAtMouse();
X
XFloat SSquareComputeLeafVar();
X
XRay	TopRay;				/* Top-level ray. */
Xint	Rectmode = FALSE,
X	Rectx0, Recty0, Rectx1, Recty1;
Xint	SuperSampleMode = 0;
X#define SSCLOSED (SuperSampleMode + 1)
X
XRender(argc, argv)
Xint argc;
Xchar **argv;
X{
X	/*
X	 * Do an adaptive trace, displaying results in a
X	 * window as we go.
X	 */
X	SSquare *cursq;
X	Pixel *pixelbuf;
X	int y, x;
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 = Camera.pos;
X	TopRay.media = (Medium *)0;
X	TopRay.depth = 0;
X	/*
X	 * Doesn't handle motion blur yet.
X	 */
X	TopRay.time = Options.framestart;
X
X	GraphicsInit(Screen.xsize, Screen.ysize, "rayview");
X	/*
X	 * Allocate array of samples.
X	 */
X	Image=(RGB **)Malloc(Screen.ysize*sizeof(RGB *));
X	SampleMap = (char **)Malloc(Screen.ysize * sizeof(char *));
X	for (y = 0; y < Screen.ysize; y++) {
X		Image[y]=(RGB *)Calloc(Screen.xsize, sizeof(RGB));
X		SampleMap[y] = (char *)Calloc(Screen.xsize, sizeof(char));
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	 * Take initial four samples
X	 */
X	SSquareSample(0, 0, FALSE);
X	SSquareSample(Screen.xsize -1, 0, FALSE);
X	SSquareSample(Screen.xsize -1, Screen.ysize -1, FALSE);
X	SSquareSample(0, Screen.ysize -1, FALSE);
X	SSquares = SSquareInstall(0, 0, Screen.xsize -1, Screen.ysize -1,
X				  0, (SSquare *) NULL);
X
X	for (y = 1; y <= 5; y++) {
X		/*
X		 * Create and draw every region at depth y.
X		 */
X		SSquareDivideToDepth(SSquares, y);
X	}
X		
X
X	while ((cursq = SSquareSelect(SSquares)) != (SSquare *)NULL) {
X		SSquareDivide(cursq);
X		if (GraphicsRedraw())
X			SSquareDraw(SSquares);
X		if (GraphicsMiddleMouseEvent()) 
X			SSetRectMode();
X	}
X
X	/*
X	 * Finished the image; write to image file.
X	 */
X	pixelbuf = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
X	PictureStart(argv);
X	for (y = 0; y < Screen.ysize; y++) {
X		/*
X		 * This is really disgusting.
X		 */
X		for (x = 0; x < Screen.xsize; x++) {
X			pixelbuf[x].r = Image[y][x][0] / 255.;
X			pixelbuf[x].g = Image[y][x][1] / 255.;
X			pixelbuf[x].b = Image[y][x][2] / 255.;
X			pixelbuf[x].alpha = SampleMap[y][x];
X		}
X		PictureWriteLine(pixelbuf);
X	}
X	PictureEnd();
X	free((voidstar)pixelbuf);
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
XSSetRectMode()
X{
X	int x1,y1,x2,y2;
X
X	if (Rectmode) {
X		Rectmode = FALSE;
X		RecomputePriority(SSquares);
X	}
X	GraphicsGetMousePos(&x1, &y1);
X	while (GraphicsMiddleMouseEvent())
X		;
X	GraphicsGetMousePos(&x2, &y2);
X	if (x1 < x2) {
X		Rectx0 = (x1 < 0) ? 0 : x1;
X		Rectx1 = (x2 >= Screen.xsize) ? Screen.xsize - 1 : x2;
X	} else {
X		Rectx0 = (x2 < 0) ? 0 : x2;
X		Rectx1 = (x1 >= Screen.xsize) ? Screen.xsize - 1 : x1;
X	} if (y1 < y2) {
X		Recty0 = (y1 < 0) ? 0 : y1;
X		Recty1 = (y2 >= Screen.ysize) ? Screen.ysize - 1 : y2;
X	} else {
X		Recty0 = (y2 < 0) ? 0 : y2;
X		Recty1 = (y1 >= Screen.ysize) ? Screen.ysize - 1 : y1;
X	}
X	Rectmode = TRUE;
X	/* setup current rect */
X	(void)RecomputePriority(SSquares);
X}
X
XRecomputePriority(sq)
XSSquare *sq;
X{
X	Float maxp;
X
X	if (!OVERLAPS_RECT(sq)) {
X		sq->closed = SSCLOSED;
X		return FALSE;
X	}
X
X	if (sq->leaf) {
X		if (SQ_AREA(sq) >= MINAREA)
X			sq->closed = NOT_CLOSED;
X		return TRUE;
X	}
X	maxp = 0.;
X	if (RecomputePriority(sq->child[0]))
X		maxp = max(maxp, sq->child[0]->prio);
X	if (RecomputePriority(sq->child[1]))
X		maxp = max(maxp, sq->child[1]->prio);
X	if (RecomputePriority(sq->child[2]))
X		maxp = max(maxp, sq->child[2]->prio);
X	if (RecomputePriority(sq->child[3]))
X		maxp = max(maxp, sq->child[3]->prio);
X	sq->prio = maxp;
X#if 0
X	if ((sq->child[0]->closed == CLOSED_SUPER) &&
X	    (sq->child[1]->closed == CLOSED_SUPER) &&
X	    (sq->child[2]->closed == CLOSED_SUPER) &&
X	    (sq->child[3]->closed == CLOSED_SUPER))
X		sq->closed = CLOSED_SUPER;
X	else if (sq->child[0]->closed && sq->child[1]->closed &&
X		 sq->child[2]->closed && sq->child[3]->closed)
X		sq->closed = CLOSED_PARTIAL;
X	else
X		sq->closed = NOT_CLOSED;
X#else
X	if ((sq->child[0]->closed >= SSCLOSED) &&
X	    (sq->child[1]->closed >= SSCLOSED) &&
X	    (sq->child[2]->closed >= SSCLOSED) &&
X	    (sq->child[3]->closed >= SSCLOSED))
X		sq->closed = SSCLOSED;
X	else
X		sq->closed = NOT_CLOSED;
X#endif
X	return TRUE;
X}
X
XSSquareSample(x, y, supersample)
Xint x, y, supersample;
X{
X	Float upos, vpos, u, v;
X	int xx, yy, xp, sampnum, usamp, vsamp;
X	Pixel ctmp;
X	Pixel p;
X	extern unsigned char correct();
X
X	if (SampleMap[y][x] >= 128 + supersample)
X		/* already a sample there */
X		return;
X	SampleMap[y][x] = 128 + supersample;
X	if (supersample) {
X		p.r = p.g = p.b = p.alpha = 0;
X		sampnum = 0;
X		xp = x + Screen.minx;
X		vpos = Screen.miny + y - 0.5*Sampling.filterwidth;
X		for (yy = 0; yy < Sampling.sidesamples; yy++,
X		     vpos += Sampling.filterdelta) {
X			upos = xp - 0.5*Sampling.filterwidth;
X			for (xx = 0; xx < Sampling.sidesamples; xx++,
X			     upos += Sampling.filterdelta) {
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				p.r += ctmp.r*Sampling.filter[xx][yy];
X				p.g += ctmp.g*Sampling.filter[xx][yy];
X				p.b += ctmp.b*Sampling.filter[xx][yy];
X				if (++sampnum == Sampling.totsamples)
X					sampnum = 0;
X			}
X		}
X	}
X	else {
X		sampnum = nrand() * Sampling.totsamples;
X		usamp = sampnum % Sampling.sidesamples;
X		vsamp = sampnum / Sampling.sidesamples;
X
X		vpos = Screen.miny + y - 0.5*Sampling.filterwidth
X			+ 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[sampnum]);
X		SampleScreen(upos, vpos, &TopRay, &p, SampleNumbers[sampnum]);
X	}
X	Image[y][x][0] = CORRECT(p.r);
X	Image[y][x][1] = CORRECT(p.g);
X	Image[y][x][2] = CORRECT(p.b);
X}
X
XSSquare *
XSSquareCreate(xp, yp, xs, ys, d, parent)
Xint xp, yp, xs, ys, d;
XSSquare *parent;
X{
X	SSquare *new;
X	Float i1, i2, i3, i4;
X
X	new = (SSquare *)Calloc(1, sizeof(SSquare));
X	new->xpos = xp; new->ypos = yp;
X	new->xsize = xs; new->ysize = ys;
X	new->depth = d;
X	new->parent = parent;
X	i1 = INTENSITY(Image[new->ypos][new->xpos]);
X	i2 = INTENSITY(Image[new->ypos+new->ysize][new->xpos]);
X	i3 = INTENSITY(Image[new->ypos+new->ysize][new->xpos+new->xsize]);
X	i4 = INTENSITY(Image[new->ypos][new->xpos+new->xsize]);
X	new->mean = 0.25 * (i1+i2+i3+i4);
X	if (SQ_AREA(new) < MINAREA) {
X		new->prio = 0;
X		new->closed = SSCLOSED;
X	} else {
X		new->var = SSquareComputeLeafVar(new, i1, i2, i3, i4);
X		new->prio = PRIORITY(new);
X		new->closed = NOT_CLOSED;
X	}
X	new->leaf = TRUE;
X	return new;
X}
X
XFloat
XSSquareComputeLeafVar(sq, i1, i2, i3, i4)
XSSquare *sq;
XFloat i1, i2, i3, i4;
X{
X	Float res, diff;
X
X	diff = i1 - sq->mean;
X	res = diff*diff;
X	diff = i2 - sq->mean;
X	res += diff*diff;
X	diff = i3 - sq->mean;
X	res += diff*diff;
X	diff = i4 - sq->mean;
X	return res + diff*diff;
X}
X
XSSquareDivideToDepth(sq, d)
XSSquare *sq;
Xint d;
X{
X	if (sq->depth == d)
X		return;
X	if (sq->leaf)
X		SSquareDivide(sq);
X	SSquareDivideToDepth(sq->child[0], d);
X	SSquareDivideToDepth(sq->child[1], d);
X	SSquareDivideToDepth(sq->child[2], d);
X	SSquareDivideToDepth(sq->child[3], d);
X}
X
XSSquareDivide(sq)
XSSquare *sq;
X{
X	int lowx, lowy, midx, midy, hix, hiy;
X	int newxsize, newysize, ndepth, supersample;
X	/*
X	 * Divide the square into fourths by tracing 12
X	 * new samples if necessary.
X	 */
X	newxsize = sq->xsize / 2;
X	newysize = sq->ysize / 2;
X	lowx = sq->xpos; lowy = sq->ypos;
X	midx = sq->xpos + newxsize;
X	midy = sq->ypos + newysize;
X	hix  = sq->xpos + sq->xsize;
X	hiy  = sq->ypos + sq->ysize;
X	ndepth = sq->depth + 1;
X	/* create new samples */
X	supersample = FALSE;
X	SSquareSample(midx, lowy, supersample);
X	SSquareSample(lowx, midy, supersample);
X	SSquareSample(midx, midy, supersample);
X	SSquareSample(hix,  midy, supersample);
X	SSquareSample(midx, hiy, supersample);
X#ifdef SHARED_EDGES
X	/* create and draw new squares */
X	sq->child[0] = SSquareInstall(lowx,lowy,newxsize,newysize,ndepth,sq);
X	sq->child[1] = SSquareInstall(midx, lowy, sq->xsize - newxsize,
X			newysize, ndepth,sq);
X	sq->child[2] = SSquareInstall(lowx, midy, newxsize,
X			sq->ysize - newysize, ndepth,sq);
X	sq->child[3] = SSquareInstall(midx, midy, sq->xsize - newxsize,
X			 sq->ysize - newysize, ndepth,sq);
X#else
X	/*
X	 *  draw additional samples in order to subdivide such that
X	 * edges of regions do not overlap
X	 */
X	SSquareSample(midx +1, lowy, supersample);
X	SSquareSample(midx +1, midy, supersample);
X	SSquareSample(lowx, midy +1, supersample);
X	SSquareSample(midx, midy +1, supersample);
X	SSquareSample(midx +1, midy +1, supersample);
X	SSquareSample(hix, midy +1, supersample);
X	SSquareSample(midx +1, hiy, supersample);
X
X	/* create and draw new squares */
X	sq->child[0] = SSquareInstall(lowx,lowy,
X				newxsize,newysize,ndepth,sq);
X	sq->child[1] = SSquareInstall(midx+1, lowy, sq->xsize - newxsize -1,
X			newysize, ndepth,sq);
X	sq->child[2] = SSquareInstall(lowx, midy+1, newxsize,
X			sq->ysize - newysize-1, ndepth,sq);
X	sq->child[3] = SSquareInstall(midx+1, midy+1, sq->xsize - newxsize-1,
X			 sq->ysize - newysize-1, ndepth,sq);
X#endif
X	sq->leaf = FALSE;
X	/*
X	 * Recompute parent's mean and variance.
X	 */
X	if (OVERLAPS_RECT(sq))
X		SSquareRecomputeStats(sq);
X}
X
XSSquareRecomputeStats(sq)
XSSquare *sq;
X{
X	Float maxp;
X	int in[4], closeflag;
X
X	in[0] = OVERLAPS_RECT(sq->child[0]);
X	in[1] = OVERLAPS_RECT(sq->child[1]);
X	in[2] = OVERLAPS_RECT(sq->child[2]);
X	in[3] = OVERLAPS_RECT(sq->child[3]);
X
X	if ((in[0] && (sq->child[0]->closed < SSCLOSED)) ||
X	    (in[1] && (sq->child[1]->closed < SSCLOSED)) ||
X	    (in[2] && (sq->child[2]->closed < SSCLOSED)) ||
X	    (in[3] && (sq->child[3]->closed < SSCLOSED))) {
X		maxp = 0.;
X		if (in[0])
X			maxp = max(maxp, sq->child[0]->prio);
X		if (in[1])
X			maxp = max(maxp, sq->child[1]->prio);
X		if (in[2])
X			maxp = max(maxp, sq->child[2]->prio);
X		if (in[3])
X			maxp = max(maxp, sq->child[3]->prio);
X		sq->closed = NOT_CLOSED;
X		sq->prio = maxp;
X	} else if ((sq->child[0]->closed == CLOSED_SUPER) &&
X		   (sq->child[1]->closed == CLOSED_SUPER) &&
X		   (sq->child[2]->closed == CLOSED_SUPER) &&
X		   (sq->child[3]->closed == CLOSED_SUPER)) {
X		sq->closed = CLOSED_SUPER;
X		sq->prio = 0;
X#if 0
X	} else if ((!in[0] || sq->child[0]->closed >= SSCLOSED) &&
X		   (!in[1] || sq->child[1]->closed >= SSCLOSED) &&
X		   (!in[2] || sq->child[2]->closed >= SSCLOSED) &&
X		   (!in[3] || sq->child[3]->closed >= SSCLOSED)) {
X		sq->closed = SSCLOSED;
X		sq->prio = 0;
X#endif
X	} else {
X		sq->closed = SSCLOSED;
X		sq->prio = 0;
X	}
X	if (sq->parent)
X		SSquareRecomputeStats(sq->parent);
X}
X
XSSquare *
XSSquareInstall(xp, yp, xs, ys, d, parent)
Xint xp, yp, xs, ys, d;
XSSquare *parent;
X{
X	SSquare *new;
X
X	new = SSquareCreate(xp, yp, xs, ys, d, parent);
X	SSquareDraw(new);
X	return new;
X}
X
XSSquare *
XSSquareSelect(list)
XSSquare *list;
X{
X	int i;
X	SSquare *res, *which;
X
X	/*
X	 * If mousebutton is pressed,
X	 * find the square in which the mouse is located and
X	 * return that if not a leaf (pixel-sized square).
X	 */
X	if (GraphicsLeftMouseEvent()) {
X		SuperSampleMode = 0;
X		if ((res = SSquareFetchAtMouse(list)) != (SSquare *)NULL)
X			return res;
X	}
X	else if (GraphicsRightMouseEvent()) {
X		SuperSampleMode = 1;
X		if ((res = SSquareFetchAtMouse(list)) != (SSquare *)NULL) {
X			return res;
X		}
X	}
X	if (list->closed >= SSCLOSED) {
X		if (Rectmode) {
X			Rectmode = FALSE;
X			RecomputePriority(SSquares);
X			return SSquareSelect(list);
X		}
X		return (SSquare *)NULL;
X	}
X	/*
X	 * Otherwise, find the square with the greatest
X	 * 'priority'.
X	 */
X	res = list;
X	while (res && !res->leaf) {
X		which = (SSquare *)NULL;
X		for (i = 0; i < 4; i++) {
X			if ((res->child[i]->closed < SSCLOSED) &&
X			    OVERLAPS_RECT(res->child[i])) {
X				which = res->child[i];
X				break;
X			}
X		}
X		while (++i < 4) {
X			if ((res->child[i]->closed < SSCLOSED) &&
X			    which->prio < res->child[i]->prio &&
X			    OVERLAPS_RECT(res->child[i]))
X				which = res->child[i];
X		}
X		res = which;
X	}
X	return res;
X}
X
XSSquare *
XSSquareFetchAtMouse(list)
XSSquare *list;
X{
X	SSquare *res;
X	int x, y;
X
X	/*
X	 * Get mouse position.
X	 */
X	GraphicsGetMousePos(&x, &y);
X	res = list;
X	while (!res->leaf && (res->closed < SSCLOSED)) { 
X		/*
X		 * Find in which child the mouse is located.
X		 */
X		if (x < res->child[1]->xpos) {
X			if (y < res->child[2]->ypos)
X				res = res->child[0];
X			else
X				res = res->child[2];
X		} else if (y < res->child[3]->ypos)
X			res = res->child[1];
X		else
X			res = res->child[3];
X	}
X	if (res->closed >= SSCLOSED)
X		return (SSquare *)NULL;
X	return res;
X}
X
XSSquareDraw(sq)
XSSquare *sq;
X{
X	if (SQ_AREA(sq) >= MINAREA)
X		GraphicsDrawRectangle(sq->xpos, sq->ypos, sq->xsize, sq->ysize,
X			Image[sq->ypos][sq->xpos],
X			Image[sq->ypos][sq->xpos+sq->xsize],
X			Image[sq->ypos+sq->ysize][sq->xpos+sq->xsize],
X			Image[sq->ypos+sq->ysize][sq->xpos]);
X	else
X		DrawPixels(sq->xpos, sq->ypos, sq->xsize, sq->ysize);
X	if (!sq->leaf) {
X		SSquareDraw(sq->child[0]);
X		SSquareDraw(sq->child[1]);
X		SSquareDraw(sq->child[2]);
X		SSquareDraw(sq->child[3]);
X	}
X}
X
XDrawPixels(xp, yp, xs, ys)
Xint xp, yp, xs, ys;
X{
X	int x, y, xi, yi;
X
X	yi = yp;
X	for (y = 0; y <= ys; y++, yi++) {
X		xi = xp;
X		for (x = 0; x <= xs; x++, xi++) {
X			SSquareSample(xi, yi, SuperSampleMode);
X			GraphicsDrawPixel(xi, yi, Image[yi][xi]);
X		}
X	}
X}
END_OF_FILE
if test 17440 -ne `wc -c <'raypaint/render.c'`; then
    echo shar: \"'raypaint/render.c'\" unpacked with wrong size!
fi
# end of 'raypaint/render.c'
fi
echo shar: End of archive 16 \(of 19\).
cp /dev/null ark16isdone
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