craig@weedeater.math.yale.edu (Craig Kolb) (10/15/90)
Posting-number: Volume 15, Issue 25
Submitted-by: Craig Kolb <craig@weedeater.math.yale.edu>
Archive-name: ggems/part03
#! /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 3 (of 5)."
# Contents: AALines/utah.c Albers.c BoundSphere.c Dissolve.c
# DoubleLine.c GraphicsGems.h LineEdge.c MatrixInvert.c
# PolyScan/poly_clip.c PolyScan/poly_scan.c Roots3And4.c
# Wrapped by craig@weedeater on Fri Oct 12 15:53:13 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'AALines/utah.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'AALines/utah.c'\"
else
echo shar: Extracting \"'AALines/utah.c'\" \(4532 characters\)
sed "s/^X//" >'AALines/utah.c' <<'END_OF_FILE'
X/*
X file: utah.c
X description: interface to Utah RLE toolkit
X author: A. T. Campbell
X date: October 27, 1989
X*/
X
X#ifndef lint
Xstatic char sccsid[] = "%W% %G%"; /* SCCS info */
X#endif lint
X
X#include <math.h>
X#include <stdio.h>
X#ifdef sequent
X#include <strings.h>
X#else
X#include <string.h>
X#endif
X#include "utah.h"
X
X/******************************************************************************/
X
X/* return values */
Xextern void free();
Xextern char *malloc();
X
X/******************************************************************************/
X
Xutah_read_close(ufp)
XUTAH_FILE *ufp;
X{
X return(0);
X}
X
X/******************************************************************************/
X
XUTAH_FILE *
Xutah_read_init(fname, ht, wd)
X
Xchar *fname;
Xint *ht, *wd;
X{
X FILE *fp;
X UTAH_FILE *ufp;
X
X /* open output stream */
X if (!strcmp(fname, ""))
X fp = stdin;
X else {
X if ((fp = fopen(fname, "r")) == NULL)
X return(NULL);
X }
X
X /* change the default sv_globals struct to match what we need */
X ufp = (UTAH_FILE *) malloc(sizeof(UTAH_FILE));
X *ufp = sv_globals;
X ufp->svfb_fd = fp;
X
X /* read the header in the input file */
X if (rle_get_setup(ufp) != 0)
X return(NULL);
X
X /* get image size */
X *wd = ufp->sv_xmax - ufp->sv_xmin + 1;
X *ht = ufp->sv_ymax - ufp->sv_ymin + 1;
X
X /* normal termination */
X return(ufp);
X}
X
X/******************************************************************************/
X
Xutah_read_pixels(ufp, pixels)
X
XUTAH_FILE *ufp;
Xunsigned char pixels[][3];
X{
X static unsigned n = 0;
X static unsigned char *r = NULL, *g = NULL, *b = NULL;
X int i, width;
X
X /* allocate storage */
X width = ufp->sv_xmax + 1;
X if (width > n) {
X if (n > 0) {
X free((char *)r);
X free((char *)g);
X free((char *)b);
X }
X n = width;
X r = (unsigned char *) malloc(n * sizeof(unsigned char));
X g = (unsigned char *) malloc(n * sizeof(unsigned char));
X b = (unsigned char *) malloc(n * sizeof(unsigned char));
X }
X
X /* read this row */
X utah_read_rgb(ufp, r, g, b);
X
X /* convert to pixels */
X for (i = 0; i < width; i++) {
X pixels[i][0] = r[i];
X pixels[i][1] = g[i];
X pixels[i][2] = b[i];
X }
X
X return(0);
X}
X
X/******************************************************************************/
X
Xutah_read_rgb(ufp, r, g, b)
X
XUTAH_FILE *ufp;
Xunsigned char r[], g[], b[];
X{
X rle_pixel *rows[3];
X
X /* set color channels */
X rows[0] = r;
X rows[1] = g;
X rows[2] = b;
X
X /* read this row */
X rle_getrow(ufp, rows);
X return(0);
X}
X
X/******************************************************************************/
X
Xutah_write_close(ufp)
X
XUTAH_FILE *ufp;
X{
X if (!ufp) return(-1);
X sv_puteof(ufp);
X return(0);
X}
X
X/******************************************************************************/
X
XUTAH_FILE *
Xutah_write_init(fname, ht, wd)
X
Xchar *fname;
Xint ht, wd;
X{
X FILE *fp;
X UTAH_FILE *ufp;
X
X /* open output stream */
X if (!strcmp(fname, ""))
X fp = stdout;
X else {
X if ((fp = fopen(fname, "w")) == NULL)
X return(NULL);
X }
X
X /* change the default sv_globals struct to match what we need */
X ufp = (UTAH_FILE *) malloc(sizeof(UTAH_FILE));
X *ufp = sv_globals;
X ufp->svfb_fd = fp;
X ufp->sv_xmax = wd - 1;
X ufp->sv_ymax = ht - 1;
X ufp->sv_alpha = 0; /* No coverage (alpha) */
X
X /* create the header in the output file */
X sv_setup(RUN_DISPATCH, ufp);
X
X /* normal termination */
X return(ufp);
X}
X
X/******************************************************************************/
X
Xutah_write_pixels(ufp, pixels)
X
XUTAH_FILE *ufp;
Xunsigned char pixels[][3];
X{
X static unsigned n = 0;
X static unsigned char *r = NULL, *g = NULL, *b = NULL;
X int i, width;
X
X /* allocate storage */
X width = ufp->sv_xmax + 1;
X if (width > n) {
X if (n > 0) {
X free((char *)r);
X free((char *)g);
X free((char *)b);
X }
X n = width;
X r = (unsigned char *) malloc(n * sizeof(unsigned char));
X g = (unsigned char *) malloc(n * sizeof(unsigned char));
X b = (unsigned char *) malloc(n * sizeof(unsigned char));
X }
X
X /* convert to color channels */
X for (i = 0; i < width; i++) {
X r[i] = pixels[i][0];
X g[i] = pixels[i][1];
X b[i] = pixels[i][2];
X }
X
X /* write this row */
X utah_write_rgb(ufp, r, g, b);
X return(0);
X}
X
X/******************************************************************************/
X
Xutah_write_rgb(ufp, r, g, b)
X
XUTAH_FILE *ufp;
Xunsigned char r[], g[], b[];
X{
X rle_pixel *rows[3];
X int width;
X
X /* set color channels */
X rows[0] = r;
X rows[1] = g;
X rows[2] = b;
X
X /* write this row */
X width = ufp->sv_xmax - ufp->sv_xmin + 1;
X sv_putrow(rows, width, ufp);
X return(0);
X}
X
X/******************************************************************************/
END_OF_FILE
if test 4532 -ne `wc -c <'AALines/utah.c'`; then
echo shar: \"'AALines/utah.c'\" unpacked with wrong size!
fi
# end of 'AALines/utah.c'
fi
if test -f 'Albers.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Albers.c'\"
else
echo shar: Extracting \"'Albers.c'\" \(3994 characters\)
sed "s/^X//" >'Albers.c' <<'END_OF_FILE'
X/*
XAlbers Equal-Area Conic Map Projection
Xby Paul Bame
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/*
X * Albers Conic Equal-Area Projection
X * Formulae taken from "Map Projections Used by the U.S.
X * Geological Survey" Bulletin #1532
X *
X * Equation reference numbers and some variable names taken
X * from the reference.
X * To use, call albers setup() once and then albers_project()
X * for each coordinate pair of interest.
X*/
X
X#include "GraphicsGems.h"
X#include <stdio.h>
X#include <math.h>
X
X/*
X * This is the Clarke 1866 Earth spheroid data which is often
X * used by the USGS - there are other spheroids however - see the
X * book.
X */
X
X/*
X * Earth radii in different units */
X#define CONST_EradiusKM (6378.2064) /* Kilometers */
X#define CONST_EradiusMI (CONST_EradiusKM/1.609) /* Miles */
X#define CONST_Ec (0.082271854) /* Eccentricity */
X#define CONST_Ecsq (0.006768658) /* Eccentricity squared */
X
X
X
X/*
X * To keep things simple, assume Earth radius is 1.0. Projected
X * coordinates (X and Y obtained from albers project ()) are
X * dimensionless and may be multiplied by any desired radius
X * to convert to desired units (usually Kilometers or Miles).
X*/
X#define CONST_Eradius 1.0
X
X/* Pre-computed variables */
Xstatic double middlelon; /* longitude at center of map */
Xstatic double bigC, cone_const, r0; /* See the reference */
X
Xstatic
Xcalc_q_msq(lat, qp, msqp)
Xdouble lat;
Xdouble *qp;
Xdouble *msqp;
X/*
X * Given latitude, calculate 'q' [eq 3-12]
X * if msqp is != NULL, m^2 [eq. 12-15].
X*/
X{
X double s, c, es;
X
X s = sin(lat);
X es = s * CONST_Ec;
X
X *qp = (1.0 - CONST_Ecsq) * ((s / (1 - es * es))-
X (1 / (2 * CONST_Ec)) * log((1 - es) / (1 + es)));
X
X if (msqp != NULL)
X {
X c = cos(lat);
X *msqp = c * c/ (1 - es * es);
X }
X}
X
X
X
X
Xalbers_setup(southlat, northlat, originlon, originlat)
Xdouble southlat, northlat;
Xdouble originlon;
Xdouble originlat;
X/*
X * Pre-compute a bunch of variables which are used by the
X * albers_project()
X *
X * southlat Southern latitude for Albers projection
X * northlat Northern latitute for Albers projection
X * originlon Longitude for origin of projected map
X * originlat Latitude for origin of projected map -
X * often (northlat + southlat) / 2
X*/
X{
X double q1, q2, q0;
X double m1sq, m2sq;
X
X middlelon = originlon;
X
X cal_q_msq(southlat, &q1, &m1sq);
X cal_q_msq(northlat, &q2, &m2sq);
X cal_q_msq(originlat, &q0, NULL);
X
X cone_const = (m1sq - m2sq) / (q2 - q1);
X bigC = m1sq + cone_const * q1;
X r0 = CONST_Eradius * sqrt(bigC - cone_const *q0) / cone_const;
X}
X
X/***************************************************************/
X
Xalbers_project(lon, lat, xp, yp)
Xdouble lon, lat;
Xdouble *xp, *yp;
X/*
X * Project lon/lat (in radians) according to albers_setup and
X * return the results via xp, yp. Units of *xp and *yp are same
X * as the units used by CONST_Eradius.
X*/
X{
X double q, r, theta;
X calc_q_msq(lat, &q, NULL);
X theta = cone_const * (lon -middlelon);
X r = CONST_Eradius * sqrt(bigC - cone_const * q) / cone_const;
X *xp = r * sin(theta);
X *yp = r0 - r * cos(theta);
X}
X
X#ifdef TESTPROGRAM
X
X/*
X * Test value from the USGS book. Because of roundoff, the
X * actual values are printed for visual inspection rather
X * than guessing what constitutes "close enough".
X*/
X/* Convert a degress, minutes pair to radians */
X#define DEG_MIN2RAD(degrees, minutes) \
X ((double) ((degrees + minutes / 60.0) * M_PI / 180.0))
X
X#define Nlat DEG_MIN2RAD(29, 30) /* 29 degrees, 30' North Lat */
X#define Slat DEG_MIN2RAD(45, 30)
X#define Originlat DEG_MIN2RAD(23, 0)
X#define Originlon DEG_MIN2RAD(-96, 0) /* West longitude is negative */
X
X#define Testlat DEG_MIN2RAD(35, 0)
X#define Testlon DEG_MIN2RAD(-75, 0)
X
X#define TestX 1885.4727
X#define TestY 1535.9250
X
Xmain()
X{
X int i;
X double x, y;
X
X/* Setup is also from USGS book test set */
X albers_setup(Slat, Nlat, Originlon, Originlat);
X
X albers_project(Testlon, Testlat, &x, &y);
X printf("%.41f, %.41f =?= %.41f, %.41f/n",
X x * CONST_EradiusKM, y * CONST_EradiusKM,
X TestX, TestY);
X}
X#endif /* TESTPROGRAM */
END_OF_FILE
if test 3994 -ne `wc -c <'Albers.c'`; then
echo shar: \"'Albers.c'\" unpacked with wrong size!
fi
# end of 'Albers.c'
fi
if test -f 'BoundSphere.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'BoundSphere.c'\"
else
echo shar: Extracting \"'BoundSphere.c'\" \(3767 characters\)
sed "s/^X//" >'BoundSphere.c' <<'END_OF_FILE'
X/*
XAn Efficient Bounding Sphere
Xby Jack Ritter
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/* Routine to calculate tight bounding sphere over */
X/* a set of points in 3D */
X/* This contains the routine find_bounding_sphere(), */
X/* the struct definition, and the globals used for parameters. */
X/* The abs() of all coordinates must be < BIGNUMBER */
X/* Code written by Jack Ritter and Lyle Rains. */
X
X#include <stdio.h>
X#include <math.h>
X#include "GraphicsGems.h"
X
X#define BIGNUMBER 100000000.0 /* hundred million */
X
X/* GLOBALS. These are used as input and output parameters. */
X
Xstruct Point3Struct caller_p,cen;
Xdouble rad;
Xint NUM_POINTS;
X
X/* Call with no parameters. Caller must set NUM_POINTS */
X/* before calling. */
X/* Caller must supply the routine GET_iTH_POINT(i), which loads his */
X/* ith point into the global struct caller_p. (0 <= i < NUM_POINTS). */
X/* The calling order of the points is irrelevant. */
X/* The final bounding sphere will be put into the globals */
X/* cen and rad. */
X
X
Xfind_bounding_sphere()
X{
Xregister int i;
Xregister double dx,dy,dz;
Xregister double rad_sq,xspan,yspan,zspan,maxspan;
Xdouble old_to_p,old_to_p_sq,old_to_new;
Xstruct Point3Struct xmin,xmax,ymin,ymax,zmin,zmax,dia1,dia2;
X
X
X/* FIRST PASS: find 6 minima/maxima points */
Xxmin.x=ymin.y=zmin.z= BIGNUMBER; /* initialize for min/max compare */
Xxmax.x=ymax.y=zmax.z= -BIGNUMBER;
Xfor (i=0;i<NUM_POINTS;i++)
X {
X GET_iTH_POINT(i); /* load global struct caller_p with */
X /* his ith point. */
X if (caller_p.x<xmin.x)
X xmin = caller_p; /* New xminimum point */
X if (caller_p.x>xmax.x)
X xmax = caller_p;
X if (caller_p.y<ymin.y)
X ymin = caller_p;
X if (caller_p.y>ymax.y)
X ymax = caller_p;
X if (caller_p.z<zmin.z)
X zmin = caller_p;
X if (caller_p.z>zmax.z)
X zmax = caller_p;
X }
X/* Set xspan = distance between the 2 points xmin & xmax (squared) */
Xdx = xmax.x - xmin.x;
Xdy = xmax.y - xmin.y;
Xdz = xmax.z - xmin.z;
Xxspan = dx*dx + dy*dy + dz*dz;
X
X/* Same for y & z spans */
Xdx = ymax.x - ymin.x;
Xdy = ymax.y - ymin.y;
Xdz = ymax.z - ymin.z;
Xyspan = dx*dx + dy*dy + dz*dz;
X
Xdx = zmax.x - zmin.x;
Xdy = zmax.y - zmin.y;
Xdz = zmax.z - zmin.z;
Xzspan = dx*dx + dy*dy + dz*dz;
X
X/* Set points dia1 & dia2 to the maximally seperated pair */
Xdia1 = xmin; dia2 = xmax; /* assume xspan biggest */
Xmaxspan = xspan;
Xif (yspan>maxspan)
X {
X maxspan = yspan;
X dia1 = ymin; dia2 = ymax;
X }
Xif (zspan>maxspan)
X {
X dia1 = zmin; dia2 = zmax;
X }
X
X
X/* dia1,dia2 is a diameter of initial sphere */
X/* calc initial center */
Xcen.x = (dia1.x+dia2.x)/2.0;
Xcen.y = (dia1.y+dia2.y)/2.0;
Xcen.z = (dia1.z+dia2.z)/2.0;
X/* calculate initial radius**2 and radius */
Xdx = dia2.x-cen.x; /* x componant of radius vector */
Xdy = dia2.y-cen.y; /* y componant of radius vector */
Xdz = dia2.z-cen.z; /* z componant of radius vector */
Xrad_sq = dx*dx + dy*dy + dz*dz;
Xrad = sqrt(rad_sq);
X
X/* SECOND PASS: increment current sphere */
X
Xfor (i=0;i<NUM_POINTS;i++)
X {
X GET_iTH_POINT(i); /* load global struct caller_p */
X /* with his ith point. */
X dx = caller_p.x-cen.x;
X dy = caller_p.y-cen.y;
X dz = caller_p.z-cen.z;
X old_to_p_sq = dx*dx + dy*dy + dz*dz;
X if (old_to_p_sq > rad_sq) /* do r**2 test first */
X { /* this point is outside of current sphere */
X old_to_p = sqrt(old_to_p_sq);
X /* calc radius of new sphere */
X rad = (rad + old_to_p) / 2.0;
X rad_sq = rad*rad; /* for next r**2 compare */
X old_to_new = old_to_p - rad;
X /* calc center of new sphere */
X cen.x = (rad*cen.x + old_to_new*caller_p.x) / old_to_p;
X cen.y = (rad*cen.y + old_to_new*caller_p.y) / old_to_p;
X cen.z = (rad*cen.z + old_to_new*caller_p.z) / old_to_p;
X /* Suppress if desired */
X printf("\n New sphere: cen,rad = %f %f %f %f",
X cen.x,cen.y,cen.z, rad);
X }
X }
X} /* end of find_bounding_sphere() */
END_OF_FILE
if test 3767 -ne `wc -c <'BoundSphere.c'`; then
echo shar: \"'BoundSphere.c'\" unpacked with wrong size!
fi
# end of 'BoundSphere.c'
fi
if test -f 'Dissolve.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Dissolve.c'\"
else
echo shar: Extracting \"'Dissolve.c'\" \(4596 characters\)
sed "s/^X//" >'Dissolve.c' <<'END_OF_FILE'
X/* A Digital Dissolve Effect
Xby Mike Morton
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/*
X * Code fragment to advance from one element to the next.
X *
X * int reg; /* current sequence element
X * reg = 1; /* start in any non-zero state
X * if (reg & 1) /* is the bottom bit set?
X * reg = (reg >>1) ^ MASK; /* yes: toss out 1 bit; XOR in mask
X * else reg = reg >>1; /* no: toss out 0 bit
X */
X
Xint randmasks[32]; /* Gotta fill this in yourself. */
X
Xdissolve1 (height, width) /* first version of the dissolve /* algorithm */
X int height, width; /* number of rows, columns */
X{
X int pixels, lastnum; /* number of pixels; */
X /* last pixel's number */
X int regwidth; /* "width" of sequence generator */
X register long mask; /* mask to XOR with to*/
X /* create sequence */
X register unsigned long element;
X /* one element of random sequence */
X register int row, column;
X /* row and column numbers for a pixel */
X
X /* Find smallest register which produces enough pixel numbers */
X pixels = height * width; /* compute number of pixels */
X /* to dissolve */
X lastnum = pixels-1; /* find last element (they go 0..lastnum) */
X regwidth = bitwidth (lastnum); /* how wide must the */
X /* register be? */
X mask = randmasks [regwidth]; /* which mask is for that width? */
X
X /* Now cycle through all sequence elements. */
X
X element = 1; /* 1st element (could be any nonzero) */
X
X
X do {
X row = element / width; /* how many rows down is this pixel? */
X column = element % width; /* and how many columns across? */
X if (row < height) /* is this seq element in the array? */
X copy (row, column); /* yes: copy the (r,c)'th pixel */
X
X /* Compute the next sequence element */
X if (element & 1) /* is the low bit set? */
X element = (element >>1)^mask; /* yes: shift value, */
X /* XOR in mask */
X else element = (element >>1); /* no: just shift the value */
X } while (element != 1); /* loop until we return */
X /* to original element */
X copy (0, 0); /* kludge: the loop doesn't produce (0,0) */
X} /* end of dissolve1() */
X
X
X
Xint bitwidth (N) /* find "bit-width" needed to represent N */
X unsigned int N; /* number to compute the width of */
X{
X int width = 0; /* initially, no bits needed to represent N */
X while (N != 0) { /* loop 'til N has been whittled down to 0 */
X N >>= 1; /* shift N right 1 bit (NB: N is unsigned) */
X width++; /* and remember how wide N is */
X } /* end of loop shrinking N down to nothing *
X return (width); /* return bit positions counted */
X
X} /* end of bitwidth() */
X
X
X
Xdissolve2 (height, width) /* fast version of the dissolve algorithm */
X int height, width; /* number of rows, columns */
X{
X int rwidth, cwidth; /* bit width for rows, for columns */
X int regwidth; /* "width" of sequence generator */
X register long mask; /* mask to XOR with to create sequence */
X register int rowshift; /* shift distance to get row */
X /* from element */
X register int colmask; /* mask to extract column from element */
X register unsigned long element; /* one element of random */ /* sequence */
X register int row, column; /* row and column for one pixel */
X
X
X /* Find the mask to produce all rows and columns. */
X
X rwidth = bitwidth (height); /* how many bits needed for height? */
X cwidth = bitwidth (width); /* how many bits needed for width? */
X regwidth = rwidth + cwidth; /* how wide must the register be? */
X mask = randmasks [regwidth]; /* which mask is for that width? */
X
X /* Find values to extract row and col numbers from each element. */
X rowshift = cwidth; /* find dist to shift to get top bits (row) */
X colmask = (1<<cwidth)-1; /* find mask to extract */
X /* bottom bits (col) */
X
X /* Now cycle through all sequence elements. */
X
X element = 1; /* 1st element (could be any nonzero) */
X do {
X row = element >> rowshift; /* find row number for this pixel */
X column = element & colmask; /* and how many columns across? */
X if ((row < height) /* does element fall in the array? */
X && (column < width)) /* ...must check row AND column */
X copy (row, column); /* in bounds: copy the (r,c)'th pixel */
X
X /* Compute the next sequence element */
X if (element & 1) /* is the low bit set? */
X element = (element >>1)^mask; /* yes: shift value, /*
X /* XOR in mask */
X else element = (element >>1); /* no: just shift the value */
X } while (element != 1); /* loop until we return to */
X /* original element */
X
X copy (0, 0); /* kludge: element never comes up zero */
X} /* end of dissolve2() */
END_OF_FILE
if test 4596 -ne `wc -c <'Dissolve.c'`; then
echo shar: \"'Dissolve.c'\" unpacked with wrong size!
fi
# end of 'Dissolve.c'
fi
if test -f 'DoubleLine.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'DoubleLine.c'\"
else
echo shar: Extracting \"'DoubleLine.c'\" \(4647 characters\)
sed "s/^X//" >'DoubleLine.c' <<'END_OF_FILE'
X/*
XSymmetric Double Step Line Algorithm
Xby Brian Wyvill
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#define swap(a,b) {a^=b; b^=a; a^=b;}
X#define absolute(i,j,k) ( (i-j)*(k = ( (i-j)<0 ? -1 : 1)))
Xint
Xsymwuline(a1, b1, a2, b2) int a1, b1, a2, b2;
X{
X int dx, dy, incr1, incr2, D, x, y, xend, c, pixels_left;
X int x1, y1;
X int sign_x, sign_y, step, reverse, i;
X
X dx = absolute(a2, a1, sign_x);
X dy = absolute(b2, b1, sign_y);
X /* decide increment sign by the slope sign */
X if (sign_x == sign_y)
X step = 1;
X else
X step = -1;
X
X if (dy > dx) { /* chooses axis of greatest movement (make
X * dx) */
X swap(a1, b1);
X swap(a2, b2);
X swap(dx, dy);
X reverse = 1;
X } else
X reverse = 0;
X /* note error check for dx==0 should be included here */
X if (a1 > a2) { /* start from the smaller coordinate */
X x = a2;
X y = b2;
X x1 = a1;
X y1 = b1;
X } else {
X x = a1;
X y = b1;
X x1 = a2;
X y1 = b2;
X }
X
X
X /* Note dx=n implies 0 - n or (dx+1) pixels to be set */
X /* Go round loop dx/4 times then plot last 0,1,2 or 3 pixels */
X /* In fact (dx-1)/4 as 2 pixels are already plottted */
X xend = (dx - 1) / 4;
X pixels_left = (dx - 1) % 4; /* number of pixels left over at the
X * end */
X plot(x, y, reverse);
X plot(x1, y1, reverse); /* plot first two points */
X incr2 = 4 * dy - 2 * dx;
X if (incr2 < 0) { /* slope less than 1/2 */
X c = 2 * dy;
X incr1 = 2 * c;
X D = incr1 - dx;
X
X for (i = 0; i < xend; i++) { /* plotting loop */
X ++x;
X --x1;
X if (D < 0) {
X /* pattern 1 forwards */
X plot(x, y, reverse);
X plot(++x, y, reverse);
X /* pattern 1 backwards */
X plot(x1, y1, reverse);
X plot(--x1, y1, reverse);
X D += incr1;
X } else {
X if (D < c) {
X /* pattern 2 forwards */
X plot(x, y, reverse);
X plot(++x, y += step, reverse);
X /* pattern 2 backwards */
X plot(x1, y1, reverse);
X plot(--x1, y1 -= step, reverse);
X } else {
X /* pattern 3 forwards */
X plot(x, y += step, reverse);
X plot(++x, y, reverse);
X /* pattern 3 backwards */
X plot(x1, y1 -= step, reverse);
X plot(--x1, y1, reverse);
X }
X D + = incr2;
X }
X } /* end for */
X
X /* plot last pattern */
X if (pixels_left) {
X if (D < 0) {
X plot(++x, y, reverse); /* pattern 1 */
X if (pixels_left > 1)
X plot(++x, y, reverse);
X if (pixels_left > 2)
X plot(--x1, y1, reverse);
X } else {
X if (D < c) {
X plot(++x, y, reverse); /* pattern 2 */
X if (pixels_left > 1)
X plot(++x, y += step, reverse);
X if (pixels_left > 2)
X plot(--x1, y1, reverse);
X } else {
X /* pattern 3 */
X plot(++x, y += step, reverse);
X if (pixels_left > 1)
X plot(++x, y, reverse);
X if (pixels_left > 2)
X plot(--x1, y1 -= step, reverse);
X }
X }
X } /* end if pixels_left */
X }
X /* end slope < 1/2 */
X else { /* slope greater than 1/2 */
X c = 2 * (dy - dx);
X incr1 = 2 * c;
X D = incr1 + dx;
X for (i = 0; i < xend; i++) {
X ++x;
X --x1;
X if (D > 0) {
X /* pattern 4 forwards */
X plot(x, y += step, reverse);
X plot(++x, y += step, reverse);
X /* pattern 4 backwards */
X plot(x1, y1 -= step, reverse);
X plot(--x1, y1 -= step, reverse);
X D += incr1;
X } else {
X if (D < c) {
X /* pattern 2 forwards */
X plot(x, y, reverse);
X plot(++x, y += step, reverse);
X
X /* pattern 2 backwards */
X plot(x1, y1, reverse);
X plot(--x1, y1 -= step, reverse);
X } else {
X /* pattern 3 forwards */
X plot(x, y += step, reverse);
X plot(++x, y, reverse);
X /* pattern 3 backwards */
X plot(x1, y1 -= step, reverse);
X plot(--x1, y1, reverse);
X }
X D += incr2;
X }
X } /* end for */
X /* plot last pattern */
X if (pixels_left) {
X if (D > 0) {
X plot(++x, y += step, reverse); /* pattern 4 */
X if (pixels_left > 1)
X plot(++x, y += step, reverse);
X if (pixels_left > 2)
X plot(--x1, y1 -= step, reverse);
X } else {
X if (D < c) {
X plot(++x, y, reverse); /* pattern 2 */
X if (pixels_left > 1)
X plot(++x, y += step, reverse);
X if (pixels_left > 2)
X plot(--x1, y1, reverse);
X } else {
X /* pattern 3 */
X plot(++x, y += step, reverse);
X if (pixels_left > 1)
X plot(++x, y, reverse);
X if (pixels_left > 2) {
X if (D > c) /* step 3 */
X plot(--x1, y1 -= step, reverse);
X else /* step 2 */
X plot(--x1, y1, reverse);
X }
X }
X }
X }
X }
X}
X/* non-zero flag indicates the pixels needing swap back. */
Xplot(x, y, flag) int x, y, flag;
X{
X if (flag)
X setpixel(y, x);
X else
X setpixel(x, y);
X}
X
END_OF_FILE
if test 4647 -ne `wc -c <'DoubleLine.c'`; then
echo shar: \"'DoubleLine.c'\" unpacked with wrong size!
fi
# end of 'DoubleLine.c'
fi
if test -f 'GraphicsGems.h' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'GraphicsGems.h'\"
else
echo shar: Extracting \"'GraphicsGems.h'\" \(4049 characters\)
sed "s/^X//" >'GraphicsGems.h' <<'END_OF_FILE'
X/*
X * GraphicsGems.h
X * Version 1.0 - Andrew Glassner
X * from "Graphics Gems", Academic Press, 1990
X */
X
X#ifndef GG_H
X
X#define GG_H 1
X
X/*********************/
X/* 2d geometry types */
X/*********************/
X
Xtypedef struct Point2Struct { /* 2d point */
X double x, y;
X } Point2;
Xtypedef Point2 Vector2;
X
Xtypedef struct IntPoint2Struct { /* 2d integer point */
X int x, y;
X } IntPoint2;
X
Xtypedef struct Matrix3Struct { /* 3-by-3 matrix */
X double element[3][3];
X } Matrix3;
X
Xtypedef struct Box2dStruct { /* 2d box */
X Point2 min, max;
X } Box2;
X
X
X/*********************/
X/* 3d geometry types */
X/*********************/
X
Xtypedef struct Point3Struct { /* 3d point */
X double x, y, z;
X } Point3;
Xtypedef Point3 Vector3;
X
Xtypedef struct IntPoint3Struct { /* 3d integer point */
X int x, y, z;
X } IntPoint3;
X
X
Xtypedef struct Matrix4Struct { /* 4-by-4 matrix */
X double element[4][4];
X } Matrix4;
X
Xtypedef struct Box3dStruct { /* 3d box */
X Point3 min, max;
X } Box3;
X
X
X
X/***********************/
X/* one-argument macros */
X/***********************/
X
X/* absolute value of a */
X#define ABS(a) (((a)<0) ? -(a) : (a))
X
X/* round a to nearest integer towards 0 */
X#define FLOOR(a) ((a)>0 ? (int)(a) : -(int)(-a))
X
X/* round a to nearest integer away from 0 */
X#define CEILING(a) \
X((a)==(int)(a) ? (a) : (a)>0 ? 1+(int)(a) : -(1+(int)(-a)))
X
X/* round a to nearest int */
X#define ROUND(a) ((a)>0 ? (int)(a+0.5) : -(int)(0.5-a))
X
X/* take sign of a, either -1, 0, or 1 */
X#define ZSGN(a) (((a)<0) ? -1 : (a)>0 ? 1 : 0)
X
X/* take binary sign of a, either -1, or 1 if >= 0 */
X#define SGN(a) (((a)<0) ? -1 : 0)
X
X/* shout if something that should be true isn't */
X#define ASSERT(x) \
Xif (!(x)) fprintf(stderr," Assert failed: x\n");
X
X/* square a */
X#define SQR(a) ((a)*(a))
X
X
X/***********************/
X/* two-argument macros */
X/***********************/
X
X/* find minimum of a and b */
X#define MIN(a,b) (((a)<(b))?(a):(b))
X
X/* find maximum of a and b */
X#define MAX(a,b) (((a)>(b))?(a):(b))
X
X/* swap a and b (see Gem by Wyvill) */
X#define SWAP(a,b) { a^=b; b^=a; a^=b; }
X
X/* linear interpolation from l (when a=0) to h (when a=1)*/
X/* (equal to (a*h)+((1-a)*l) */
X#define LERP(a,l,h) ((l)+(((h)-(l))*(a)))
X
X/* clamp the input to the specified range */
X#define CLAMP(v,l,h) ((v)<(l) ? (l) : (v) > (h) ? (h) : v)
X
X
X/****************************/
X/* memory allocation macros */
X/****************************/
X
X/* create a new instance of a structure (see Gem by Hultquist) */
X#define NEWSTRUCT(x) (struct x *)(malloc((unsigned)sizeof(struct x)))
X
X/* create a new instance of a type */
X#define NEWTYPE(x) (x *)(malloc((unsigned)sizeof(x)))
X
X
X/********************/
X/* useful constants */
X/********************/
X
X#define PI 3.141592 /* the venerable pi */
X#define PITIMES2 6.283185 /* 2 * pi */
X#define PIOVER2 1.570796 /* pi / 2 */
X#define E 2.718282 /* the venerable e */
X#define SQRT2 1.414214 /* sqrt(2) */
X#define SQRT3 1.732051 /* sqrt(3) */
X#define GOLDEN 1.618034 /* the golden ratio */
X#define DTOR 0.017453 /* convert degrees to radians */
X#define RTOD 57.29578 /* convert radians to degrees */
X
X
X/************/
X/* booleans */
X/************/
X
X#define TRUE 1
X#define FALSE 0
X#define ON 1
X#define OFF 0
Xtypedef int boolean; /* boolean data type */
Xtypedef boolean flag; /* flag data type */
X
Xextern double V2SquaredLength(), V2Length();
Xextern double V2Dot(), V2DistanceBetween2Points();
Xextern Vector2 *V2Negate(), *V2Normalize(), *V2Scale(), *V2Add(), *V2Sub();
Xextern Vector2 *V2Lerp(), *V2Combine(), *V2Mul(), *V2MakePerpendicular();
Xextern Vector2 *V2New(), *V2Duplicate();
Xextern Point2 *V2MulPointByMatrix();
Xextern Matrix3 *V2MatMul();
X
Xextern double V3SquaredLength(), V3Length();
Xextern double V3Dot(), V3DistanceBetween2Points();
Xextern Vector3 *V3Normalize(), *V3Scale(), *V3Add(), *V3Sub();
Xextern Vector3 *V3Lerp(), *V3Combine(), *V3Mul(), *V3Cross();
Xextern Vector3 *V3New(), *V3Duplicate();
Xextern Point3 *V3MulPointByMatrix();
Xextern Matrix4 *V3MatMul();
X
Xextern double RegulaFalsi(), NewtonRaphson(), findroot();
X
X#endif
END_OF_FILE
if test 4049 -ne `wc -c <'GraphicsGems.h'`; then
echo shar: \"'GraphicsGems.h'\" unpacked with wrong size!
fi
# end of 'GraphicsGems.h'
fi
if test -f 'LineEdge.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'LineEdge.c'\"
else
echo shar: Extracting \"'LineEdge.c'\" \(3852 characters\)
sed "s/^X//" >'LineEdge.c' <<'END_OF_FILE'
X/*
XFast Line-Edge Intersections on a Uniform Grid
Xby Andrew Shapira
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#include "GraphicsGems.h"
X
X#define OCTANT(f1, f2, f3, f4, f5, i1, s1, r1, r2) \
X for (f1, f2, f3, nr = 0; f4; f5) { \
X if (nr < liconst) { \
X if (i1) \
X r1(&C); \
X else \
X vertex(&C); \
X } \
X else { \
X s1; \
X if (nr -= liconst) { \
X r2(&C); \
X r1(&C); \
X } \
X else \
X vertex(&C); \
X } \
X }
X
Xfind_intersections(Pptr, Qptr)
XIntPoint2 *Pptr, *Qptr; /* P and Q as described in gem text */
X{
X IntPoint2 P, Q; /* P and Q, dereferenced for speed */
X IntPoint2 C; /* current grid point */
X int nr; /* remainder */
X int deltax, deltay; /* Q.x - P.x, Q.y - P.y */
X int liconst; /* loop-invariant constant */
X
X P.x = Pptr->x;
X P.y = Pptr->y;
X Q.x = Qptr->x;
X Q.y = Qptr->y;
X deltax = Q.x - P.x;
X deltay = Q.y - P.y;
X
X
X /* for reference purposes, let theta be the angle from P to Q */
X
X if ((deltax >= 0) && (deltay >= 0) && (deltay < deltax))
X /* 0 <= theta < 45 */
X OCTANT(C.x = P.x + 1, C.y = P.y, liconst = deltax - deltay,
X C.x < Q.x, C.x++, nr += deltay, C.y++, up, left)
X else if ((deltax > 0) && (deltay >= 0) && (deltay >= deltax))
X /* 45 <= theta < 90 */
X OCTANT(C.y = P.y + 1, C.x = P.x, liconst = deltay - deltax,
X C.y < Q.y, C.y++, nr += deltax, C.x++, right, down)
X else if ((deltax <= 0) && (deltay >= 0) && (deltay > -deltax))
X /* 90 <= theta < 135 */
X OCTANT(C.y = P.y + 1, C.x = P.x, liconst = deltay + deltax,
X C.y < Q.y, C.y++, nr -= deltax, C.x--, left, down)
X else if ((deltax <= 0) && (deltay > 0) && (deltay <= -deltax))
X /* 135 <= theta < 180 */
X OCTANT(C.x = P.x - 1, C.y = P.y, liconst = -deltax - deltay,
X C.x > Q.x, C.x--, nr += deltay, C.y++, up, right)
X else if ((deltax <= 0) && (deltay <= 0) && (deltay > deltax))
X /* 180 <= theta < 225 */
X OCTANT(C.x = P.x - 1, C.y = P.y, liconst = -deltax + deltay,
X C.x > Q.x, C.x--, nr -= deltay, C.y--, down, right)
X else if ((deltax < 0) && (deltay <= 0) && (deltay <= deltax))
X /* 225 <= theta < 270 */
X OCTANT(C.y = P.y - 1, C.x = P.x, liconst = -deltay + deltax,
X C.y > Q.y, C.y--, nr -= deltax, C.x--, left, up)
X else if ((deltax >= 0) && (deltay <= 0) && (-deltay > deltax))
X /* 270 <= theta < 315 */
X OCTANT(C.y = P.y - 1, C.x = P.x, liconst = -deltay - deltax,
X C.y > Q.y, C.y--, nr += deltax, C.x++, right, up)
X else if ((deltax >= 0) && (deltay < 0) && (-deltay <= deltax))
X /* 315 <= theta < 360 */
X OCTANT(C.x = P.x + 1, C.y = P.y, liconst = deltax + deltay,
X C.x < Q.x, C.x++, nr -= deltay, C.y--, down, left)
X else {} /* P = Q */
X}
X
Xvertex(I)
XIntPoint2 *I;
X{
X /* Note: replace printf with code to process vertex, if desired */
X (void) printf("vertex at %d %d\n", I->x, I->y);
X}
X
Xleft(I)
XIntPoint2 *I;
X{
X
X /* Note: replace printf with code to process leftward */
X /* intersection, if desired */
X (void) printf("left from %d %d\n", I->x, I->y);
X}
X
Xup(I)
XIntPoint2 *I;
X{
X /* Note: replace printf with code to process upward */
X /* intersection, if desired */
X (void) printf("up from %d %d\n", I->x, I->y);
X}
X
Xright(I)
XIntPoint2 *I;
X{
X /* Note: replace printf with code to process rightward */
X /* intersection, if desired */
X (void) printf("right from %d %d\n", I->x, I->y);
X}
X
Xdown(I)
XIntPoint2 *I;
X{
X /* Note: replace printf with code to process downward */
X /* intersection, if desired */
X (void) printf("down from %d %d\n", I->x, I->y);
X}
END_OF_FILE
if test 3852 -ne `wc -c <'LineEdge.c'`; then
echo shar: \"'LineEdge.c'\" unpacked with wrong size!
fi
# end of 'LineEdge.c'
fi
if test -f 'MatrixInvert.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'MatrixInvert.c'\"
else
echo shar: Extracting \"'MatrixInvert.c'\" \(4893 characters\)
sed "s/^X//" >'MatrixInvert.c' <<'END_OF_FILE'
X/*
XMatrix Inversion
Xby Richard Carling
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#define SMALL_NUMBER 1.e-8
X/*
X * inverse( original_matrix, inverse_matrix )
X *
X * calculate the inverse of a 4x4 matrix
X *
X * -1
X * A = ___1__ adjoint A
X * det A
X */
X
X#include "GraphicsGems.h"
Xinverse( in, out ) Matrix4 *in, *out;
X{
X int i, j;
X double det, det4x4();
X
X /* calculate the adjoint matrix */
X
X adjoint( in, out );
X
X /* calculate the 4x4 determinent
X * if the determinent is zero,
X * then the inverse matrix is not unique.
X */
X
X det = det4x4( in );
X
X if ( fabs( det ) < SMALL_NUMBER) {
X printf("Non-singular matrix, no inverse!\n");
X exit();
X }
X
X /* scale the adjoint matrix to get the inverse */
X
X for (i=0; i<4; i++)
X for(j=0; j<4; j++)
X out->element[i][j] = out->element[i][j] / det;
X}
X
X
X/*
X * adjoint( original_matrix, inverse_matrix )
X *
X * calculate the adjoint of a 4x4 matrix
X *
X * Let a denote the minor determinant of matrix A obtained by
X * ij
X *
X * deleting the ith row and jth column from A.
X *
X * i+j
X * Let b = (-1) a
X * ij ji
X *
X * The matrix B = (b ) is the adjoint of A
X * ij
X */
X
Xadjoint( in, out ) Matrix4 *in; Matrix4 *out;
X{
X double a1, a2, a3, a4, b1, b2, b3, b4;
X double c1, c2, c3, c4, d1, d2, d3, d4;
X double det3x3();
X
X /* assign to individual variable names to aid */
X /* selecting correct values */
X
X a1 = in->element[0][0]; b1 = in->element[0][1];
X c1 = in->element[0][2]; d1 = in->element[0][3];
X
X a2 = in->element[1][0]; b2 = in->element[1][1];
X c2 = in->element[1][2]; d2 = in->element[1][3];
X
X a3 = in->element[2][0]; b3 = in->element[2][1];
X c3 = in->element[2][2]; d3 = in->element[2][3];
X
X a4 = in->element[3][0]; b4 = in->element[3][1];
X c4 = in->element[3][2]; d4 = in->element[3][3];
X
X
X /* row column labeling reversed since we transpose rows & columns */
X
X out->element[0][0] = det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4);
X out->element[1][0] = - det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4);
X out->element[2][0] = det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4);
X out->element[3][0] = - det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4);
X
X out->element[0][1] = - det3x3( b1, b3, b4, c1, c3, c4, d1, d3, d4);
X out->element[1][1] = det3x3( a1, a3, a4, c1, c3, c4, d1, d3, d4);
X out->element[2][1] = - det3x3( a1, a3, a4, b1, b3, b4, d1, d3, d4);
X out->element[3][1] = det3x3( a1, a3, a4, b1, b3, b4, c1, c3, c4);
X
X out->element[0][2] = det3x3( b1, b2, b4, c1, c2, c4, d1, d2, d4);
X out->element[1][2] = - det3x3( a1, a2, a4, c1, c2, c4, d1, d2, d4);
X out->element[2][2] = det3x3( a1, a2, a4, b1, b2, b4, d1, d2, d4);
X out->element[3][2] = - det3x3( a1, a2, a4, b1, b2, b4, c1, c2, c4);
X
X out->element[0][3] = - det3x3( b1, b2, b3, c1, c2, c3, d1, d2, d3);
X out->element[1][3] = det3x3( a1, a2, a3, c1, c2, c3, d1, d2, d3);
X out->element[2][3] = - det3x3( a1, a2, a3, b1, b2, b3, d1, d2, d3);
X out->element[3][3] = det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3);
X}
X/*
X * double = det4x4( matrix )
X *
X * calculate the determinent of a 4x4 matrix.
X */
Xdouble det4x4( m ) Matrix4 *m;
X{
X double ans;
X double a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
X double det3x3();
X
X /* assign to individual variable names to aid selecting */
X /* correct elements */
X
X a1 = m->element[0][0]; b1 = m->element[0][1];
X c1 = m->element[0][2]; d1 = m->element[0][3];
X
X a2 = m->element[1][0]; b2 = m->element[1][1];
X c2 = m->element[1][2]; d2 = m->element[1][3];
X
X a3 = m->element[2][0]; b3 = m->element[2][1];
X c3 = m->element[2][2]; d3 = m->element[2][3];
X
X a4 = m->element[3][0]; b4 = m->element[3][1];
X c4 = m->element[3][2]; d4 = m->element[3][3];
X
X ans = a1 * det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4)
X - b1 * det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4)
X + c1 * det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4)
X - d1 * det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4);
X return ans;
X}
X
X
X
X/*
X * double = det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3 )
X *
X * calculate the determinent of a 3x3 matrix
X * in the form
X *
X * | a1, b1, c1 |
X * | a2, b2, c2 |
X * | a3, b3, c3 |
X */
X
Xdouble det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3 )
Xdouble a1, a2, a3, b1, b2, b3, c1, c2, c3;
X{
X double ans;
X double det2x2();
X
X ans = a1 * det2x2( b2, b3, c2, c3 )
X - b1 * det2x2( a2, a3, c2, c3 )
X + c1 * det2x2( a2, a3, b2, b3 );
X return ans;
X}
X
X/*
X * double = det2x2( double a, double b, double c, double d )
X *
X * calculate the determinent of a 2x2 matrix.
X */
X
Xdouble det2x2( a, b, c, d)
Xdouble a, b, c, d;
X{
X double ans;
X ans = a * d - b * c;
X return ans;
X}
X
X
END_OF_FILE
if test 4893 -ne `wc -c <'MatrixInvert.c'`; then
echo shar: \"'MatrixInvert.c'\" unpacked with wrong size!
fi
# end of 'MatrixInvert.c'
fi
if test -f 'PolyScan/poly_clip.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'PolyScan/poly_clip.c'\"
else
echo shar: Extracting \"'PolyScan/poly_clip.c'\" \(4217 characters\)
sed "s/^X//" >'PolyScan/poly_clip.c' <<'END_OF_FILE'
X/*
X * Generic Convex Polygon Scan Conversion and Clipping
X * by Paul Heckbert
X * from "Graphics Gems", Academic Press, 1990
X */
X
X/*
X * poly_clip.c: homogeneous 3-D convex polygon clipper
X *
X * Paul Heckbert 1985, Dec 1989
X */
X
X#include "poly.h"
X
X#define SWAP(a, b, temp) {temp = a; a = b; b = temp;}
X#define COORD(vert, i) ((double *)(vert))[i]
X
X#define CLIP_AND_SWAP(elem, sign, k, p, q, r) { \
X poly_clip_to_halfspace(p, q, &v->elem-(double *)v, sign, sign*k); \
X if (q->n==0) {p1->n = 0; return POLY_CLIP_OUT;} \
X SWAP(p, q, r); \
X}
X
X/*
X * poly_clip_to_box: Clip the convex polygon p1 to the screen space box
X * using the homogeneous screen coordinates (sx, sy, sz, sw) of each vertex,
X * testing if v->sx/v->sw > box->x0 and v->sx/v->sw < box->x1,
X * and similar tests for y and z, for each vertex v of the polygon.
X * If polygon is entirely inside box, then POLY_CLIP_IN is returned.
X * If polygon is entirely outside box, then POLY_CLIP_OUT is returned.
X * Otherwise, if the polygon is cut by the box, p1 is modified and
X * POLY_CLIP_PARTIAL is returned.
X */
X
Xint poly_clip_to_box(p1, box)
Xregister Poly *p1;
Xregister Poly_box *box;
X{
X int x0out = 0, x1out = 0, y0out = 0, y1out = 0, z0out = 0, z1out = 0;
X register int i;
X register Poly_vert *v;
X Poly p2, *p, *q, *r;
X
X /* count vertices "outside" with respect to each of the six planes */
X for (v=p1->vert, i=p1->n; i>0; i--, v++) {
X if (v->sx < box->x0*v->sw) x0out++; /* out on left */
X if (v->sx > box->x1*v->sw) x1out++; /* out on right */
X if (v->sy < box->y0*v->sw) y0out++; /* out on top */
X if (v->sy > box->y1*v->sw) y1out++; /* out on bottom */
X if (v->sz < box->z0*v->sw) z0out++; /* out on near */
X if (v->sz > box->z1*v->sw) z1out++; /* out on far */
X }
X
X /* check if all vertices inside */
X if (x0out+x1out+y0out+y1out+z0out+z1out == 0) return POLY_CLIP_IN;
X
X /* check if all vertices are "outside" any of the six planes */
X if (x0out==p1->n || x1out==p1->n || y0out==p1->n ||
X y1out==p1->n || z0out==p1->n || z1out==p1->n) {
X p1->n = 0;
X return POLY_CLIP_OUT;
X }
X
X /*
X * now clip against each of the planes that might cut the polygon,
X * at each step toggling between polygons p1 and p2
X */
X p = p1;
X q = &p2;
X if (x0out) CLIP_AND_SWAP(sx, -1., box->x0, p, q, r);
X if (x1out) CLIP_AND_SWAP(sx, 1., box->x1, p, q, r);
X if (y0out) CLIP_AND_SWAP(sy, -1., box->y0, p, q, r);
X if (y1out) CLIP_AND_SWAP(sy, 1., box->y1, p, q, r);
X if (z0out) CLIP_AND_SWAP(sz, -1., box->z0, p, q, r);
X if (z1out) CLIP_AND_SWAP(sz, 1., box->z1, p, q, r);
X
X /* if result ended up in p2 then copy it to p1 */
X if (p==&p2)
X bcopy(&p2, p1, sizeof(Poly)-(POLY_NMAX-p2.n)*sizeof(Poly_vert));
X return POLY_CLIP_PARTIAL;
X}
X
X/*
X * poly_clip_to_halfspace: clip convex polygon p against a plane,
X * copying the portion satisfying sign*s[index] < k*sw into q,
X * where s is a Poly_vert* cast as a double*.
X * index is an index into the array of doubles at each vertex, such that
X * s[index] is sx, sy, or sz (screen space x, y, or z).
X * Thus, to clip against xmin, use
X * poly_clip_to_halfspace(p, q, XINDEX, -1., -xmin);
X * and to clip against xmax, use
X * poly_clip_to_halfspace(p, q, XINDEX, 1., xmax);
X */
X
Xvoid poly_clip_to_halfspace(p, q, index, sign, k)
XPoly *p, *q;
Xregister int index;
Xdouble sign, k;
X{
X register int m;
X register double *up, *vp, *wp;
X register Poly_vert *v;
X int i;
X Poly_vert *u;
X double t, tu, tv;
X
X q->n = 0;
X q->mask = p->mask;
X
X /* start with u=vert[n-1], v=vert[0] */
X u = &p->vert[p->n-1];
X tu = sign*COORD(u, index) - u->sw*k;
X for (v= &p->vert[0], i=p->n; i>0; i--, u=v, tu=tv, v++) {
X /* on old polygon (p), u is previous vertex, v is current vertex */
X /* tv is negative if vertex v is in */
X tv = sign*COORD(v, index) - v->sw*k;
X if (tu<=0. ^ tv<=0.) {
X /* edge crosses plane; add intersection point to q */
X t = tu/(tu-tv);
X up = (double *)u;
X vp = (double *)v;
X wp = (double *)&q->vert[q->n];
X for (m=p->mask; m!=0; m>>=1, up++, vp++, wp++)
X if (m&1) *wp = *up+t*(*vp-*up);
X q->n++;
X }
X if (tv<=0.) /* vertex v is in, copy it to q */
X q->vert[q->n++] = *v;
X }
X}
END_OF_FILE
if test 4217 -ne `wc -c <'PolyScan/poly_clip.c'`; then
echo shar: \"'PolyScan/poly_clip.c'\" unpacked with wrong size!
fi
# end of 'PolyScan/poly_clip.c'
fi
if test -f 'PolyScan/poly_scan.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'PolyScan/poly_scan.c'\"
else
echo shar: Extracting \"'PolyScan/poly_scan.c'\" \(4638 characters\)
sed "s/^X//" >'PolyScan/poly_scan.c' <<'END_OF_FILE'
X/*
X * Generic Convex Polygon Scan Conversion and Clipping
X * by Paul Heckbert
X * from "Graphics Gems", Academic Press, 1990
X */
X
X/*
X * poly_scan.c: point-sampled scan conversion of convex polygons
X *
X * Paul Heckbert 1985, Dec 1989
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "poly.h"
X
X/*
X * poly_scan: Scan convert a polygon, calling pixelproc at each pixel with an
X * interpolated Poly_vert structure. Polygon can be clockwise or ccw.
X * Polygon is clipped in 2-D to win, the screen space window.
X *
X * Scan conversion is done on the basis of Poly_vert fields sx and sy.
X * These two must always be interpolated, and only they have special meaning
X * to this code; any other fields are blindly interpolated regardless of
X * their semantics.
X *
X * The pixelproc subroutine takes the arguments:
X *
X * pixelproc(x, y, point)
X * int x, y;
X * Poly_vert *point;
X *
X * All the fields of point indicated by p->mask will be valid inside pixelproc
X * except sx and sy. If they were computed, they would have values
X * sx=x+.5 and sy=y+.5, since sampling is done at pixel centers.
X */
X
Xvoid poly_scan(p, win, pixelproc)
Xregister Poly *p; /* polygon */
XWindow *win; /* 2-D screen space clipping window */
Xvoid (*pixelproc)(); /* procedure called at each pixel */
X{
X register int i, li, ri, y, ly, ry, top, rem, mask;
X double ymin;
X Poly_vert l, r, dl, dr;
X
X if (p->n>POLY_NMAX) {
X fprintf(stderr, "poly_scan: too many vertices: %d\n", p->n);
X return;
X }
X
X ymin = HUGE;
X for (i=0; i<p->n; i++) /* find top vertex (y points down) */
X if (p->vert[i].sy < ymin) {
X ymin = p->vert[i].sy;
X top = i;
X }
X
X li = ri = top; /* left and right vertex indices */
X rem = p->n; /* number of vertices remaining */
X y = ceil(ymin-.5); /* current scan line */
X ly = ry = y-1; /* lower end of left & right edges */
X mask = p->mask & ~POLY_MASK(sy); /* stop interpolating screen y */
X
X while (rem>0) { /* scan in y, activating new edges on left & right */
X /* as scan line passes over new vertices */
X
X while (ly<=y && rem>0) { /* advance left edge? */
X rem--;
X i = li-1; /* step ccw down left side */
X if (i<0) i = p->n-1;
X incrementalize_y(&p->vert[li], &p->vert[i], &l, &dl, y, mask);
X ly = floor(p->vert[i].sy+.5);
X li = i;
X }
X while (ry<=y && rem>0) { /* advance right edge? */
X rem--;
X i = ri+1; /* step cw down right edge */
X if (i>=p->n) i = 0;
X incrementalize_y(&p->vert[ri], &p->vert[i], &r, &dr, y, mask);
X ry = floor(p->vert[i].sy+.5);
X ri = i;
X }
X
X while (y<ly && y<ry) { /* do scanlines till end of l or r edge */
X if (y>=win->y0 && y<=win->y1)
X if (l.sx<=r.sx) scanline(y, &l, &r, win, pixelproc, mask);
X else scanline(y, &r, &l, win, pixelproc, mask);
X y++;
X increment(&l, &dl, mask);
X increment(&r, &dr, mask);
X }
X }
X}
X
X/* scanline: output scanline by sampling polygon at Y=y+.5 */
X
Xstatic scanline(y, l, r, win, pixelproc, mask)
Xint y, mask;
XPoly_vert *l, *r;
XWindow *win;
Xvoid (*pixelproc)();
X{
X int x, lx, rx;
X Poly_vert p, dp;
X
X mask &= ~POLY_MASK(sx); /* stop interpolating screen x */
X lx = ceil(l->sx-.5);
X if (lx<win->x0) lx = win->x0;
X rx = floor(r->sx-.5);
X if (rx>win->x1) rx = win->x1;
X if (lx>rx) return;
X incrementalize_x(l, r, &p, &dp, lx, mask);
X for (x=lx; x<=rx; x++) { /* scan in x, generating pixels */
X (*pixelproc)(x, y, &p);
X increment(&p, &dp, mask);
X }
X}
X
X/*
X * incrementalize_y: put intersection of line Y=y+.5 with edge between points
X * p1 and p2 in p, put change with respect to y in dp
X */
X
Xstatic incrementalize_y(p1, p2, p, dp, y, mask)
Xregister double *p1, *p2, *p, *dp;
Xregister int mask;
Xint y;
X{
X double dy, frac;
X
X dy = ((Poly_vert *)p2)->sy - ((Poly_vert *)p1)->sy;
X if (dy==0.) dy = 1.;
X frac = y+.5 - ((Poly_vert *)p1)->sy;
X
X for (; mask!=0; mask>>=1, p1++, p2++, p++, dp++)
X if (mask&1) {
X *dp = (*p2-*p1)/dy;
X *p = *p1+*dp*frac;
X }
X}
X
X/*
X * incrementalize_x: put intersection of line X=x+.5 with edge between points
X * p1 and p2 in p, put change with respect to x in dp
X */
X
Xstatic incrementalize_x(p1, p2, p, dp, x, mask)
Xregister double *p1, *p2, *p, *dp;
Xregister int mask;
Xint x;
X{
X double dx, frac;
X
X dx = ((Poly_vert *)p2)->sx - ((Poly_vert *)p1)->sx;
X if (dx==0.) dx = 1.;
X frac = x+.5 - ((Poly_vert *)p1)->sx;
X
X for (; mask!=0; mask>>=1, p1++, p2++, p++, dp++)
X if (mask&1) {
X *dp = (*p2-*p1)/dx;
X *p = *p1+*dp*frac;
X }
X}
X
Xstatic increment(p, dp, mask)
Xregister double *p, *dp;
Xregister int mask;
X{
X for (; mask!=0; mask>>=1, p++, dp++)
X if (mask&1)
X *p += *dp;
X}
END_OF_FILE
if test 4638 -ne `wc -c <'PolyScan/poly_scan.c'`; then
echo shar: \"'PolyScan/poly_scan.c'\" unpacked with wrong size!
fi
# end of 'PolyScan/poly_scan.c'
fi
if test -f 'Roots3And4.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Roots3And4.c'\"
else
echo shar: Extracting \"'Roots3And4.c'\" \(4127 characters\)
sed "s/^X//" >'Roots3And4.c' <<'END_OF_FILE'
X/*
X * Roots3And4.c
X *
X * Utility functions to find cubic and quartic roots,
X * coefficients are passed like this:
X *
X * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
X *
X * The functions return the number of non-complex roots and
X * put the values into the s array.
X *
X * Author: Jochen Schwarze (schwarze@isa.de)
X *
X * Jan 26, 1990 Version for Graphics Gems
X * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic
X * (reported by Mark Podlipec),
X * Old-style function definitions,
X * IsZero() as a macro
X */
X
X#include <math.h>
X
X/* epsilon surrounding for near zero values */
X
X#define EQN_EPS 1e-9
X#define IsZero(x) ((x) > -EQN_EPS && (x) < EQN_EPS)
X
X
Xint SolveQuadric(c, s)
X double c[ 3 ];
X double s[ 2 ];
X{
X double p, q, D;
X
X /* normal form: x^2 + px + q = 0 */
X
X p = c[ 1 ] / (2 * c[ 2 ]);
X q = c[ 0 ] / c[ 2 ];
X
X D = p * p - q;
X
X if (IsZero(D))
X {
X s[ 0 ] = - p;
X return 1;
X }
X else if (D < 0)
X {
X return 0;
X }
X else if (D > 0)
X {
X double sqrt_D = sqrt(D);
X
X s[ 0 ] = sqrt_D - p;
X s[ 1 ] = - sqrt_D - p;
X return 2;
X }
X}
X
X
Xint SolveCubic(c, s)
X double c[ 4 ];
X double s[ 3 ];
X{
X int i, num;
X double sub;
X double A, B, C;
X double sq_A, p, q;
X double cb_p, D;
X
X /* normal form: x^3 + Ax^2 + Bx + C = 0 */
X
X A = c[ 2 ] / c[ 3 ];
X B = c[ 1 ] / c[ 3 ];
X C = c[ 0 ] / c[ 3 ];
X
X /* substitute x = y - A/3 to eliminate quadric term:
X x^3 +px + q = 0 */
X
X sq_A = A * A;
X p = 1.0/3 * (- 1.0/3 * sq_A + B);
X q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
X
X /* use Cardano's formula */
X
X cb_p = p * p * p;
X D = q * q + cb_p;
X
X if (IsZero(D))
X {
X if (IsZero(q)) /* one triple solution */
X {
X s[ 0 ] = 0;
X num = 1;
X }
X else /* one single and one double solution */
X {
X double u = cbrt(-q);
X s[ 0 ] = 2 * u;
X s[ 1 ] = - u;
X num = 2;
X }
X }
X else if (D < 0) /* Casus irreducibilis: three real solutions */
X {
X double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
X double t = 2 * sqrt(-p);
X
X s[ 0 ] = t * cos(phi);
X s[ 1 ] = - t * cos(phi + M_PI / 3);
X s[ 2 ] = - t * cos(phi - M_PI / 3);
X num = 3;
X }
X else /* one real solution */
X {
X double sqrt_D = sqrt(D);
X double u = cbrt(sqrt_D - q);
X double v = - cbrt(sqrt_D + q);
X
X s[ 0 ] = u + v;
X num = 1;
X }
X
X /* resubstitute */
X
X sub = 1.0/3 * A;
X
X for (i = 0; i < num; ++i)
X s[ i ] -= sub;
X
X return num;
X}
X
X
Xint SolveQuartic(c, s)
X double c[ 5 ];
X double s[ 4 ];
X{
X double coeffs[ 4 ];
X double z, u, v, sub;
X double A, B, C, D;
X double sq_A, p, q, r;
X int i, num;
X
X /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
X
X A = c[ 3 ] / c[ 4 ];
X B = c[ 2 ] / c[ 4 ];
X C = c[ 1 ] / c[ 4 ];
X D = c[ 0 ] / c[ 4 ];
X
X /* substitute x = y - A/4 to eliminate cubic term:
X x^4 + px^2 + qx + r = 0 */
X
X sq_A = A * A;
X p = - 3.0/8 * sq_A + B;
X q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
X r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
X
X if (IsZero(r))
X {
X /* no absolute term: y(y^3 + py + q) = 0 */
X
X coeffs[ 0 ] = q;
X coeffs[ 1 ] = p;
X coeffs[ 2 ] = 0;
X coeffs[ 3 ] = 1;
X
X num = SolveCubic(coeffs, s);
X
X s[ num++ ] = 0;
X }
X else
X {
X /* solve the resolvent cubic ... */
X
X coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
X coeffs[ 1 ] = - r;
X coeffs[ 2 ] = - 1.0/2 * p;
X coeffs[ 3 ] = 1;
X
X (void) SolveCubic(coeffs, s);
X
X /* ... and take the one real solution ... */
X
X z = s[ 0 ];
X
X /* ... to build two quadric equations */
X
X u = z * z - r;
X v = 2 * z - p;
X
X if (IsZero(u))
X u = 0;
X else if (u > 0)
X u = sqrt(u);
X else
X return 0;
X
X if (IsZero(v))
X v = 0;
X else if (v > 0)
X v = sqrt(v);
X else
X return 0;
X
X coeffs[ 0 ] = z - u;
X coeffs[ 1 ] = q < 0 ? -v : v;
X coeffs[ 2 ] = 1;
X
X num = SolveQuadric(coeffs, s);
X
X coeffs[ 0 ]= z + u;
X coeffs[ 1 ] = q < 0 ? v : -v;
X coeffs[ 2 ] = 1;
X
X num += SolveQuadric(coeffs, s + num);
X }
X
X /* resubstitute */
X
X sub = 1.0/4 * A;
X
X for (i = 0; i < num; ++i)
X s[ i ] -= sub;
X
X return num;
X}
X
END_OF_FILE
if test 4127 -ne `wc -c <'Roots3And4.c'`; then
echo shar: \"'Roots3And4.c'\" unpacked with wrong size!
fi
# end of 'Roots3And4.c'
fi
echo shar: End of archive 3 \(of 5\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 4 5 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 5 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