Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/misc/volume21/rayshade/part17

Text File | 1991-07-21 | 39.2 KB | 1,690 lines

Newsgroups: comp.sources.misc
From: Rayshade Construction Co. <[e-mail]>
Subject:  v21i020:  rayshade - A raytracing package for UNIX, Part17/19
Message-ID: <[e-mail]>
X-Md4-Signature: a4d125964e11af669e118c42f54e0e56
Date: Sun, 21 Jul 1991 22:34:16 GMT
Approved: [e-mail]

Submitted-by: Rayshade Construction Co. <[e-mail]>
Posting-number: Volume 21, Issue 20
Archive-name: rayshade/part17
Environment: UNIX, !16BIT

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 17 (of 19)."
# Contents:  libray/libobj/blob.c rayview/glmethods.c
# Wrapped by kolb@woody on Wed Jul 17 17:57:00 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'libray/libobj/blob.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'libray/libobj/blob.c'\"
else
echo shar: Extracting \"'libray/libobj/blob.c'\" \(18207 characters\)
sed "s/^X//" >'libray/libobj/blob.c' <<'END_OF_FILE'
X/*
X * blob.c
X *
X * Copyright (C) 1990, 1991, Mark Podlipec, Craig E. Kolb
X * All rights reserved.
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: blob.c,v 4.0 91/07/17 14:36:07 kolb Exp Locker: kolb $
X *
X * $Log:    blob.c,v $
X * Revision 4.0  91/07/17  14:36:07  kolb
X * Initial version.
X *
X */
X#include "geom.h"
X#include "blob.h"
X
Xstatic Methods *iBlobMethods = NULL;
Xstatic char blobName[] = "blob";
X
Xunsigned long BlobTests, BlobHits;
X
X/*
X * Blob/Metaball Description
X *
X * In this implementation a Blob is made up of a threshold T, and a
X * group of 1 or more Metaballs.  Each metaball 'i'  is defined by
X * its position (xi,yi,zi), its radius ri, and its strength si.
X * Around each Metaball is a density function di(Ri) defined by:
X *
X *         di(Ri) =  c4i * Ri^4  +  c2i * Ri^2  + c0i     0 <= Ri <= ri
X *         di(Ri) =  0                                    ri < Ri
X *
X * where
X *       c4i  = si / (ri^4)
X *       c2i  = -(2 * si) / (ri^2)
X *       c0i  = si
X *       Ri^2 is the distance from a point (x,y,z) to the center of
X *            the Metaball - (xi,yi,zi).
X *
X * The density function looks sort of like a W (Float-u) with the
X * center hump crossing the d-axis at  si  and the two humps on the
X * side being tangent to the R-axis at  -ri  and  +ri. Only the
X * range  [0,ri] is being used.
X * I chose this function so that derivative = 0  at the points  0 and +ri.
X * This keeps the surface smooth when the densities are added. I also
X * wanted no  Ri^3  and  Ri  terms, since the equations are easier to
X * solve without them. This led to the symmetry about the d-axis and
X * the derivative being equal to zero at -ri as well.
X *
X * The surface of the Blob is defined by:
X *
X *
X *                  N
X *                 ---
X *      F(x,y,z) = \    di(Ri)  = T
X *                 /
X *                 ---
X *                 i=0
X *
X * where
X *
X *     di(Ri)    is   given above
X *     Ri^2      =    (x-xi)^2  +  (y-yi)^2  + (z-zi)^2
X *      N        is   number of Metaballs in Blob
X *      T        is   the threshold
X *    (xi,yi,zi) is   the center of Metaball i
X *
X */
X
X/*****************************************************************************
X * Create & return reference to a metaball.
X */
XBlob *
XBlobCreate(T, mlist, npoints)
XFloat T;
XMetaList *mlist;
Xint npoints;
X{
X    Blob *blob;
X    int i;
X    MetaList *cur;
X
X/*
X * There has to be at least one metaball in the blob.
X * Note: if there's only one metaball, the blob
X * will be a sphere.
X */
X    if (npoints < 1)
X    {
X        RLerror(RL_WARN, "blob field not correct.\n");
X        return (Blob *)NULL;
X    }
X
X/*
X * Initialize primitive and Geom structures
X */
X    blob = (Blob *)Malloc(sizeof(Blob));
X    blob->T = T;
X    blob->list=(MetaVector *)
X        Malloc( (unsigned)(npoints*sizeof(MetaVector)) );
X    blob->num = npoints;
X
X/*
X * Initialize Metaball list
X */
X    for(i=0;i<npoints;i++)
X    {
X        cur = mlist;
X        if ( (cur->mvec.c0 < EPSILON) || (cur->mvec.rs < EPSILON) )
X        {
X            RLerror(RL_WARN,"Degenerate metaball in blob (sr).\n");
X            return (Blob *)NULL;
X        }
X        /* store radius squared */
X        blob->list[i].rs = cur->mvec.rs * cur->mvec.rs;
X        /* Calculate and store coefficients for each metaball */
X        blob->list[i].c0 = cur->mvec.c0;
X        blob->list[i].c2 = -(2.0 * cur->mvec.c0) / blob->list[i].rs;
X        blob->list[i].c4 = cur->mvec.c0 /
X                (blob->list[i].rs * blob->list[i].rs);
X        blob->list[i].x = cur->mvec.x;
X        blob->list[i].y = cur->mvec.y;
X        blob->list[i].z = cur->mvec.z;
X        mlist = mlist->next;
X        free((voidstar)cur);
X    }
X    /*
X     * Allocate room for Intersection Structures and
X     * Allocate room for an array of pointers to point to
X     * the Intersection Structures.
X     */
X    blob->ilist = (MetaInt *)Malloc( 2 * blob->num * sizeof(MetaInt));
X    blob->iarr = (MetaInt **)Malloc( 2 * blob->num * sizeof(MetaInt *));
X    return blob;
X}
X
XMethods *
XBlobMethods()
X{
X    if (iBlobMethods == (Methods *)NULL) {
X        iBlobMethods = MethodsCreate();
X        iBlobMethods->create = (GeomCreateFunc *)BlobCreate;
X        iBlobMethods->methods = BlobMethods;
X        iBlobMethods->name = BlobName;
X        iBlobMethods->intersect = BlobIntersect;
X        iBlobMethods->normal = BlobNormal;
X        iBlobMethods->bounds = BlobBounds;
X        iBlobMethods->stats = BlobStats;
X        iBlobMethods->checkbounds = TRUE;
X        iBlobMethods->closed = TRUE;
X    }
X    return iBlobMethods;
X}
X
X/*****************************************************************************
X * Function used by qsort() when sorting the Ray/Blob intersection list
X */
Xint
XMetaCompare(A,B)
Xchar *A,*B;
X{
X    MetaInt **AA,**BB;
X
X    AA = (MetaInt **) A;
X    BB = (MetaInt **) B;
X    if (AA[0]->bound == BB[0]->bound) return(0);
X    if (AA[0]->bound <  BB[0]->bound) return(-1);
X    return(1);  /* AA[0]->bound is > BB[0]->bound */
X}
X
X/*****************************************************************************
X * Ray/metaball intersection test.
X */
Xint
XBlobIntersect(blob, ray, mindist, maxdist)
XBlob *blob;
XRay *ray;
XFloat mindist, *maxdist;
X{
X    double c[5], s[4];
X    Float dist;
X    MetaInt *ilist,**iarr;
X    register int i,j,inum;
X    extern void qsort();
X
X    BlobTests++;
X
X    ilist = blob->ilist;
X    iarr = blob->iarr;
X
X/*
X * The first step in calculating the Ray/Blob intersection is to
X * divide the Ray into intervals such that only a fixed set of
X * Metaballs contribute to the density function along that interval.
X * This is done by finding the set of intersections between the Ray
X * and each Metaball's Sphere/Region of influence, which has a
X * radius  ri  and is centered at (xi,yi,zi).
X * Intersection information is kept track of in the MetaInt
X * structure and consists of:
X *
X *   type    indicates whether this intersection is the start(R_START)
X *           of a Region or the end(R_END) of one.
X *   pnt     the Metaball of this intersection
X *   bound   the distance from Ray origin to this intersection
X *
X * This list is then sorted by  bound  and used later to find the Ray/Blob
X * intersection.
X */
X
X    inum = 0;
X    for(i=0; i < blob->num; i++)
X    {
X        register Float xadj, yadj, zadj;
X        register Float b, t, rs;
X        register Float dmin,dmax;
X
X        rs   = blob->list[i].rs;
X        xadj = blob->list[i].x - ray->pos.x;
X        yadj = blob->list[i].y - ray->pos.y;
X        zadj = blob->list[i].z - ray->pos.z;
X
X    /*
X     * Ray/Sphere of Influence intersection
X     */
X        b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
X        t = b * b - xadj * xadj - yadj * yadj - zadj * zadj + rs;
X
X    /*
X     * don't except imaginary or single roots. A single root is a ray
X     * tangent to the Metaball's Sphere/Region. The Metaball's
X     * contribution to the overall density function at this point is
X     * zero anyway.
X     */
X        if (t > 0.0) /* we have two roots */
X        {
X            t = sqrt(t);
X            dmin = b - t;
X    /*
X     * only interested in stuff in front of ray origin
X     */
X            if (dmin < mindist) dmin = mindist;
X            dmax = b + t;
X            if (dmax > dmin) /* we have a valid Region */
X            {
X        /*
X             * Initialize min/start and max/end Intersections Structures
X             * for this Metaball
X             */
X                ilist[inum].type = R_START;
X                ilist[inum].pnt = i;
X                ilist[inum].bound = dmin;
X                for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
X                iarr[inum] = &(ilist[inum]);
X                inum++;
X
X                ilist[inum].type = R_END;
X                ilist[inum].pnt = i;
X                ilist[inum].bound = dmax;
X                for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
X                iarr[inum] = &(ilist[inum]);
X                inum++;
X            } /* End of valid Region */
X        } /* End of two roots */
X    } /* End of loop through metaballs */
X
X    /*
X     * If there are no Ray/Metaball intersections there will
X     * not be a Ray/Blob intersection. Exit now.
X     */
X    if (inum == 0)
X    {
X        return FALSE;
X    }
X
X    /*
X     * Sort Intersection list. No sense using qsort if there's only
X     * two intersections.
X     *
X     * Note: we actually aren't sorting the Intersection structures, but
X     * an array of pointers to the Intersection structures.
X     * This is faster than sorting the Intersection structures themselves.
X     */
X    if (inum==2)
X    {
X        MetaInt *t;
X        if (ilist[0].bound > ilist[1].bound)
X        {
X            t = iarr[0];
X            iarr[0] = iarr[1];
X            iarr[1] = t;
X        }
X    }
X    else qsort((voidstar)(iarr), (unsigned)inum, sizeof(MetaInt *),
X            MetaCompare);
X
X/*
X* Finding the Ray/Blob Intersection
X*
X* The non-zero part of the density function for each Metaball is
X*
X*   di(Ri) = c4i * Ri^4  +  c2i * Ri^2  +  c0i
X*
X* To find find the Ray/Blob intersection for one metaball
X* substitute for distance
X*
X*     Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2
X*
X* and then substitute the Ray equations:
X*
X*     x  = x0 + x1 * t
X*     y  = y0 + y1 * t
X*     z  = z0 + z1 * t
X*
X* to get a big mess :^). Actually, it's a Quartic in t and it's fully
X* listed farther down. Here's a short version:
X*
X*   c[4] * t^4  +  c[3] * t^3  +  c[2] * t^2  +  c[1] * t  +  c[0]  =  T
X*
X* Don't forget that the Blob is defined by the density being equal to
X* the threshold T.
X* To calculate the intersection of a Ray and two or more Metaballs,
X* the coefficients are calculated for each Metaball and then added
X* together. We can do this since we're working with polynomials.
X* The points of intersection are the roots of the resultant equation.
X*
X* The algorithm loops through the intersection list, calculating
X* the coefficients if an intersection is the start of a Region and
X* adding them to all intersections in that region.
X* When it detects a valid interval, it calculates the
X* roots from the starting intersection's coefficients and if any of
X* the roots are in the interval, the smallest one is returned.
X*
X*/
X
X    {
X        register Float *tmpc;
X        MetaInt *strt,*tmp;
X        register int istrt,itmp;
X        register int num,exitflag,inside;
X
X    /*
X     * Start the algorithm outside the first interval
X     */
X        inside = 0;
X        istrt = 0;
X        strt = iarr[istrt];
X        if (strt->type != R_START)
X            RLerror(RL_WARN,"MetaInt sanity check FAILED!\n");
X
X        /*
X         * Loop through intersection. If a root is found the code
X         * will return at that point.
X         */
X        while( istrt < inum )
X        {
X            /*
X             * Check for multiple intersections  at the same point.
X             * This is also where the coefficients are calculated
X             * and spread throughout that Metaball's sphere of
X             * influence.
X             */
X            do
X            {
X                /* find out which metaball */
X                i = strt->pnt;
X                /* only at starting boundaries do this */
X                if (strt->type == R_START)
X                {
X                    register MetaVector *ml;
X                    register Float a1,a0;
X                    register Float xd,yd,zd;
X                    register Float c4,c2,c0;
X
X                    /* we're inside */
X                    inside++;
X
X    /*  As promised, the full equations
X     *
X     *   c[4] = c4*a2*a2;
X     *   c[3] = 4.0*c4*a1*a2;
X     *   c[2] = 4.0*c4*a1*a1 + 2.0*c4*a2*a0 + c2*a2;
X     *   c[1] = 4.0*c4*a1*a0 + 2.0*c2*a1;
X     *   c[0] = c4*a0*a0 + c2*a0 + c0;
X     *
X     * where
X     *        a2 = (x1*x1 + y1*y1 + z1*z1) = 1.0 because the ray
X     *                                           is normalized
X     *        a1 = (xd*x1 + yd*y1 + zd*z1)
X     *        a0 = (xd*xd + yd*yd + zd*zd)
X     *        xd = (x0 - xi)
X     *        yd = (y0 - yi)
X     *        zd = (z0 - zi)
X     *        (xi,yi,zi) is center of Metaball
X     *        (x0,y0,z0) is Ray origin
X     *        (x1,y1,z1) is normalized Ray direction
X     *        c4,c2,c0   are the coefficients for the
X     *                       Metaball's density function
X     *
X     */
X
X                    /*
X                     * Some compilers would recalculate
X                     * this each time its used.
X                     */
X                    ml = &(blob->list[i]);
X
X                    xd = ray->pos.x - ml->x;
X                    yd = ray->pos.y - ml->y;
X                    zd = ray->pos.z - ml->z;
X                    a1 = (xd * ray->dir.x + yd * ray->dir.y
X                        + zd * ray->dir.z);
X                    a0 = (xd * xd + yd * yd + zd * zd);
X
X                    c0 = ml->c0;
X                    c2 = ml->c2;
X                    c4 = ml->c4;
X
X                    c[4] = c4;
X                    c[3] = 4.0*c4*a1;
X                    c[2] = 2.0*c4*(2.0*a1*a1 + a0) + c2;
X                    c[1] = 2.0*a1*(2.0*c4*a0 + c2);
X                    c[0] = c4*a0*a0 + c2*a0 + c0;
X
X        /*
X         * Since we've gone through the trouble of calculating the
X         * coefficients, put them where they'll be used.
X         * Starting at the current intersection(It's also the start of
X         * a region), add the coefficients to each intersection until
X         * we reach the end of this region.
X         */
X                    itmp = istrt;
X                    tmp = strt;
X                    while( (tmp->pnt != i) ||
X                           (tmp->type != R_END) )
X                    {
X                        tmpc = tmp->c;
X                        for(j=0;j<5;j++)
X                            *tmpc++ += c[j];
X                        itmp++;
X                        tmp = iarr[itmp];
X                    }
X
X                } /* End of start of a Region */
X                /*
X                 * Since the intersection wasn't the start
X                 * of a region, it must the end of one.
X                 */
X                else inside--;
X
X        /*
X         * If we are inside a region(or regions) and the next
X         * intersection is not at the same place, then we have
X         * the start of an interval. Set the exitflag.
X         */
X                if ((inside > 0)
X                    && (strt->bound != iarr[istrt+1]->bound) )
X                    exitflag = 1;
X                else
X                /* move to next intersection */
X                {
X                    istrt++;
X                    strt = iarr[istrt];
X                    exitflag = 0;
X                }
X        /*
X         * Check to see if we're at the end. If so then exit with
X         * no intersection found.
X         */
X                if (istrt >= inum)
X                {
X                    return FALSE;
X                }
X            } while(!exitflag);
X                /* End of Search for valid interval */
X
X        /*
X         * Find Roots along this interval
X         */
X
X            /* get coefficients from Intersection structure */
X            tmpc = strt->c;
X            for(j=0;j<5;j++) c[j] = *tmpc++;
X
X            /* Don't forget to put in threshold */
X            c[0] -= blob->T;
X
X            /* Use Graphic Gem's Quartic Root Finder */
X            num = SolveQuartic(c,s);
X
X        /*
X         * If there are roots, search for roots within the interval.
X         */
X            dist = 0.0;
X            if (num>0)
X            {
X                for(j=0;j<num;j++)
X                {
X        /*
X         * Not sure if EPSILON is truly needed, but it might cause
X         * small cracks between intervals in some cases. In any case
X         * I don't believe it hurts.
X         */
X                    if ((s[j] >= (strt->bound - EPSILON))
X                        && (s[j] <= (iarr[istrt+1]->bound +
X                            EPSILON) ) )
X                    {
X                        if (dist == 0.0)
X                        /* first valid root */
X                            dist = s[j];
X                        else if (s[j] < dist)
X                         /* else only if smaller */
X                            dist = s[j];
X                    }
X                }
X                /*
X                 * Found a valid root
X                 */
X                if (dist > mindist && dist < *maxdist)
X                {
X                    *maxdist = dist;
X                    BlobHits++;
X                    return TRUE;
X                    /* Yeah! Return valid root */
X                }
X            }
X
X        /*
X         * Move to next intersection
X         */
X            istrt++;
X            strt = iarr[istrt];
X
X        } /* End of loop through Intersection List */
X    } /* End of Solving for Ray/Blob Intersection */
X
X    /*
X     * return negative
X     */
X    return FALSE;
X}
X
X
X/***********************************************
X * Find the Normal of a Blob at a given point
X *
X * The normal of a surface at a point (x0,y0,z0) is the gradient of that
X * surface at that point. The gradient is the vector
X *
X *            Fx(x0,y0,z0) , Fy(x0,y0,z0) , Fz(x0,y0,z0)
X *
X * where Fx(),Fy(),Fz() are the partial dervitives of F() with respect
X * to x,y and z respectively. Since the surface of a Blob is made up
X * of the sum of one or more polynomials, the normal of a Blob at a point
X * is the sum of the gradients of the individual polynomials at that point.
X * The individual polynomials in this case are di(Ri) i = 0,...,num .
X *
X * To find the gradient of the contributing polynomials
X * take di(Ri) and substitute U = Ri^2;
X *
X *      di(U)    = c4i * U^2  +  c2i * U  + c0i
X *
X *      dix(U)   = 2 * c4i * Ux * U  +  c2i * Ux
X *
X *        U      = (x-xi)^2 + (y-yi)^2 + (z-zi)^2
X *
X *        Ux     = 2 * (x-xi)
X *
X * Rearranging we get
X *
X *    dix(x0,y0,z0) = 4 * c4i * xdi * disti + 2 * c2i * xdi
X *    diy(x0,y0,z0) = 4 * c4i * ydi * disti + 2 * c2i * ydi
X *    diz(x0,y0,z0) = 4 * c4i * zdi * disti + 2 * c2i * zdi
X *
X * where
X *         xdi       =   x0-xi
X *         ydi       =   y0-yi
X *         zdi       =   y0-yi
X *        disti      =   xdi*xdi + ydi*ydi + zdi*zdi
X *       (xi,yi,zi)  is  the center of Metaball i
X *       c4i,c2i,c0i are the coefficients of Metaball i's density function
X */
Xint
XBlobNormal(blob, pos, nrm, gnrm)
XBlob *blob;
XVector *pos, *nrm, *gnrm;
X{
X    register int i;
X
X    /*
X     * Initialize normals to zero
X     */
X    nrm->x = nrm->y = nrm->z = 0.0;
X    /*
X     * Loop through Metaballs. If the point is within a Metaball's
X     * Sphere of influence, calculate the gradient and add it to the
X     * normals
X     */
X    for(i=0;i < blob->num; i++)
X    {
X        register MetaVector *sl;
X        register Float dist,xd,yd,zd;
X
X        sl = &(blob->list[i]);
X        xd = pos->x - sl->x;
X        yd = pos->y - sl->y;
X        zd = pos->z - sl->z;
X
X        dist  = xd*xd + yd*yd + zd*zd;
X        if (dist <= sl->rs )
X        {
X            register Float temp;
X
X            /* temp is negative so normal points out of blob */
X            temp = -2.0 * (2.0 * sl->c4 * dist  +  sl->c2);
X            nrm->x += xd * temp;
X            nrm->y += yd * temp;
X            nrm->z += zd * temp;
X
X        }
X    }
X    (void)VecNormalize(nrm);
X    *gnrm = *nrm;
X    return FALSE;
X}
X
X
X/*****************************************************************************
X * Calculate the extent of the Blob
X */
Xvoid
XBlobBounds(blob, bounds)
XBlob *blob;
XFloat bounds[2][3];
X{
X    Float r;
X    Float minx,miny,minz,maxx,maxy,maxz;
X    int i;
X
X    /*
X     * Loop through Metaballs to find the minimun and maximum in each
X     * direction.
X     */
X    for( i=0;  i< blob->num;  i++)
X    {
X        register Float d;
X
X        r = sqrt(blob->list[i].rs);
X        if (i==0)
X        {
X            minx = blob->list[i].x - r;
X            miny = blob->list[i].y - r;
X            minz = blob->list[i].z - r;
X            maxx = blob->list[i].x + r;
X            maxy = blob->list[i].y + r;
X            maxz = blob->list[i].z + r;
X        }
X        else
X        {
X            d = blob->list[i].x - r;
X            if (d<minx) minx = d;
X            d = blob->list[i].y - r;
X            if (d<miny) miny = d;
X            d = blob->list[i].z - r;
X            if (d<minz) minz = d;
X            d = blob->list[i].x + r;
X            if (d>maxx) maxx = d;
X            d = blob->list[i].y + r;
X            if (d>maxy) maxy = d;
X            d = blob->list[i].z + r;
X            if (d>maxz) maxz = d;
X        }
X    }
X    bounds[LOW][X]  = minx;
X    bounds[HIGH][X] = maxx;
X    bounds[LOW][Y]  = miny;
X    bounds[HIGH][Y] = maxy;
X    bounds[LOW][Z]  = minz;
X    bounds[HIGH][Z] = maxz;
X}
X
Xchar *
XBlobName()
X{
X    return blobName;
X}
X
Xvoid
XBlobStats(tests, hits)
Xunsigned long *tests, *hits;
X{
X    *tests = BlobTests;
X    *hits = BlobHits;
X}
X
Xvoid
XBlobMethodRegister(meth)
XUserMethodType meth;
X{
X    if (iBlobMethods)
X        iBlobMethods->user = meth;
X}
END_OF_FILE
if test 18207 -ne `wc -c <'libray/libobj/blob.c'`; then
    echo shar: \"'libray/libobj/blob.c'\" unpacked with wrong size!
fi
# end of 'libray/libobj/blob.c'
fi
if test -f 'rayview/glmethods.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'rayview/glmethods.c'\"
else
echo shar: Extracting \"'rayview/glmethods.c'\" \(17863 characters\)
sed "s/^X//" >'rayview/glmethods.c' <<'END_OF_FILE'
X/*
X * glmethods.c
X *
X * Support routines for SGI/RS6000 machines.
X *
X * Copyright (C) 1989, 1991, Craig E. Kolb, Allan Snyder
X *
X * This software may be freely copied, modified, and redistributed
X * provided that this copyright notice is preserved on all copies.
X *
X * You may not distribute this software, in whole or in part, as part of
X * any commercial product without the express consent of the authors.
X *
X * There is no warranty or other guarantee of fitness of this software
X * for any purpose.  It is provided solely "as is".
X *
X * $Id: glmethods.c,v 4.0 91/07/17 17:38:43 kolb Exp Locker: kolb $
X *
X * $Log:    glmethods.c,v $
X * Revision 4.0  91/07/17  17:38:43  kolb
X * Initial version
X *
X */
X#include <gl.h>
X#include <device.h>
X#include "rayshade.h"
X#include "options.h"
X#include "viewing.h"
X
X#include "libsurf/surface.h"
X
X#include "libobj/box.h"
X#include "libobj/cone.h"
X#include "libobj/csg.h"
X#include "libobj/cylinder.h"
X#include "libobj/disc.h"
X#include "libobj/grid.h"
X#include "libobj/instance.h"
X#include "libobj/list.h"
X#include "libobj/plane.h"
X#include "libobj/poly.h"
X#include "libobj/sphere.h"
X#include "libobj/triangle.h"
X
X#include "liblight/light.h"
X#include "liblight/extended.h"
X#include "liblight/infinite.h"
X#include "liblight/point.h"
X#include "liblight/spot.h"
X
X#define CIRCLE_SAMPLES    20    /* # of samples around circle (cone/cyl/etc) */
X#define DEFAULT_NEAR    .2    /* default value for near clipping plane */
X#define DEFAULT_FAR    350.    /* default far clipping plane */
X#define PLANE_RAD    450.    /* radius of disc used to represent planes */
X
X#ifndef SPHERELIB
X#define SPHERE_LEVEL    3
X#endif
X
X/*
X * Pass a normal stored in a Vector to the geometry pipeline.
X */
X#define GLNormal(w)    (glnrm[0] = (w)->x, glnrm[1] = (w)->y, \
X             glnrm[2] = (w)->z, n3f(glnrm))
X
Xstatic void    GLBoxDraw(), GLConeDraw(), GLCsgDraw(), GLCylinderDraw(),
X        GLDiscDraw(),
X        GLGridDraw(), GLHfDraw(), GLInstanceDraw(),
X        GLListDraw(), GLPlaneDraw(), GLPolygonDraw(),
X        GLSphereDraw(), GLTorusDraw(), GLTriangleDraw(),
X        GLBoundsDraw(),
X        GLInfiniteLight(), GLPointLight(), GLSpotLight(),
X        GLExtendedLight();
X
Xstatic short    cursurf = 1,
X        curlight = 1;
X
Xstatic Object    GLBoxObject,
X        GLCylObject,
X        GLSphereObject,
X        GLBoxObjectDefine(),
X        GLCylObjectDefine();
X
Xextern Object    GLSphereObjectDefine();
X
Xstatic int    Doublebuffered;
X
Xstatic float    Ident[4][4] =    {1., 0., 0., 0.,
X                 0., 1., 0., 0.,
X                 0., 0., 1., 0.,
X                 0., 0., 0., 1.},
X        glnrm[3];
X
Xstatic float **CirclePoints;
X
Xstatic void LightDraw(), GeomDraw(), GLMultMatrix(),
X        GLPushSurface(), GLPopSurface(), LightDrawInit(),
X        ObjectInit(), GLDrawFrame(), ScreenDrawInit(), DrawInit(),
X        GLUnitCircle(), GLDisc(), ComputeClippingPlanes();
X
XFloat CurrentTime;
X
Xstatic unsigned long BackPack;    /* Packed background color */
Xstatic long Zinit;        /* maximum Zbuffer value */
X
XMethodsRegister()
X{
X    BoxMethodRegister(GLBoxDraw);
X    ConeMethodRegister(GLConeDraw);
X    CsgMethodRegister(GLCsgDraw);
X    CylinderMethodRegister(GLCylinderDraw);
X    DiscMethodRegister(GLDiscDraw);
X    GridMethodRegister(GLGridDraw);
X    /*HfMethodRegister(GLHfDraw);*/
X    InstanceMethodRegister(GLInstanceDraw);
X    ListMethodRegister(GLListDraw);
X    PlaneMethodRegister(GLPlaneDraw);
X    PolygonMethodRegister(GLPolygonDraw);
X    SphereMethodRegister(GLSphereDraw);
X    /* Have to write sphere routine for RS6000... */
X    /*TorusMethodRegister(GLTorusDraw);*/
X    TriangleMethodRegister(GLTriangleDraw);
X
X    InfiniteMethodRegister(GLInfiniteLight);
X    PointMethodRegister(GLPointLight);
X    ExtendedMethodRegister(GLExtendedLight);
X    SpotMethodRegister(GLSpotLight);
X}
X
Xvoid
XRender()
X{
X    short val;
X    float tmp;
X
X    DrawInit();
X    qdevice(ESCKEY);
X    qdevice(SPACEKEY);
X    qdevice(REDRAW);
X
X    DrawFrames();
X
X    while (TRUE) {
X        switch (qread(&val)) {
X            case ESCKEY:
X                exit(0);
X                break;
X            case REDRAW:
X                reshapeviewport();
X                DrawFrames();
X                break;
X            case SPACEKEY:
X                if (!val)
X                    DrawFrames();
X                break;
X        }
X    }
X}
X
XDrawFrames()
X{
X    int i;
X    for (i = 0; i < Options.totalframes; i++)
X        GLDrawFrame(i);
X}
X
Xstatic void
XDrawInit()
X{
X    extern Surface DefaultSurface;
X
X    ScreenDrawInit();
X    ObjectInit();
X    LightDrawInit();
X    /*
X     * Push the default surface.
X     */
X    GLPushSurface(&DefaultSurface);
X}
X
Xstatic void
XScreenDrawInit()
X{
X    /*
X     * Open window, initialize graphics, etc.
X     */
X#ifdef sgi
X    foreground();
X#endif
X    prefsize(Screen.xsize, Screen.ysize);
X    winopen("rayview");
X    RGBmode();
X    mmode(MVIEWING);
X
X    if (Options.totalframes > 1) {
X        Doublebuffered = TRUE;
X        doublebuffer();
X    }
X
X    gconfig();
X    /*
X     * Initialize viewing matrix.
X     */
X    GLViewingInit();
X
X    zbuffer(TRUE);
X
X    BackPack = (unsigned char)(255*Screen.background.r) |
X        ((unsigned char)(255*Screen.background.g) << 8) |
X        ((unsigned char)(255*Screen.background.b) << 16);
X    Zinit = getgdesc(GD_ZMAX);
X    lsetdepth(0, Zinit);
X}
X
X/*
X * Draw the World object.
X */
Xstatic void
XGLDrawFrame(frame)
Xint frame;
X{
X    extern Geom *World;
X
X    RSStartFrame(frame);
X    CurrentTime = Options.framestart;
X    TimeSet(CurrentTime);
X
X    czclear(BackPack, Zinit);
X
X    /*
X     * Draw the World object
X     */
X    GeomDraw(World);
X    if (Doublebuffered)
X        swapbuffers();
X}
X
XGLViewingInit()
X{
X    float near, far, aspect, dist, T[4][4];
X    short ang;
X    extern Geom *World;
X
X    ang = Camera.vfov * 10 + 0.5;
X    aspect = Camera.hfov / Camera.vfov;
X
X    T[0][0]=Screen.scrni.x; T[0][1]=Screen.scrnj.x; T[0][2]= -Camera.dir.x;
X    T[1][0]=Screen.scrni.y; T[1][1]=Screen.scrnj.y; T[1][2]= -Camera.dir.y;
X    T[2][0]=Screen.scrni.z; T[2][1]=Screen.scrnj.z; T[2][2]= -Camera.dir.z;
X
X    T[0][3] = T[1][3] = T[2][3] = 0.;
X
X    T[3][0] = -dotp(&Camera.lookp, &Screen.scrni);
X    T[3][1] = -dotp(&Camera.lookp, &Screen.scrnj);
X    T[3][2] = dotp(&Camera.lookp, &Camera.dir);
X    T[3][3] = 1.;
X
X    ComputeClippingPlanes(&near, &far, World->bounds);
X
X    loadmatrix(Ident);
X    perspective(ang, aspect, near, far);
X    polarview((float)Camera.lookdist, 0, 0, 0);
X    multmatrix(T);
X}
X
X
Xstatic void
XObjectInit()
X{
X    CircleDefine();
X    GLBoxObject = GLBoxObjectDefine();
X
X#ifdef SPHERELIB
X    GLSphereObject = GLSphereObjectDefine();
X#else
X    GLSphereObject = GLSphereObjectDefine(SPHERE_LEVEL);
X#endif
X
X    GLCylObject = GLCylObjectDefine();
X}
X
Xstatic void
XGeomDraw(obj)
XGeom *obj;
X{
X    Trans *ct;
X    /*
X     * If object has a surface associated with it,
X     * start using it.
X     */
X    if (obj->surf)
X        GLPushSurface(obj->surf);
X    if (obj->trans) {
X        /*
X         * Take care of any animated transformations that
X         * exist.
X         */
X        if (obj->animtrans && !equal(obj->timenow, CurrentTime)) {
X            TransResolveAssoc(obj->trans);
X            obj->timenow = CurrentTime;
X        }
X        pushmatrix();
X        /*
X          * Apply in reverse order.
X         */
X        for (ct = obj->transtail; ct; ct = ct->prev)
X            GLMultMatrix(&ct->trans);
X    }
X
X    if (obj->methods->user) {
X        /*
X          * Call proper method
X          */
X        (*obj->methods->user)(obj->obj);
X    } else {
X        /*
X         * Draw bounding box.
X         */
X        GLBoundsDraw(obj->bounds);
X    }
X
X    if (obj->surf)
X        GLPopSurface();
X    if (obj->trans)
X        popmatrix();
X}
X
Xstatic float surfprops[] =    {AMBIENT, 0, 0, 0,
X                 DIFFUSE, 0, 0, 0,
X                 SPECULAR, 0, 0, 0,
X                 SHININESS, 0,
X                 LMNULL};
Xstatic float    *amb = &surfprops[1],
X        *diff = &surfprops[5],
X        *spec = &surfprops[9],
X        *shine = &surfprops[13];
X
Xstatic void
XGLPushSurface(surf)
XSurface *surf;
X{
X    static Surface *lastsurf = NULL;
X
X    /*
X     * Start using the given surface.
X     * In the case of the use of an "applysurf",
X     * the same surface will be applied to many
X     * primitives individually.  By saving the
X     * pointer to the last surface, we keep from
X     * lmdef'ing the surface when we don't need to.
X     */
X    if (surf != lastsurf) {
X        amb[0] = surf->amb.r; amb[1] = surf->amb.g;
X            amb[2] = surf->amb.b;
X        diff[0] = surf->diff.r; diff[1] = surf->diff.g;
X            diff[2] = surf->diff.b;
X        spec[0] = surf->spec.r; spec[1] = surf->spec.g;
X            spec[2] = surf->spec.b;
X        shine[0] = surf->srexp;
X        lmdef(DEFMATERIAL, cursurf, 15, surfprops);
X        lastsurf = surf;
X    }
X    lmbind(MATERIAL, cursurf);
X    cursurf++;
X}
X
Xstatic void
XGLMultMatrix(trans)
XRSMatrix *trans;
X{
X    float newmat[4][4];
X    /*
X     * Multiply in the given transformation.
X     */
X    newmat[0][0] = trans->matrix[0][0]; newmat[0][1] = trans->matrix[0][1];
X    newmat[0][2] = trans->matrix[0][2]; newmat[0][3] = 0.;
X    newmat[1][0] = trans->matrix[1][0]; newmat[1][1] = trans->matrix[1][1];
X    newmat[1][2] = trans->matrix[1][2]; newmat[1][3] = 0.;
X    newmat[2][0] = trans->matrix[2][0]; newmat[2][1] = trans->matrix[2][1];
X    newmat[2][2] = trans->matrix[2][2]; newmat[2][3] = 0.;
X    newmat[3][0] = trans->translate.x; newmat[3][1] = trans->translate.y;
X    newmat[3][2] = trans->translate.z ; newmat[3][3] = 1.;
X    multmatrix(newmat);
X}
X
Xstatic void
XGLPopSurface()
X{
X    cursurf--;
X    lmbind(MATERIAL, cursurf-1);
X}
X
Xstatic void
XGLBoundsDraw(bounds)
XFloat bounds[2][3];
X{
X    float sx, sy, sz;
X
X    pushmatrix();
X
X    translate(bounds[LOW][X], bounds[LOW][Y], bounds[LOW][Z]);
X    scale(    bounds[HIGH][X] - bounds[LOW][X],
X        bounds[HIGH][Y] - bounds[LOW][Y],
X        bounds[HIGH][Z] - bounds[LOW][Z]);
X    pushmatrix();
X    callobj(GLBoxObject);
X    popmatrix();
X    popmatrix();
X}
X
Xstatic void
XGLBoxDraw(box)
XBox *box;
X{
X    GLBoundsDraw(box->bounds);
X}
X
X
Xstatic void
XGLConeDraw(cone)
XCone *cone;
X{
X    int i, j;
X    float norm[3], normconst, ZeroVect[3], tmpv[3];
X
X    ZeroVect[0] = ZeroVect[1] = ZeroVect[2] = 0.;
X    /*
X     * Sides.
X     * For the normal, assume we are finding the normal at
X     * (C_P[i].x, C_P[i].y, 1.) at this point, the unnormalized
X     * normal = (C_P[i].x, C_P[i].y, -tantheta^2).  When we normalize,
X     * then length of the vector is simply equal to:
X     * sqrt(CP[i].x^2 + CP[i].y^2 + tantheta^4)  ==
X     * sqrt(1. + tantheta^4)
X     * In our case, tantheta = 1., so...
X     */
X
X    normconst = 1. / sqrt(2.);
X    norm[2] = -normconst;
X
X    pushmatrix();
X    GLMultMatrix(&cone->trans.trans);
X
X    for (i = 0; i < CIRCLE_SAMPLES; i++) {
X        j = (i + 1) % CIRCLE_SAMPLES;
X        norm[0] = CirclePoints[i][0] * normconst;
X        norm[1] = CirclePoints[i][1] * normconst;
X        n3f(norm);
X        bgnpolygon();
X        if (cone->start_dist > EPSILON) {
X            tmpv[2] = cone->start_dist;
X            tmpv[0] = CirclePoints[i][0] * cone->start_dist;
X            tmpv[1] = CirclePoints[i][1] * cone->start_dist;
X            v3f(tmpv);
X            tmpv[0] = CirclePoints[j][0] * cone->start_dist;
X            tmpv[1] = CirclePoints[j][1] * cone->start_dist;
X            v3f(tmpv);
X        } else
X            v3f(ZeroVect);
X        tmpv[2] = 1.;
X        tmpv[0] = CirclePoints[j][0];
X        tmpv[1] = CirclePoints[j][1];
X        v3f(tmpv);
X        tmpv[0] = CirclePoints[i][0];
X        tmpv[1] = CirclePoints[i][1];
X        v3f(tmpv);
X        endpolygon();
X    }
X    popmatrix();
X}
X
Xstatic void
XGLCsgDraw(csg)
XCsg *csg;
X{
X    /*
X     * Punt by drawing both objects.
X     */
X    GeomDraw(csg->obj1);
X    GeomDraw(csg->obj2);
X}
X
Xstatic void
XGLCylinderDraw(cyl)
XCylinder *cyl;
X{
X    pushmatrix();
X    GLMultMatrix(&cyl->trans.trans);
X    callobj(GLCylObject);
X    popmatrix();
X}
X
Xstatic void
XGLDiscDraw(disc)
XDisc *disc;
X{
X    GLDisc((Float)sqrt(disc->radius), &disc->pos, &disc->norm);
X}
X
Xstatic void
XGLDisc(rad, pos, norm)
XFloat rad;
XVector *pos, *norm;
X{
X    RSMatrix m, tmp;
X    Vector atmp;
X
X    /*
X     * This is kinda disgusting...
X     */
X    ScaleMatrix(rad, rad, 1., &m);
X    if (fabs(norm->z) == 1.) {
X        atmp.x = 1.;
X        atmp.y = atmp.z = 0.;
X    } else {
X        atmp.x = norm->y;
X        atmp.y = -norm->x;
X        atmp.z= 0.;
X    }
X
X    RotationMatrix(atmp.x, atmp.y, atmp.z, -acos(norm->z), &tmp);
X    MatrixMult(&m, &tmp, &m);
X    TranslationMatrix(pos->x, pos->y, pos->z, &tmp);
X    MatrixMult(&m, &tmp, &m);
X    pushmatrix();
X    GLMultMatrix(&m);
X    /*
X     * Draw unit circle.
X     */
X    GLUnitCircle(0., TRUE);
X    popmatrix();
X}
X
Xstatic void
XGLGridDraw(grid)
XGrid *grid;
X{
X    Geom *ltmp;
X
X    for (ltmp = grid->unbounded; ltmp; ltmp = ltmp->next)
X        GeomDraw(ltmp);
X    for (ltmp = grid->objects; ltmp; ltmp = ltmp->next)
X        GeomDraw(ltmp);
X}
X
Xstatic void
XGLHfDraw(){}
X
Xstatic void
XGLInstanceDraw(inst)
XInstance *inst;
X{
X    GeomDraw(inst->obj);
X}
X
Xstatic void
XGLListDraw(list)
XList *list;
X{
X    Geom *ltmp;
X
X    for (ltmp = list->unbounded; ltmp; ltmp = ltmp->next)
X        GeomDraw(ltmp);
X    for (ltmp = list->list; ltmp; ltmp = ltmp->next)
X        GeomDraw(ltmp);
X}
X
Xstatic void
XGLPlaneDraw(plane)
XPlane *plane;
X{
X    /*
X     * Draw a really big disc.
X     */
X    GLDisc((Float)PLANE_RAD, &plane->pos, &plane->norm);
X}
X
Xstatic void
XGLPolygonDraw(poly)
Xregister Polygon *poly;
X{
X    register int i;
X
X    GLNormal(&poly->norm);
X
X    bgnpolygon();
X    for (i = 0; i < poly->npoints; i++)
X        v3d(&poly->points[i]);
X    endpolygon();
X}
X
Xstatic void
XGLSphereDraw(sph)
XSphere *sph;
X{
X    pushmatrix();
X    translate(sph->x, sph->y, sph->z);
X    scale(sph->r, sph->r, sph->r);
X    callobj(GLSphereObject);
X    popmatrix();
X}
X
Xstatic void
XGLTorusDraw(){}
X
Xstatic void
XGLTriangleDraw(tri)
Xregister Triangle *tri;
X{
X    /*
X     * If Float is a double, use v3d,
X     * otherwise use v3f.
X     */
X    if (tri->type != PHONGTRI) {
X        GLNormal(&tri->nrm);
X        bgnpolygon();
X        v3d(&tri->p[0]); v3d(&tri->p[1]); v3d(&tri->p[2]);
X        endpolygon();
X    } else {
X        bgnpolygon();
X        GLNormal(&tri->vnorm[0]);
X        v3d(&tri->p[0]);
X        GLNormal(&tri->vnorm[1]);
X        v3d(&tri->p[1]);
X        GLNormal(&tri->vnorm[2]);
X        v3d(&tri->p[2]);
X        endpolygon();
X    }
X}
X
Xfloat lightprops[] =    {POSITION, 0., 0., 0., 0.,
X             LCOLOR, 0, 0, 0,
X             LMNULL};
X
Xfloat     *lpos = &lightprops[1],
X    *lcolor = &lightprops[6];
X
Xfloat lmodel[] =    {AMBIENT, 1., 1., 1.,
X             ATTENUATION, 1., 0.,
X             LOCALVIEWER, 0.,
X#ifdef sgi
X             ATTENUATION2, 0.,
X             TWOSIDE, 1.,
X#endif
X             LMNULL};
X
Xstatic void
XLightDrawInit()
X{
X    Light *light;
X    extern Light *Lights;
X
X    for (light = Lights; light; light = light->next)
X        LightDraw(light);
X
X    switch (curlight-1) {
X        case 8:
X            lmbind(LIGHT7, 8);
X        case 7:
X            lmbind(LIGHT6, 7);
X        case 6:
X            lmbind(LIGHT5, 6);
X        case 5:
X            lmbind(LIGHT4, 5);
X        case 4:
X            lmbind(LIGHT3, 4);
X        case 3:
X            lmbind(LIGHT2, 3);
X        case 2:
X            lmbind(LIGHT1, 2);
X        case 1:
X            lmbind(LIGHT0, 1);
X    }
X
X    lmodel[1] = Options.ambient.r;
X    lmodel[2] = Options.ambient.g;
X    lmodel[3] = Options.ambient.b;
X
X#ifdef sgi
X    lmdef(DEFLMODEL, 1, 14, lmodel);
X#else
X    lmdef(DEFLMODEL, 1, 10, lmodel);
X#endif
X    lmbind(LMODEL, 1);
X}
X
Xstatic void
XLightDraw(light)
XLight *light;
X{
X    if (!light || !light->methods->user)
X        return;
X    lcolor[0] = light->color.r;
X    lcolor[1] = light->color.g;
X    lcolor[2] = light->color.b;
X    (*light->methods->user)(light->light);
X}
X
Xstatic void
XGLExtendedLight(ext)
XExtended *ext;
X{
X    lpos[0] = ext->pos.x;
X    lpos[1] = ext->pos.y;
X    lpos[2] = ext->pos.z;
X    lpos[3] = 1.;
X    lmdef(DEFLIGHT, curlight++, 10, lightprops);
X}
X
X
Xstatic void
XGLInfiniteLight(inf)
XInfinite *inf;
X{
X    lpos[0] = inf->dir.x;
X    lpos[1] = inf->dir.y;
X    lpos[2] = inf->dir.z;
X    lpos[3] = 0.;
X    lmdef(DEFLIGHT, curlight++, 10, lightprops);
X}
X
Xstatic void
XGLPointLight(pt)
XPointlight *pt;
X{
X    lpos[0] = pt->pos.x;
X    lpos[1] = pt->pos.y;
X    lpos[2] = pt->pos.z;
X    lpos[3] = 1.;
X    lmdef(DEFLIGHT, curlight++, 10, lightprops);
X}
X
Xstatic void
XGLSpotLight(spot)
XSpotlight *spot;
X{
X}
X
Xstatic float boxfaces[6][4][3] = {
X{     {1., 1., 1.},
X    {0., 1., 1.},
X    {0., 0., 1.},
X    {1., 0., 1.}    },
X{    {1., 0., 1.},
X    {0., 0., 1.},
X    {0., 0., 0.},
X    {1., 0., 0.}    },
X{    {1., 0., 0.},
X    {0., 0., 0.},
X    {0., 1., 0.},
X    {1., 1., 0.}    },
X{    {0., 1., 0.},
X    {0., 1., 1.},
X    {1., 1., 1.},
X    {1., 1., 0.}    },
X{    {1., 1., 1.},
X    {1., 0., 1.},
X    {1., 0., 0.},
X    {1., 1., 0.}    },
X{    {0., 0., 1.},
X    {0., 1., 1.},
X    {0., 1., 0.},
X    {0., 0., 0.}}};
X
Xfloat boxnorms[6][3] = {
X    {0, 0, 1},
X    {0, -1, 0},
X    {0, 0, -1},
X    {0, 1, 0},
X    {1, 0, 0,},
X    {-1, 0, 0}};
X
X/*
X * Define a unit cube centered at (0.5, 0.5, 0.5)
X */
Xstatic Object
XGLBoxObjectDefine()
X{
X    int i, box;
X
X    makeobj(box = genobj());
X    for (i = 0; i < 6; i++)    {
X        n3f(boxnorms[i]);
X        polf(4, boxfaces[i]);
X    }
X    closeobj();
X    return box;
X}
X
Xstatic void
XGLUnitCircle(z, up)
Xfloat z;
Xint up;
X{
X    int i;
X    float norm[3];
X
X    norm[0] = norm[1] = 0.;
X    bgnpolygon();
X
X    if (up) {
X        norm[2] = 1.;
X        n3f(norm);
X        for (i = 0; i < CIRCLE_SAMPLES; i++) {
X            CirclePoints[i][2] = z;
X            v3f(CirclePoints[i]);
X        }
X    } else {
X        norm[2] = -1.;
X        n3f(norm);
X        for (i = CIRCLE_SAMPLES -1; i; i--) {
X            CirclePoints[i][2] = z;
X            v3f(CirclePoints[i]);
X        }
X    }
X
X    endpolygon();
X}
X
Xstatic Object
XGLCylObjectDefine()
X{
X    int i, j, cyl;
X    float norm[3];
X
X    makeobj(cyl = genobj());
X    norm[2] = 0;
X    for (i = 0; i < CIRCLE_SAMPLES; i++) {
X        j = (i+1)%CIRCLE_SAMPLES;
X        norm[0] = CirclePoints[i][0];
X        norm[1] = CirclePoints[i][1];
X#ifdef sgi
X        n3f(norm);
X        bgnpolygon();
X        CirclePoints[i][2] = 0;
X        v3f(CirclePoints[i]);
X        CirclePoints[j][2] = 0;
X        v3f(CirclePoints[j]);
X        CirclePoints[j][2] = 1;
X        v3f(CirclePoints[j]);
X        CirclePoints[i][2] = 1;
X        v3f(CirclePoints[i]);
X        endpolygon();
X#else
X        n3f(norm);
X        pmv(CirclePoints[i][0], CirclePoints[i][1], 0.);
X        pdr(CirclePoints[j][0], CirclePoints[j][1], 0.);
X        pdr(CirclePoints[j][0], CirclePoints[j][1], 1.);
X        pdr(CirclePoints[i][0], CirclePoints[i][1], 1.);
X        pclos();
X#endif
X    }
X    closeobj();
X    return cyl;
X}
X
X/*
X * Fill CirclePoints[t] with X and Y values cooresponding to points
X * along the unit circle.  The size of CirclePoints is equal to
X * CIRCLE_SAMPLES.
X */
XCircleDefine()
X{
X    int i;
X    double theta, dtheta;
X
X    dtheta = 2.*M_PI / (double)CIRCLE_SAMPLES;
X    CirclePoints = (float **)malloc(CIRCLE_SAMPLES * sizeof(float *));
X    for (i = 0, theta = 0; i < CIRCLE_SAMPLES; i++, theta += dtheta) {
X        CirclePoints[i] = (float *)malloc(3 * sizeof(float));
X        CirclePoints[i][0] = cos(theta);
X        CirclePoints[i][1] = sin(theta);
X        /* Z is left unset. */
X    }
X}
X
X/*
X * Given world bounds, determine near and far
X * clipping planes.  Only problem occurs when there are
X * unbbounded objects (planes)...
X */
Xstatic void
XComputeClippingPlanes(near, far, bounds)
Xfloat *near, *far;
XFloat bounds[2][3];
X{
X    Vector e, c;
X    double rad, d;
X
X
X    /*
X     * Compute 'radius' of scene.
X     */
X    rad = max(max((bounds[HIGH][X] - bounds[LOW][X])*0.5,
X              (bounds[HIGH][Y] - bounds[LOW][Y])*0.5),
X              (bounds[HIGH][Z] - bounds[LOW][Z])*0.5);
X
X    c.x = (bounds[LOW][X] + bounds[HIGH][X])*0.5;
X    c.y = (bounds[LOW][Y] + bounds[HIGH][Y])*0.5;
X    c.z = (bounds[LOW][Z] + bounds[HIGH][Z])*0.5;
X
X    /*
X     * Compute position of eye relative to the center of the sphere.
X     */
X    VecSub(Camera.pos, c, &e);
X    d = VecNormalize(&e);
X    if (d <= rad)
X        /*
X         * Eye is inside sphere.
X         */
X        *near = DEFAULT_NEAR;
X    else
X        *near = (d - rad)*0.2;
X    *far = (d + rad)*2.;
X}
END_OF_FILE
if test 17863 -ne `wc -c <'rayview/glmethods.c'`; then
    echo shar: \"'rayview/glmethods.c'\" unpacked with wrong size!
fi
# end of 'rayview/glmethods.c'
fi
echo shar: End of archive 17 \(of 19\).
cp /dev/null ark17isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 19 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0

exit 0 # Just in case...