Method IO Graphics PART03.SH

# ! /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 Wed Dec 12 20:45:14 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'\" \(4724 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 *  Nov 23, 1990    Some systems do not declare acos() and cbrt() in
X *                  <math.h>, though the functions exist in the library.
X *                  If large coefficients are used, EQN_EPS should be
X *                  reduced considerably (e.g. to 1E-30), results will be
X *                  correct but multiple roots might be reported more
X *                  than once.
X */
X
X#include    <math.h>
Xextern double   sqrt(), cbrt(), cos(), acos();
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#ifdef NOCBRT
X#define     cbrt(x)     ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
X                          ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
X#endif
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 4724 -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