#! /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 17 (of 19)."
# Contents:  libray/libobj/blob.c raypaint/render.c
# Wrapped by cek@oyoy on Sun Feb  9 21:56:52 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'libray/libobj/blob.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'libray/libobj/blob.c'\"
else
echo shar: Extracting \"'libray/libobj/blob.c'\" \(18207 characters\)
sed "s/^X//" >'libray/libobj/blob.c' <<'END_OF_FILE'
X/*
X * blob.c
X *
X * Copyright (C) 1990, 1991, Mark Podlipec, 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: blob.c,v 4.0 91/07/17 14:36:07 kolb Exp Locker: kolb $
X *
X * $Log:	blob.c,v $
X * Revision 4.0  91/07/17  14:36:07  kolb
X * Initial version.
X * 
X */
X#include "geom.h"
X#include "blob.h"
X
Xstatic Methods *iBlobMethods = NULL;
Xstatic char blobName[] = "blob";
X
Xunsigned long BlobTests, BlobHits;
X
X/*
X * Blob/Metaball Description
X *
X * In this implementation a Blob is made up of a threshold T, and a 
X * group of 1 or more Metaballs.  Each metaball 'i'  is defined by 
X * its position (xi,yi,zi), its radius ri, and its strength si.
X * Around each Metaball is a density function di(Ri) defined by:
X *
X *         di(Ri) =  c4i * Ri^4  +  c2i * Ri^2  + c0i     0 <= Ri <= ri
X *         di(Ri) =  0                                    ri < Ri 
X *
X * where
X *       c4i  = si / (ri^4)
X *       c2i  = -(2 * si) / (ri^2)
X *       c0i  = si
X *       Ri^2 is the distance from a point (x,y,z) to the center of
X *            the Metaball - (xi,yi,zi).
X *
X * The density function looks sort of like a W (Float-u) with the 
X * center hump crossing the d-axis at  si  and the two humps on the
X * side being tangent to the R-axis at  -ri  and  +ri. Only the 
X * range  [0,ri] is being used. 
X * I chose this function so that derivative = 0  at the points  0 and +ri. 
X * This keeps the surface smooth when the densities are added. I also 
X * wanted no  Ri^3  and  Ri  terms, since the equations are easier to 
X * solve without them. This led to the symmetry about the d-axis and
X * the derivative being equal to zero at -ri as well.
X *
X * The surface of the Blob is defined by:
X *
X *
X *                  N
X *                 ---  
X *      F(x,y,z) = \    di(Ri)  = T
X *                 /
X *                 ---
X *                 i=0
X *
X * where
X *
X *     di(Ri)    is   given above
X *     Ri^2      =    (x-xi)^2  +  (y-yi)^2  + (z-zi)^2
X *      N        is   number of Metaballs in Blob
X *      T        is   the threshold
X *    (xi,yi,zi) is   the center of Metaball i
X *
X */
X
X/*****************************************************************************
X * Create & return reference to a metaball.
X */
XBlob *
XBlobCreate(T, mlist, npoints)
XFloat T;
XMetaList *mlist;
Xint npoints;
X{
X	Blob *blob;
X	int i;
X	MetaList *cur;
X
X/* 
X * There has to be at least one metaball in the blob.
X * Note: if there's only one metaball, the blob 
X * will be a sphere.
X */
X	if (npoints < 1)
X	{
X		RLerror(RL_WARN, "blob field not correct.\n");
X		return (Blob *)NULL;
X	}
X
X/*  
X * Initialize primitive and Geom structures
X */
X	blob = (Blob *)Malloc(sizeof(Blob));
X	blob->T = T;
X	blob->list=(MetaVector *)
X	    Malloc( (unsigned)(npoints*sizeof(MetaVector)) );
X	blob->num = npoints;
X
X/*
X * Initialize Metaball list
X */
X	for(i=0;i<npoints;i++)
X	{
X		cur = mlist;
X		if ( (cur->mvec.c0 < EPSILON) || (cur->mvec.rs < EPSILON) )
X		{
X			RLerror(RL_WARN,"Degenerate metaball in blob (sr).\n");
X			return (Blob *)NULL;
X		}
X		/* store radius squared */
X		blob->list[i].rs = cur->mvec.rs * cur->mvec.rs;
X		/* Calculate and store coefficients for each metaball */
X		blob->list[i].c0 = cur->mvec.c0;
X		blob->list[i].c2 = -(2.0 * cur->mvec.c0) / blob->list[i].rs;
X		blob->list[i].c4 = cur->mvec.c0 /
X				(blob->list[i].rs * blob->list[i].rs);
X		blob->list[i].x = cur->mvec.x;
X		blob->list[i].y = cur->mvec.y;
X		blob->list[i].z = cur->mvec.z;
X		mlist = mlist->next;
X		free((voidstar)cur);
X	}
X	/* 
X	 * Allocate room for Intersection Structures and
X	 * Allocate room for an array of pointers to point to
X	 * the Intersection Structures.
X	 */
X	blob->ilist = (MetaInt *)Malloc( 2 * blob->num * sizeof(MetaInt));
X	blob->iarr = (MetaInt **)Malloc( 2 * blob->num * sizeof(MetaInt *));
X	return blob;
X}
X
XMethods *
XBlobMethods()
X{
X	if (iBlobMethods == (Methods *)NULL) {
X		iBlobMethods = MethodsCreate();
X		iBlobMethods->create = (GeomCreateFunc *)BlobCreate;
X		iBlobMethods->methods = BlobMethods;
X		iBlobMethods->name = BlobName;
X		iBlobMethods->intersect = BlobIntersect;
X		iBlobMethods->normal = BlobNormal;
X		iBlobMethods->bounds = BlobBounds;
X		iBlobMethods->stats = BlobStats;
X		iBlobMethods->checkbounds = TRUE;
X		iBlobMethods->closed = TRUE;
X	}
X	return iBlobMethods;
X}
X
X/*****************************************************************************
X * Function used by qsort() when sorting the Ray/Blob intersection list
X */
Xint
XMetaCompare(A,B)
Xchar *A,*B;
X{
X	MetaInt **AA,**BB;
X
X	AA = (MetaInt **) A;
X	BB = (MetaInt **) B;
X	if (AA[0]->bound == BB[0]->bound) return(0);
X	if (AA[0]->bound <  BB[0]->bound) return(-1);
X	return(1);  /* AA[0]->bound is > BB[0]->bound */
X}
X
X/*****************************************************************************
X * Ray/metaball intersection test.
X */
Xint
XBlobIntersect(blob, ray, mindist, maxdist)
XBlob *blob;
XRay *ray;
XFloat mindist, *maxdist;
X{
X	double c[5], s[4];
X	Float dist;
X	MetaInt *ilist,**iarr;
X	register int i,j,inum;
X	extern void qsort();
X
X	BlobTests++;
X
X	ilist = blob->ilist;
X	iarr = blob->iarr;
X
X/*
X * The first step in calculating the Ray/Blob intersection is to 
X * divide the Ray into intervals such that only a fixed set of 
X * Metaballs contribute to the density function along that interval. 
X * This is done by finding the set of intersections between the Ray 
X * and each Metaball's Sphere/Region of influence, which has a 
X * radius  ri  and is centered at (xi,yi,zi).
X * Intersection information is kept track of in the MetaInt 
X * structure and consists of:
X *
X *   type    indicates whether this intersection is the start(R_START)
X *           of a Region or the end(R_END) of one. 
X *   pnt     the Metaball of this intersection
X *   bound   the distance from Ray origin to this intersection
X *
X * This list is then sorted by  bound  and used later to find the Ray/Blob 
X * intersection.
X */
X
X	inum = 0;
X	for(i=0; i < blob->num; i++)
X	{
X		register Float xadj, yadj, zadj;
X		register Float b, t, rs;
X		register Float dmin,dmax;
X
X		rs   = blob->list[i].rs;
X		xadj = blob->list[i].x - ray->pos.x;
X		yadj = blob->list[i].y - ray->pos.y;
X		zadj = blob->list[i].z - ray->pos.z;
X
X	/*
X	 * Ray/Sphere of Influence intersection
X	 */
X		b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
X		t = b * b - xadj * xadj - yadj * yadj - zadj * zadj + rs;
X
X	/* 
X	 * don't except imaginary or single roots. A single root is a ray
X	 * tangent to the Metaball's Sphere/Region. The Metaball's
X	 * contribution to the overall density function at this point is
X	 * zero anyway.
X	 */
X		if (t > 0.0) /* we have two roots */
X		{
X			t = sqrt(t);
X			dmin = b - t;
X	/* 
X	 * only interested in stuff in front of ray origin 
X	 */
X			if (dmin < mindist) dmin = mindist;
X			dmax = b + t;
X			if (dmax > dmin) /* we have a valid Region */
X			{
X		/*
X	         * Initialize min/start and max/end Intersections Structures
X	         * for this Metaball
X	         */
X				ilist[inum].type = R_START;
X				ilist[inum].pnt = i;
X				ilist[inum].bound = dmin;
X				for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
X				iarr[inum] = &(ilist[inum]);
X				inum++;
X
X				ilist[inum].type = R_END;
X				ilist[inum].pnt = i;
X				ilist[inum].bound = dmax;
X				for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
X				iarr[inum] = &(ilist[inum]);
X				inum++;
X			} /* End of valid Region */
X		} /* End of two roots */
X	} /* End of loop through metaballs */
X
X	/*
X	 * If there are no Ray/Metaball intersections there will 
X	 * not be a Ray/Blob intersection. Exit now.
X	 */
X	if (inum == 0)
X	{
X		return FALSE;
X	}
X
X	/* 
X	 * Sort Intersection list. No sense using qsort if there's only
X	 * two intersections.
X	 *
X	 * Note: we actually aren't sorting the Intersection structures, but
X	 * an array of pointers to the Intersection structures. 
X	 * This is faster than sorting the Intersection structures themselves.
X	 */
X	if (inum==2)
X	{
X		MetaInt *t;
X		if (ilist[0].bound > ilist[1].bound)
X		{
X			t = iarr[0];
X			iarr[0] = iarr[1];
X			iarr[1] = t;
X		}
X	}
X	else qsort((voidstar)(iarr), (unsigned)inum, sizeof(MetaInt *),
X			MetaCompare);
X
X/*
X* Finding the Ray/Blob Intersection
X* 
X* The non-zero part of the density function for each Metaball is
X* 
X*   di(Ri) = c4i * Ri^4  +  c2i * Ri^2  +  c0i 
X*
X* To find find the Ray/Blob intersection for one metaball
X* substitute for distance   
X*
X*     Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2
X* 
X* and then substitute the Ray equations:
X*   
X*     x  = x0 + x1 * t
X*     y  = y0 + y1 * t
X*     z  = z0 + z1 * t
X*
X* to get a big mess :^). Actually, it's a Quartic in t and it's fully 
X* listed farther down. Here's a short version:
X*
X*   c[4] * t^4  +  c[3] * t^3  +  c[2] * t^2  +  c[1] * t  +  c[0]  =  T
X*
X* Don't forget that the Blob is defined by the density being equal to 
X* the threshold T.
X* To calculate the intersection of a Ray and two or more Metaballs, 
X* the coefficients are calculated for each Metaball and then added 
X* together. We can do this since we're working with polynomials.
X* The points of intersection are the roots of the resultant equation.
X*
X* The algorithm loops through the intersection list, calculating
X* the coefficients if an intersection is the start of a Region and 
X* adding them to all intersections in that region.
X* When it detects a valid interval, it calculates the 
X* roots from the starting intersection's coefficients and if any of 
X* the roots are in the interval, the smallest one is returned.
X*
X*/
X
X	{
X		register Float *tmpc;
X		MetaInt *strt,*tmp;
X		register int istrt,itmp;
X		register int num,exitflag,inside;
X
X	/*
X	 * Start the algorithm outside the first interval
X	 */
X		inside = 0;
X		istrt = 0; 
X		strt = iarr[istrt];
X		if (strt->type != R_START)
X			RLerror(RL_WARN,"MetaInt sanity check FAILED!\n");
X
X		/*
X		 * Loop through intersection. If a root is found the code
X		 * will return at that point.
X		 */
X		while( istrt < inum )
X		{
X			/* 
X			 * Check for multiple intersections  at the same point.
X			 * This is also where the coefficients are calculated
X			 * and spread throughout that Metaball's sphere of
X			 * influence.
X			 */
X			do
X			{
X				/* find out which metaball */
X				i = strt->pnt;
X				/* only at starting boundaries do this */
X				if (strt->type == R_START)
X				{
X					register MetaVector *ml;
X					register Float a1,a0;
X					register Float xd,yd,zd;
X					register Float c4,c2,c0;
X
X					/* we're inside */
X					inside++;
X
X	/*  As promised, the full equations
X	 *
X	 *   c[4] = c4*a2*a2; 
X	 *   c[3] = 4.0*c4*a1*a2;
X	 *   c[2] = 4.0*c4*a1*a1 + 2.0*c4*a2*a0 + c2*a2;
X	 *   c[1] = 4.0*c4*a1*a0 + 2.0*c2*a1;
X	 *   c[0] = c4*a0*a0 + c2*a0 + c0;
X	 *
X	 * where
X	 *        a2 = (x1*x1 + y1*y1 + z1*z1) = 1.0 because the ray 
X	 *                                           is normalized
X	 *        a1 = (xd*x1 + yd*y1 + zd*z1)
X	 *        a0 = (xd*xd + yd*yd + zd*zd)
X	 *        xd = (x0 - xi)
X	 *        yd = (y0 - yi)
X	 *        zd = (z0 - zi)
X	 *        (xi,yi,zi) is center of Metaball
X	 *        (x0,y0,z0) is Ray origin
X	 *        (x1,y1,z1) is normalized Ray direction
X	 *        c4,c2,c0   are the coefficients for the
X	 *                       Metaball's density function
X	 *
X	 */
X
X					/*
X					 * Some compilers would recalculate
X					 * this each time its used.
X					 */
X					ml = &(blob->list[i]);
X
X					xd = ray->pos.x - ml->x;
X					yd = ray->pos.y - ml->y;
X					zd = ray->pos.z - ml->z;
X					a1 = (xd * ray->dir.x + yd * ray->dir.y
X						+ zd * ray->dir.z);
X					a0 = (xd * xd + yd * yd + zd * zd);
X
X					c0 = ml->c0;
X					c2 = ml->c2;
X					c4 = ml->c4;
X
X					c[4] = c4;
X					c[3] = 4.0*c4*a1;
X					c[2] = 2.0*c4*(2.0*a1*a1 + a0) + c2;
X					c[1] = 2.0*a1*(2.0*c4*a0 + c2);
X					c[0] = c4*a0*a0 + c2*a0 + c0;
X
X		/* 
X		 * Since we've gone through the trouble of calculating the 
X		 * coefficients, put them where they'll be used.
X		 * Starting at the current intersection(It's also the start of
X		 * a region), add the coefficients to each intersection until
X		 * we reach the end of this region.
X		 */
X					itmp = istrt; 
X					tmp = strt;
X					while( (tmp->pnt != i) ||
X					       (tmp->type != R_END) )
X					{
X						tmpc = tmp->c;
X						for(j=0;j<5;j++)
X							*tmpc++ += c[j];
X						itmp++; 
X						tmp = iarr[itmp];
X					}
X
X				} /* End of start of a Region */ 
X				/* 
X				 * Since the intersection wasn't the start
X				 * of a region, it must the end of one.
X				 */
X				else inside--;
X
X		/*
X		 * If we are inside a region(or regions) and the next
X		 * intersection is not at the same place, then we have
X		 * the start of an interval. Set the exitflag.
X		 */
X				if ((inside > 0)
X				    && (strt->bound != iarr[istrt+1]->bound) )
X					exitflag = 1;
X				else 
X				/* move to next intersection */
X				{
X					istrt++; 
X					strt = iarr[istrt];
X					exitflag = 0;
X				}
X		/* 
X		 * Check to see if we're at the end. If so then exit with
X		 * no intersection found.
X		 */
X				if (istrt >= inum)
X				{
X					return FALSE;
X				}
X			} while(!exitflag);
X				/* End of Search for valid interval */
X
X		/*
X		 * Find Roots along this interval
X		 */
X
X			/* get coefficients from Intersection structure */
X			tmpc = strt->c;
X			for(j=0;j<5;j++) c[j] = *tmpc++;
X
X			/* Don't forget to put in threshold */
X			c[0] -= blob->T;
X
X			/* Use Graphic Gem's Quartic Root Finder */
X			num = SolveQuartic(c,s);
X
X		/*
X		 * If there are roots, search for roots within the interval.
X		 */
X			dist = 0.0;
X			if (num>0)
X			{
X				for(j=0;j<num;j++)
X				{
X		/* 
X		 * Not sure if EPSILON is truly needed, but it might cause
X		 * small cracks between intervals in some cases. In any case 
X		 * I don't believe it hurts.
X		 */
X					if ((s[j] >= (strt->bound - EPSILON))
X					    && (s[j] <= (iarr[istrt+1]->bound +
X					        EPSILON) ) )
X					{
X						if (dist == 0.0)
X						/* first valid root */
X							dist = s[j];
X						else if (s[j] < dist)
X	 					/* else only if smaller */
X							dist = s[j];
X					}
X				}
X				/*
X				 * Found a valid root 
X				 */
X				if (dist > mindist && dist < *maxdist)
X				{
X					*maxdist = dist;
X					BlobHits++;
X					return TRUE;
X					/* Yeah! Return valid root */
X				}
X			}
X
X		/* 
X		 * Move to next intersection
X		 */
X			istrt++;
X			strt = iarr[istrt];
X
X		} /* End of loop through Intersection List */
X	} /* End of Solving for Ray/Blob Intersection */
X
X	/* 
X	 * return negative
X	 */
X	return FALSE;
X}
X
X
X/***********************************************
X * Find the Normal of a Blob at a given point
X *
X * The normal of a surface at a point (x0,y0,z0) is the gradient of that
X * surface at that point. The gradient is the vector 
X *
X *            Fx(x0,y0,z0) , Fy(x0,y0,z0) , Fz(x0,y0,z0)
X *
X * where Fx(),Fy(),Fz() are the partial dervitives of F() with respect
X * to x,y and z respectively. Since the surface of a Blob is made up
X * of the sum of one or more polynomials, the normal of a Blob at a point
X * is the sum of the gradients of the individual polynomials at that point.
X * The individual polynomials in this case are di(Ri) i = 0,...,num .
X *
X * To find the gradient of the contributing polynomials
X * take di(Ri) and substitute U = Ri^2;
X *
X *      di(U)    = c4i * U^2  +  c2i * U  + c0i
X *
X *      dix(U)   = 2 * c4i * Ux * U  +  c2i * Ux
X *
X *        U      = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 
X *
X *        Ux     = 2 * (x-xi) 
X *
X * Rearranging we get
X *  
X *    dix(x0,y0,z0) = 4 * c4i * xdi * disti + 2 * c2i * xdi
X *    diy(x0,y0,z0) = 4 * c4i * ydi * disti + 2 * c2i * ydi
X *    diz(x0,y0,z0) = 4 * c4i * zdi * disti + 2 * c2i * zdi
X * 
X * where
X *         xdi       =   x0-xi
X *         ydi       =   y0-yi
X *         zdi       =   y0-yi
X *        disti      =   xdi*xdi + ydi*ydi + zdi*zdi   
X *       (xi,yi,zi)  is  the center of Metaball i
X *       c4i,c2i,c0i are the coefficients of Metaball i's density function
X */
Xint
XBlobNormal(blob, pos, nrm, gnrm)
XBlob *blob;
XVector *pos, *nrm, *gnrm;
X{
X	register int i;
X
X	/*  
X	 * Initialize normals to zero 
X	 */
X	nrm->x = nrm->y = nrm->z = 0.0;
X	/*
X	 * Loop through Metaballs. If the point is within a Metaball's
X	 * Sphere of influence, calculate the gradient and add it to the
X	 * normals
X	 */
X	for(i=0;i < blob->num; i++)
X	{
X		register MetaVector *sl;
X		register Float dist,xd,yd,zd;
X
X		sl = &(blob->list[i]);
X		xd = pos->x - sl->x;
X		yd = pos->y - sl->y;
X		zd = pos->z - sl->z;
X
X		dist  = xd*xd + yd*yd + zd*zd;
X		if (dist <= sl->rs )
X		{
X			register Float temp;
X
X			/* temp is negative so normal points out of blob */
X			temp = -2.0 * (2.0 * sl->c4 * dist  +  sl->c2);
X			nrm->x += xd * temp;
X			nrm->y += yd * temp;
X			nrm->z += zd * temp;
X
X		}
X	}
X	(void)VecNormalize(nrm);
X	*gnrm = *nrm;
X	return FALSE;
X}
X
X
X/*****************************************************************************
X * Calculate the extent of the Blob
X */
Xvoid
XBlobBounds(blob, bounds)
XBlob *blob;
XFloat bounds[2][3];
X{
X	Float r;
X	Float minx,miny,minz,maxx,maxy,maxz;
X	int i;
X
X	/*
X	 * Loop through Metaballs to find the minimun and maximum in each
X	 * direction.
X	 */
X	for( i=0;  i< blob->num;  i++)
X	{
X		register Float d;
X
X		r = sqrt(blob->list[i].rs);
X		if (i==0)
X		{
X			minx = blob->list[i].x - r;
X			miny = blob->list[i].y - r;
X			minz = blob->list[i].z - r;
X			maxx = blob->list[i].x + r;
X			maxy = blob->list[i].y + r;
X			maxz = blob->list[i].z + r;
X		}
X		else
X		{
X			d = blob->list[i].x - r; 
X			if (d<minx) minx = d;
X			d = blob->list[i].y - r; 
X			if (d<miny) miny = d;
X			d = blob->list[i].z - r; 
X			if (d<minz) minz = d;
X			d = blob->list[i].x + r; 
X			if (d>maxx) maxx = d;
X			d = blob->list[i].y + r; 
X			if (d>maxy) maxy = d;
X			d = blob->list[i].z + r; 
X			if (d>maxz) maxz = d;
X		}
X	}
X	bounds[LOW][X]  = minx;
X	bounds[HIGH][X] = maxx;
X	bounds[LOW][Y]  = miny;
X	bounds[HIGH][Y] = maxy;
X	bounds[LOW][Z]  = minz;
X	bounds[HIGH][Z] = maxz;
X}
X
Xchar *
XBlobName()
X{
X	return blobName;
X}
X
Xvoid
XBlobStats(tests, hits)
Xunsigned long *tests, *hits;
X{
X	*tests = BlobTests;
X	*hits = BlobHits;
X}
X
Xvoid
XBlobMethodRegister(meth)
XUserMethodType meth;
X{
X	if (iBlobMethods)
X		iBlobMethods->user = meth;
X}
END_OF_FILE
if test 18207 -ne `wc -c <'libray/libobj/blob.c'`; then
    echo shar: \"'libray/libobj/blob.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/blob.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 17 \(of 19\).
cp /dev/null ark17isdone
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