From bacchus.pa.dec.com!decwrl!uunet!allbery Mon Oct 15 17:35:24 PDT 1990 Article: 1856 of comp.sources.misc Path: bacchus.pa.dec.com!decwrl!uunet!allbery From: craig@weedeater.math.yale.edu (Craig Kolb) Newsgroups: comp.sources.misc Subject: v15i024: Graphics Gems, Part 2/5 Message-ID: <108571@uunet.UU.NET> Date: 15 Oct 90 01:12:16 GMT Sender: allbery@uunet.UU.NET Lines: 2113 Approved: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) X-UNIX-From: craig@weedeater.math.yale.edu Posting-number: Volume 15, Issue 24 Submitted-by: Craig Kolb Archive-name: ggems/part02 #! /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 '2DClip/bio.c' <<'END_OF_FILE' X X/* X * file bio.c X * contains the basic operations X * X */ X#include X#include "GraphicsGems.h" X#include "line.h" X X/* X * def_contour X * X * Purpose: X * add a contour to the list X * NOTE: coordinates must already be converted into longs! X * X * x x coordinates of the end points of segments X * y y coordinates of the end points of segments X * n number of coordinate pairs X * no contour number (id), does no have to be unique! X * type type of clip operation desired CLIP_NORMAL means X * clip everything inside the contour X */ Xdef_contour(x, y, n, no, type) Xlong x[], y[]; Xint n, no, type; X{ X short i; X long dx1, dx2, dy1, dy2; X long minx, miny, maxx, maxy; X CONTOUR *cp; X SEGMENT *sp, *spp; X X if((cp = CL) == (CONTOUR *)NULL) { X cp = CL = NEWTYPE(CONTOUR); X } X else { X while(cp->_next != (CONTOUR *)NULL) X cp = cp->_next; X i = cp->_no; X cp = cp->_next = NEWTYPE(CONTOUR); X } X X cp->_next = (CONTOUR *)NULL; X cp->_no = no; X SET_ON(cp); X if(type == CLIP_NORMAL) X SET_INVERSE(cp); X else X SET_NORMAL(cp); X minx = miny = 1000000; X maxx = maxy = -1; X for(i=0; i_s = sp = NEWTYPE(SEGMENT); X sp->_from._x = x[0]; X sp->_from._y = y[0]; X sp->_to._x = x[1]; X sp->_to._y = y[1]; X } X else { X /* X * if necessary stretch the contour X * and skip the point X */ X dx1 = sp->_to._x - sp->_from._x; X dy1 = sp->_to._y - sp->_from._y; X dx2 = x[(i == n-1) ? 0 : i+1] - sp->_to._x; X dy2 = y[(i == n-1) ? 0 : i+1] - sp->_to._y; X if(dy2*dx1 == dy1*dx2) { X sp->_to._x = x[(i == n-1) ? 0 : i+1]; X sp->_to._y = y[(i == n-1) ? 0 : i+1]; X } X else { X spp = sp; X sp = sp->_next = NEWTYPE(SEGMENT); X sp->_prev = spp; X sp->_from._x = x[i]; X sp->_from._y = y[i]; X sp->_to._x = x[(i == n-1) ? 0 : i+1]; X sp->_to._y = y[(i == n-1) ? 0 : i+1]; X } X } X X/* X * calculate the enclosing box X */ X if(x[i] < minx) X minx = x[i]; X if(x[i] > maxx) X maxx = x[i]; X if(y[i] < miny) X miny = y[i]; X if(y[i] > maxy) X maxy = y[i]; X X } X cp->_minx = minx; X cp->_maxx = maxx; X cp->_miny = miny; X cp->_maxy = maxy; X sp->_next = cp->_s; X cp->_s->_prev = sp; X cp->_next = (CONTOUR *)NULL; X} X X/* X * get_contour_ptr X * X * PURPOSE X * get the pointer to a contour given its id X * with multiple id's first fit algorithm is X * used. Returns NULL in case of error. X * X * no id of contour X */ XCONTOUR *get_contour_ptr(no) Xint no; X{ X CONTOUR *cp; X X if((cp = CL) == (CONTOUR *)NULL) X return((CONTOUR *)NULL); X else { X while(cp != (CONTOUR *)NULL) { X if(cp->_no == no) X return(cp); X else X cp = cp->_next; X } X return((CONTOUR *)NULL); X } X} X X X/* X * del_contour X * X * PURPOSE X * delete contour(s) from the list with id X * no X */ Xdel_contour(no) Xint no; X{ X CONTOUR *cp, *cpp; X CONTOUR *qp = (CONTOUR *)NULL; X SEGMENT *sp, *spp; X X if((cp = CL) == (CONTOUR *)NULL) X return; X while(cp != (CONTOUR *)NULL) { X if(cp->_no == no) { X sp = cp->_s; X do { X spp = sp->_next; X free(sp); X sp = spp; X } while(sp != cp->_s); X cpp = cp->_next; X free(cp); X if(qp) X qp->_next = cpp; X else X CL = cpp; X cp = cpp; X } X else { X qp = cp; X cp = cp->_next; X } X } X} X X END_OF_FILE if test 3146 -ne `wc -c <'2DClip/bio.c'`; then echo shar: \"'2DClip/bio.c'\" unpacked with wrong size! fi # end of '2DClip/bio.c' fi if test -f '2DClip/clip.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'2DClip/clip.c'\" else echo shar: Extracting \"'2DClip/clip.c'\" \(2904 characters\) sed "s/^X//" >'2DClip/clip.c' <<'END_OF_FILE' X/* X * file clip.c X * contains the actual clipping routines X */ X#include X#include "GraphicsGems.h" X#include "line.h" X X/* X * vis_vector X * X * PURPOSE X * actual user interface. Draws a clipped line X * NOTE: coordinates are given in converted LONGS! X * X * xf, yf from coordinates of vector to be drawn X * xt, yt to coordinates of vector to be drawn X */ Xvis_vector(xf, yf, xt, yt) Xlong xf, yf, xt, yt; X{ X SEGMENT l; X X if(xf == xt && yf == yt) X return; X l._from._x = xf; X l._from._y = yf; X l._to._x = xt; X l._to._y = yt; X/* X * start at top of list X */ X clip(CL, &l); X} X X/* X * clip X * X * PURPOSE X * X * p pointer to polygon X * l pointer to line segment X */ Xclip(p, l) XCONTOUR *p; XSEGMENT *l; X{ X SEGMENT *sp, ss; X CLIST *sol; X POINT pt; X boolean up, delay, inside, p_inside(), disjunct(); X int i; X short nsol, nsmax = 2; X X X/* X * list exhausted do what you like X * we want to plot X */ X if(p == (CONTOUR *)NULL) { X move((l->_from._x), (l->_from._y)); X cont((l->_to._x), (l->_to._y)); X return; X } X/* X * polygon is switched off X * take next one X */ X if(!IS_ON(p)) { X clip(p->_next, l); X return; X } X/* X * comparison on basis of the X * enclosing rectangle X */ X if(disjunct(p, l)) { X if(!IS_NORMAL(p)) { X clip(p->_next, l); X } X return; X } X/* X * calculate possible intersections X */ X sol = (CLIST *) calloc(2, sizeof(CLIST)); X sol[0]._p._x = l->_from._x; X sol[0]._p._y = l->_from._y; X sol[0]._type = STD; X sol[1]._p._x = l->_to._x; X sol[1]._p._y = l->_to._y; X sol[1]._type = STD; X nsol = 2; X cross_calc(p, l, &sol, &nsol, nsmax); X pt._x = sol[0]._p._x; X pt._y = sol[0]._p._y; X X/* X * determine status of first point X */ X inside = p_inside(p, &pt); X if((!inside && IS_NORMAL(p)) || (inside && !IS_NORMAL(p))) X up = TRUE; X else X up = FALSE; X delay = FALSE; X/* X * process list of intersections X */ X for(i=1; i_next, &ss); X } X if(!delay) { X if(sol[i]._type != DELAY) X up = (up) ? FALSE : TRUE; X else X delay = TRUE; X } X else { X up = (up) ? FALSE : TRUE; X delay = FALSE; X } X } X free(sol); X} X X/* X * disjunct X * X * PURPOSE X * determine if the box enclosing the polygon X * stored in p and the box enclosing the line X * segment stored in l are disjunct. X * Return TRUE if disjunct else FALSE X * X * p points to the polygon structure X * l points to the linesegment structure X * X */ Xboolean disjunct(p, l) XCONTOUR *p; XSEGMENT *l; X X{ X if((max(l->_from._x, l->_to._x) < p->_minx) || X (min(l->_from._x, l->_to._x) > p->_maxx) || X (max(l->_from._y, l->_to._y) < p->_miny) || X (min(l->_from._y, l->_to._y) > p->_maxy) ) X return(TRUE); X else X return(FALSE); X} X X#define DEBUG X#ifdef DEBUG Xmove(x, y) Xlong x, y; X{ X printf("(%d,%d) ->", x, y); X} X Xcont(x, y) Xlong x, y; X{ X printf("(%d,%d)\n", x, y); X} X X#endif X X X X END_OF_FILE if test 2904 -ne `wc -c <'2DClip/clip.c'`; then echo shar: \"'2DClip/clip.c'\" unpacked with wrong size! fi # end of '2DClip/clip.c' fi if test -f 'BoxSphere.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'BoxSphere.c'\" else echo shar: Extracting \"'BoxSphere.c'\" \(3760 characters\) sed "s/^X//" >'BoxSphere.c' <<'END_OF_FILE' X/* XA Simple Method for Box-Sphere Intersection Testing Xby Jim Arvo Xfrom "Graphics Gems", Academic Press, 1990 X*/ X X#include "GraphicsGems.h" X X/* X * This routine tests for intersection between an n-dimensional X * axis-aligned box and an n-dimensional sphere. The mode argument X * indicates whether the objects are to be regarded as surfaces or X * solids. The values are: X * X * mode X * X * 0 Hollow Box, Hollow Sphere X * 1 Hollow Box, Solid Sphere X * 2 Solid Box, Hollow Sphere X * 3 Solid Box, Solid Sphere X * X*/ Xint Box_Sphere_Intersect( n, Bmin, Bmax, C, r, mode ) Xint n; /* The dimension of the space. */ Xfloat Bmin[]; /* The minium of the box for each axis. */ Xfloat Bmax[]; /* The maximum of the box for each axis. */ Xfloat C[]; /* The sphere center in n-space. */ Xfloat r; /* The radius of the sphere. */ Xint mode; /* Selects hollow or solid. */ X { X float a, b; X float dmin, dmax; X float r2 = SQR( r ); X int i, face; X X X switch( mode ) X { X case 0: /* Hollow Box - Hollow Sphere */ X dmin = 0; X dmax = 0; X face = FALSE; X for( i = 0; i < n; i++ ) { X a = SQR( C[i] - Bmin[i] ); X b = SQR( C[i] - Bmax[i] ); X dmax += MAX( a, b ); X if( C[i] < Bmin[i] ) { X face = TRUE; X dmin += a; X } X else if( C[i] > Bmax[i] ) { X face = TRUE; X dmin += b; X } X else if( MIN( a, b ) <= r2 ) face = TRUE; X } X if(face && ( dmin <= r2 ) && ( r2 <= dmax)) return(TRUE); X break; X X case 1: /* Hollow Box - Solid Sphere */ X dmin = 0; X face = FALSE; X for( i = 0; i < n; i++ ) { X if( C[i] < Bmin[i] ) { X face = TRUE; X dmin += SQR( C[i] - Bmin[i] ); X } X else if( C[i] > Bmax[i] ) { X face = TRUE; X dmin += SQR( C[i] - Bmax[i] ); X } X else if( C[i] - Bmin[i] <= r ) face = TRUE; X else if( Bmax[i] - C[i] <= r ) face = TRUE; X } X if( face && ( dmin <= r2 ) ) return( TRUE ); X break; X X X case 2: /* Solid Box - Hollow Sphere */ X dmax = 0; X dmin = 0; X for( i = 0; i < n; i++ ) { X a = SQR( C[i] - Bmin[i] ); X b = SQR( C[i] - Bmax[i] ); X dmax += MAX( a, b ); X if( C[i] < Bmin[i] ) dmin += a; else X if( C[i] > Bmax[i] ) dmin += b; X } X if( dmin <= r2 && r2 <= dmax ) return( TRUE ); X break; X X case 3: /* Solid Box - Solid Sphere */ X dmin = 0; X for( i = 0; i < n; i++ ) { X if( C[i] < Bmin[i] ) dmin += SQR(C[i] - Bmin[i] ); else X if( C[i] > Bmax[i] ) dmin += SQR( C[i] - Bmax[i] ); X } X if( dmin <= r2 ) return( TRUE ); X break; X X } /* end switch */ X X return( FALSE ); X } END_OF_FILE if test 3760 -ne `wc -c <'BoxSphere.c'`; then echo shar: \"'BoxSphere.c'\" unpacked with wrong size! fi # end of 'BoxSphere.c' fi if test -f 'Label.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'Label.c'\" else echo shar: Extracting \"'Label.c'\" \(2291 characters\) sed "s/^X//" >'Label.c' <<'END_OF_FILE' X/* X * Nice Numbers for Graph Labels X * by Paul Heckbert X * from "Graphics Gems", Academic Press, 1990 X */ X X/* X * label.c: demonstrate nice graph labeling X * X * Paul Heckbert 2 Dec 88 X */ X X#include X#include X#include "GraphicsGems.h" X Xdouble nicenum(); X X/* expt(a,n)=a^n for integer n */ X X#ifdef POW_NOT_TRUSTWORTHY X/* if roundoff errors in pow cause problems, use this: */ X Xdouble expt(a, n) Xdouble a; Xregister int n; X{ X double x; X X x = 1.; X if (n>0) for (; n>0; n--) x *= a; X else for (; n<0; n++) x /= a; X return x; X} X X#else X# define expt(a, n) pow(a, (double)(n)) X#endif X X#define NTICK 5 /* desired number of tick marks */ X Xmain(ac, av) Xint ac; Xchar **av; X{ X double min, max; X X if (ac!=3) { X fprintf(stderr, "Usage: label \n"); X exit(1); X } X min = atof(av[1]); X max = atof(av[2]); X X loose_label(min, max); X} X X/* X * loose_label: demonstrate loose labeling of data range from min to max. X * (tight method is similar) X */ X Xloose_label(min, max) Xdouble min, max; X{ X char str[6], temp[20]; X int nfrac; X double d; /* tick mark spacing */ X double graphmin, graphmax; /* graph range min and max */ X double range, x; X X /* we expect min!=max */ X range = nicenum(max-min, 0); X d = nicenum(range/(NTICK-1), 1); X graphmin = floor(min/d)*d; X graphmax = ceil(max/d)*d; X nfrac = MAX(-floor(log10(d)), 0); /* # of fractional digits to show */ X sprintf(str, "%%.%df", nfrac); /* simplest axis labels */ X X printf("graphmin=%g graphmax=%g increment=%g\n", graphmin, graphmax, d); X for (x=graphmin; x'Makefile' <<'END_OF_FILE' X# X# Makefile for Graphics Gems source X# X# Craig Kolb, 8/90 X# X# This make file will build "gemslib.a" and a number of executables. X# Gemslib is built solely for debugging purposes -- it is not intended X# to be used as a library. X# X# Note that some of the gems need additional macros, functions, tables X# driving routines, etc. before they will compile or run properly. X# These include: X# X# AALines X# Needs variables set in the Makefile. X# Dissolve X# Needs a table filled. X# FitCurves X# Needs a functioning DrawBezierCurve(). X# MatrixOrtho X# Needs a number of matrix routines. X# PolyScan X# Needs driving routines and other additional code. X# RayPolygon X# Needs surrounding function structure, declaration of X# variable types, etc. X X# X# C compiler flags X# XCFLAGS = -g X# X# Location of Graphics Gems library X# XLIBFILE = gemslib.a X X# X# Graphics Gems Vector Library X# XVECLIB = GGVecLib.o X XMFLAGS = "LIBFILE = ../$(LIBFILE)" "GENCFLAGS = $(CFLAGS)" X XFTPDIR = /usr/spool/uucppublic/ftp/pub/GraphicsGems/src XZOOFILE = Gems.zoo X XSHELL = /bin/sh X XOFILES = PntOnLine.o ViewTrans.o AAPolyScan.o Albers.o \ X Interleave.o BoundSphere.o BoxSphere.o CircleRect.o \ X ConcaveScan.o Roots3And4.o Dissolve.o DigitalLine.o \ X FastJitter.o FixedTrig.o HSLtoRGB.o HypotApprox.o LineEdge.o \ X MatrixInvert.o MatrixPost.o Median.o PixelInteger.o \ X TriPoints.o Quaternions.o RGBTo4Bits.o RayBox.o \ X SeedFill.o SquareRoot.o DoubleLine.o TransBox.o X XDIRS = 2DClip PolyScan Sturm X XALL = Hash3D FitCurves Forms NearestPoint Label \ X OrderDither BinRec $(LIBFILE) X Xall: $(ALL) X @for d in $(DIRS) ; do \ X (cd $$d ; $(MAKE) $(MFLAGS)) ;\ X done X X$(LIBFILE): $(OFILES) $(VECLIB) X ar rcs $(LIBFILE) $(OFILES) $(VECLIB) X XHash3D: Hash3D.o X $(CC) $(CFLAGS) -o $@ Hash3D.o X XFitCurves: FitCurves.o $(VECLIB) X $(CC) $(CFLAGS) -o $@ FitCurves.o $(VECLIB) -lm X XForms: Forms.o X $(CC) $(CFLAGS) -o $@ Forms.o X XNearestPoint: NearestPoint.o $(VECLIB) X $(CC) $(CFLAGS) -o $@ NearestPoint.o $(VECLIB) -lm X XLabel: Label.o X $(CC) $(CFLAGS) -o $@ Label.o -lm X XOrderDither: OrderDither.o X $(CC) $(CFLAGS) -o $@ OrderDither.o X XBinRec: BinRec.o X $(CC) $(CFLAGS) -o $@ BinRec.o X Xftpreadme: X /bin/rm -f $(FTPDIR)/README X sed -e "s/__DATE__/`date`/g" < README.ftp > $(FTPDIR)/README X Xreadme: X /bin/rm -rf README X sed -e "s/__DATE__/`date`/g" < README.dist > README X X Xclean: X @for d in $(DIRS) ; do \ X (cd $$d ; $(MAKE) $(MFLAGS) clean) ;\ X done X /bin/rm -f $(OFILES) $(VECLIB) X /bin/rm -f Hash3D.o FitCurves.o \ X Forms.o NearestPoint.o Label.o OrderDither.o BinRec.o \ X Hash3D FitCurves Forms NearestPoint Label OrderDither BinRec \ X bugs a.out core Part??.Z Part?? $(ZOOFILE) X Xkit: readme X /bin/rm -f Part??.Z Part?? X (makekit -iPACKING_LIST -oMANIFEST) X Xzoo: readme X /bin/rm -f Gems.zoo X cut -d" " -f2 PACKING_LIST | zoo aI $(ZOOFILE) X Xftp: kit zoo ftpreadme X compress Part?? X /bin/cp Part??.Z $(ZOOFILE) $(FTPDIR) X X$(ALL): GraphicsGems.h END_OF_FILE if test 2894 -ne `wc -c <'Makefile'`; then echo shar: \"'Makefile'\" unpacked with wrong size! fi # end of 'Makefile' fi if test -f 'MatrixPost.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'MatrixPost.c'\" else echo shar: Extracting \"'MatrixPost.c'\" \(3390 characters\) sed "s/^X//" >'MatrixPost.c' <<'END_OF_FILE' X/* XEfficient Post-Concatenation of Transformation Matrics Xby Joseph M. Cychosz Xfrom "Graphics Gems", Academic Press, 1990 X*/ X X#include X#include "GraphicsGems.h" X X X/* M4xform.c - Basic 4x4 transformation package. X * X * Description: X * M4xform.c contains a collection of routines used to perform X * direct post-concatenated transformation operations. X * X * X * Contents: X * M4RotateX Post-concatenate a x-axis rotation. X * M4RotateY Post-concatenate a y-axis rotation. X * M4RotateZ Post-concatenate a z-axis rotation. X * M4Scale Post-concatenate a scaling. X * M4Translate Post-concatenate a translation. X * M4ZPerspective Post-concatenate a z-axis perspective X * transformation. X * X * Externals: X * cos, sin. X */ X X X X/* M4RotateX - Post-concatenate a x-axis rotation matrix. */ X XMatrix4 *M4RotateX (m,a) X Matrix4 *m; /* Current transformation matrix*/ X double a; /* Rotation angle */ X{ X double c, s; X double t; X int i; X X X c = cos (a); s = sin (a); X X for (i = 0 ; i < 4 ; i++) { X t = m->element[i][1]; X m->element[i][1] = t*c - m->element[i][2]*s; X m->element[i][2] = t*s + m->element[i][2]*c; X } X return (m); X} X X X/* M4RotateY - Post-concatenate a y-axis rotation matrix. */ X XMatrix4 *M4RotateY (m,a) X Matrix4 *m; /* Current transformation matrix*/ X double a; /* Rotation angle */ X{ X double c, s; X double t; X int i; X X c = cos (a); s = sin (a); X X for (i = 0 ; i < 4 ; i++) { X t = m->element[i][0]; X m->element[i][0] = t*c + m->element[i][2]*s; X m->element[i][2] = m->element[i][2]*c - t*s; X } X return (m); X} X X X/* M4RotateZ - Post-concatenate a z-axis rotation matrix. */ X XMatrix4 *M4RotateZ (m,a) X Matrix4 *m; /* Current transformation matrix*/ X double a; /* Rotation angle */ X{ X double c, s; X double t; X int i; X X c = cos (a); s = sin (a); X X for (i = 0 ; i < 4 ; i++) { X t = m->element[i][0]; X m->element[i][0] = t*c - m->element[i][1]*s; X m->element[i][1] = t*s + m->element[i][1]*c; X } X return (m); X} X X X X/* M4Scale - Post-concatenate a scaling. */ X XMatrix4 *M4Scale (m,sx,sy,sz) X Matrix4 *m; /* Current transformation matrix */ X double sx, sy, sz; /* Scale factors about x,y,z */ X{ X int i; X X for (i = 0 ; i < 4 ; i++) { X m->element[i][0] *= sx; X m->element[i][1] *= sy; X m->element[i][2] *= sz; X } X return (m); X} X X X/* M4Translate - Post-concatenate a translation. */ X XMatrix4 *M4Translate (m,tx,ty,tz) X Matrix4 *m; /* Current transformation matrix */ X double tx, ty, tz; /* Translation distance */ X{ X int i; X X for (i = 0 ; i < 4 ; i++) { X m->element[i][0] += m->element[i][3]*tx; X m->element[i][1] += m->element[i][3]*ty; X m->element[i][2] += m->element[i][3]*tz; X } X return (m); X} X X X /* M4ZPerspective Post-concatenate a perspective */ /*transformation. */ X /* */ X /* Perspective is along the z-axis with the eye at +z. */ X X XMatrix4 *M4ZPerspective (m,d) X Matrix4 *m; /* Current transformation matrix */ X double d; /* Perspective distance */ X{ X int i; X double f = 1. / d; X X for (i = 0 ; i < 4 ; i++) { X m->element[i][3] += m->element[i][2]*f; X m->element[i][2] = 0.; X } X return (m); X} X X X X END_OF_FILE if test 3390 -ne `wc -c <'MatrixPost.c'`; then echo shar: \"'MatrixPost.c'\" unpacked with wrong size! fi # end of 'MatrixPost.c' fi if test -f 'Median.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'Median.c'\" else echo shar: Extracting \"'Median.c'\" \(2941 characters\) sed "s/^X//" >'Median.c' <<'END_OF_FILE' X/* XMedian Finding on a 3-by-3 Grid Xby Alan Paeth Xfrom "Graphics Gems", Academic Press, 1990 X*/ X X#define s2(a,b) {register int t; if ((t=b-a)<0) {a+=t; b-=t;}} X#define mn3(a,b,c) s2(a,b); s2(a,c); X#define mx3(a,b,c) s2(b,c); s2(a,c); X#define mnmx3(a,b,c) mx3(a,b,c); s2(a,b); X#define mnmx4(a,b,c,d) s2(a,b); s2(c,d); s2(a,c); s2(b,d); X#define mnmx5(a,b,c,d,e) s2(a,b); s2(c,d); mn3(a,c,e); mx3(b,d,e); X#define mnmx6(a,b,c,d,e,f) s2(a,d); s2(b,e); s2(c,f);\ X mn3(a,b,c); mx3(d,e,f); Xmed3x3(b1, b2, b3) X int *b1, *b2, *b3; X/* X * Find median on a 3x3 input box of integers. X * b1, b2, b3 are pointers to the left-hand edge of three X * parallel scan-lines to form a 3x3 spatial median. X * Rewriting b2 and b3 as b1 yields code which forms median X * on input presented as a linear array of nine elements. X */ X { X register int r1, r2, r3, r4, r5, r6; X r1 = *b1++; r2 = *b1++; r3 = *b1++; X r4 = *b2++; r5 = *b2++; r6 = *b2++; X mnmx6(r1, r2, r3, r4, r5, r6); X r1 = *b3++; X mnmx5(r1, r2, r3, r4, r5); X r1 = *b3++; X mnmx4(r1, r2, r3, r4); X r1 = *b3++; X mnmx3(r1, r2, r3); X return(r2); X } X X X/* t2(i,j) transposes elements in A[] such that A[i] <= A[j] */ X X#define t2(i, j) s2(A[i-1], A[j-1]) X X Xint median25(A) X int A[25]; X { X/* X * median25(A) partitions the array A[0..24] such that element X * A[12] is the median and subarrays A[0..11] and A[13..24] X * are partitions containing elements of smaller and larger X * value (rank), respectively. X * X * The exchange table lists element indices on the range 1..25, X * this accounts for the "-1" offsets in the macro t2 and in X * the final return value used to adjust subscripts to C-code X * conventions (array indices begin at zero). X */ X t2( 1, 2); t2( 4, 5); t2( 3, 5); t2( 3, 4); t2( 7, 8); X t2( 6, 8); t2( 6, 7); t2(10,11); t2( 9,11); t2( 9,10); X t2(13,14); t2(12,14); t2(12,13); t2(16,17); t2(15,17); X t2(15,16); t2(19,20); t2(18,20); t2(18,19); t2(22,23); X t2(21,23); t2(21,22); t2(24,25); t2( 3, 6); t2( 4, 7); X t2( 1, 7); t2( 1, 4); t2( 5, 8); t2( 2, 8); t2( 2, 5); X t2(12,15); t2( 9,15); t2( 9,12); t2(13,16); t2(10,16); X t2(10,13); t2(14,17); t2(11,17); t2(11,14); t2(21,24); X t2(18,24); t2(18,21); t2(22,25); t2(19,25); t2(19,22); X t2(20,23); t2( 9,18); t2(10,19); t2( 1,19); t2( 1,10); X t2(11,20); t2( 2,20); t2( 2,11); t2(12,21); t2( 3,21); X t2( 3,12); t2(13,22); t2( 4,22); t2( 4,13); t2(14,23); X t2( 5,23); t2( 5,14); t2(15,24); t2( 6,24); t2( 6,15); X t2(16,25); t2( 7,25); t2( 7,16); t2( 8,17); t2( 8,20); X t2(14,22); t2(16,24); t2( 8,14); t2( 8,16); t2( 2,10); X t2( 4,12); t2( 6,18); t2(12,18); t2(10,18); t2( 5,11); X t2( 7,13); t2( 8,15); t2( 5, 7); t2( 5, 8); t2(13,15); X t2(11,15); t2( 7, 8); t2(11,13); t2( 7,11); t2( 7,18); X t2(13,18); t2( 8,18); t2( 8,11); t2(13,19); t2( 8,13); X t2(11,19); t2(13,21); t2(11,21); t2(11,13); X return(A[13-1]); X } END_OF_FILE if test 2941 -ne `wc -c <'Median.c'`; then echo shar: \"'Median.c'\" unpacked with wrong size! fi # end of 'Median.c' fi if test -f 'OrderDither.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'OrderDither.c'\" else echo shar: Extracting \"'OrderDither.c'\" \(2483 characters\) sed "s/^X//" >'OrderDither.c' <<'END_OF_FILE' X/* XOrdered Dithering Xby Stephen Hawley Xfrom "Graphics Gems", Academic Press, 1990 X*/ X X/* Program to generate dithering matrices. X * written by Jim Blandy, Oberlin College, jimb@occs.oberlin.edu X * Gifted to, documented and revised by Stephen Hawley, X * sdh@flash.bellcore.com X * X * Generates a dithering matrix from the command line arguments. X * The first argument, size, determines the dimensions of the X * matrix: 2^size by 2^size X * The optional range argument is the range of values to be X * dithered over. By default, it is (2^size)^2, or simply the X * total number of elements in the matrix. X * The final output is suitable for inclusion in a C program. X * A typical dithering function is something like this: X * extern int dm[], size; X * X * int X * dither(x,y, level) X * register int x,y, level; X * { X * return(level > dm[(x % size) + size * (y % size)]; X * } X */ X Xmain(argc, argv) Xint argc; Xchar **argv; X{ X register int size, range; X X if (argc >= 2) size = atoi(argv[1]); X else size = 2; X X if (argc == 3) range = atoi(argv[2]); X else range = (1 << size) * (1 << size); X X printdither (size, range); X} X X Xprintdither (size, range) Xregister int size, range; X{ X register int l = (1 << size), i; X /* X * print a dithering matrix. X * l is the length on a side. X */ X range = range / (l * l); X puts("int dm[] = {"); X for (i=0; i < l*l; i++) { X if (i % l == 0) /* tab in 4 spaces per row */ X printf(" "); X /* print the dither value for this location X * scaled to the given range X */ X printf("%4d", range * dithervalue(i / l, i % l, size)); X X /* commas after all but the last */ X if (i + 1 < l * l) X putchar(','); X /* newline at the end of the row */ X if ((i + 1) % l == 0) X putchar('\n'); X } X puts("\n}; "); X} Xdithervalue(x, y, size) Xregister int x, y, size; X{ X register int d; X /* X * calculate the dither value at a particular X * (x, y) over the size of the matrix. X */ X d=0; X while (size-->0) { X /* Think of d as the density. At every iteration, X * d is shifted left one and a new bit is put in the X * low bit based on x and y. If x is odd and y is even, X * or x is even and y is odd, a bit is put in. This X * generates the checkerboard seen in dithering. X * This quantity is shifted left again and the low bit of X * y is added in. X * This whole thing interleaves a checkerboard bit pattern X * and y's bits, which is the value you want. X */ X d = (d <<1 | (x&1 ^ y&1))<<1 | y&1; X x >>= 1; X y >>= 1; X } X return(d); X} X END_OF_FILE if test 2483 -ne `wc -c <'OrderDither.c'`; then echo shar: \"'OrderDither.c'\" unpacked with wrong size! fi # end of 'OrderDither.c' fi if test -f 'PntOnLine.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'PntOnLine.c'\" else echo shar: Extracting \"'PntOnLine.c'\" \(2024 characters\) sed "s/^X//" >'PntOnLine.c' <<'END_OF_FILE' X/* XA Fast 2D Point-On-Line Test Xby Alan Paeth Xfrom "Graphics Gems", Academic Press, 1990 X*/ X X#include "GraphicsGems.h" X Xint PntOnLine(px,py,qx,qy,tx,ty) X int px, py, qx, qy, tx, ty; X { X/* X * given a line through P:(px,py) Q:(qx,qy) and T:(tx,ty) X * return 0 if T is not on the line through <--P--Q--> X * 1 if T is on the open ray ending at P: <--P X * 2 if T is on the closed interior along: P--Q X * 3 if T is on the open ray beginning at Q: Q--> X * X * Example: consider the line P = (3,2), Q = (17,7). A plot X * of the test points T(x,y) (with 0 mapped onto '.') yields: X * X * 8| . . . . . . . . . . . . . . . . . 3 3 X * Y 7| . . . . . . . . . . . . . . 2 2 Q 3 3 Q = 2 X * 6| . . . . . . . . . . . 2 2 2 2 2 . . . X * a 5| . . . . . . . . 2 2 2 2 2 2 . . . . . X * x 4| . . . . . 2 2 2 2 2 2 . . . . . . . . X * i 3| . . . 2 2 2 2 2 . . . . . . . . . . . X * s 2| 1 1 P 2 2 . . . . . . . . . . . . . . P = 2 X * 1| 1 1 . . . . . . . . . . . . . . . . . X * +-------------------------------------- X * 1 2 3 4 5 X-axis 10 15 19 X * X * Point-Line distance is normalized with the Infinity Norm X * avoiding square-root code and tightening the test vs the X * Manhattan Norm. All math is done on the field of integers. X * The latter replaces the initial ">= MAX(...)" test with X * "> (ABS(qx-px) + ABS(qy-py))" loosening both inequality X * and norm, yielding a broader target line for selection. X * The tightest test is employed here for best discrimination X * in merging collinear (to grid coordinates) vertex chains X * into a larger, spanning vectors within the Lemming editor. X */ X X X if ( ABS((qy-py)*(tx-px)-(ty-py)*(qx-px)) >= X (MAX(ABS(qx-px), ABS(qy-py)))) return(0); X if (((qx'PolyScan/fancytest.c' <<'END_OF_FILE' X/* X * fancytest.c: subroutine illustrating the use of poly_clip and poly_scan X * for Phong-shading and texture mapping. X * X * Note: lines enclosed in angle brackets '<', '>' should be replaced X * with the code described. X * Makes calls to hypothetical packages "shade", "image", "texture", "zbuffer". X * X * Paul Heckbert Dec 1989 X */ X X#include X#include X#include "poly.h" X X#define XMAX 1280 /* hypothetical image width */ X#define YMAX 1024 /* hypothetical image height */ X#define LIGHT_INTEN 255 /* light source intensity */ X Xvoid pixelproc(); X Xfancytest() X{ X int i; X double WS[4][4]; /* world space to screen space transform */ X Poly p; /* a polygon */ X Poly_vert *v; X X static Poly_box box = {0, XMAX, 0, YMAX, -32678, 32767.99}; X /* 3-D screen clipping box, continuous coordinates */ X X static Window win = {0, 0, XMAX-1, YMAX-1}; X /* 2-D screen clipping window, discrete coordinates */ X X X X X /* transform vertices from world space to homogeneous screen space */ X for (i=0; ix, v->y, v->z, 1., WS, &v->sx, &v->sy, &v->sz, &v->sw); X } X X /* interpolate sx, sy, sz, sw, nx, ny, nz, u, v in poly_clip */ X p.mask = POLY_MASK(sx) | POLY_MASK(sy) | POLY_MASK(sz) | POLY_MASK(sw) | X POLY_MASK(nx) | POLY_MASK(ny) | POLY_MASK(nz) | X POLY_MASK(u) | POLY_MASK(v); X X poly_print("before clipping", &p); X if (poly_clip_to_box(&p, &box) == POLY_CLIP_OUT) /* clip polygon */ X return; /* quit if off-screen */ X X /* do homogeneous division of screen position and texture position */ X for (i=0; isx /= v->sw; X v->sy /= v->sw; X v->sz /= v->sw; X v->u /= v->sw; X v->v /= v->sw; X v->q = 1./v->sw; X } X /* X * previously we ignored q (assumed q=1), henceforth ignore sw (assume sw=1) X * Interpolate sx, sy, sz, nx, ny, nz, u, v, q in poly_scan X */ X p.mask &= ~POLY_MASK(sw); X p.mask |= POLY_MASK(q); X X poly_print("scan converting the polygon", &p); X poly_scan(&p, &win, pixelproc); /* scan convert! */ X} X Xstatic void pixelproc(x, y, pt) /* called at each pixel by poly_scan */ Xint x, y; XPoly_vert *pt; X{ X int sz, u, v, inten; X double len, nor[3], diff, spec; X X sz = pt->sz; X if (sz < zbuffer_read(x, y)) { X len = sqrt(pt->nx*pt->nx + pt->ny*pt->ny + pt->nz*pt->nz); X nor[0] = pt->nx/len; /* unitize the normal vector */ X nor[1] = pt->ny/len; X nor[2] = pt->nz/len; X shade(nor, &diff, &spec); /* compute specular and diffuse coeffs*/ X u = pt->u/pt->q; /* do homogeneous div. of texture pos */ X v = pt->v/pt->q; X inten = texture_read(u, v)*diff + LIGHT_INTEN*spec; X image_write(x, y, inten); X zbuffer_write(x, y, sz); X } X} X X/* mx4_transform: transform 4-vector p by matrix m yielding q: q = p*m */ X Xmx4_transform(px, py, pz, pw, m, qxp, qyp, qzp, qwp) Xdouble px, py, pz, pw, m[4][4], *qxp, *qyp, *qzp, *qwp; X{ X *qxp = px*m[0][0] + py*m[1][0] + pz*m[2][0] + pw*m[3][0]; X *qyp = px*m[0][1] + py*m[1][1] + pz*m[2][1] + pw*m[3][1]; X *qzp = px*m[0][2] + py*m[1][2] + pz*m[2][2] + pw*m[3][2]; X *qwp = px*m[0][3] + py*m[1][3] + pz*m[2][3] + pw*m[3][3]; X} END_OF_FILE if test 3298 -ne `wc -c <'PolyScan/fancytest.c'`; then echo shar: \"'PolyScan/fancytest.c'\" unpacked with wrong size! fi # end of 'PolyScan/fancytest.c' fi if test -f 'PolyScan/poly.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'PolyScan/poly.c'\" else echo shar: Extracting \"'PolyScan/poly.c'\" \(2310 characters\) sed "s/^X//" >'PolyScan/poly.c' <<'END_OF_FILE' X/* X * poly.c: simple utilities for polygon data structures X */ X X#include "poly.h" X XPoly_vert *poly_dummy; /* used superficially by POLY_MASK macro */ X X/* X * poly_print: print Poly p to stdout, prefixed by the label str X */ X Xvoid poly_print(str, p) Xchar *str; XPoly *p; X{ X int i; X X printf("%s: %d sides\n", str, p->n); X poly_vert_label(" ", p->mask); X for (i=0; in; i++) { X printf(" v[%d] ", i); X poly_vert_print("", &p->vert[i], p->mask); X } X} X Xvoid poly_vert_label(str, mask) Xchar *str; Xint mask; X{ X printf("%s", str); X if (mask&POLY_MASK(sx)) printf(" sx "); X if (mask&POLY_MASK(sy)) printf(" sy "); X if (mask&POLY_MASK(sz)) printf(" sz "); X if (mask&POLY_MASK(sw)) printf(" sw "); X if (mask&POLY_MASK(x)) printf(" x "); X if (mask&POLY_MASK(y)) printf(" y "); X if (mask&POLY_MASK(z)) printf(" z "); X if (mask&POLY_MASK(u)) printf(" u "); X if (mask&POLY_MASK(v)) printf(" v "); X if (mask&POLY_MASK(q)) printf(" q "); X if (mask&POLY_MASK(r)) printf(" r "); X if (mask&POLY_MASK(g)) printf(" g "); X if (mask&POLY_MASK(b)) printf(" b "); X if (mask&POLY_MASK(nx)) printf(" nx "); X if (mask&POLY_MASK(ny)) printf(" ny "); X if (mask&POLY_MASK(nz)) printf(" nz "); X printf("\n"); X} X Xvoid poly_vert_print(str, v, mask) Xchar *str; XPoly_vert *v; Xint mask; X{ X printf("%s", str); X if (mask&POLY_MASK(sx)) printf(" %6.1f", v->sx); X if (mask&POLY_MASK(sy)) printf(" %6.1f", v->sy); X if (mask&POLY_MASK(sz)) printf(" %6.2f", v->sz); X if (mask&POLY_MASK(sw)) printf(" %6.2f", v->sw); X if (mask&POLY_MASK(x)) printf(" %6.2f", v->x); X if (mask&POLY_MASK(y)) printf(" %6.2f", v->y); X if (mask&POLY_MASK(z)) printf(" %6.2f", v->z); X if (mask&POLY_MASK(u)) printf(" %6.2f", v->u); X if (mask&POLY_MASK(v)) printf(" %6.2f", v->v); X if (mask&POLY_MASK(q)) printf(" %6.2f", v->q); X if (mask&POLY_MASK(r)) printf(" %6.4f", v->r); X if (mask&POLY_MASK(g)) printf(" %6.4f", v->g); X if (mask&POLY_MASK(b)) printf(" %6.4f", v->b); X if (mask&POLY_MASK(nx)) printf(" %6.3f", v->nx); X if (mask&POLY_MASK(ny)) printf(" %6.3f", v->ny); X if (mask&POLY_MASK(nz)) printf(" %6.3f", v->nz); X printf("\n"); X} END_OF_FILE if test 2310 -ne `wc -c <'PolyScan/poly.c'`; then echo shar: \"'PolyScan/poly.c'\" unpacked with wrong size! fi # end of 'PolyScan/poly.c' fi if test -f 'PolyScan/poly.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'PolyScan/poly.h'\" else echo shar: Extracting \"'PolyScan/poly.h'\" \(1914 characters\) sed "s/^X//" >'PolyScan/poly.h' <<'END_OF_FILE' X/* poly.h: definitions for polygon package */ X X#ifndef POLY_HDR X#define POLY_HDR X X#define POLY_NMAX 8 /* max #sides to a polygon; change if needed */ X Xtypedef struct { /* A POLYGON VERTEX */ X double sx, sy, sz, sw; /* screen space position (sometimes homo.) */ X double x, y, z; /* world space position */ X double u, v, q; /* texture position (sometimes homogeneous) */ X double r, g, b; /* (red,green,blue) color */ X double nx, ny, nz; /* world space normal vector */ X} Poly_vert; X/* update poly.c if you change this structure */ X Xtypedef struct { /* A POLYGON */ X int n; /* number of sides */ X int mask; /* interpolation mask for vertex elems */ X Poly_vert vert[POLY_NMAX]; /* vertices */ X} Poly; X/* X * mask is an interpolation mask whose kth bit indicates whether the kth X * double in a Poly_vert is relevant. X * For example, if the valid attributes are sx, sy, and sz, then set X * mask = POLY_MASK(sx) | POLY_MASK(sy) | POLY_MASK(sz); X */ X Xtypedef struct { /* A BOX (TYPICALLY IN SCREEN SPACE) */ X double x0, x1; /* left and right */ X double y0, y1; /* top and bottom */ X double z0, z1; /* near and far */ X} Poly_box; X Xtypedef struct { /* WINDOW: A DISCRETE 2-D RECTANGLE */ X int x0, y0; /* xmin and ymin */ X int x1, y1; /* xmax and ymax (inclusive) */ X} Window; X X#define POLY_MASK(elem) (1 << (&poly_dummy->elem - (double *)poly_dummy)) X X#define POLY_CLIP_OUT 0 /* polygon entirely outside box */ X#define POLY_CLIP_PARTIAL 1 /* polygon partially inside */ X#define POLY_CLIP_IN 2 /* polygon entirely inside box */ X Xextern Poly_vert *poly_dummy; /* used superficially by POLY_MASK macro */ X Xvoid poly_print(/* str, p */); Xvoid poly_vert_label(/* str, mask */); Xvoid poly_vert_print(/* str, v, mask */); Xint poly_clip_to_box(/* p1, box */); Xvoid poly_clip_to_halfspace(/* p, q, index, sign, k, name */); Xvoid poly_scan(/* p, win, pixelproc */); X X#endif END_OF_FILE if test 1914 -ne `wc -c <'PolyScan/poly.h'`; then echo shar: \"'PolyScan/poly.h'\" unpacked with wrong size! fi # end of 'PolyScan/poly.h' fi if test -f 'PolyScan/scantest.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'PolyScan/scantest.c'\" else echo shar: Extracting \"'PolyScan/scantest.c'\" \(1854 characters\) sed "s/^X//" >'PolyScan/scantest.c' <<'END_OF_FILE' X/* X * scantest.c: use poly_scan() for Gouraud shading and z-buffer demo. X * Given the screen space X, Y, and Z of N-gon on command line, X * print out all pixels during scan conversion. X * This code could easily be modified to actually read and write pixels. X * X * Paul Heckbert Dec 1989 X */ X X#include X#include X#include "poly.h" X X#define XMAX 1280 /* hypothetical image width */ X#define YMAX 1024 /* hypothetical image height */ X X#define FRANDOM() ((rand()&32767)/32767.) /* random number between 0 and 1 */ X Xstatic void pixelproc(); X Xmain(ac, av) Xint ac; Xchar **av; X{ X int i; X Poly p; X static Window win = {0, 0, XMAX-1, YMAX-1}; /* screen clipping window */ X X if (ac<2 || ac != 2+3*(p.n = atoi(av[1]))) { X fprintf(stderr, "Usage: scantest N X1 Y1 Z1 X2 Y2 Z2 ... XN YN ZN\n"); X exit(1); X } X for (i=0; isz, point->r, point->g, point->b); X X /* X * in real graphics program you could read and write pixels, e.g.: X * X * if (point->sz < zbuffer_read(x, y)) { X * image_write_rgb(x, y, point->r, point->g, point->b); X * zbuffer_write(x, y, point->sz); X * } X */ X} END_OF_FILE if test 1854 -ne `wc -c <'PolyScan/scantest.c'`; then echo shar: \"'PolyScan/scantest.c'\" unpacked with wrong size! fi # end of 'PolyScan/scantest.c' fi if test -f 'Quaternions.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'Quaternions.c'\" else echo shar: Extracting \"'Quaternions.c'\" \(2216 characters\) sed "s/^X//" >'Quaternions.c' <<'END_OF_FILE' X/* XUsing Quaternions for Coding 3D Transformations XPatrick-Gilles Maillot Xfrom "Graphics Gems", Academic Press, 1990 X*/ X Xextern double P[3], Q[4], M[4][4]; X Xset_obs_position(x,y,z) Xfloat x, y, z; X{ Xint i; X X/* X * Set the values of the eye's position. X * The position here represents the position of the orthonormal base X * in respect to the observer. X */ X P[0] = -x; X P[1] = -y; X P[2] = -z; X/* X * Set the visualization to be in the decreasing x axis X */ X Q[0] = 1.; X for (i = 1; i < 4; i++) Q[i] = 0.; X} X Xtranslate_quaternion(x,i,w) Xfloat x; Xint i, w; X{ Xint j, k; Xfloat A, B, D, E, F; X X if (w < 0) { X/* X * The observer moves in respect to the scene. X */ X P[i - 1] -= x; X } else { X X/* X * The scene moves in respect to the observer. X * Compute the successor axis of i [1,2,3]; X * and then the successor axis of j [1,2,3]; X */ X if ((j = i + 1) > 3) j = 1; X if ((k = j + 1) > 3) k = 1; X A = Q[j]; B = Q[k]; F = Q[0]; E = Q[i]; X P[i - 1] += x * (E * E + F * F - A * A - B * B); X D = x + x; X P[j - 1] += D * (E * A + F * B); X P[k - 1] += D * (E * B + F * A); X } X} X Xrotate_quaternion(x,y,i,w) Xfloat x, y; Xint i, w; X{ Xint j, k; Xfloat E, F, R1; X/* X * Compute the successor axis of i [1,2,3] and j [1,2,3]; X */ X if ((j = i + 1) > 3) j = 1; X if ((k = j + 1) > 3) k = 1; X E = Q[i]; X Q[i] = E * x + w * y * Q[0]; X Q[0] = Q[0] * x - w * y * E; X E = Q[j]; X Q[j] = E * x + y * Q[k]; X Q[k] = Q[k] * x - y * E; X if (w < 0) { X/* Compute a new position if the observer moves in respect to the scene. */ X j -= 1; k -= 1; X R1 = x * x - y * y; X F = 2. * x * y; X E = P[j]; X P[j] = E * R1 + F * P[k]; X P[k] = P[k] * R1 - F * E; X } X} X X XEvaluate_matrix() X{ Xfloat e, f, r[4]; Xint i, j, k; X/* X * We will need some square values! X */ X for (i = 0; i < 4; i++) r[i] = Q[i] * Q[i]; X/* X * Compute each element of the matrix. X * j is the successor of i (in 1,2,3), while k is the successor of j. X */ X for (i = 1; i < 4; i++) { X if ((j = i + 1) > 3) j = 1; X if ((k = j + 1) > 3) k = 1; X e = 2. * Q[i] * Q[j]; X f = 2. * Q[k] * Q[0]; X M[j][i] = e - f; X M[i][j] = e + f; X M[i][i] = r[i] + r[0] - r[j] - r[k]; X M[0][i] = P[i - 1]; X M[i][0] = 0.; X } X M[0][0] = 1.; X} X X X X END_OF_FILE if test 2216 -ne `wc -c <'Quaternions.c'`; then echo shar: \"'Quaternions.c'\" unpacked with wrong size! fi # end of 'Quaternions.c' fi if test -f 'SeedFill.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'SeedFill.c'\" else echo shar: Extracting \"'SeedFill.c'\" \(2371 characters\) sed "s/^X//" >'SeedFill.c' <<'END_OF_FILE' X/* X * A Seed Fill Algorithm X * by Paul Heckbert X * from "Graphics Gems", Academic Press, 1990 X */ X X/* X * fill.c : simple seed fill program X * Calls pixelread() to read pixels, pixelwrite() to write pixels. X * X * Paul Heckbert 13 Sept 1982, 28 Jan 1987 X */ X Xtypedef struct { /* window: a discrete 2-D rectangle */ X int x0, y0; /* xmin and ymin */ X int x1, y1; /* xmax and ymax (inclusive) */ X} Window; X Xtypedef int Pixel; /* 1-channel frame buffer assumed */ X XPixel pixelread(); X Xtypedef struct {short y, xl, xr, dy;} Segment; X/* X * Filled horizontal segment of scanline y for xl<=x<=xr. X * Parent segment was on line y-dy. dy=1 or -1 X */ X X#define MAX 10000 /* max depth of stack */ X X#define PUSH(Y, XL, XR, DY) /* push new segment on stack */ \ X if (sp=win->y0 && Y+(DY)<=win->y1) \ X {sp->y = Y; sp->xl = XL; sp->xr = XR; sp->dy = DY; sp++;} X X#define POP(Y, XL, XR, DY) /* pop segment off stack */ \ X {sp--; Y = sp->y+(DY = sp->dy); XL = sp->xl; XR = sp->xr;} X X/* X * fill: set the pixel at (x,y) and all of its 4-connected neighbors X * with the same pixel value to the new pixel value nv. X * A 4-connected neighbor is a pixel above, below, left, or right of a pixel. X */ X Xfill(x, y, win, nv) Xint x, y; /* seed point */ XWindow *win; /* screen window */ XPixel nv; /* new pixel value */ X{ X int l, x1, x2, dy; X Pixel ov; /* old pixel value */ X Segment stack[MAX], *sp = stack; /* stack of filled segments */ X X ov = pixelread(x, y); /* read pv at seed point */ X if (ov==nv || xx0 || x>win->x1 || yy0 || y>win->y1) return; X PUSH(y, x, x, 1); /* needed in some cases */ X PUSH(y+1, x, x, -1); /* seed segment (popped 1st) */ X X while (sp>stack) { X /* pop segment off stack and fill a neighboring scan line */ X POP(y, x1, x2, dy); X /* X * segment of scan line y-dy for x1<=x<=x2 was previously filled, X * now explore adjacent pixels in scan line y X */ X for (x=x1; x>=win->x0 && pixelread(x, y)==ov; x--) X pixelwrite(x, y, nv); X if (x>=x1) goto skip; X l = x+1; X if (lx1 && pixelread(x, y)==ov; x++) X pixelwrite(x, y, nv); X PUSH(y, l, x-1, dy); X if (x>x2+1) PUSH(y, x2+1, x-1, -dy); /* leak on right? */ Xskip: for (x++; x<=x2 && pixelread(x, y)!=ov; x++); X l = x; X } while (x<=x2); X } X} END_OF_FILE if test 2371 -ne `wc -c <'SeedFill.c'`; then echo shar: \"'SeedFill.c'\" unpacked with wrong size! fi # end of 'SeedFill.c' fi if test -f 'SquareRoot.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'SquareRoot.c'\" else echo shar: Extracting \"'SquareRoot.c'\" \(2048 characters\) sed "s/^X//" >'SquareRoot.c' <<'END_OF_FILE' X/* X * A High Speed, Low Precision Square Root X * by Paul Lalonde and Robert Dawson X * from "Graphics Gems", Academic Press, 1990 X */ X X X/* X * SPARC implementation of a fast square root by table X * lookup. X * SPARC floating point format is as follows: X * X * BIT 31 30 23 22 0 X * sign exponent mantissa X */ X X#include X Xstatic short sqrttab[0x100]; /* declare table of square roots */ X Xvoid Xbuild_table() X{ X unsigned short i; X float f; X unsigned int *fi = &f; /* to access the bits of a float in */ X /* C quickly we must misuse pointers */ X X for(i=0; i<= 0x7f; i++) { X *fi = 0; X X /* X * Build a float with the bit pattern i as mantissa X * and an exponent of 0, stored as 127 X */ X X *fi = (i << 16) | (127 << 23); X f = sqrt(f); X X /* X * Take the square root then strip the first 7 bits of X * the mantissa into the table X */ X X sqrttab[i] = (*fi & 0x7fffff) >> 16; X X /* X * Repeat the process, this time with an exponent of X * 1, stored as 128 X */ X X *fi = 0; X *fi = (i << 16) | (128 << 23); X f = sqrt(f); X sqrttab[i+0x80] = (*fi & 0x7fffff) >> 16; X } X} X X/* X * fsqrt - fast square root by table lookup X */ Xfloat Xfsqrt(float n) X{ X unsigned int *num = &n; /* to access the bits of a float in C X * we must misuse pointers */ X X short e; /* the exponent */ X if (n == 0) return (0); /* check for square root of 0 */ X e = (*num >> 23) - 127; /* get the exponent - on a SPARC the */ X /* exponent is stored with 127 added */ X *num & = 0x7fffff; /* leave only the mantissa */ X if (e & 0x01) *num | = 0x800000; X /* the exponent is odd so we have to */ X /* look it up in the second half of */ X /* the lookup table, so we set the */ X /* high bit */ X e >>= 1; /* divide the exponent by two */ X /* note that in C the shift */ X /* operators are sign preserving */ X /* for signed operands */ X /* Do the table lookup, based on the quaternary mantissa, X * then reconstruct the result back into a float X */ X *num = ((sqrttab[*num >> 16]) << 16) | ((e + 127) << 23); X return(n); X} END_OF_FILE if test 2048 -ne `wc -c <'SquareRoot.c'`; then echo shar: \"'SquareRoot.c'\" unpacked with wrong size! fi # end of 'SquareRoot.c' fi if test -f 'Sturm/main.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'Sturm/main.c'\" else echo shar: Extracting \"'Sturm/main.c'\" \(2173 characters\) sed "s/^X//" >'Sturm/main.c' <<'END_OF_FILE' X/* XUsing Sturm Sequences to Bracket Real Roots of Polynomial Equations Xby D.G. Hook and P.R. McAree Xfrom "Graphics Gems", Academic Press, 1990 X*/ X X/* X * main.c X * X * a sample driver program. X */ X#include X#include X#include "solve.h" X X/* X * a driver program for a root solver. X */ Xmain() X{ X poly sseq[MAX_ORDER]; X double min, max, roots[MAX_ORDER]; X int i, j, order, nroots, nchanges, np, atmin, atmax; X X /* X * get the details... X */ X X printf("Please enter order of polynomial: "); X scanf("%d", &order); X X printf("\n"); X X for (i = order; i >= 0; i--) { X printf("Please enter coefficient number %d: ", i); X scanf("%lf", &sseq[0].coef[i]); X } X X printf("\n"); X X /* X * build the Sturm sequence X */ X np = buildsturm(order, sseq); X X printf("Sturm sequence for:\n"); X X for (i = order; i >= 0; i--) X printf("%lf ", sseq[0].coef[i]); X X printf("\n\n"); X X for (i = 0; i <= np; i++) { X for (j = sseq[i].ord; j >= 0; j--) X printf("%lf ", sseq[i].coef[j]); X printf("\n"); X } X X printf("\n"); X X X /* X * get the number of real roots X */ X nroots = numroots(np, sseq, &atmin, &atmax); X X if (nroots == 0) { X printf("solve: no real roots\n"); X exit(0); X } X X /* X * calculate the bracket that the roots live in X */ X min = -1.0; X nchanges = numchanges(np, sseq, min); X for (i = 0; nchanges != atmin && i != MAXPOW; i++) { X min *= 10.0; X nchanges = numchanges(np, sseq, min); X } X X if (nchanges != atmin) { X printf("solve: unable to bracket all negetive roots\n"); X atmin = nchanges; X } X X max = 1.0; X nchanges = numchanges(np, sseq, max); X for (i = 0; nchanges != atmax && i != MAXPOW; i++) { X max *= 10.0; X nchanges = numchanges(np, sseq, max); X } X X if (nchanges != atmax) { X printf("solve: unable to bracket all positive roots\n"); X atmax = nchanges; X } X X nroots = atmin - atmax; X X /* X * perform the bisection. X */ X sbisect(np, sseq, min, max, atmin, atmax, roots); X X X /* X * write out the roots... X */ X if (nroots == 1) { X printf("\n1 distinct real root at x = %f\n", roots[0]); X } else { X printf("\n%d distinct real roots for x: ", nroots); X X for (i = 0; i != nroots; i++) X printf("%f ", roots[i]); X printf("\n"); X } X} END_OF_FILE if test 2173 -ne `wc -c <'Sturm/main.c'`; then echo shar: \"'Sturm/main.c'\" unpacked with wrong size! fi # end of 'Sturm/main.c' fi if test -f 'TriPoints.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'TriPoints.c'\" else echo shar: Extracting \"'TriPoints.c'\" \(2746 characters\) sed "s/^X//" >'TriPoints.c' <<'END_OF_FILE' X/* XGenerating Random Points in Triangles Xby Greg Turk Xfrom "Graphics Gems", Academic Press, 1990 X*/ X X#include X#include "GraphicsGems.h" X X/***************************************************************** XCompute relative areas of sub-triangles that form a convex polygon. XThere are vcount-2 sub-triangles, each defined by the first point Xin the polygon and two other adjacent points. X XThis procedure should be called once before using Xsquare_to_polygon(). X XEntry: X vertices - list of the vertices of a convex polygon X vcount - number of vertices of polygon XExit: X areas - relative areas of sub-triangles of polygon X*******************************************************************/ X Xtriangle_areas (vertices, vcount, areas) X Point3 vertices[]; X int vcount; X float areas[]; X{ X int i; X float area_sum = 0; X Vector3 v1,v2,v3; X X /* compute relative areas of the sub-triangles of polygon */ X X for (i = 0; i < vcount - 2; i++) { X V3Sub(&vertices[i+1], &vertices[0], &v1); X V3Sub(&vertices[i+2], &vertices[0], &v2); X V3Cross(&v1, &v2, &v3); X areas[i] = LengthVector3(&v3); X area_sum += areas[i]; X } X X /* normalize areas so that the sum of all sub-triangles is one */ X X for (i = 0; i < vcount - 2; i++) X areas[i] /= area_sum; X} X X X/****************************************************************** XMap a point from the square [0,1] x [0,1] into a convex polygon. XUniform random points in the square will generate uniform random Xpoints in the polygon. X XThe procedure triangle_areas() must be called once to compute X'areas', and then this procedure can be called repeatedly. X XEntry: X vertices - list of the vertices of a convex polygon X vcount - number of vertices of polygon X areas - relative areas of sub-triangles of polygon X s,t - position in the square [0,1] x [0,1] XExit: X p - position in polygon X*********************************************************************/ X Xsquare_to_polygon (vertices, vcount, areas, s, t, p) X Point3 vertices[]; X int vcount; X float areas[]; X float s,t; X Point3 *p; X{ X int i; X float area_sum = 0; X float a,b,c; X X /* use 's' to pick one sub-triangle, weighted by relative */ X /* area of triangles */ X X for (i = 0; i < vcount - 2; i++) { X area_sum += areas[i]; X if (area_sum >= s) X break; X } X X /* map 's' into the interval [0,1] */ X X s = (s - area_sum + areas[i]) / areas[i]; X X /* map (s,t) to a point in that sub-triangle */ X X t = sqrt(t); X X a = 1 - t; X b = (1 - s) * t; X c = s * t; X X p->x = a * vertices[0].x + b * vertices[i+1].x + c * vertices[i+2].x; X p->y = a * vertices[0].y + b * vertices[i+1].y + c * vertices[i+2].y; X p->z = a * vertices[0].z + b * vertices[i+1].z + c * vertices[i+2].z; X} END_OF_FILE if test 2746 -ne `wc -c <'TriPoints.c'`; then echo shar: \"'TriPoints.c'\" unpacked with wrong size! fi # end of 'TriPoints.c' fi echo shar: End of archive 2 \(of 5\). cp /dev/null ark2isdone 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