Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/unix/volume21/rayshade/part05

Internet Message Format | 1990-03-21 | 47.2 KB

Subject:  v21i012:  A ray tracing program, Part05/08
Newsgroups: comp.sources.unix
Sender: sources
Approved: [e-mail]

Submitted-by: Craig Kolb <[e-mail]>
Posting-number: Volume 21, Issue 12
Archive-name: rayshade/part05

#! /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 5 (of 8)."
# Contents:  doc/texture.ms src/grid.c src/matrix.c src/shade.c
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'doc/texture.ms' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'doc/texture.ms'\"
else
echo shar: Extracting \"'doc/texture.ms'\" \(13859 characters\)
sed "s/^X//" >'doc/texture.ms' <<'END_OF_FILE'
X.\"
X.\" Brief tutorial on adding textures to rayshade.
X.\" Craig Kolb 10/89
X.\"
X.\" $Id: texture.ms,v 3.0 89/10/23 16:42:57 craig Exp $
X.\"
X.\" $Log:    texture.ms,v $
X.\" Revision 3.0  89/10/23  16:42:57  craig
X.\" Baseline for first official release.
X.\"
X.de D(
X.DS
X.nr PS 8
X.ps 8
X.nr VS 12p
X.vs 12p
X.cs 1 20
X..
X.de D)
X.DE
X.nr PS 10
X.ps 10
X.nr VS 18p
X.vs 18p
X.cs 1
X..
X.DS
X.ND October, 1989
X.RP
X.de ]H
X.nr PS 10
X.ps 10
X'sp |0.5i-1
X.tl 'Rayshade Texture Tutorial'-%-'Draft'
X.ps
X.ft
X'sp |1.0i
X.ns
X..
X.wh0 ]H
X.pn 2
X.LP
X.TL
XAdding Textures to Rayshade
X.TL
X\fR\fIDraft\fR
X.AU
XCraig Kolb
X.AI
XDepartment of Mathematics
XYale University
XNew Haven, CT  06520
X.sp .5i
X.nr VS 18pts
X.PP
XThis tutorial describes the process of adding new textures to
Xrayshade.  Although the texturing interface is relatively
Xstraight-forward, it is difficult to see the "big picture" by
Xstudying the source code, as changes must be made in a
Xnumber of different files.  While this tutorial is primarily
Xmeant for those interested getting their hands dirty, and
Xassumes familiarity with solid texturing,
Xit provides a overview of the design of at least one portion
Xof rayshade and thus may be of interest even if you are not
Xplanning on making modifications.
X.LP
XAdding new textures to rayshade involves modifying at least
Xfour source files:
X.NH texture.h
Xtexture.h
X.LP
XA numerical type is given to the texture.  The routines to
Xcreate a reference to the new texture and the texturing function
Xitself are declared.
X.NH
Xtexture.c
X.LP
XAt least two new routines are added.  The first creates and
Xreturns a reference to a texture of the given type, and the
Xsecond performs the actual texturing.  In addition, an
Xarray of pointers to functions is modified to include the
Xnew texturing function.
X.NH
Xinput_lex.l
X.LP
XThe keyword used to identify the new texture is added to
Xthe list of recognized keywords.
X.NH
Xinput_yacc.y
X.LP
XThe new texture is added to the list of textures parsed by
Xthe yacc grammar.  The added code will call the routine
Xwhich creates and returns a reference to the texture type
Xwhich was defined in texture.c
X.PP
XIn this tutorial, a new texture,
X.I mountain,
Xis added to rayshade.  This
Xtexture will modify the diffuse component of a surface as a function of the
XZ component of the point of intersection.  If the name of a colormap
Xis given, an index into
Xthe colormap is used to a color to be used as the diffuse component the
Xsurface.
XOtherwise, the diffuse component of the surface is simply scaled.
XTo avoid strictly horizontal boundries
Xbetween colors when using a colormap, we add a bit of "noise" to the
XZ component of the point of intersection.  The magnitude and nature of the
Xnoise
Xis defined by the user.
X.br
X.NR PS 12
X.ps 12
X.sp .5
X\fBThe Texture type\fR
X.nr PS 10
X.ps 10
X.sp .5
X.LP
XAll textures in rayshade are referenced using a single Texture structure.
XThis structure is defined as:
X.D(
Xtypedef struct Texture {
X        char type;              /* Texture type */
X        Surface *surf1;         /* Alternate surface */
X        double size;            /* Scale/size factor */
X        double *args;           /* Random arguments. */
X        Color *colormap;        /* Colormap */
X        Trans *trans;           /* Transformation matrices. */
X        struct Texture *next;   /* Pointer to next texture. */
X} Texture;
X.D)
X.LP
XThe
X.I type
Xfield is used by apply_textures() to determine which texturing
Xfunction to call.  The
X.I trans
Xfield and
X.I next
Xfield are used internally
Xby rayshade.
X.LP
XThe rest of the fields are for use by the texturing functions.
X.LP
XThe
X.I args
Xfield provides a pointer to space for storing an arbitrary
Xnumber of floating-point parameters.  The
X.I size
Xfield is a handy
Xgeneral-purpose floating-point value (the idea being that you get
Xone parameter "free"
Xwith every Texture structure).
XThe
X.I colormap
Xfield is generally used to store an array of Colors read from
Xan ascii colormap file.
XThe
X.I surf1
Xfield is often set to point to a secondary surface.  This is
Xuseful if the texture performs some sort of interpolation between two
Xsurfaces -- the first being the surface assigned to the object being textured
Xand the second specified as an argument to the texture
X(e.g., the blotch texture).
Xfor an example).
X.NH 0
XModifying texture.h
X.LP
XThe file texture.h contains a list of #defines which look something like:
X.D(
X#define CHECKER         0       /* Checkerboard */
X#define BLOTCH          1       /* Color blotches */
X#define BUMP            2       /* Bump mapping */
X#define MARBLE          3       /* marble texture */
X#define FBM             4       /* fBm texture */
X#define FBMBUMP         5       /* fBm bump map */
X#define WOOD            6       /* wood-like texture */
X.D)
X.LP
XThese numerical types are used to identify the type of a given
Xtexture structure (via the
X.I type
Xfield in the Texture structure).
XWe need to add a similar definition for our new texture.
XAfter the WOOD definition, we add:
X.D(
X#define MOUNTAIN        7       /* bad mountain texture */
X.D)
X.LP
XIn addition, we must declare the two new functions which we will
Xadd to texture.c.  The first function, NewMountainText(), will
Xreturn a pointer to a Texture structure:
X.D(
XTexture *NewMountainText();
X.D)
XThe second routine, MountainText, returns nothing, but needs to be
Xdeclared:
X.D(
Xint MountainText();
X.D)
X.NH
XModifying texture.c
X.LP
XFirstly, we must include the new texturing function in the array of
Xtexturing functions used by rayshade.  This array, indexed by texture
Xtype, is used by apply_textures() to call the correct texturing function.
X.LP
XSo, we modify the textures[] definition to read:
X.D(
Xint (*textures[])() =
X        {CheckerText, BlotchText, BumpText, MarbleText, fBmText, fBmBumpText,
X         WoodText, MountainText};
X.D)
X.LP
XNote that MOUNTAIN was defined to be 7, and that MountainText is texture
Xnumber 7 (starting from 0) in the array.
X.LP
XNext, we need to write NewMountainText(), which will create and return a
Xreference to the texture.  Our new texture will be a function of 5 parameters:
X.D(
X        scale           amount to scale \fINoise()\fR by
X        omega, lambda   fBm parameters
X        octaves         number of octaves of \fINoise()\fR in fBm
X        colormap        name of colormap file, if any
X.D)
XThus, we add to the end of texture.c:
X.D(
XTexture *
XNewMountainText(scale, omega, lambda, octaves, mapname)
Xdouble scale, omega, lambda;
Xint octaves;
Xchar *mapname;
X{
X        /*
X         * Pointer to new texture.
X         */
X        Texture *text;
X
X        /*
X         * Allocate new texture of type MOUNTAIN
X         */
X        text = new_texture(MOUNTAIN);
X        /*
X         * Allocate space to store omega, lambda and octaves.
X         */
X        text->args = (double *)Malloc(3 * sizeof(double));
X        text->args[0] = omega;
X        text->args[1] = lambda;
X        text->args[2] = (double)octaves;
X        /*
X         * scale is stored in 'size'.
X         */
X        text->size = scale;
X        /*
X         * If a colormap name was specified, read it into 'colormap'.
X         */
X        if (mapname != (char *)0)
X                text->colormap = read_colormap(mapname);
X        /*
X         * All done -- return new texture.
X         */
X        return text;
X}
X.D)
X.LP
XThus, NewMountainText is called with the desired parameters and a
Xnew Texture is returned.
X.LP
XFinally, we must write the routine which actually performs the texturing.
XEach texturing function is called by apply_textures() with the following
Xarguments:
X.D(
X        text
X                a pointer to the Texture being applied
X        pos
X                a pointer to the coordinates of the point of intersection
X        norm
X                a pointer to the surface normal at the point of intersection
X        surf
X                the pointer to a copy of the surface of the object being
X                textured  (a copy is used so that the surface can be
X                modified for a given shading calculation without affecting
X                subsequent calculations).
X.D)
X.LP
XThus, we write:
X.D(
XMountainText(text, pos, norm, surf)
XTexture *text;
XVector *pos, *norm;
XSurface *surf;
X{
X        double val;
X        int index;
X
X        /*
X         * Compute value of fBm (fractional Brownian motion) for
X         * the given point.
X         */
X        val = fBm(pos, text->args[0], text->args[1], (int)text->args[2]);
X        /*
X         * Scale the result as requested and add in the Z component of
X         * the point of intersection.  Note that in a better texture
X         * we'd probably have additional parameters to afford
X         * greater control of val.
X         */
X        val = pos->z + text->size * val;
X
X        if (text->colormap) {
X                /*
X                 * If we're using a colormap, compute an index into
X                 * the colormap and use the appropriate color as the
X                 * diffuse component of the surface.
X                 */
X                index = 255. * val;
X                if (index > 255)
X                        index = 255;
X                if (index < 0)
X                        index = 0;
X                surf->diff = text->colormap[index];
X        } else {
X                /*
X                 * If there's no colormap, simply scale the diffuse
X                 * component.
X                 */
X                surf->diff = ScaleColor(val, surf->diff);
X        }
X}
X.D)
X.LP
X.NH
Xinput_lex.l
X.LP
XNow that we have the hard parts written, all that is left is making
Xthe parser recognize the new texture.  To do this, we first need to
Xadd a keyword for our texture to the list of keywords recognized by
Xrayshade.  This is done by editing input_lex.l.
X.LP
XThe file input_lex.l contains, among other things, an alphabetical list of
Xrayshade's keywords.  To add a new keyword, one simply follows the
Xexample of the other keywords.  Thus, we add the line:
X.D(
Xmountain                     {return(tMOUNTAIN);}
X.D)
Xbetween the lines defining 'mist' and 'object'.  This line directs
Xlex to return the token tMOUNTAIN whenever the string "mountain" is
Xencountered in an input file.
X.NH
Xinput_yacc.y
X.LP
XFinally, we need to write the code which will actually create an instance
Xof the new texture by calling NewMountainText().  This is done in
Xinput_yacc.y.
X.LP
XIn input_yacc.y, there are a series of lines which look something like:
X.D(
X%token tBACKGROUND tBLOTCH tBOX tBUMP tCONE tCYL tDIRECTIONAL
X%token tENDDEF tEXTENDED tEYEP tFBM tFBMBUMP tFOG tFOV tGRID
X%token tHEIGHTFIELD tLIGHT tLIST tLOOKP tMARBLE tMAXDEPTH tMIST
X%token tOBJECT tOUTFILE
X%token tPLANE tPOINT tPOLY tROTATE
X%token tSCALE tSCREEN tSPHERE tSTARTDEF tSUPERQ tSURFACE tRESOLUTION
X%token tTHRESH tTRANSLATE tTRANSFORM tTRIANGLE tUP tENDFILE
X%token tTEXTURE tCHECKER tWOOD
X.D)
X.LP
XThese lines declare the tokens returned by lex.  We need to declare
XtMOUNTAIN in a similar manner.  So, we change the last line to
Xread:
X.D(
X%token tTEXTURE tCHECKER tWOOD tMOUNTAIN
X.D)
XNext, we need to call NewMountainText() in the proper place.  In input_yacc.y,
Xthere is a production which reads something like:
X.D(
XTexturetype     : tCHECKER String
X                {
X                        $$ = NewCheckText($2);
X                }
X                | ...
X                ...
X                | tWOOD
X                {
X                        $$ = NewWoodText();
X                }
X                ;
X.D)
X.LP
XThese productions invoke the proper texture creation routine when appropriate.
XFor example, when the keyword corresponding to tCHECKER is followed by
Xa String, yacc will invoke NewCheckText() with the string (the name of
Xa surface, in this case) as an argument.  The Yacc grammar understands
Xthe following datatypes, among others:
X.D(
X        String          a series of alphanumerics surrounded by
X                        white space (i.e., the string need not be quoted)
X        Fnumber         a floating-point number
X        tINT            an integer
X        Vector          a vector (x, y, z)
X        Color           a color (r, g, b)
X.D)
XTo add a texture to the list of recognized textures, we change:
X.D(
X                ...
X                | tWOOD
X                {
X                        $$ = NewWoodText();
X                }
X                ;
X.D)
Xto:
X.D(
X                | tWOOD
X                {
X                        $$ = NewWoodText();
X                }
X                | tMOUNTAIN Fnumber Fnumber Fnumber tINT
X                {
X                        $$ = NewMountainText($2, $3, $4, $5, (char *)0);
X                }
X                | tMOUNTAIN Fnumber Fnumber Fnumber tINT String
X                {
X                        $$ = NewMountainText($2, $3, $4, $5, $6);
X                }
X                ;
X.D)
X.LP
XThe first new production invokes NewMountainText() when the keyword
Xassociated with tMOUNTAIN ("mountain") appears in an appropriate place
Xin the input file followed by
X.I four
Xparameters (scale, omega,
Xlambda, and octaves).  In this case, NewMountainText is passed a
XNULL pointer as the name of the colormap to use.
XThis code creates a reference to a mountain texture that does
X.I not
Xmake
Xuse of a colormap.  So, this production is invoked whenever a line
Xsuch as the following is encountered in the input file:
X.D(
X        texture mountain 0.2 0.5 2.0 6
X.D)
X.LP
XThe second production works in a similar manner, except that it passes
Xa colormap name to NewMountainText().  It handles lines such as:
X.D(
X        texture mountain 0.2 0.5 2.0 6 mountain.map
X.D)
X.NH
XCompiling
X.LP
XThat, in theory, is all there is to it.  Run 'make' to recompile
Xinput_lex.o, input_yacc.o, and texture.o.
X.NH
XTesting
X.LP
XA good test input file for the new texture might be something like:
X.D(
Xscreen 512 512
Xeyep 0 -10 0
Xlookp 0 0 0
X
Xfov 20.
X
Xlight 1.0 directional 1. -1. 1.
Xsurface boring .1 .1 .1 .8 .8 .8 0 0 0 0 0 0 0
X
Xsphere boring 1. 0. 0. 0. texture mountain 0.2 0.5 2.0 6 planet.map
X.D)
END_OF_FILE
if test 13859 -ne `wc -c <'doc/texture.ms'`; then
    echo shar: \"'doc/texture.ms'\" unpacked with wrong size!
fi
# end of 'doc/texture.ms'
fi
if test -f 'src/grid.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'src/grid.c'\"
else
echo shar: Extracting \"'src/grid.c'\" \(9938 characters\)
sed "s/^X//" >'src/grid.c' <<'END_OF_FILE'
X/*
X * grid.c
X *
X * Copyright (C) 1989, Craig E. Kolb
X *
X * This software may be freely copied, modified, and redistributed,
X * provided that this copyright notice is preserved on all copies.
X *
X * There is no warranty or other guarantee of fitness for this software,
X * it is provided solely .  Bug reports or fixes may be sent
X * to the author, who may or may not act on them as he desires.
X *
X * You may not include this software in a program or other software product
X * without supplying the source, or without informing the end-user that the
X * source is available for no extra charge.
X *
X * If you modify this software, you should include a notice giving the
X * name of the person performing the modification, the date of modification,
X * and the reason for such modification.
X *
X * $Id: grid.c,v 3.0 89/10/27 02:05:50 craig Exp $
X *
X * $Log:    grid.c,v $
X * Revision 3.0  89/10/27  02:05:50  craig
X * Baseline for first official release.
X *
X */
X#include <stdio.h>
X#include "constants.h"
X#include "typedefs.h"
X#include "funcdefs.h"
X
Xextern double CheckVoxel(), intersect();
Xunsigned long    raynumber = 1;        /* Current "ray number". */
X                    /* (should be "grid number") */
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 */
Xdouble
Xint_grid(grid, source, ray, hitinfo)
XGrid *grid;
XPrimitive *source;
XRay *ray;
XHitInfo *hitinfo;
X{
X    ObjList *list;
X    Object *tmpobj;
X    double dist, mindist, offset, tMaxX, tMaxY, tMaxZ;
X    double tDeltaX, tDeltaY, tDeltaZ, *raybounds[2][3];
X    int stepX, stepY, stepZ, outX, outY, outZ, x, y, z;
X    HitInfo infotmp;
X    Vector curpos, nXp, nYp, nZp, np, pDeltaX, pDeltaY, pDeltaZ;
X    unsigned long counter;
X
X    mindist = FAR_AWAY;
X
X    /*
X     * The "raynumber" scheme prevents us from performing unnecessary
X     * ray/object intersection tests. Since an object may span more than
X     * one voxel, it is possible that such an object will be tested more
X     * than once as we move from voxel-to-voxel. The idea here is to
X     * increment "raynumber" once when we begin a ray/grid intersection
X     * test. Whenever we check for intersection with an object in
X     * the grid, we assign that object's "counter" field the current value
X     * of "raynumber".  On subsequent calls to CheckVoxel, if the object's
X     * "counter" is greater than or equal to the value "raynumber" had
X     * when we started tracing the current grid, we know that we've already
X     * checked that object.  This isn't as straight-forward as one
X     * might think due to instantiation -- one may have several
X     * instances of the same object in a grid.  By incrementing raynumber
X     * whenever we being a new ray/grid int. test, we guarantee that
X     * we will check for intersection with the objects in the grid
X     * at most once.
X     */
X    infotmp.totaltrans = hitinfo->totaltrans;
X    /*
X     * Check unbounded objects.
X     */
X    for (list = grid->unbounded ; list; list = list->next) {
X        tmpobj = list->data;
X        dist = intersect(tmpobj, source, ray, &infotmp);
X        if (dist > 0. && dist < mindist) {
X            *hitinfo = infotmp;
X            mindist = dist;
X        }
X    }
X
X    /*
X     * If outside of the bounding box, check to see if we
X     * hit it.
X     */
X    if (OutOfBounds(&ray->pos, grid->bounds)) {
X        offset = IntBounds(ray, grid->bounds);
X        if (offset <= 0.)
X            /*
X             * Ray never hit grid space.
X             */
X            return (mindist == FAR_AWAY ? 0. : mindist);
X        else if (mindist < offset)
X            /*
X             * Hit an unbounded object closer than grid.
X             */
X            return mindist;
X        /*
X         * else
X         *    The ray enters voxels space before it hits
X         *     an unbounded object.
X         */
X        addscaledvec(ray->pos, offset, ray->dir, &curpos);
X    } else {
X        offset = 0.;
X        curpos = ray->pos;
X    }
X
X    counter = raynumber++;
X    /*
X     * Kludge:  Voxel space is mapped as [start, end),
X     * and we want it to be [start, end].
X     */
X    x = x2voxel(grid, curpos.x);
X    if (x == grid->xsize)
X        x = grid->xsize -1;
X    y = y2voxel(grid, curpos.y);
X    if (y == grid->ysize)
X        y = grid->ysize -1;
X    z = z2voxel(grid, curpos.z);
X    if (z == grid->zsize)
X        z = grid->zsize -1;
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    if (ray->dir.x < 0.) {
X        stepX = outX = -1;
X        tMaxX = (voxel2x(grid, x) - curpos.x) / ray->dir.x;
X        tDeltaX = grid->voxsize[X] / - ray->dir.x;
X        raybounds[LOW][X] = &np.x;
X        raybounds[HIGH][X] = &curpos.x;
X    } else if (ray->dir.x > 0.) {
X        stepX = 1;
X        outX = grid->xsize;
X        tMaxX = (voxel2x(grid, x+1) - curpos.x) / ray->dir.x;
X        tDeltaX = grid->voxsize[X] / ray->dir.x;
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    if (ray->dir.y < 0.) {
X        stepY = outY = -1;
X        tMaxY = (voxel2y(grid, y) - curpos.y) / ray->dir.y;
X        tDeltaY = grid->voxsize[Y] / - ray->dir.y;
X        raybounds[LOW][Y] = &np.y;
X        raybounds[HIGH][Y] = &curpos.y;
X    } else if (ray->dir.y > 0.) {
X        stepY = 1;
X        outY = grid->ysize;
X        tMaxY = (voxel2y(grid, y+1) - curpos.y) / ray->dir.y;
X        tDeltaY = grid->voxsize[Y] / ray->dir.y;
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    if (ray->dir.z < 0.) {
X        stepZ = outZ = -1;
X        tMaxZ = (voxel2z(grid, z) - curpos.z) / ray->dir.z;
X        tDeltaZ = grid->voxsize[Z] / - ray->dir.z;
X        raybounds[LOW][Z] = &np.z;
X        raybounds[HIGH][Z] = &curpos.z;
X    } else if (ray->dir.z > 0.) {
X        stepZ = 1;
X        outZ = grid->zsize;
X        tMaxZ = (voxel2z(grid, z+1) - curpos.z) / ray->dir.z;
X        tDeltaZ = grid->voxsize[Z] / ray->dir.z;
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    /*
X     * We've already moved from the ray origin along the ray to "hitpoint"
X     * by "offset" units in order to reach the grid.
X     */
X    tMaxX += offset;
X    tMaxY += offset;
X    tMaxZ += offset;
X
X    pDeltaX.x = tDeltaX * ray->dir.x;
X    pDeltaX.y = tDeltaX * ray->dir.y;
X    pDeltaX.z = tDeltaX * ray->dir.z;
X    pDeltaY.x = tDeltaY * ray->dir.x;
X    pDeltaY.y = tDeltaY * ray->dir.y;
X    pDeltaY.z = tDeltaY * ray->dir.z;
X    pDeltaZ.x = tDeltaZ * ray->dir.x;
X    pDeltaZ.y = tDeltaZ * ray->dir.y;
X    pDeltaZ.z = tDeltaZ * ray->dir.z;
X
X    nXp.x = ray->pos.x + tMaxX * ray->dir.x;
X    nXp.y = ray->pos.y + tMaxX * ray->dir.y;
X    nXp.z = ray->pos.z + tMaxX * ray->dir.z;
X    nYp.x = ray->pos.x + tMaxY * ray->dir.x;
X    nYp.y = ray->pos.y + tMaxY * ray->dir.y;
X    nYp.z = ray->pos.z + tMaxY * ray->dir.z;
X    nZp.x = ray->pos.x + tMaxZ * ray->dir.x;
X    nZp.y = ray->pos.y + tMaxZ * ray->dir.y;
X    nZp.z = ray->pos.z + tMaxZ * ray->dir.z;
X
X    while (1) {
X        if (tMaxX < tMaxY) {
X            if (tMaxX < tMaxZ) {
X                np = nXp;
X                if ((list = grid->cells[x][y][z]))
X                    mindist = CheckVoxel(list,source,ray,
X                    raybounds,hitinfo,counter,mindist);
X                x += stepX;
X                if (mindist < 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 {
X                np = nZp;
X                if ((list = grid->cells[x][y][z]))
X                    mindist = CheckVoxel(list,source,ray,
X                    raybounds, hitinfo,counter,mindist);
X                z += stepZ;
X                if (mindist < 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            }
X        } else {
X            if (tMaxY < tMaxZ) {
X                np = nYp;
X                if ((list = grid->cells[x][y][z]))
X                    mindist = CheckVoxel(list,source,ray,
X                    raybounds,hitinfo,counter,mindist);
X                y += stepY;
X                if (mindist < 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            } else {
X                np = nZp;
X                if ((list = grid->cells[x][y][z]))
X                    mindist = CheckVoxel(list,source,ray,
X                    raybounds, hitinfo,counter,mindist);
X                z += stepZ;
X                if (mindist < 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            }
X        }
X
X    }
X    if (mindist == FAR_AWAY)
X        return 0.;
X    return mindist;
X
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 */
Xdouble
XCheckVoxel(list,source,ray,raybounds,hitinfo,counter,mindist)
XObjList *list;
XPrimitive *source;
XRay *ray;
Xdouble *raybounds[2][3];
XHitInfo *hitinfo;
Xunsigned long counter;
Xdouble mindist;
X{
X    Object *obj;
X    HitInfo tmpinfo;
X    double dist, lx, hx, ly, hy, lz, hz;
X
X    tmpinfo.totaltrans = hitinfo->totaltrans;
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    for (; list; list = list->next) {
X        obj = list->data;
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 LINDA
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            continue;
X#ifdef LINDA
X        *obj->counter = counter;
X#else
X        obj->counter = counter;
X#endif
X        dist = intersect(obj, source, ray, &tmpinfo);
X        if (dist > 0. && dist < mindist) {
X            *hitinfo = tmpinfo;
X            mindist = dist;
X        }
X    }
X
X    return mindist;
X}
END_OF_FILE
if test 9938 -ne `wc -c <'src/grid.c'`; then
    echo shar: \"'src/grid.c'\" unpacked with wrong size!
fi
# end of 'src/grid.c'
fi
if test -f 'src/matrix.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'src/matrix.c'\"
else
echo shar: Extracting \"'src/matrix.c'\" \(9704 characters\)
sed "s/^X//" >'src/matrix.c' <<'END_OF_FILE'
X/*
X * matrix.c
X *
X * Copyright (C) 1989, Craig E. Kolb
X *
X * This software may be freely copied, modified, and redistributed,
X * provided that this copyright notice is preserved on all copies.
X *
X * There is no warranty or other guarantee of fitness for this software,
X * it is provided solely .  Bug reports or fixes may be sent
X * to the author, who may or may not act on them as he desires.
X *
X * You may not include this software in a program or other software product
X * without supplying the source, or without informing the end-user that the
X * source is available for no extra charge.
X *
X * If you modify this software, you should include a notice giving the
X * name of the person performing the modification, the date of modification,
X * and the reason for such modification.
X *
X * $Id: matrix.c,v 3.0 89/10/27 02:05:56 craig Exp $
X *
X * $Log:    matrix.c,v $
X * Revision 3.0  89/10/27  02:05:56  craig
X * Baseline for first official release.
X *
X */
X#include <math.h>
X#include <stdio.h>
X#include "typedefs.h"
X#include "constants.h"
X#include "funcdefs.h"
X
X/*
X * Matrices are indexed row-first; that is:
X * matrix[ROW][COLUMN]
X */
X/*
X * Allocate new structure which holds both object-to-world and
X * world-to-object space transformation structures.  It probably
X * should hold pointers to these structures.
X */
XTrans *
Xnew_trans(t, it)
XTransInfo *t, *it;
X{
X    Trans *res;
X
X    res = (Trans *)Malloc(sizeof(Trans));
X    res->obj2world = *t;
X    res->world2obj = *it;
X    return res;
X}
X
X/*
X * Allocate new transformation structure.
X */
XTransInfo *
Xnew_transinfo()
X{
X    TransInfo *res;
X
X    res = (TransInfo *)Malloc(sizeof(TransInfo));
X    init_trans(res);
X    return res;
X}
X
X/*
X * Set up matrix for counter-clockwise rotation about vector (x, y, z) by
X * theta radians.
X */
Xrotate_trans(trans, x, y, z, theta)
XTransInfo *trans;
Xdouble x, y, z, theta;
X{
X    double n1, n2, n3, sintheta, costheta;
X    Vector vector;
X
X    vector.x = x;
X    vector.y = y;
X    vector.z = z;
X    (void)normalize(&vector);
X
X    sintheta = sin(theta);
X    costheta = cos(theta);
X
X    n1 = vector.x; n2 = vector.y; n3 = vector.z;
X    trans->matrix[0][0] = (double)(n1*n1 + (1. - n1*n1)*costheta);
X    trans->matrix[1][0] = (double)(n1*n2*(1 - costheta) + n3*sintheta);
X    trans->matrix[2][0] = (double)(n1*n3*(1 - costheta) - n2*sintheta);
X    trans->matrix[0][1] = (double)(n1*n2*(1 - costheta) - n3*sintheta);
X    trans->matrix[1][1] = (double)(n2*n2 + (1 - n2*n2)*costheta);
X    trans->matrix[2][1] = (double)(n2*n3*(1 - costheta) + n1*sintheta);
X    trans->matrix[0][2] = (double)(n1*n3*(1 - costheta) + n2*sintheta);
X    trans->matrix[1][2] = (double)(n2*n3*(1 - costheta) - n1*sintheta);
X    trans->matrix[2][2] = (double)(n3*n3 + (1 - n3*n3)*costheta);
X    trans->translate.x = trans->translate.y = trans->translate.z = 0.;
X}
X
X/*
X * Apply translation by (vec) to trans.
X */
Xtranslate(trans, vec)
XTransInfo *trans;
XVector *vec;
X{
X    trans->translate.x += vec->x;
X    trans->translate.y += vec->y;
X    trans->translate.z += vec->z;
X}
X
X/*
X * Apply rotation about (dir) to matrix.
X */
Xrotate(trans, dir, radians)
XTransInfo *trans;
Xdouble radians;
XVector *dir;
X{
X    TransInfo ttmp;
X
X    rotate_trans(&ttmp, dir->x, dir->y, dir->z, radians);
X    mmult(trans, &ttmp, trans);
X}
X
X/*
X * Apply scale of (x, y, z) to trans.
X */
Xscale(trans, x, y, z)
XTransInfo *trans;
Xdouble x, y, z;
X{
X    trans->matrix[0][0] *= x;
X    trans->matrix[1][0] *= x;
X    trans->matrix[2][0] *= x;
X    trans->matrix[0][1] *= y;
X    trans->matrix[1][1] *= y;
X    trans->matrix[2][1] *= y;
X    trans->matrix[0][2] *= z;
X    trans->matrix[1][2] *= z;
X    trans->matrix[2][2] *= z;
X
X    trans->translate.x *= x;
X    trans->translate.y *= y;
X    trans->translate.z *= z;
X
X}
X
X/*
X * Multiply m1 and m2, copy result into "res".
X */
Xmmult(t1, t2, res)
XTransInfo *t1, *t2, *res;
X{
X    register int i, j, k;
X    TransInfo tmp;
X
X    for(i = 0; i < 3; i++) {
X        for(j = 0; j < 3;j++) {
X            tmp.matrix[i][j] = 0.;
X            for(k = 0; k < 3;k++)
X                tmp.matrix[i][j] += t1->matrix[i][k] *
X                            t2->matrix[k][j];
X
X        }
X    }
X    tmp.translate.x = t1->translate.x * t2->matrix[0][0] +
X                t1->translate.y * t2->matrix[1][0] +
X                t1->translate.z * t2->matrix[2][0] +
X                    t2->translate.x;
X    tmp.translate.y = t1->translate.x * t2->matrix[0][1] +
X                t1->translate.y * t2->matrix[1][1] +
X                t1->translate.z * t2->matrix[2][1] +
X                    t2->translate.y;
X    tmp.translate.z = t1->translate.x * t2->matrix[0][2] +
X                t1->translate.y * t2->matrix[1][2] +
X                t1->translate.z * t2->matrix[2][2] +
X                    t2->translate.z;
X    copy_trans(res, &tmp);
X}
X
X/*
X * Copy a given transformation structure.
X */
Xcopy_trans(into, from)
XTransInfo *into, *from;
X{
X    register int i;
X
X    for(i = 0;i < 3;i++) {
X        into->matrix[i][0] = from->matrix[i][0];
X        into->matrix[i][1] = from->matrix[i][1];
X        into->matrix[i][2] = from->matrix[i][2];
X    }
X
X    into->translate = from->translate;
X}
X
X/*
X * Initialize transformation structure.
X */
Xinit_trans(trans)
XTransInfo *trans;
X{
X    register int i, j;
X
X    for(i = 0;i < 3;i++)
X        for(j = 0;j < 3;j++)
X            trans->matrix[i][j] = (double)(i == j);
X
X    trans->translate.x = trans->translate.y = trans->translate.z = 0.;
X}
X
X/*
X * Calculate inverse of the given transformation structure.
X */
Xinvert_trans(inverse, trans)
XTransInfo *inverse, *trans;
X{
X    TransInfo ttmp;
X    int i;
X    double d;
X    extern int yylineno;
X
X    ttmp.matrix[0][0] = trans->matrix[1][1]*trans->matrix[2][2] -
X                trans->matrix[1][2]*trans->matrix[2][1];
X    ttmp.matrix[1][0] = trans->matrix[1][0]*trans->matrix[2][2] -
X                trans->matrix[1][2]*trans->matrix[2][0];
X    ttmp.matrix[2][0] = trans->matrix[1][0]*trans->matrix[2][1] -
X                trans->matrix[1][1]*trans->matrix[2][0];
X
X    ttmp.matrix[0][1] = trans->matrix[0][1]*trans->matrix[2][2] -
X                trans->matrix[0][2]*trans->matrix[2][1];
X    ttmp.matrix[1][1] = trans->matrix[0][0]*trans->matrix[2][2] -
X                trans->matrix[0][2]*trans->matrix[2][0];
X    ttmp.matrix[2][1] = trans->matrix[0][0]*trans->matrix[2][1] -
X                trans->matrix[0][1]*trans->matrix[2][0];
X
X    ttmp.matrix[0][2] = trans->matrix[0][1]*trans->matrix[1][2] -
X                trans->matrix[0][2]*trans->matrix[1][1];
X    ttmp.matrix[1][2] = trans->matrix[0][0]*trans->matrix[1][2] -
X                trans->matrix[0][2]*trans->matrix[1][0];
X    ttmp.matrix[2][2] = trans->matrix[0][0]*trans->matrix[1][1] -
X                trans->matrix[0][1]*trans->matrix[1][0];
X
X    d = trans->matrix[0][0]*ttmp.matrix[0][0] -
X        trans->matrix[0][1]*ttmp.matrix[1][0] +
X        trans->matrix[0][2]*ttmp.matrix[2][0];
X
X    if (abs(d) < EPSILON*EPSILON) {
X        fprintf(stderr,"Singular matrix (line %d):\n",yylineno);
X        for (i = 0; i < 3; i++)
X            fprintf(stderr,"%g %g %g\n", trans->matrix[i][0],
X                trans->matrix[i][1], trans->matrix[i][2]);
X        exit(0);
X    }
X
X    ttmp.matrix[0][0] /= d;
X    ttmp.matrix[0][2] /= d;
X    ttmp.matrix[1][1] /= d;
X    ttmp.matrix[2][0] /= d;
X    ttmp.matrix[2][2] /= d;
X
X    d = -d;
X
X    ttmp.matrix[0][1] /= d;
X    ttmp.matrix[1][0] /= d;
X    ttmp.matrix[1][2] /= d;
X    ttmp.matrix[2][1] /= d;
X
X    ttmp.translate.x = -(ttmp.matrix[0][0]*trans->translate.x +
X                 ttmp.matrix[1][0]*trans->translate.y +
X                 ttmp.matrix[2][0]*trans->translate.z);
X    ttmp.translate.y = -(ttmp.matrix[0][1]*trans->translate.x +
X                 ttmp.matrix[1][1]*trans->translate.y +
X                 ttmp.matrix[2][1]*trans->translate.z);
X    ttmp.translate.z = -(ttmp.matrix[0][2]*trans->translate.x +
X                 ttmp.matrix[1][2]*trans->translate.y +
X                 ttmp.matrix[2][2]*trans->translate.z);
X
X    copy_trans(inverse, &ttmp);
X}
X
X/*
X * Apply a transformation to a point (translation affects the point).
X */
Xtransform_point(vec, trans)
XVector *vec;
XTransInfo *trans;
X{
X    Vector tmp;
X
X    tmp.x = vec->x * trans->matrix[0][0] + vec->y * trans->matrix[1][0] +
X            vec->z * trans->matrix[2][0] + trans->translate.x;
X    tmp.y = vec->x * trans->matrix[0][1] + vec->y * trans->matrix[1][1] +
X            vec->z * trans->matrix[2][1] + trans->translate.y;
X    tmp.z = vec->x * trans->matrix[0][2] + vec->y * trans->matrix[1][2] +
X            vec->z * trans->matrix[2][2] + trans->translate.z;
X    *vec = tmp;
X}
X
X/*
X * Apply an explicit transformation to the given transformation structure.
X * 'c1x' is the X (0th) component of the first column, and so on.
X */
Xexplicit_trans(trans, c1x,c1y,c1z, c2x,c2y,c2z, c3x,c3y,c3z, tx, ty, tz, res)
XTransInfo *trans, *res;
Xdouble c1x, c1y, c1z, c2x, c2y, c2z, c3x, c3y, c3z, tx, ty, tz;
X{
X    TransInfo tmp;
X
X    tmp.matrix[0][0] = c1x;
X    tmp.matrix[1][0] = c1y;
X    tmp.matrix[2][0] = c1z;
X
X    tmp.matrix[0][1] = c2x;
X    tmp.matrix[1][1] = c2y;
X    tmp.matrix[2][1] = c2z;
X
X    tmp.matrix[0][2] = c3x;
X    tmp.matrix[1][2] = c3y;
X    tmp.matrix[2][2] = c3z;
X
X    tmp.translate.x = tx;
X    tmp.translate.y = ty;
X    tmp.translate.z = tz;
X
X    mmult(trans, &tmp, res);
X}
X
X/*
X * Apply transformation to a vector (translations have no effect).
X */
Xtransform_vector(vec, trans)
XVector *vec;
XTransInfo *trans;
X{
X    Vector tmp;
X
X    tmp.x = vec->x*trans->matrix[0][0] +
X        vec->y*trans->matrix[1][0] + vec->z*trans->matrix[2][0];
X    tmp.y = vec->x*trans->matrix[0][1] +
X        vec->y*trans->matrix[1][1] + vec->z*trans->matrix[2][1];
X    tmp.z = vec->x*trans->matrix[0][2] +
X        vec->y*trans->matrix[1][2] + vec->z*trans->matrix[2][2];
X
X    *vec = tmp;
X}
X
X/*
X * Transform "ray" by transforming the origin point and direction vector.
X */
Xdouble
XTransformRay(ray, trans)
XRay *ray;
XTransInfo *trans;
X{
X    transform_point(&ray->pos, trans);
X    transform_vector(&ray->dir, trans);
X    return normalize(&ray->dir);
X}
X
X/*
X * Transform normal -- multiply by the transpose of the given
X * matrix (which is the inverse of the desired transformation).
X */
XTransformNormal(norm, it)
XVector *norm;
XTransInfo *it;
X{
X    Vector onorm;
X
X    onorm = *norm;
X
X    norm->x = onorm.x*it->matrix[0][0] + onorm.y*it->matrix[0][1] +
X                onorm.z*it->matrix[0][2];
X    norm->y = onorm.x*it->matrix[1][0] + onorm.y*it->matrix[1][1] +
X                onorm.z*it->matrix[1][2];
X    norm->z = onorm.x*it->matrix[2][0] + onorm.y*it->matrix[2][1] +
X                onorm.z*it->matrix[2][2];
X}
END_OF_FILE
if test 9704 -ne `wc -c <'src/matrix.c'`; then
    echo shar: \"'src/matrix.c'\" unpacked with wrong size!
fi
# end of 'src/matrix.c'
fi
if test -f 'src/shade.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'src/shade.c'\"
else
echo shar: Extracting \"'src/shade.c'\" \(10618 characters\)
sed "s/^X//" >'src/shade.c' <<'END_OF_FILE'
X/*
X * shade.c
X *
X * Copyright (C) 1989, Craig E. Kolb
X *
X * This software may be freely copied, modified, and redistributed,
X * provided that this copyright notice is preserved on all copies.
X *
X * There is no warranty or other guarantee of fitness for this software,
X * it is provided solely .  Bug reports or fixes may be sent
X * to the author, who may or may not act on them as he desires.
X *
X * You may not include this software in a program or other software product
X * without supplying the source, or without informing the end-user that the
X * source is available for no extra charge.
X *
X * If you modify this software, you should include a notice giving the
X * name of the person performing the modification, the date of modification,
X * and the reason for such modification.
X *
X * $Id: shade.c,v 3.0 89/10/27 02:06:03 craig Exp $
X *
X * $Log:    shade.c,v $
X * Revision 3.0  89/10/27  02:06:03  craig
X * Baseline for first official release.
X *
X */
X#include <math.h>
X#include <stdio.h>
X#include "constants.h"
X#include "typedefs.h"
X#include "funcdefs.h"
X
Xint    level, maxlevel;    /* Current tree depth, max depth */
Xdouble    DefIndex = 1.0;        /* Default index of refraction. */
Xdouble    TreeCutoff = UNSET;    /* Minimum contribution of any ray. */
X
X/*
X * Calculate color of ray.
X */
XShadeRay(hitinfo, ray, dist, back, color, contrib)
XHitInfo *hitinfo;        /* Information about point of intersection. */
XRay *ray;            /* Direction and origin of ray. */
Xdouble dist;            /* Distance from origin of intersection. */
XColor *back;            /* "Background" color */
XColor *color;            /* Color to assign current ray. */
Xdouble contrib;            /* Contribution of this ray to final color */
X{
X    Vector hit;
X    extern unsigned long HitRays;
X    extern Fog *GlobalFog;
X    extern Mist *GlobalMist;
X
X    /*
X     * If we got here, then a ray hit something, so...
X     */
X    HitRays++;
X
X    (void)normalize(&hitinfo->norm);
X    /*
X      * "hit" is the location of intersection in world space.
X     * hitinfo->pos is the intersection point in object space.
X     */
X    addscaledvec(ray->pos, dist, ray->dir, &hit);
X    /*
X     * Calculate ray color.
X     */
X    shade(&hit, ray, &hitinfo->norm, hitinfo->prim, &hitinfo->surf,
X            back, color, contrib);
X    /*
X     * If fog or mist is present, modify computed color.
X     */
X    if (GlobalFog)
X        ComputeFog(GlobalFog, dist, color);
X    if (GlobalMist)
X        ComputeMist(GlobalMist, &ray->pos, &hit, dist, color);
X}
X
Xshade(pos, ray, nrm, prim, surf, back, color, contrib)
XVector *pos, *nrm;
XRay *ray;
XPrimitive *prim;
XSurface *surf;
XColor *back, *color;
Xdouble contrib;
X{
X    int lnum, entering;
X    double dist, k;
X    Color newcol;
X    Ray NewRay;
X    HitInfo hitinfo;
X    Vector refl;
X    Light *lp;
X    extern int nlight;
X    extern Light light[];
X    extern unsigned long ReflectRays, RefractRays;
X    extern double TraceRay();
X
X    /*
X     * Ambient color is always included.
X     */
X    *color = surf->amb;
X
X    /*
X     * Calculate direction of reflected ray.
X     */
X    k = -dotp(&ray->dir, nrm);
X    addscaledvec(ray->dir, 2.*k, *nrm, &refl);
X
X    /*
X     * Calculate intensity contributed by each light source.
X     */
X    for(lnum = 0, lp = light;lnum < nlight; lnum++, lp++) {
X        if (lp->type == EXTENDED)
X            extended_lightray(pos, nrm, &refl, lp,
X                      prim, surf, color);
X        else
X            generic_lightray(pos, nrm, &refl, lp,
X                     prim, surf, color);
X    }
X
X
X    if (level >= maxlevel)
X        /*
X         * Don't spawn any refracted/reflected rays.
X         */
X        return;
X    /*
X     * Specular reflection.
X     */
X    if(surf->refl > 0. && contrib * surf->refl > TreeCutoff) {
X        level++;
X        NewRay.shadow = FALSE;
X        NewRay.pos = *pos;        /* Origin == hit point */
X        NewRay.dir = refl;        /* Direction == reflection */
X        NewRay.media = ray->media;    /* Medium == old medium */
X        ReflectRays++;
X        dist = TraceRay(prim, &NewRay, &hitinfo);
X        if (dist > EPSILON) {
X            ShadeRay(&hitinfo, &NewRay, dist, back, &newcol,
X                    contrib*surf->refl);
X            AddScaledColor(*color, surf->refl, newcol, color);
X        } else
X            AddScaledColor(*color, surf->refl, *back, color);
X        level--;
X    }
X    /*
X     * Specular transmission (refraction).
X     */
X    if(surf->transp > 0. && contrib * surf->transp > TreeCutoff) {
X        NewRay.shadow = FALSE;
X        NewRay.pos = *pos;        /* Origin == hit point */
X        NewRay.media = ray->media;    /* Media == old media */
X        if (k < 0.) {
X            /*
X             * Normal points "away" from incoming ray.
X             * Hit "inside" surface -- assume we're exiting.
X             * Pop medium from stack.
X             */
X            if (NewRay.media == (SurfaceList *)0)
X                /*
X                 * We had a funky intersection at some
X                 * point -- e.g. we hit at the intersection
X                 * of two refracting surfaces. Skip it.
X                 */
X                return;
X            free((char *)NewRay.media);
X            NewRay.media = NewRay.media->next;
X            if (refract(&NewRay.dir, surf->kref,
X                NewRay.media ? NewRay.media->surf->kref :
X                DefIndex, ray->dir, *nrm, k))
X                return;
X            entering = FALSE;
X        } else {
X            /*
X             * Entering surface.
X             */
X            if (refract(&NewRay.dir,
X                NewRay.media ? NewRay.media->surf->kref :
X                DefIndex, surf->kref, ray->dir, *nrm, k))
X                return;
X            NewRay.media = add_surface(surf, NewRay.media);
X            entering = TRUE;
X        }
X        level++;
X        RefractRays++;
X        dist = TraceRay((Primitive *)NULL, &NewRay, &hitinfo);
X        if (dist > EPSILON) {
X            ShadeRay(&hitinfo, &NewRay, dist, back, &newcol,
X                contrib * surf->transp);
X            AddScaledColor(*color, surf->transp, newcol, color);
X        } else
X            AddScaledColor(*color, surf->transp, *back, color);
X        if (entering)
X            free((char *)NewRay.media);
X        level--;
X    }
X}
X
X/*
X * Sample an extended (area) light source.
X */
Xextended_lightray(pos, norm, refl, lp, obj, surf, color)
XVector *pos, *norm, *refl;
XLight *lp;
XPrimitive *obj;
XSurface *surf;
XColor *color;
X{
X    int uSample, vSample;
X    double jit, vbase, ubase, vpos, upos;
X    Color newcol;
X    Vector Uaxis, Vaxis, toLight, SampleDir;
X    extern double lightdist, JitterWeight, SampleSpacing;
X    extern int JitSamples, Jittered, SampleNumber;
X
X    /*
X     * Determinte two orthoganal vectors which line in the plane
X     * whose normal is defined by the vector from the center
X     * of the light source to the point of intersection and
X     * which passes through the center of the light source.
X     */
X    LightCoordSys(lp, pos, &toLight, &Uaxis, &Vaxis);
X    jit = 2. * lp->radius * SampleSpacing;
X
X    if (Jittered) {
X        /*
X         * Sample a single point, determined by SampleNumber,
X         * on the extended source.
X         */
X        vpos = -lp->radius + (SampleNumber % JitSamples) * jit;
X        upos = -lp->radius + (SampleNumber / JitSamples) * jit;
X        vpos += nrand() * jit;
X        upos += nrand() * jit;
X        veccomb(upos, Uaxis, vpos, Vaxis, &SampleDir);
X        vecadd(toLight, SampleDir, &SampleDir);
X        /*
X         * Lightdist, the distance to the light source, is
X         * used by inshadow(), called from Lighting().
X         */
X        lightdist = normalize(&SampleDir);
X        Lighting(pos,norm,refl,lp,obj,surf,&SampleDir,color);
X        return;
X    }
X
X    newcol.r = newcol.g = newcol.b = 0.;
X
X    /*
X     * Sample JitSamples^2 -4 points arranged in a square lying on the
X     * plane calculated above.  The size of the square is equal to
X     * the diameter of the light source.  We sample the square at equal
X     * intervals in the U and V direction, with "jitter" thrown in to mask
X     * aliasing.  The corners of the square are skipped to speed up
X     * the calculation and to more closely model a circular source.
X     */
X    ubase = -lp->radius;
X    for(uSample = 1;uSample <= JitSamples;uSample++, ubase += jit) {
X        vbase = -lp->radius;
X        for(vSample=1;vSample <= JitSamples;vSample++,vbase += jit) {
X            /*
X             * Skip corners.
X             */
X            if ((uSample == 1 || uSample == JitSamples) &&
X                (vSample == 1 || vSample == JitSamples))
X                continue;
X            vpos = vbase + nrand() * jit;
X            upos = ubase + nrand() * jit;
X            veccomb(upos, Uaxis, vpos, Vaxis, &SampleDir);
X            vecadd(toLight, SampleDir, &SampleDir);
X            lightdist = normalize(&SampleDir);
X            Lighting(pos, norm, refl, lp, obj, surf,
X                &SampleDir, &newcol);
X        }
X    }
X
X    AddScaledColor(*color, JitterWeight, newcol, color);
X}
X
X/*
X * Lighting calculations for a point or directional light source.
X */
Xgeneric_lightray(pos, norm, reflect, lp, obj, surf, color)
XVector *pos, *norm, *reflect;
XLight *lp;
XPrimitive *obj;
XSurface *surf;
XColor *color;
X{
X    Vector Nlightray;
X
X    lightray(lp, pos, &Nlightray);
X
X    Lighting(pos, norm, reflect, lp, obj, surf, &Nlightray, color);
X}
X
X/*
X * Compute shading function (diffuse reflection and specular highlight)
X * given the point of intersection, incident ray, normal to the surface,
X * direction of the reflected ray,
X * the current light source, the object hit, the surface hit, and
X * the ray to the current light source.
X *
X * This function *adds* the computed color to "color".
X */
XLighting(pos, norm, refl, lp, obj, surf, tolight, color)
XVector *pos, *norm, *refl, *tolight;
XLight *lp;
XPrimitive *obj;
XSurface *surf;
XColor *color;
X{
X    Color bright;
X    double intens;
X
X    /*
X     * Diffuse reflection.
X     * Falls off as the cosine of the angle between
X     * the normal and the ray to the light.
X     */
X    intens = dotp(norm, tolight);
X    if(intens <= 0.) {
X        /*
X         * Object does not face light source.
X         * If the surface is translucent, add in
X         * diffuse and specular components due to
X         * light hitting the other side of the surface.
X         * (This is a "poor man's" diffuse transmission
X         * and specularly transmitted highlights.  Note
X         * that we ignore index of refraction here...)
X         */
X        if (surf->translucency == 0.)
X            return;
X        if (inshadow(&bright, obj, lp, pos, tolight))
X            return;
X        /*
X         * We use -translucency to counter the fact
X         * that intens is < 0.
X         */
X        intens *= -surf->translucency;
X        color->r += surf->diff.r * intens * bright.r;
X        color->g += surf->diff.g * intens * bright.g;
X        color->b += surf->diff.b * intens * bright.b;
X        intens = -dotp(refl, tolight);
X        if (intens <= 0.)
X            return;
X        intens = surf->translucency * pow(intens, surf->stcoef);
X        color->r += surf->spec.r * intens * bright.r;
X        color->g += surf->spec.g * intens * bright.g;
X        color->b += surf->spec.b * intens * bright.b;
X        return;
X    }
X    /*
X     * Add diffuse reflection component.
X     */
X    if (inshadow(&bright, obj, lp, pos, tolight))
X        return;        /* Object in shadow. */
X    color->r += surf->diff.r * intens * bright.r;
X    color->g += surf->diff.g * intens * bright.g;
X    color->b += surf->diff.b * intens * bright.b;
X    /*
X     * Specularly reflected highlights.
X     * Fall off as the cosine of the angle
X     * between the reflected ray and the ray to the light source.
X     */
X    if (surf->coef < EPSILON)
X        return;
X    intens = dotp(refl, tolight);
X    if(intens <= 0.)
X        return;
X    /*
X     * Specular highlight = cosine of the angle raised to the
X     * appropriate power.
X     */
X    intens = pow(intens, surf->coef);
X    color->r += surf->spec.r * intens * bright.r;
X    color->g += surf->spec.g * intens * bright.g;
X    color->b += surf->spec.b * intens * bright.b;
X}
END_OF_FILE
if test 10618 -ne `wc -c <'src/shade.c'`; then
    echo shar: \"'src/shade.c'\" unpacked with wrong size!
fi
# end of 'src/shade.c'
fi
echo shar: End of archive 5 \(of 8\).
cp /dev/null ark5isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 8 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0