#! /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 14 (of 19)."
# Contents:  Doc/Guide/texture.tex libray/libobj/grid.c
#   libray/libobj/triangle.c rayshade/raytrace.c
# Wrapped by kolb@woody on Wed Jul 17 17:56:53 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'Doc/Guide/texture.tex' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Doc/Guide/texture.tex'\"
else
echo shar: Extracting \"'Doc/Guide/texture.tex'\" \(13009 characters\)
sed "s/^X//" >'Doc/Guide/texture.tex' <<'END_OF_FILE'
X\chapter {Texturing}
X
XTextures are used to modify the appearance of an object through the
Xuse of procedural functions.  A texture may modify any surface characteristic,
Xsuch as diffuse color, reflectivity, or transparency, or it may 
Xmodify the surface normal (``bump mapping'') in order to give the
Xappearance of a rough surface.
X
XAny number of textures may be associated with an object.  If more than
Xone texture is specified, they are applied in the order given.  This allows
Xone to compose texturing functions and create, for example
Xa tiled marble ground plane using the {\em checker} and {\em marble}
Xtextures.
X
XTextures are associated with objects by following the object specification
Xby a number of lines of the form:
X
X\begin{center}
X{\tt texture} {\em name} $<${\em Texturing Arguments}$>$ [{\em Transformations}]
X\end{center}
X
XTransformations may be applied to the texture in order to, for example,
Xshrink or grow feature size, change the orientation of features, and
Xchange the position of features.
X
XSeveral of the texturing functions take the name of a colormap as an
Xargument.  A colormap is 256-line ASCII file, with each line containing
Xthree space-separated values ranging from 0 to 255.  Each line gives
Xthe red, green, and blue values for a single entry in the colormap.
X
X\section {Texturing Functions}
X
X\begin{defkey}{blotch}{{\em BlendFactor surface}}
XProduces a mildly interesting blotchy-looking surface.
X{\em BlendFactor} is used to control the interpolation between
Xthe default surface characteristics and the characteristics of
Xthe given surface.  A value of 0 results in a roughly 50-50 mix
Xof the two surfaces.  Higher values result in a great portion of
Xthe default surface characteristics.
X\end{defkey}
X
X\begin{defkey}{bump}{{\em scale}}
XApply a random bump map.  The point of intersection is passed to
X{\em DNoise()}.
XThe returned normalized vector is weighted by {\em scale}
Xand the result is added to the normal vector at the point of intersection.
X\end{defkey}
X
X\begin{defkey}{checker}{$<${\em Surface}$>$}
XApplies a 3D checkerboard texture.  Every point that falls within an
X``even'' unit cube will be assigned the characteristics of the named surface
Xapplied to it, while points that fall within ``odd'' cubes will have
Xits usual surface characteristics.  Be wary of strange effects due
Xto roundoff error that occur when a planar checkered surface lies
Xin a plane of constant integral value (e.g., $z=0$) in texture space.
XIn such cases,
Xsimply translate the texture to ensure that the planar surface is not
Xcoincident with an integral plane in texture space
X(e.g., {\tt translate 0 0 0.1}).
X\end{defkey}
X
X\begin{defkey}{cloud}{{\em scale H $\lambda$ octaves cthresh lthresh tscale}}
X	This texture is a variant on Geoff Gardner's ellipsoid-texturing
X	algorithm.  It should be applied to unit spheres centered
X	at the origin.  These spheres may, of course, be transformed
X	at will to form the appropriately-shaped cloud or tree.
X
X	A sample of normalized {\em fBm} (see the {\em fbm} texture) is
X	generated
X	at the point of intersection.  This sample is used to
X	modulate the surface transparency.  The final transparency
X	if a function of the sample value, the
X	the proximity of the point of intersection to the edge of
X	the sphere (as seen from the ray origin), and three parameters
X	to control the overall ``density.''  The proximity of the point
X	to the sphere edge is determined by evaluating a {\em limb} function,
X	which varies from 0 on the limb to 1 at the center of the sphere.
X\[
Xtransp = 1. - \frac{fBm - cthresh - (lthresh - cthresh)(1 - limb)}{tscale}
X\]
X\end{defkey}
X
X\begin{defkey}{fbm}{{\em offset scale H $\lambda$ octaves thresh}
X[{\em colormap}]}
XGenerate a sample of discretized fractional Brownian motion (fBm) and
Xuses it to scale the diffuse and ambient component of an object's surface.
X{\em Scale} is used to scale the value
Xreturned by the fBm function.  {\em Offset} allows one to control the minimum
Xvalue of the fBm function.  {\em H} is the {\em Holder exponent}
Xused in the fBm function (a value of 0.5 works well).  $\lambda$ is
Xused to control {\em lacunarity}, and specifies the the frequency
Xdifference between successive samples of the fBm basis function (a
Xvalue of 2.0 will suffice).  {\em Octaves} specifies the number of
Xoctaves (samples) to take of the fBm basis function (in this case, Noise()).
XBetween five and seven octaves usually works well.  {\em Thresh} is used
Xto specify a lower bound onthe output of the fBm function.  Any
Xvalue lower than {\em thresh} is set to zero.
X
XIf a {\em colormap} is named, a 256-entry colormap is read from the named
Xfile, and the sample of fBm is scaled by 255 and is used as an index into
Xthe colormap.  The resulting colormap entry
Xis used to scale the ambient and diffuse components of the
Xobject's surface.
X\end{defkey}
X
X\begin{defkey}{fbmbump}{{\em offset scale H $\lambda$ octaves}}
XSimilar to the {\em fbm} texture.  Rather than modifying the color of
Xa surface, this texture acts as a bump map.
X\end{defkey}
X
X\begin{defkey}{gloss}{{\em glossiness}}
XGives reflective surfaces a glossy appearance. This texture perturbs
Xthe object's surface normal such that the normal ``samples'' a cone of
Xunit height with radius $1. - glossiness$.  A value of 1 results
Xin perfect mirror-like reflections, while a value of 0 results
Xin extremely fuzzy reflections.  For best results, jittered sampling
Xshould be used to render scenes that make use of this texture.
X\end{defkey}
X
X\begin{defkey}{marble}{[{\em colormap}]}
XGives a surface a marble-like appearance.  The texture is implemented as
Xroughly parallel alternating veins of marble, each of which is
Xseparated by 1/7 of a unit and runs perpendicular to the Z axis.
XIf a colormap is named, the surface's ambient and diffuse colors
Xwill be scaled using the RGB values in the colormap.  If no colormap is
Xgiven, the diffuse and ambient components are simply scaled by the
Xvalue of the marble function.  One may transform the texture to
Xcontrol the density and orientation of the marble veins.
X\end{defkey}
X
X\begin{defkey}{sky}{{\em scale H $\lambda$ octaves cthresh ltresh}}
X	Similar to the {\em fbm} texture.  Rather than modifying the
X	color of a surface, this texture modulates its transparency.
X	{\em cthresh} is the value of the {\em fBm} function above
X	which the surface is totally opaque.  Below {\em lthresh},
X	the surface is totally transparent.
X\end{defkey}
X
X\begin{defkey}{stripe}{$<${\em Surface}$>$ {\em size bump} $<$Mapping$>$}
X	Apply a ``raised'' stripe pattern to the surface.
X	The surface properties used to color the stripe are those
X	of the given surface.  The width of the stripe, as compared
X	to the unit interval, is given by {\em size}.  The magnitude
X	of {\em bump} controls the extent to which the bump appears
X	to be displaced from the rest of the surface.  If negative,
X	the stripe will appear to
X	sink into the surface; if negative, it will appear to stand
X	out of the surface.
X\end{defkey}
XMapping functions are described below.
X
X\begin{defkey}{wood}{}
XGives a surface a wood-like appearance.  The feature size of this texture
Xis approximately $0.01$ of a unit, making it often necessary to
Xscale the texture in order to achieve the desired appearance.
X\end{defkey}
X
X\section {Image Texturing}
X
X{\em Rayshade} also supports an {\tt image} texture.  This texture
Xallows you to use images to modify the characteristics of a surface.  You
Xcan use three-channel images to modify the any or all of
Xthe ambient, diffuse, and specular colors of a surface.
XIf you are using the Utah Raster Toolkit,
Xyou can also use single-channel images to modify surface reflectance,
Xtransparency, and the specular exponent.  You can also use a single-channel
Ximage to apply a bump map to a surface.
X
XIn all but the bump-mapping case,
Xa component is modified by multiplying the given value by the value
Xcomputed by the texturing function.  When using the Utah Raster Toolkit,
Xsurface characteristics are modified in proportion to the value of
Xthe {\em alpha} channel in the image.  If there is no {\em alpha} channel,
Xor you are not using the Utah Raster Toolkit, {\em alpha} is assumed to be
Xeverywhere
Xequal to $1$.
X
X\begin{defkey}{component}{$<${\em Component}$>$}
X	The named component will be modified.
X\end{defkey}
XPossible surface components are:
X{\tt ambient} (modify ambient color),
X{\tt diffuse} (modify diffuse color),
X{\tt specular} (modify specular color),
X{\tt specpow}, (modify specular exponent),
X{\tt reflect},	(modify reflectivity),
X{\tt transp} (modify transparency),
X{\tt bump}, (modify surface normal).
XThe {\tt specpow}, {\tt reflect}, {\tt transp}, and {\tt bump}
Xcomponents require the use of a single-channel image.
X
X\begin{defkey}{range}{{\em high low}}
X	Specify the range of values to which the values in the
X	image should be mapped.  An value of $1$ will
X	be mapped {\em high}, $0$ to {\em low}.  Intermediate
X	values will be linearly interpolated.
X\end{defkey}
X
X\begin{defkey}{smooth}{}
X	When given, pixel averaging will be performed in order
X	to smooth the sampled image.  If not specified, no averaging
X	will occur.
X\end{defkey}
X
X\begin{defkey}{textsurf}{$<${\em Surface Specification}$>$}
X	For use when modifying surface colors, this keyword specifies
X	that the given surface should be used as the base
X	to be modified when the {\em alpha} value in the image
X	is non-zero.  When {\em alpha} is zero, the object's
X	unmodified default surface characteristics are retained.
X\end{defkey}
XThe usual behavior is for the object's default surface properties to
Xbe used.
X
X\begin{defkey}{tile}{{\em un vn}}
X	Specify how the image should be tiled (repeated) along the
X	$u$ and $v$ axes.
X	If positive, the value of {\em un} gives the number of
X	times the image should be repeated along the $u$ axis, starting
X	from the origin of the texture, and positive {\em vn} gives the
X	number of times it
X	should be repeated along the $v$ axis.  If either value is zero,
X	the image is repeated infinitely along the appropriate axis.
X\end{defkey}
XTiling is usually only a concern when linear mapping is being used,
Xthough it may also be used if image textures are being scaled.  By default
X{\em un} and {\em vn} are both zero.
X 
X
XA mapping function may also be associated with an image texture.
X
X\section {Mapping Functions}
X
XMapping functions are used to apply two-dimensional textures to
Xsurfaces.  Each mapping functions defines a different method of transforming
Xa three dimensional point of intersection to a two dimensional $u-v$ pair
Xtermed texturing coordinates.
XTypically, the arguments to a mapping method define a center of
Xa projection and two non-parallel axes that define a local coordinate system.
X
XThe default mapping method is termed $u-v$ mapping or {\em inverse mapping}.
XNormally, there is a different inverse mapping method for each primitive type
X(see chapter 5).
XWhen inverse mapping is used, the point of intersection is passed to
Xthe $uv$ method for the primitive that was hit.
X
X\begin{defkey}{map}{{\tt uv}}
X	Use the $uv$ (inverse mapping) method associated with the
X	object that was intersected in order to map from 3D to determine
X	texturing coordinates.
X\end{defkey}
XThe inverse mapping method for each primitive is described in Chapter 5.
X
X\begin{defkey}{map}{{\tt linear} [\evec{origin} \evec{vaxis} \evec{uaxis}]}
X	Use a linear mapping method. The 2D texture is transformed
X	so that its $u$ axis is given by \evec{uaxis} and its $v$
X	axis by $vaxis$.  The texture is projected along the vector
X	defined by the cross product of the $u$ and $v$ axes, with
X	the (0,0) in texture space mapped to \evec{origin}.
X\end{defkey}
X
X\begin{defkey}{map}{{\tt cylindrical} [\evec{origin} \evec{vaxis} \evec{uaxis}]}
X	Use a cylindrical mapping method.  The point of intersection
X	is projected onto an imaginary cylinder, and the location
X	of the projected point is used to determine the texture coordinates.
X	If given, \evec{origin} and
X	\evec{vaxis} define the cylinder's axis, and \evec{uaxis} defines
X	where $u=0$ is located.
X\end{defkey}
XSee the description of the inverse mapping method for the 
Xcylinder in Chapter 5.  By default, the point of intersection is
Xprojected onto a cylinder that runs through the origin along the $z$
Xaxis, with \evec{uaxis} equal to the $x$ axis.
X
X\begin{defkey}{map}{{\tt spherical} [\evec{origin} \evec{vaxis} \evec{uaxis}]}
X	Use a spherical mapping method.  The intersection point is
X	projected onto an imaginary sphere, and the location of the
X	projected point	 is used to determine the texturing coordinates
X	in a manner identical to that used in the inverse mapping method
X	for the sphere primitive.
X	If given, the center of
X	the projection is \evec{origin}, \evec{vaxis} defines
X	the sphere axis, and the point where the
X	non-parallel \evec{uaxis} intersects the sphere
X	defines where $u=0$ is located.
X\end{defkey}
XBy default, a spherical mapping projects points towards the origin,
Xwith \evec{vaxis} defined to be the $z$ axis and
X\evec{uaxis} defined to be the $x$ axis.
END_OF_FILE
if test 13009 -ne `wc -c <'Doc/Guide/texture.tex'`; then
    echo shar: \"'Doc/Guide/texture.tex'\" unpacked with wrong size!
fi
# end of 'Doc/Guide/texture.tex'
fi
if test -f 'libray/libobj/grid.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'libray/libobj/grid.c'\"
else
echo shar: Extracting \"'libray/libobj/grid.c'\" \(11666 characters\)
sed "s/^X//" >'libray/libobj/grid.c' <<'END_OF_FILE'
X/*
X * grid.c
X *
X * Copyright (C) 1989, 1991, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: grid.c,v 4.0 91/07/17 14:38:02 kolb Exp Locker: kolb $
X *
X * $Log:	grid.c,v $
X * Revision 4.0  91/07/17  14:38:02  kolb
X * Initial version.
X * 
X */
X#include "geom.h"
X#include "grid.h"
X
Xstatic Methods *iGridMethods = NULL;
Xstatic char gridName[] = "grid";
X
Xstatic unsigned long raynumber = 1;		/* Current "ray number". */
X						/* (should be "grid number") */
Xstatic void engrid(), GridFreeVoxels();
Xstatic int pos2grid(), CheckVoxel();
X
XGrid *
XGridCreate(x, y, z)
Xint x, y, z;
X{
X	Grid *grid;
X
X	if (x < 1 || y < 1 || z < 1) {
X		RLerror(RL_WARN, "Invalid grid specification.\n");
X		return (Grid *)NULL;
X	}
X	grid = (Grid *)share_calloc(1, sizeof(Grid));
X	grid->xsize = x;
X	grid->ysize = y;
X	grid->zsize = z;
X	return grid;	
X}
X
Xchar *
XGridName()
X{
X	return gridName;
X}
X
X/*
X * Intersect ray with grid, returning distance from "pos" to
X * nearest intersection with an object in the grid.  Returns 0.
X * if no intersection.
X */
Xint
XGridIntersect(grid, ray, hitlist, mindist, maxdist)
XGrid *grid;
XRay *ray;
XHitList *hitlist;
XFloat mindist, *maxdist;
X{
X	GeomList *list;
X	Geom *obj;
X	int hit;
X	Float offset, tMaxX, tMaxY, tMaxZ;
X	Float tDeltaX, tDeltaY, tDeltaZ, *raybounds[2][3];
X	int stepX, stepY, stepZ, outX, outY, outZ, x, y, z;
X	Vector curpos, nXp, nYp, nZp, np, pDeltaX, pDeltaY, pDeltaZ;
X	unsigned long counter;
X
X	hit = FALSE;
X	/*
X	 * Check unbounded objects.
X	 */
X	for (obj = grid->unbounded ; obj; obj = obj->next) {
X		if (intersect(obj, ray, hitlist, mindist, maxdist))
X			hit = TRUE;
X	}
X
X	/*
X	 * If outside of the bounding box, check to see if we
X	 * hit it.
X	 */
X	VecAddScaled(ray->pos, mindist, ray->dir, &curpos);
X	if (OutOfBounds(&curpos, grid->bounds)) {
X		offset = *maxdist;
X		if (!BoundsIntersect(ray, grid->bounds, mindist, &offset))
X			/*
X			 * Ray never hit grid space.
X			 */
X			return hit;
X		/*
X		 * else
X		 *	The ray enters voxel space before it hits
X		 * 	an unbounded object.
X		 */
X		VecAddScaled(ray->pos, offset, ray->dir, &curpos);
X	} else
X		offset = mindist;
X
X	counter = raynumber++;
X
X	/*
X	 * tMaxX is the absolute distance from the ray origin we must move
X	 *		to get to the next voxel in the X
X	 *		direction.  It is incrementally updated
X	 * 		by DDA as we move from voxel-to-voxel.
X	 * tDeltaX is the relative amount along the ray we move to
X	 *		get to the next voxel in the X direction. Thus,
X	 *		when we decide to move in the X direction,
X	 * 		we increment tMaxX by tDeltaX.
X	 */
X	x = x2voxel(grid, curpos.x);
X	if (x == grid->xsize)
X		x--;
X	if (ray->dir.x < 0.) {
X		tMaxX = offset + (voxel2x(grid, x) - curpos.x) / ray->dir.x;
X		tDeltaX = grid->voxsize[X] / - ray->dir.x;
X		stepX = outX = -1;
X		raybounds[LOW][X] = &np.x;
X		raybounds[HIGH][X] = &curpos.x;
X	} else if (ray->dir.x > 0.) {
X		tMaxX = offset + (voxel2x(grid, x+1) - curpos.x) / ray->dir.x;
X		tDeltaX = grid->voxsize[X] / ray->dir.x;
X		stepX = 1;
X		outX = grid->xsize;
X		raybounds[LOW][X] = &curpos.x;
X		raybounds[HIGH][X] = &np.x;
X	} else {
X		tMaxX = FAR_AWAY;
X		raybounds[LOW][X] = &curpos.x;
X		raybounds[HIGH][X] = &np.x;
X		tDeltaX = 0.;
X	}
X
X	y = y2voxel(grid, curpos.y);
X	if (y == grid->ysize)
X		y--;
X	if (ray->dir.y < 0.) {
X		tMaxY = offset + (voxel2y(grid, y) - curpos.y) / ray->dir.y;
X		tDeltaY = grid->voxsize[Y] / - ray->dir.y;
X		stepY = outY = -1;
X		raybounds[LOW][Y] = &np.y;
X		raybounds[HIGH][Y] = &curpos.y;
X	} else if (ray->dir.y > 0.) {
X		tMaxY = offset + (voxel2y(grid, y+1) - curpos.y) / ray->dir.y;
X		tDeltaY = grid->voxsize[Y] / ray->dir.y;
X		stepY = 1;
X		outY = grid->ysize;
X		raybounds[LOW][Y] = &curpos.y;
X		raybounds[HIGH][Y] = &np.y;
X	} else {
X		tMaxY = FAR_AWAY;
X		raybounds[LOW][Y] = &curpos.y;
X		raybounds[HIGH][Y] = &np.y;
X		tDeltaY = 0.;
X	}
X
X	z = z2voxel(grid, curpos.z);
X	if (z == grid->zsize)
X		z--;
X	if (ray->dir.z < 0.) {
X		tMaxZ = offset + (voxel2z(grid, z) - curpos.z) / ray->dir.z;
X		tDeltaZ = grid->voxsize[Z] / - ray->dir.z;
X		stepZ = outZ = -1;
X		raybounds[LOW][Z] = &np.z;
X		raybounds[HIGH][Z] = &curpos.z;
X	} else if (ray->dir.z > 0.) {
X		tMaxZ = offset + (voxel2z(grid, z+1) - curpos.z) / ray->dir.z;
X		tDeltaZ = grid->voxsize[Z] / ray->dir.z;
X		stepZ = 1;
X		outZ = grid->zsize;
X		raybounds[LOW][Z] = &curpos.z;
X		raybounds[HIGH][Z] = &np.z;
X	} else {
X		tMaxZ = FAR_AWAY;
X		raybounds[LOW][Z] = &curpos.z;
X		raybounds[HIGH][Z] = &np.z;
X		tDeltaZ = 0.;
X	}
X
X	VecScale(tDeltaX, ray->dir, &pDeltaX);
X	VecScale(tDeltaY, ray->dir, &pDeltaY);
X	VecScale(tDeltaZ, ray->dir, &pDeltaZ);
X
X	VecAddScaled(ray->pos, tMaxX, ray->dir, &nXp);
X	VecAddScaled(ray->pos, tMaxY, ray->dir, &nYp);
X	VecAddScaled(ray->pos, tMaxZ, ray->dir, &nZp);
X
X	while (TRUE) {
X		list = grid->cells[x][y][z];
X		if (tMaxX < tMaxY && tMaxX < tMaxZ) {
X			if (list) {
X				np = nXp;
X			    	if (CheckVoxel(list,ray,raybounds,
X			    	    hitlist,counter,offset,maxdist))
X					hit = TRUE;
X			}
X			x += stepX;
X			if (*maxdist < tMaxX || x == outX)
X				break;
X			tMaxX += tDeltaX;
X			curpos = nXp;
X			nXp.x += pDeltaX.x;
X			nXp.y += pDeltaX.y;
X			nXp.z += pDeltaX.z;
X		} else if (tMaxZ < tMaxY) {
X			if (list) {
X				np = nZp;
X			    	if (CheckVoxel(list,ray, raybounds,
X			    	    hitlist,counter,offset,maxdist))
X					hit = TRUE;
X			}
X			z += stepZ;
X			if (*maxdist < tMaxZ || z == outZ)
X				break;
X			tMaxZ += tDeltaZ;
X			curpos = nZp;
X			nZp.x += pDeltaZ.x;
X			nZp.y += pDeltaZ.y;
X			nZp.z += pDeltaZ.z;
X		} else {
X			if (list) {
X				np = nYp;
X			    	if (CheckVoxel(list,ray,raybounds,
X			    	    hitlist,counter,offset,maxdist))
X					hit = TRUE;
X			}
X			y += stepY;
X			if (*maxdist < tMaxY || y == outY)
X				break;
X			tMaxY += tDeltaY;
X			curpos = nYp;
X			nYp.x += pDeltaY.x;
X			nYp.y += pDeltaY.y;
X			nYp.z += pDeltaY.z;
X		}
X	}
X	return hit;
X}
X
X/*
X * Intersect ray with objects in grid cell.  Note that there are a many ways
X * to speed up this routine, all of which uglify the code to a large extent.
X */
Xstatic int
XCheckVoxel(list,ray,raybounds,hitlist,counter,mindist,maxdist)
XGeomList *list;
XRay *ray;
XFloat *raybounds[2][3];
XHitList *hitlist;
Xunsigned long counter;
XFloat mindist, *maxdist;
X{
X	Geom *obj;
X	int hit;
X	Float lx, hx, ly, hy, lz, hz;
X
X	lx = *raybounds[LOW][X];
X	hx = *raybounds[HIGH][X];
X	ly = *raybounds[LOW][Y];
X	hy = *raybounds[HIGH][Y];
X	lz = *raybounds[LOW][Z];
X	hz = *raybounds[HIGH][Z];
X
X	hit = FALSE;
X
X	do {
X		obj = list->obj;
X		/*
X		 * If object's counter is greater than or equal to the
X		 * number associated with the current grid,
X		 * don't bother checking again.  In addition, if the
X		 * bounding box of the ray's extent in the voxel does
X		 * not intersect the bounding box of the object, don't bother.
X		 */
X#ifdef SHAREDMEM
X		if (*obj->counter < counter &&
X#else
X		if (obj->counter < counter &&
X#endif
X		    obj->bounds[LOW][X] <= hx  &&
X		    obj->bounds[HIGH][X] >= lx &&
X		    obj->bounds[LOW][Y] <= hy  &&
X		    obj->bounds[HIGH][Y] >= ly &&
X		    obj->bounds[LOW][Z] <= hz  &&
X		    obj->bounds[HIGH][Z] >= lz) {
X#ifdef SHAREDMEM
X			*obj->counter = counter;
X#else
X			obj->counter = counter;
X#endif
X			if (intersect(obj, ray, hitlist, mindist, maxdist))
X				hit = TRUE;
X		}
X	} while ((list = list->next) != (GeomList *)0);
X
X	return hit;
X}
X
Xint
XGridConvert(grid, objlist)
XGrid *grid;
XGeom *objlist;
X{
X	int num;
X
X	/*
X	 * Keep linked list of all bounded objects in grid; it may come
X	 * in handy.
X	 */
X	grid->objects = objlist;
X	for (num = 0; objlist; objlist = objlist->next)
X		num += objlist->prims;
X
X	return num;
X}
X
Xvoid
XGridBounds(grid, bounds)
XGrid *grid;
XFloat bounds[2][3];
X{
X	Geom *obj, *ltmp;
X	int x, y, i;
X
X	BoundsInit(bounds);
X	/*
X	 * For each object on the list,
X	 * compute its bounds...
X	 */
X	/*
X	 * Find bounding box of bounded objects and get list of
X	 * unbounded objects.
X	 */
X	grid->unbounded = GeomComputeAggregateBounds(&grid->objects,
X				grid->unbounded, grid->bounds);
X
X	BoundsCopy(grid->bounds, bounds);
X
X	grid->voxsize[X] = (grid->bounds[HIGH][X]-grid->bounds[LOW][X])/
X				grid->xsize;
X	grid->voxsize[Y] = (grid->bounds[HIGH][Y]-grid->bounds[LOW][Y])/
X				grid->ysize;
X	grid->voxsize[Z] = (grid->bounds[HIGH][Z]-grid->bounds[LOW][Z])/
X				grid->zsize;
X
X	if (grid->cells == (GeomList ****)NULL) {
X		/*
X	 	 * Allocate voxels.
X	 	 */
X		grid->cells = (GeomList ****)share_malloc(grid->xsize *
X					sizeof(GeomList ***));
X		for (x = 0; x < grid->xsize; x++) {
X			grid->cells[x] = (GeomList ***)share_malloc(grid->ysize *
X				sizeof(GeomList **));
X			for (y = 0; y < grid->ysize; y++)
X				grid->cells[x][y] = (GeomList **)share_calloc(
X					(unsigned)grid->zsize,sizeof(GeomList *));
X		}
X	} else {
X		/*
X		 * New frame...
X		 * Free up the objlists in each voxel.
X		 */
X		GridFreeVoxels(grid);
X	}
X
X	/*
X	 * objlist now holds a linked list of bounded objects.
X	 */
X	for (ltmp = grid->objects; ltmp != (Geom *)0; ltmp = ltmp->next)
X		engrid(ltmp, grid);
X}
X
Xstatic void
XGridFreeVoxels(grid)
XGrid *grid;
X{
X	int x, y, z;
X	GeomList *cell, *next;
X
X	for (x = 0; x < grid->xsize; x++) {
X		for (y = 0; y < grid->ysize; y++) {
X			for (z = 0; z < grid->zsize; z++) {
X				for (cell = grid->cells[x][y][z]; cell; cell = next) {
X					next = cell->next;
X					free((voidstar)cell);
X				}
X				grid->cells[x][y][z] = (GeomList *)NULL;
X			}
X		}
X	}
X}
X
XMethods *
XGridMethods()
X{
X	if (iGridMethods == (Methods *)NULL) {
X		iGridMethods = MethodsCreate();
X		iGridMethods->methods = GridMethods;
X		iGridMethods->create = (GeomCreateFunc *)GridCreate;
X		iGridMethods->intersect = GridIntersect;
X		iGridMethods->name = GridName;
X		iGridMethods->convert = GridConvert;
X		iGridMethods->bounds = GridBounds;
X		iGridMethods->checkbounds = FALSE;
X		iGridMethods->closed = TRUE;
X	}
X	return iGridMethods;
X}
X
X/*
X * Place an object in a grid.
X */
Xstatic void
Xengrid(obj, grid)
XGeom *obj;
XGrid *grid;
X{
X	int x, y, z, low[3], high[3];
X	GeomList *ltmp;
X
X	/*
X	 * This routine should *never* be passed an unbounded object, but...
X	 */
X	if (!pos2grid(grid, obj->bounds[LOW], low) ||
X	    !pos2grid(grid, obj->bounds[HIGH], high) ||
X	    obj->bounds[LOW][X] > obj->bounds[HIGH][X]) {
X		/*
X		 * Geom is partially on wholly outside of
X		 * grid -- this should never happen, but just
X		 * in case...
X		 */
X		RLerror(RL_ABORT, "Engrid got an unbounded object?!\n");
X		return;
X	    }
X
X	/*
X	 * For each voxel that intersects the object's bounding
X	 * box, add pointer to this object to voxel's linked list.
X 	 */
X	for (x = low[X]; x <= high[X]; x++) {
X		for (y = low[Y]; y <= high[Y]; y++) {
X			for (z = low[Z]; z <= high[Z]; z++) {
X				ltmp = (GeomList *)share_malloc(sizeof(GeomList));
X				ltmp->obj = obj;
X				ltmp->next = grid->cells[x][y][z];
X				grid->cells[x][y][z] = ltmp;
X			}
X		}
X	}
X}
X
X/*
X * Convert 3D point to index into grid's voxels.
X */
Xstatic int
Xpos2grid(grid, pos, index)
XGrid *grid;
XFloat pos[3];
Xint index[3];
X{
X	index[X] = (int)(x2voxel(grid, pos[0]));
X	index[Y] = (int)(y2voxel(grid, pos[1]));
X	index[Z] = (int)(z2voxel(grid, pos[2]));
X
X	if (index[X] == grid->xsize)
X		index[X]--;
X	if (index[Y] == grid->ysize)
X		index[Y]--;
X	if (index[Z] == grid->zsize)
X		index[Z]--;
X
X	if (index[X] < 0 || index[X] >= grid->xsize ||
X	    index[Y] < 0 || index[Y] >= grid->ysize ||
X	    index[Z] < 0 || index[Z] >= grid->zsize)
X		return FALSE;
X	return TRUE;
X}
X
Xvoid
XGridMethodRegister(meth)
XUserMethodType meth;
X{
X	if (iGridMethods)
X		iGridMethods->user = meth;
X}
END_OF_FILE
if test 11666 -ne `wc -c <'libray/libobj/grid.c'`; then
    echo shar: \"'libray/libobj/grid.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/grid.c'
fi
if test -f 'libray/libobj/triangle.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'libray/libobj/triangle.c'\"
else
echo shar: Extracting \"'libray/libobj/triangle.c'\" \(11924 characters\)
sed "s/^X//" >'libray/libobj/triangle.c' <<'END_OF_FILE'
X/*
X * triangle.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: triangle.c,v 4.0 91/07/17 14:39:38 kolb Exp Locker: kolb $
X *
X * $Log:	triangle.c,v $
X * Revision 4.0  91/07/17  14:39:38  kolb
X * Initial version.
X * 
X */
X#include "geom.h"
X#include "triangle.h"
X
Xstatic Methods *iTriangleMethods = NULL;
Xstatic char triName[] = "triangle";
X
Xunsigned long TriTests, TriHits;
X
Xstatic void TriangleSetdPdUV();
X
X/*
X * Create and return reference to a triangle.
X */
XTriangle *
XTriangleCreate(type, p1, p2, p3, n1, n2, n3, u1, u2, u3, flipflag)
Xint	type;
XVector	*p1, *p2, *p3, *n1, *n2, *n3;
XVec2d	*u1, *u2, *u3;
Xint flipflag;
X{
X	Triangle *triangle;
X	Vector ptmp, anorm;
X	Float d;
X
X	/*
X	 * Allocate new triangle and primitive to point to it.
X	 */
X	triangle = (Triangle *)share_malloc(sizeof(Triangle));
X	triangle->type = type;	/* so inttri can tell the difference */
X
X	VecSub(*p2, *p1, &triangle->e[0]);
X	VecSub(*p3, *p2, &triangle->e[1]);
X	VecSub(*p1, *p3, &triangle->e[2]);
X
X	/* Find plane normal. */
X	VecCross(&triangle->e[0], &triangle->e[1], &ptmp);
X	triangle->nrm = ptmp;
X	if (VecNormalize(&triangle->nrm) == 0.) {
X		RLerror(RL_ADVISE, "Degenerate triangle.\n");
X		return (Triangle *)NULL;
X	}
X
X	if (flipflag)
X		VecScale(-1, triangle->nrm, &triangle->nrm);
X
X	triangle->d = dotp(&triangle->nrm, p1);
X
X	triangle->p[0] = *p1;
X	triangle->p[1] = *p2;
X	triangle->p[2] = *p3;
X
X	if (type == PHONGTRI) {
X		if (VecNormalize(n1) == 0. || VecNormalize(n2) == 0. ||
X		    VecNormalize(n3) == 0.) {
X			RLerror(RL_WARN, "Degenerate vertex normal.\n");
X			return (Triangle *)NULL;
X		}
X		triangle->vnorm = (Vector *)Malloc(3 * sizeof(Vector));
X		triangle->vnorm[0] = *n1;
X		triangle->vnorm[1] = *n2;
X		triangle->vnorm[2] = *n3;
X		if (flipflag) {
X			/* Flip the vertex normals */
X			VecScale(-1, triangle->vnorm[0],
X				&triangle->vnorm[0]);
X			VecScale(-1, triangle->vnorm[1],
X				&triangle->vnorm[1]);
X			VecScale(-1, triangle->vnorm[2],
X				&triangle->vnorm[2]);
X		} else if (dotp(&triangle->vnorm[0], &triangle->nrm) < 0.) {
X			/*
X			 * Reverse direction of surface normal on Phong
X			 * triangle if the surface normal points "away"
X			 * from the first vertex normal.
X			 * Note that this means that we trust the vertex
X			 * normals rather than trust that the user gave the
X			 * vertices in the correct order.
X			 */
X			RLerror(RL_ADVISE, "Inconsistant triangle normals.\n");
X			VecScale(-1., triangle->nrm, &triangle->nrm);
X			VecScale(-1., ptmp, &ptmp);
X			triangle->d = -triangle->d;
X			VecScale(-1., triangle->e[0], &triangle->e[0]);
X			VecScale(-1., triangle->e[1], &triangle->e[1]);
X			VecScale(-1., triangle->e[2], &triangle->e[2]);
X		}
X	}
X
X	/*
X	 * If UV coordinates are given for the vertices, allocate and
X	 * store them.
X	 */
X	if (u1 && u2 && u3) {
X		triangle->uv = (Vec2d *)Malloc(3 * sizeof(Vec2d));
X		triangle->uv[0] = *u1;
X		triangle->uv[1] = *u2;
X		triangle->uv[2] = *u3;
X		/* Calculate dpdu and dpdv vectors */
X		triangle->dpdu = (Vector *)Malloc(sizeof(Vector));
X		triangle->dpdv = (Vector *)Malloc(sizeof(Vector));
X		TriangleSetdPdUV(triangle->p, triangle->uv,
X				 triangle->dpdu, triangle->dpdv);
X	} else {
X		triangle->uv = (Vec2d *)NULL;
X	}
X
X	/*
X	 * Find "dominant" part of normal vector.
X	 */
X	anorm.x = fabs(ptmp.x);
X	anorm.y = fabs(ptmp.y);
X	anorm.z = fabs(ptmp.z);
X
X	/*
X	 * Scale edges by dominant part of normal.  This makes intersection
X	 * testing a bit faster.
X	 */ 
X	if (anorm.x > anorm.y && anorm.x > anorm.z) {
X		triangle->index = XNORMAL;
X		d = 1. / ptmp.x;
X	} else if (anorm.y > anorm.z) {
X		triangle->index = YNORMAL;
X		d = 1. / ptmp.y;
X	} else {
X		triangle->index = ZNORMAL;
X		d = 1. /ptmp.z;
X	}
X
X	VecScale(d, triangle->e[0], &triangle->e[0]);
X	VecScale(d, triangle->e[1], &triangle->e[1]);
X	VecScale(d, triangle->e[2], &triangle->e[2]);
X
X	return triangle;
X}
X
XMethods *
XTriangleMethods()
X{
X	if (iTriangleMethods == (Methods *)NULL) {
X		iTriangleMethods = MethodsCreate();
X		iTriangleMethods->create = (GeomCreateFunc *)TriangleCreate;
X		iTriangleMethods->methods = TriangleMethods;
X		iTriangleMethods->name = TriangleName;
X		iTriangleMethods->intersect = TriangleIntersect;
X		iTriangleMethods->normal = TriangleNormal;
X		iTriangleMethods->uv = TriangleUV;
X		iTriangleMethods->bounds = TriangleBounds;
X		iTriangleMethods->stats = TriangleStats;
X		iTriangleMethods->checkbounds = TRUE;
X		iTriangleMethods->closed = FALSE;
X	}
X	return iTriangleMethods;
X}
X
X/*
X * Intersect ray with triangle.  This is an optimized version of the
X * intersection routine from Snyder and Barr's '87 SIGGRAPH paper.
X */
Xint
XTriangleIntersect(tri, ray, mindist, maxdist)
XTriangle *tri;
XRay *ray;
XFloat mindist, *maxdist;
X{
X	Float qi1, qi2, s, k, b0, b1, b2;
X	Vector pos, dir;
X
X	TriTests++;
X	pos = ray->pos;
X	dir = ray->dir;
X	/*
X	 * Plane intersection.
X	 */
X	k = dotp(&tri->nrm, &dir);
X	if (fabs(k) < EPSILON)
X		return FALSE;
X	s = (tri->d - dotp(&tri->nrm, &pos)) / k;
X	if (s < mindist || s > *maxdist)
X		return FALSE;
X	
X	if (tri->index == XNORMAL) {
X		qi1 = pos.y + s * dir.y;
X		qi2 = pos.z + s * dir.z;
X		b0 = tri->e[1].y * (qi2 - tri->p[1].z) -
X				tri->e[1].z * (qi1 - tri->p[1].y);
X		if (b0 < 0. || b0 > 1.)
X			return FALSE;
X		b1 = tri->e[2].y * (qi2 - tri->p[2].z) -
X				tri->e[2].z * (qi1 - tri->p[2].y);
X		if (b1 < 0. || b1 > 1.)
X			return FALSE;
X		b2 = tri->e[0].y * (qi2 - tri->p[0].z) -
X				tri->e[0].z * (qi1 - tri->p[0].y);
X		if (b2 < 0. || b2 > 1.)
X			return FALSE;
X	} else if (tri->index == YNORMAL) {
X		qi1 = pos.x + s * dir.x;
X		qi2 = pos.z + s * dir.z;
X		b0 = tri->e[1].z * (qi1 - tri->p[1].x) -
X			tri->e[1].x * (qi2 - tri->p[1].z);
X		if (b0 < 0. || b0 > 1.)
X			return FALSE;
X		b1 = tri->e[2].z * (qi1 - tri->p[2].x) -
X			tri->e[2].x * (qi2 - tri->p[2].z);
X		if (b1 < 0. || b1 > 1.)
X			return FALSE;
X		b2 = tri->e[0].z * (qi1 - tri->p[0].x) -
X			tri->e[0].x * (qi2 - tri->p[0].z);
X		if (b2 < 0. || b2 > 1.)
X			return FALSE;
X	} else {
X		qi1 = pos.x + s * dir.x;
X		qi2 = pos.y + s * dir.y;
X		b0 = tri->e[1].x * (qi2 - tri->p[1].y) -
X			tri->e[1].y * (qi1 - tri->p[1].x);
X		if (b0 < 0. || b0 > 1.)
X			return FALSE;
X		b1 = tri->e[2].x * (qi2 - tri->p[2].y) -
X				tri->e[2].y * (qi1 - tri->p[2].x);
X		if (b1 < 0. || b1 > 1.)
X			return FALSE;
X		b2 = tri->e[0].x * (qi2 - tri->p[0].y) -
X				tri->e[0].y * (qi1 - tri->p[0].x);
X		if (b2 < 0. || b2 > 1.)
X			return FALSE;
X	}
X
X	tri->b[0] = b0;
X	tri->b[1] = b1;
X	tri->b[2] = b2;
X
X	TriHits++;
X	*maxdist = s;
X	return TRUE;
X}
X
Xint
XTriangleNormal(tri, pos, nrm, gnrm)
XTriangle *tri;
XVector *pos, *nrm, *gnrm;
X{
X	*gnrm = tri->nrm;
X
X	if (tri->type == FLATTRI) {
X		*nrm = tri->nrm;
X		return FALSE;
X	}
X
X	/*
X	 * Interpolate normals of Phong-shaded triangles.
X	 */
X	nrm->x = tri->b[0]*tri->vnorm[0].x+tri->b[1]*tri->vnorm[1].x+
X		tri->b[2]*tri->vnorm[2].x;
X	nrm->y = tri->b[0]*tri->vnorm[0].y+tri->b[1]*tri->vnorm[1].y+
X		tri->b[2]*tri->vnorm[2].y;
X	nrm->z = tri->b[0]*tri->vnorm[0].z+tri->b[1]*tri->vnorm[1].z+
X		tri->b[2]*tri->vnorm[2].z;
X	(void)VecNormalize(nrm);
X	return TRUE;
X}
X
X/*ARGSUSED*/
Xvoid
XTriangleUV(tri, pos, norm, uv, dpdu, dpdv)
XTriangle *tri;
XVector *pos, *norm, *dpdu, *dpdv;
XVec2d *uv;
X{
X	Float d;
X
X	/*
X	 * Normalize barycentric coordinates.
X	 */
X	d = tri->b[0]+tri->b[1]+tri->b[2];
X
X	tri->b[0] /= d;
X	tri->b[1] /= d; 
X	tri->b[2] /= d;
X
X	if (dpdu) {
X		if (tri->uv == (Vec2d *)NULL) {
X			*dpdu = tri->e[0];
X			(void)VecNormalize(dpdu);
X			VecSub(tri->p[0], *pos, dpdv);
X			(void)VecNormalize(dpdv);
X		} else {
X			*dpdu = *tri->dpdu;
X			*dpdv = *tri->dpdv;
X		}
X	}
X
X	if (tri->uv == (Vec2d *)NULL) {
X		uv->v = tri->b[2];
X		if (equal(uv->v, 1.))
X			uv->u = 0.;
X		else
X			uv->u = tri->b[1] / (tri->b[0] + tri->b[1]);
X	} else {
X		/*
X		 * Compute UV by taking weighted sum of UV coordinates.
X		 */
X		uv->u = tri->b[0]*tri->uv[0].u + tri->b[1]*tri->uv[1].u +
X			tri->b[2]*tri->uv[2].u;
X		uv->v = tri->b[0]*tri->uv[0].v + tri->b[1]*tri->uv[1].v +
X			tri->b[2]*tri->uv[2].v;
X	}
X}
X
Xvoid
XTriangleBounds(tri, bounds)
XTriangle *tri;
XFloat bounds[2][3];
X{
X	bounds[LOW][X] = bounds[HIGH][X] = tri->p[0].x;
X	bounds[LOW][Y] = bounds[HIGH][Y] = tri->p[0].y;
X	bounds[LOW][Z] = bounds[HIGH][Z] = tri->p[0].z;
X
X	if (tri->p[1].x < bounds[LOW][X]) bounds[LOW][X] = tri->p[1].x;
X	if (tri->p[1].x > bounds[HIGH][X]) bounds[HIGH][X] = tri->p[1].x;
X	if (tri->p[2].x < bounds[LOW][X]) bounds[LOW][X] = tri->p[2].x;
X	if (tri->p[2].x > bounds[HIGH][X]) bounds[HIGH][X] = tri->p[2].x;
X
X	if (tri->p[1].y < bounds[LOW][Y]) bounds[LOW][Y] = tri->p[1].y;
X	if (tri->p[1].y > bounds[HIGH][Y]) bounds[HIGH][Y] = tri->p[1].y;
X	if (tri->p[2].y < bounds[LOW][Y]) bounds[LOW][Y] = tri->p[2].y;
X	if (tri->p[2].y > bounds[HIGH][Y]) bounds[HIGH][Y] = tri->p[2].y;
X
X	if (tri->p[1].z < bounds[LOW][Z]) bounds[LOW][Z] = tri->p[1].z;
X	if (tri->p[1].z > bounds[HIGH][Z]) bounds[HIGH][Z] = tri->p[1].z;
X	if (tri->p[2].z < bounds[LOW][Z]) bounds[LOW][Z] = tri->p[2].z;
X	if (tri->p[2].z > bounds[HIGH][Z]) bounds[HIGH][Z] = tri->p[2].z;
X}
X
Xchar *
XTriangleName()
X{
X	return triName;
X}
X
Xvoid
XTriangleStats(tests, hits)
Xunsigned long *tests, *hits;
X{
X	*tests = TriTests;
X	*hits = TriHits;
X}
X
X/*
X * Given three vertices of a triangle and the uv coordinates associated
X * with each, compute directions of u and v axes.
X */
Xstatic void
XTriangleSetdPdUV(p, t, dpdu, dpdv)
XVector p[3];			/* Triangle vertices */
XVec2d t[3];			/* uv coordinates for each vertex */
XVector *dpdu, *dpdv;		/* u and v axes (return values) */
X{
X	Float scale;
X	int hi, mid, lo;
X	Vector base;
X
X	/* sort u coordinates */
X	if (t[2].u > t[1].u) {
X		if (t[1].u > t[0].u) {
X			hi = 2; mid = 1; lo = 0;
X		} else if (t[2].u > t[0].u) {
X			hi = 2; mid = 0; lo = 1;
X		} else {
X			hi = 0; mid = 2; lo = 1;
X		}
X	} else {
X		if (t[2].u > t[0].u) {
X			hi = 1; mid = 2; lo = 0;
X		} else if (t[1].u > t[0].u) {
X			hi = 1; mid = 0; lo = 2;
X		} else {
X			hi = 0; mid = 1; lo = 2;
X		}
X	}
X	if (t[hi].u - t[lo].u == 0.) {
X		/* degenerate axis */
X		dpdv->x = dpdv->y = dpdv->z = 0.;
X	} else {
X		/*
X		 * Given u coordinates of vertices forming the
X		 * 'long' edge, find where 'middle'
X		 * vertex falls on that edge given its u coordinate.
X		 */
X		scale = (t[mid].u - t[lo].u) / (t[hi].u - t[lo].u);
X		VecComb(1.0 - scale, p[lo], scale, p[hi], &base);
X		/*
X		 * v axis extends from computed basepoint to
X		 * middle vertex -- but in which direction?
X		 */
X		if (t[mid].v < ((1.0 - scale)*t[lo].v + scale*t[hi].v))
X			VecSub(base, p[mid], dpdv);
X		else
X			VecSub(p[mid], base, dpdv);
X		(void)VecNormalize(dpdv);
X	}
X
X	/* sort v coordinates */
X	if (t[2].v > t[1].v) {
X		if (t[1].v > t[0].v) {
X			hi = 2; mid = 1; lo = 0;
X		} else if (t[2].v > t[0].v) {
X			hi = 2; mid = 0; lo = 1;
X		} else {
X			hi = 0; mid = 2; lo = 1;
X		}
X	} else {
X		if (t[2].v > t[0].v) {
X			hi = 1; mid = 2; lo = 0;
X		} else if (t[1].v > t[0].v) {
X			hi = 1; mid = 0; lo = 2;
X		} else {
X			hi = 0; mid = 1; lo = 2;
X		}
X	}
X	if (t[hi].v - t[lo].v == 0.) {
X		/* degenerate axis */
X		dpdu->x = dpdu->y = dpdu->z = 0.;
X	} else {
X		/*
X		 * Given v coordinates of vertices forming the
X		 * 'long' edge, find where 'middle'
X		 * vertex falls on that edge given its v coordinate.
X		 */
X		scale = (t[mid].v - t[lo].v) / (t[hi].v - t[lo].v);
X		VecComb(1.0 - scale, p[lo], scale, p[hi], &base);
X		/*
X		 * u axis extends from computed basepoint to
X		 * middle vertex -- but in which direction?
X		 */
X		if (t[mid].u < ((1.0 - scale)*t[lo].u + scale*t[hi].u))
X			VecSub(base, p[mid], dpdu);
X		else
X			VecSub(p[mid], base, dpdu);
X		(void)VecNormalize(dpdu);
X	}
X}
X
Xvoid
XTriangleMethodRegister(meth)
XUserMethodType meth;
X{
X	if (iTriangleMethods)
X		iTriangleMethods->user = meth;
X}
END_OF_FILE
if test 11924 -ne `wc -c <'libray/libobj/triangle.c'`; then
    echo shar: \"'libray/libobj/triangle.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/triangle.c'
fi
if test -f 'rayshade/raytrace.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'rayshade/raytrace.c'\"
else
echo shar: Extracting \"'rayshade/raytrace.c'\" \(11786 characters\)
sed "s/^X//" >'rayshade/raytrace.c' <<'END_OF_FILE'
X/*
X * raytrace.c
X *
X * Copyright (C) 1989, 1991, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: raytrace.c,v 4.0 91/07/17 14:50:49 kolb Exp Locker: kolb $
X *
X * $Log:	raytrace.c,v $
X * Revision 4.0  91/07/17  14:50:49  kolb
X * Initial version.
X * 
X */
X
X#include "rayshade.h"
X#include "libsurf/atmosphere.h"
X#include "libsurf/surface.h"
X#include "libcommon/sampling.h"
X#include "options.h"
X#include "stats.h"
X#include "raytrace.h"
X#include "viewing.h"
X
X#define UNSAMPLED	-1
X#define SUPERSAMPLED	-2
X
Xtypedef struct {
X	Pixel	*pix;	/* Pixel values */
X	int	*samp;	/* Sample number */
X} Scanline;
X
Xstatic int		*SampleNumbers;
Xstatic void	RaytraceInit();
X
Xstatic Ray	TopRay;				/* Top-level ray. */
XFloat		SampleTime();
X
XPixel		WhitePix = {1., 1., 1., 1.},
X		BlackPix = {0., 0., 0., 0.};
X
X/*
X * "Dither matrices" used to encode the 'number' of a ray that samples a
X * particular portion of a pixel.  Hand-coding is ugly, but...
X */
Xstatic int OneSample[1] = 	{0};
Xstatic int TwoSamples[4] =	{0, 2,
X				 3, 1};
Xstatic int ThreeSamples[9] =	{0, 2, 7,
X				 6, 5, 1,
X				 3, 8, 4};
Xstatic int FourSamples[16] =	{ 0,  8,  2, 10,
X				 12,  4, 14,  6,
X				  3, 11,  1,  9,
X				 15,  7, 13,  5};
Xstatic int FiveSamples[25] =	{ 0,  8, 23, 17,  2,
X				 19, 12,  4, 20, 15,
X				  3, 21, 16,  9,  6,
X				 14, 10, 24,  1, 13,
X				 22,  7, 18, 11,  5};
Xstatic int SixSamples[36] =	{ 6, 32,  3, 34, 35,  1,
X				  7, 11, 27, 28,  8, 30,
X				 24, 14, 16, 15, 23, 19,
X				 13, 20, 22, 21, 17, 18,
X				 25, 29, 10,  9, 26, 12,
X				 36,  5, 33,  4,  2, 31};
Xstatic int SevenSamples[49] =	{22, 47, 16, 41, 10, 35,  4,
X				  5, 23, 48, 17, 42, 11, 29,
X				 30,  6, 24, 49, 18, 36, 12,
X				 13, 31,  7, 25, 43, 19, 37,
X				 38, 14, 32,  1, 26, 44, 20,
X				 21, 39,  8, 33,  2, 27, 45,
X				 46, 15, 40,  9, 34,  3, 28};
Xstatic int EightSamples[64] =	{ 8, 58, 59,  5,  4, 62, 63,  1,
X				 49, 15, 14, 52, 53, 11, 10, 56,
X				 41, 23, 22, 44, 45, 19, 18, 48,
X				 32, 34, 35, 29, 28, 38, 39, 25,
X				 40, 26, 27, 37, 36, 30, 31, 33,
X				 17, 47, 46, 20, 21, 43, 42, 24,
X				  9, 55, 54, 12, 13, 51, 50, 16,
X				 64,  2,  3, 61, 60,  6,  7, 57};
X
Xvoid	AdaptiveRefineScanline(), FullySamplePixel(), FullySampleScanline(),
X	SingleSampleScanline();
Xstatic int	ExcessiveContrast();
Xstatic Scanline scan0, scan1, scan2;
X
X
Xvoid
Xraytrace(argc, argv)
Xint argc;
Xchar **argv;
X{
X	int y, *tmpsamp;
X	Pixel *tmppix;
X	Float usertime, systime, lasttime;
X
X	/*
X	 * If this is the first frame,
X	 * allocate scanlines, etc.
X	 */
X	if (Options.framenum == Options.startframe)
X		RaytraceInit();
X	/*
X	 * The top-level ray TopRay always has as its origin the
X	 * eye position and as its medium NULL, indicating that it
X	 * is passing through a medium with index of refraction
X	 * equal to DefIndex.
X	 */
X	TopRay.pos = Camera.pos;
X	TopRay.media = (Medium *)0;
X	TopRay.depth = 0;
X
X	/*
X	 * Always fully sample the bottom and top rows and the left
X	 * and right column of pixels.  This minimizes artifacts that
X	 * may arise when piecing together images.
X	 */
X	FullySampleScanline(0, &scan0);
X
X	SingleSampleScanline(1, &scan1);
X	FullySamplePixel(0, 1, &scan1.pix[0], &scan1.samp[0]);
X	FullySamplePixel(Screen.xsize -1, 1, &scan1.pix[Screen.xsize -1],
X		&scan1.samp[Screen.xsize -1]);
X
X	lasttime = 0;
X	for (y = 1; y < Screen.ysize; y++) {
X		SingleSampleScanline(y+1, &scan2);
X		FullySamplePixel(0, y+1, &scan2.pix[0], &scan2.samp[0]);
X		FullySamplePixel(Screen.xsize -1, y+1,
X			&scan2.pix[Screen.xsize -1],
X			&scan2.samp[Screen.xsize -1]);
X
X		if (Sampling.sidesamples > 1)
X			AdaptiveRefineScanline(y,&scan0,&scan1,&scan2);
X
X		PictureWriteLine(scan0.pix);
X
X		tmppix = scan0.pix;
X		tmpsamp = scan0.samp;
X		scan0.pix = scan1.pix;
X		scan0.samp = scan1.samp;
X		scan1.pix = scan2.pix;
X		scan1.samp = scan2.samp;
X		scan2.pix = tmppix;
X		scan2.samp = tmpsamp;
X
X		if ((y-1) % Options.report_freq == 0) {
X			fprintf(Stats.fstats,"Finished line %d (%lu rays",y-1,
X							Stats.EyeRays);
X			if (Options.verbose) {
X				/*
X				 * Report total CPU and split times.
X				 */
X				RSGetCpuTime(&usertime, &systime);
X				fprintf(Stats.fstats,", %2.2f sec,",
X						usertime+systime);
X				fprintf(Stats.fstats," %2.2f split",
X						usertime+systime-lasttime);
X				lasttime = usertime+systime;
X			}
X			fprintf(Stats.fstats,")\n");
X			(void)fflush(Stats.fstats);
X		}
X
X	}
X	/*
X	 * Supersample last scanline.
X	 */
X	for (y = 1; y < Screen.xsize -1; y++) {
X		if (scan0.samp[y] != SUPERSAMPLED)
X			FullySamplePixel(y, Screen.ysize -1,
X				&scan0.pix[y],
X				&scan0.samp[y]);
X	}
X	PictureWriteLine(scan0.pix);
X}
X
Xvoid
XSingleSampleScanline(line, data)
Xint line;
XScanline *data;
X{
X	Float upos, vpos, yp;
X	int x, usamp, vsamp;
X	Pixel tmp;
X
X	yp = line + Screen.miny - 0.5*Sampling.filterwidth;
X	for (x = 0; x < Screen.xsize; x++) {
X		/*
X		 * Pick a sample number...
X		 */
X		data->samp[x] = nrand() * Sampling.totsamples;
X		/*
X		 * Take sample corresponding to sample #.
X		 */
X		usamp = data->samp[x] % Sampling.sidesamples;
X		vsamp = data->samp[x] / Sampling.sidesamples;
X
X		vpos = yp + vsamp * Sampling.filterdelta;
X		upos = x + Screen.minx - 0.5*Sampling.filterwidth +
X				usamp*Sampling.filterdelta;
X		if (Options.jitter) {
X			vpos += nrand()*Sampling.filterdelta;
X			upos += nrand()*Sampling.filterdelta;
X		}
X		TopRay.time = SampleTime(SampleNumbers[data->samp[x]]);
X		SampleScreen(upos, vpos, &TopRay,
X			&data->pix[x], SampleNumbers[data->samp[x]]);
X		if (Options.samplemap)
X			data->pix[x].alpha = 0;
X	}
X}
X
Xvoid
XFullySampleScanline(line, data)
Xint line;
XScanline *data;
X{
X	int x;
X
X	for (x = 0; x < Screen.xsize; x++) {
X		data->samp[x] = UNSAMPLED;
X		FullySamplePixel(x, line, &data->pix[x], &data->samp[x]);
X	}
X}
X
Xvoid
XFullySamplePixel(xp, yp, pix, prevsamp)
Xint xp, yp;
XPixel *pix;
Xint *prevsamp;
X{
X	Float upos, vpos, u, v;
X	int x, y, sampnum;
X	Pixel ctmp;
X
X	if (*prevsamp == SUPERSAMPLED)
X		return;	/* already done */
X
X	Stats.SuperSampled++;
X	if (*prevsamp == UNSAMPLED) {
X		/*
X		 * No previous sample; initialize to black.
X		 */
X		pix->r = pix->g = pix->b = pix->alpha = 0.;
X	} else {
X		if (Sampling.sidesamples == 1) {
X			*prevsamp = SUPERSAMPLED;
X			return;
X		}
X		x = *prevsamp % Sampling.sidesamples;
X		y = *prevsamp / Sampling.sidesamples;
X		pix->r *= Sampling.filter[x][y];
X		pix->g *= Sampling.filter[x][y];
X		pix->b *= Sampling.filter[x][y];
X		pix->alpha *= Sampling.filter[x][y];
X	}
X
X	sampnum = 0;
X	xp += Screen.minx;
X	vpos = Screen.miny + yp - 0.5*Sampling.filterwidth;
X	for (y = 0; y < Sampling.sidesamples; y++,
X	     vpos += Sampling.filterdelta) {
X		upos = xp - 0.5*Sampling.filterwidth;
X		for (x = 0; x < Sampling.sidesamples; x++,
X		     upos += Sampling.filterdelta) {
X			if (sampnum != *prevsamp) {
X				if (Options.jitter) {
X					u = upos + nrand()*Sampling.filterdelta;
X					v = vpos + nrand()*Sampling.filterdelta;
X				} else {
X					u = upos;
X					v = vpos;
X				}
X				TopRay.time = SampleTime(SampleNumbers[sampnum]);
X				SampleScreen(u, v, &TopRay, &ctmp,
X					SampleNumbers[sampnum]);
X				pix->r += ctmp.r*Sampling.filter[x][y];
X				pix->g += ctmp.g*Sampling.filter[x][y];
X				pix->b += ctmp.b*Sampling.filter[x][y];
X				pix->alpha += ctmp.alpha*Sampling.filter[x][y];
X			}
X			if (++sampnum == Sampling.totsamples)
X				sampnum = 0;
X		}
X	}
X
X	if (Options.samplemap)
X		pix->alpha = 255;
X
X	*prevsamp = SUPERSAMPLED;
X}
X
Xvoid
XAdaptiveRefineScanline(y, scan0, scan1, scan2)
Xint y;
XScanline *scan0, *scan1, *scan2;
X{
X	int x, done;
X
X	/*
X	 * Walk down scan1, looking at 4-neighbors for excessive contrast.
X	 * If found, supersample *all* neighbors not already supersampled.
X	 * The process is repeated until either there are no
X	 * high-contrast regions or all such regions are already supersampled.
X	 */
X
X	do {
X		done = TRUE;
X		for (x = 1; x < Screen.xsize -1; x++) {
X			/*
X		 	 * Find min and max RGB for area we care about
X			 */
X			if (ExcessiveContrast(x, scan0->pix, scan1->pix,
X			    scan2->pix)) {
X				if (scan1->samp[x-1] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x-1, y,
X						&scan1->pix[x-1],
X						&scan1->samp[x-1]);
X				}
X				if (scan0->samp[x] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x, y-1,
X						&scan0->pix[x],
X						&scan0->samp[x]);
X				}
X				if (scan1->samp[x+1] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x+1, y,
X						&scan1->pix[x+1],
X						&scan1->samp[x+1]);
X				}
X				if (scan2->samp[x] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x, y+1,
X						&scan2->pix[x],
X						&scan2->samp[x]);
X				}
X				if (scan1->samp[x] != SUPERSAMPLED) {
X					done = FALSE;
X					FullySamplePixel(x, y,
X						&scan1->pix[x],
X						&scan1->samp[x]);
X				}
X			}
X		}
X	} while (!done);
X}
X
Xstatic int
XExcessiveContrast(x, pix0, pix1, pix2)
Xint x;
XPixel *pix0, *pix1, *pix2;
X{
X	Float mini, maxi, sum, diff;
X
X	maxi = max(pix0[x].r, pix1[x-1].r);
X	if (pix1[x].r > maxi) maxi = pix1[x].r;
X	if (pix1[x+1].r > maxi) maxi = pix1[x+1].r;
X	if (pix2[x].r > maxi) maxi = pix2[x].r;
X
X	mini = min(pix0[x].r, pix1[x-1].r);
X	if (pix1[x].r < mini) mini = pix1[x].r;
X	if (pix1[x+1].r < mini) mini = pix1[x+1].r;
X	if (pix2[x].r < mini) mini = pix2[x].r;
X
X	diff = maxi - mini;
X	sum = maxi + mini;
X	if (sum > EPSILON && diff/sum > Options.contrast.r)
X		return TRUE;
X
X	maxi = max(pix0[x].g, pix1[x-1].g);
X	if (pix1[x].g > maxi) maxi = pix1[x].g;
X	if (pix1[x+1].g > maxi) maxi = pix1[x+1].g;
X	if (pix2[x].g > maxi) maxi = pix2[x].g;
X
X	mini = min(pix0[x].g, pix1[x-1].g);
X	if (pix1[x].g < mini) mini = pix1[x].g;
X	if (pix1[x+1].g < mini) mini = pix1[x+1].g;
X	if (pix2[x].g < mini) mini = pix2[x].g;
X
X	diff = maxi - mini;
X	sum = maxi + mini;
X
X	if (sum > EPSILON && diff/sum > Options.contrast.g)
X		return TRUE;
X
X	maxi = max(pix0[x].b, pix1[x-1].b);
X	if (pix1[x].b > maxi) maxi = pix1[x].b;
X	if (pix1[x+1].b > maxi) maxi = pix1[x+1].b;
X	if (pix2[x].b > maxi) maxi = pix2[x].b;
X
X	mini = min(pix0[x].b, pix1[x-1].b);
X	if (pix1[x].b < mini) mini = pix1[x].b;
X	if (pix1[x+1].b < mini) mini = pix1[x+1].b;
X	if (pix2[x].b < mini) mini = pix2[x].b;
X
X	diff = maxi - mini;
X	sum = maxi + mini;
X	if (sum > EPSILON && diff/sum > Options.contrast.b)
X		return TRUE;
X
X	return FALSE;
X}
X
XFloat
XSampleTime(sampnum)
Xint sampnum;
X{
X	Float window, jitter = 0.0, res;
X
X	if (Options.shutterspeed <= 0.)
X		return Options.framestart;
X	if (Options.jitter)
X		jitter = nrand();
X	window = Options.shutterspeed / Sampling.totsamples;
X	res = Options.framestart + window * (sampnum + jitter);
X	TimeSet(res);
X	return res;
X}
X
Xstatic void
XRaytraceInit()
X{
X
X	switch (Sampling.sidesamples) {
X		case 1:
X			SampleNumbers = OneSample;
X			break;
X		case 2:
X			SampleNumbers = TwoSamples;
X			break;
X		case 3:
X			SampleNumbers = ThreeSamples;
X			break;
X		case 4:
X			SampleNumbers = FourSamples;
X			break;
X		case 5:
X			SampleNumbers = FiveSamples;
X			break;
X		case 6:
X			SampleNumbers = SixSamples;
X			break;
X		case 7:
X			SampleNumbers = SevenSamples;
X			break;
X		case 8:
X			SampleNumbers = EightSamples;
X			break;
X		default:
X			RLerror(RL_PANIC,
X				"Sorry, %d rays/pixel not supported.\n",
X					Sampling.totsamples);
X	}
X
X	/*
X 	 * Allocate pixel arrays and arrays to store sampling info.
X 	 */
X	scan0.pix = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
X	scan1.pix = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
X	scan2.pix = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
X
X	scan0.samp = (int *)Malloc(Screen.xsize * sizeof(int));
X	scan1.samp = (int *)Malloc(Screen.xsize * sizeof(int));
X	scan2.samp = (int *)Malloc(Screen.xsize * sizeof(int));
X}
END_OF_FILE
if test 11786 -ne `wc -c <'rayshade/raytrace.c'`; then
    echo shar: \"'rayshade/raytrace.c'\" unpacked with wrong size!
fi
# end of 'rayshade/raytrace.c'
fi
echo shar: End of archive 14 \(of 19\).
cp /dev/null ark14isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 19 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0