v21i013: A ray tracing program, Part06/08
Rich Salz
rsalz at uunet.uu.net
Thu Feb 8 07:50:17 AEST 1990
Submitted-by: Craig Kolb <craig at weedeater.math.yale.edu>
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
--
Please send comp.sources.unix-related mail to rsalz at uunet.uu.net.
Use a domain-based address or give alternate paths, or you may lose out.
More information about the Comp.sources.unix
mailing list