please dont rip this site

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


file: /Techref/method/io/graphics/Part03.sh, 53KB, , updated: 1999/6/25 11:01, local time: 2022/6/30 21:12,
TOP NEW HELP FIND: 
44.201.94.72:LOG IN

 ©2022 These pages are served without commercial sponsorship. (No popup ads, etc...).Bandwidth abuse increases hosting cost forcing sponsorship or shutdown. This server aggressively defends against automated copying for any reason including offline viewing, duplication, etc... Please respect this requirement and DO NOT RIP THIS SITE. Questions?
Please DO link to this page! Digg it! / MAKE!

<A HREF="http://www.sxlist.com/techref/method/io/graphics/Part03.sh"> method io graphics Part03</A>

Did you find what you needed?

 

Welcome to sxlist.com!


Site supported by
sales, advertizing,
& kind contributors
just like you!

Please don't rip/copy
(here's why

Copies of the site on CD
are available at minimal cost.
 

Welcome to www.sxlist.com!

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  .