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

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

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

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

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