#! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # arm.c # chull.c # cube.c # dt.c # extreme.c # graham.c # i.poly.18 # inpoly.c # int.c # macros.h # quadratic.c # sphere.c # stack.c # tri.c # vector.c # This archive created: Tue Oct 4 16:20:41 1994 export PATH; PATH=/bin:$PATH echo shar: extracting "'arm.c'" '(7257 characters)' if test -f 'arm.c' then echo shar: will not over-write existing file "'arm.c'" else cat << \SHAR_EOF > 'arm.c' /* arm.c Prints out one arm configuration to reach given target. Assumes number of links >= 3. */ #include #include #define ERROR 1 #define X 0 #define Y 1 #define NLINKS 20 typedef enum { FALSE, TRUE } bool; #define DIM 2 /* Dimension of points */ typedef int tPointi[DIM]; /* type integer point */ typedef double tPointd[DIM]; /* type double point */ #define FUZZ 1.0e-10 void ReadTarget( tPointi target ); int ReadLinks( void ); void PrintLinks( void ); void PrintPointi( tPointi p ); void PrintPointd( tPointd p ); bool Solve3( int l1, int l2, int l3, tPointi target ); int TwoCircles( tPointi c1, int r1, tPointi c2, int r2, tPointd p ); int TwoCircles0a( int r1, tPointi c2, int r2, tPointd p ); int TwoCircles0b( int r1, tPointi c2, int r2, tPointd p ); void TwoCircles00( int r1, double a2, int r2, tPointd p ); void SubVec( tPointi a, tPointi b, tPointi c ); bool EqPointi( tPointi a, tPointi b ); int Length2( tPointi v ); bool Solve2( int l1, int l2, tPointi target, tPointd J ); bool Nested( int l1, int l2, int l3, tPointi target ); bool Solven( int nlinks ); double RadDeg( double x ); int linklen[NLINKS]; /* link lengths */ int nlinks; /* number of links */ tPointi target; /* target point */ main(argc, argv) int argc; char *argv[]; { ReadTarget( target ); nlinks = ReadLinks( ); if ( !Solven( nlinks ) ) printf("Solven: no solutions!\n"); } bool Solven( int nlinks ) { int i; int m; /* index of median link */ int l1, l2, l3; /* length of links between kinks */ int totlength; /* total length of all links */ int halflength; /* floor of half of total */ /* Compute total & half length. */ totlength = 0; for ( i = 0; i < nlinks; i++ ) totlength += linklen[i]; halflength = totlength / 2; /* Find median link. */ l1 = 0; for ( m = 0; m < nlinks; m++ ) { if ( (l1 + linklen[m]) > halflength) break; l1 += linklen[m]; } l2 = linklen[m]; l3 = totlength - l1 - l2; if ( Solve3( l1, l2, l3, target ) ) return TRUE; else return FALSE; } bool Solve3( int l1, int l2, int l3, tPointi target ) { tPointd Jk; /* coords of kinked joint returned by Solve2 */ tPointi J1; /* Joint1 on x-axis */ tPointi Ttarget;/* translated target */ printf("==>Solve3: links = %d,%d,%d\n", l1, l2, l3); if ( Solve2( l1 + l2, l3, target, Jk ) ) { printf("Solve3: link1=%d, link2=%d, joint=", l1 + l2, l3); PrintPointd( Jk ); return TRUE; } else if ( Solve2( l1, l2 + l3, target, Jk ) ) { printf("Solve3: link1=%d, link2=%d, joint=", l1, l2 + l3); PrintPointd( Jk ); return TRUE; } else { /* pin J0 to 0. */ /* Shift so J1 is origin. */ J1[X] = l1; J1[Y] = 0; SubVec( target, J1, Ttarget ); if ( Solve2( l2, l3, Ttarget, Jk ) ) { /* Shift solution back to origin. */ Jk[X] += l1; printf("Solve3: link1=%d, link2=%d, link3=%d,\ joints=", l1, l2, l3); PrintPointi( J1 ); PrintPointd( Jk ); return TRUE; } else return FALSE; } } bool Solve2( int l1, int l2, tPointi target, tPointd J ) { tPointi c1 = {0,0}; /* center of circle 1 */ int nsoln; /* # of solns: 0,1,2,3(infinite) */ nsoln = TwoCircles( c1, l1, target, l2, J ); printf("<==Solve2: l1=%d, l2=%d; nsoln=%d\n", l1, l2, nsoln); return nsoln != 0; } /* TwoCircles finds the intersection points between two circles. This is the general routine, no assumptions. One intersection point is placed in p. */ int TwoCircles( tPointi c1, int r1, tPointi c2, int r2, tPointd p) { tPointi c; SubVec( c2, c1, c ); return TwoCircles0a( r1, c, r2, p ); } /* TwoCircles0a assumes that the first circle is centered on the origin. It handles all the special cases, returning the number of intersection points: 0, 1, 2, 3 (inf). */ int TwoCircles0a( int r1, tPointi c2, int r2, tPointd p ) { int dc2; /* dist to center 2 squared */ int rplus2, rminus2; /* (r1 +/- r2)^2 */ double f; /* fraction along c2 for nsoln=1 */ /* Handle special cases. */ dc2 = Length2( c2 ); rplus2 = (r1 + r2) * (r1 + r2); rminus2 = (r1 - r2) * (r1 - r2); /* No solution if c2 out of reach + or -. */ if ( ( dc2 > rplus2 ) || ( dc2 < rminus2 ) ) return 0; /* One solution if c2 just reached. */ /* Then solution is r1-of-the-way (f) to c2. */ if ( dc2 == rplus2 ) { f = r1 / (double)(r1 + r2); p[X] = f * c2[X]; p[Y] = f * c2[Y]; return 1; } if ( dc2 == rminus2 ) { if ( rminus2 == 0 ) { /* Circles coincide. */ p[X] = r1; p[Y] = 0; return 3; } f = r1 / (double)(r1 - r2); p[X] = f * c2[X]; p[Y] = f * c2[Y]; return 1; } /* Two intersections. */ return TwoCircles0b( r1, c2, r2, p ); } /* TwoCircles0b also assumes that the first circle is centered on the origin. It rotates c2 to lie on the x-axis. */ int TwoCircles0b( int r1, tPointi c2, int r2, tPointd p ) { double a2; /* center of 2nd circle when rotated to x-axis */ tPointd q; /* one solution when c2 on x-axis */ double cost, sint; /* sine and cosine of angle of c2 */ /* Rotate c2 to a2 on x-axis. */ a2 = sqrt( (double) Length2( c2 ) ); cost = c2[X] / a2; sint = c2[Y] / a2; TwoCircles00( r1, a2, r2, q ); /* Rotate back */ p[X] = cost * q[X] + -sint * q[Y]; p[Y] = sint * q[X] + cost * q[Y]; return 2; } /* TwoCircles00 assume circle one is origin-centered, and that the center of the second circle is at (a2,0). */ void TwoCircles00( int r1, double a2, int r2, tPointd p ) { /* Two intersections (only one returned in p). */ p[X] = ( a2 + ( r1*r1 - r2*r2 ) / a2 ) / 2; p[Y] = sqrt( r1*r1 - p[X]*p[X] ); } void ReadTarget( tPointi target ) { int i; for( i = 0; i < DIM; i++) scanf("%d", &target[i] ); printf("Target (target) point = "); PrintPointi( target ); putchar('\n'); } void PrintPointi( tPointi p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%d", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } void PrintPointd( tPointd p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%g", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); putchar('\n'); } int ReadLinks( void ) { int length; nlinks = 0; while( scanf("%d", &length) != EOF ) { linklen[nlinks] = length; nlinks++; } PrintLinks( ); return nlinks; } void PrintLinks( ) { int i; for ( i=0; i < nlinks; i++) printf("%d:%d\t", i, linklen[i]); putchar('\n'); } bool EqPointi( tPointi a, tPointi b ) { int i; for ( i=0; i < DIM; i++ ) if ( a[i] != b[i] ) return FALSE; return TRUE; } int Length2( tPointi v ) { int i; int ss; ss = 0; for ( i=0; i < DIM; i++ ) ss += v[i] * v[i]; return ss; } double RadDeg( double x ) { return 180 * x / M_PI; } void SubVec( tPointi a, tPointi b, tPointi c ) { int i; for ( i=0; i < DIM; i++ ) { c[i] = a[i] - b[i]; } } SHAR_EOF fi # end of overwriting check echo shar: extracting "'chull.c'" '(27374 characters)' if test -f 'chull.c' then echo shar: will not over-write existing file "'chull.c'" else cat << \SHAR_EOF > 'chull.c' /* chull.c Written by Joseph O'Rourke and John Kutcher, Catherine Schevon, Susan Weller. Last modified: 18 Februrary 1993. */ #include #include /* Define Boolean type */ typedef enum { FALSE, TRUE } bool; /* Define vertex indices. */ #define X 0 #define Y 1 #define Z 2 /* Define structures for vertices, edges and faces */ typedef struct tVertexStructure tsVertex; typedef tsVertex *tVertex; typedef struct tEdgeStructure tsEdge; typedef tsEdge *tEdge; typedef struct tFaceStructure tsFace; typedef tsFace *tFace; struct tVertexStructure { int v[3]; int vnum; tEdge duplicate; /* pointer to incident cone edge (or NULL) */ bool onhull; /* T iff point on hull. */ bool mark; /* T iff point already processed. */ tVertex next, prev; }; struct tEdgeStructure { tFace adjface[2]; tVertex endpts[2]; tFace newface; /* pointer to incident cone face. */ bool delete; /* T iff edge should be delete. */ tEdge next, prev; }; struct tFaceStructure { tEdge edge[3]; tVertex vertex[3]; bool visible; /* T iff face visible from new point. */ tFace next, prev; }; /* Define flags */ #define ONHULL TRUE #define REMOVED TRUE #define VISIBLE TRUE #define PROCESSED TRUE /* Global variable definitions */ tVertex vertices = NULL; tEdge edges = NULL; tFace faces = NULL; bool debug = FALSE; tVertex vertices; tEdge edges; tFace faces; bool debug; /* Function declarations */ tVertex MakeVertex( void ); void ReadVertices( void ); void Print( void ); void Tetrahedron( void ); void ConstructHull( void ); bool AddOne( tVertex p ); int Volume6(tFace f, tVertex p); tFace MakeStructs( tEdge e, tVertex p ); void MakeCcw( tFace f, tEdge e, tVertex p ); tEdge MakeEdge( void ); tFace MakeFace( void ); void CleanUp( void ); void CleanEdges( void ); void CleanFaces( void ); void CleanVertices( void ); bool Collinear( tVertex a, tVertex b, tVertex c ); void CheckEuler(int V, int E, int F ); double Volumed( tFace f, tVertex p ); void PrintPoint( tVertex p ); void Checks( void ); void Consistency( void ); void Convexity( void ); void PrintOut( tVertex v ); void PrintVertices( void ); void PrintEdges( void ); void PrintFaces( void ); #include "macros.h" /*--------------------------------------------------------------------------*/ main( int argc, char *argv[] ) { char *s; if ( argc > 1 && argv[1][0] == '-' ) for ( s = argv[1] + 1; *s; ++s ) switch ( *s ) { case 'd': debug = TRUE; fprintf( stderr, "Debug mode\n"); break; default : debug = FALSE; } else if ( argc > 1 && argv[1][0] != '-' ) { printf ("Usage: %s -d[ebug] \n", *argv ); printf ("x y z coords of vertices from stdin\n"); exit(1); } ReadVertices(); Tetrahedron(); ConstructHull(); Print(); } /*------------------------------------------------------------------ MakeVertex: Makes a vertex, nulls out fields. --------------------------------------------------------------------*/ tVertex MakeVertex( void ) { tVertex v; NEW( v, tsVertex ); v->duplicate = NULL; v->onhull = !ONHULL; v->mark = !PROCESSED; ADD( vertices, v ); return v; } /*------------------------------------------------------------------ ReadVertices: Reads in the vertices, and links them into a circular list with MakeVertex. There is no need for the # of vertices to be the first line: the function looks for EOF instead. --------------------------------------------------------------------*/ void ReadVertices( void ) { tVertex v; int x, y, z; int vnum = 0; while ( scanf ("%d %d %d", &x, &y, &z ) != EOF ) { v = MakeVertex(); v->v[X] = x; v->v[Y] = y; v->v[Z] = z; v->vnum = vnum++; if ( abs(x) > 512 ) printf("coord %d might be too large\n", x); if ( abs(y) > 512 ) printf("coord %d might be too large\n", y); if ( abs(z) > 512 ) printf("coord %d might be too large\n", z); } } /*------------------------------------------------------------------ Print: Prints out the vertices and the faces. Uses the vnum indices corresponding to the order in which the vertices were input. --------------------------------------------------------------------*/ void Print( void ) { /* Pointers to vertices, edges, faces. */ tVertex v; tEdge e; tFace f; /* Counters for Euler's formula. */ int V = 0, E = 0 , F = 0; /* Note: lowercase==pointer, uppercase==counter. */ /* Vertices. */ printf("\n"); v = vertices; do { V++; v = v->next; } while ( v != vertices ); printf("Vertices:\tV = %d\n", V ); printf("index:\tx\ty\tz\n"); do { printf("%5d:\t%d\t%d\t%d\n", v->vnum, v->v[X], v->v[Y], v->v[Z] ); v = v->next; } while ( v != vertices ); /* Faces. */ f = faces; do { ++F; f = f ->next; } while ( f != faces ); printf("Faces:\t\tF = %d\n", F ); printf("\tv0\tv1\tv2\t(vertex indices)\n"); do { printf("\t%d\t%d\t%d\n", f->vertex[0]->vnum, f->vertex[1]->vnum, f->vertex[2]->vnum ); f = f->next; } while ( f != faces ); /* Edges. */ e = edges; do { E++; e = e->next; } while ( e != edges ); printf("Edges:\tE = %d\n", E ); /* Edges not printed out (but easily added). */ debug = TRUE; CheckEuler( V, E, F ); } /*----------------------------------------------------------------------- Tetrahedron builds the initial tetrahedron. It first finds 3 noncollinear points and makes a face out of them, and then finds a fourth point that is not coplanar with that face. The vertices are stored in the face structure in counterclockwise order so that the volume between the face and the point is negative. Lastly, the 3 newfaces to the fourth point are constructed and the data structures are cleaned up. -----------------------------------------------------------------------*/ void Tetrahedron( void ) { tVertex v1, v4, t; tFace f; tEdge e1, e2, e3, s; int vol; /* Find 3 non-Collinear points. */ v1 = vertices; while ( Collinear( v1, v1->next, v1->next->next ) ) if ( ( v1 = v1->next ) == vertices ) { printf("All points are Collinear!\n"); exit(0); } /* Mark the vertices as processed. */ v1->mark = PROCESSED; v1->next->mark = PROCESSED; v1->next->next->mark = PROCESSED; /* Create edges of the initial triangle. */ e1 = MakeEdge(); e2 = MakeEdge(); e3 = MakeEdge(); e1->endpts[0] = v1; e1->endpts[1] = v1->next; e2->endpts[0] = v1->next; e2->endpts[1] = v1->next->next; e3->endpts[0] = v1->next->next; e3->endpts[1] = v1; /* Create face for triangle. */ f = MakeFace(); f->edge[0] = e1; f->edge[1] = e2; f->edge[2] = e3; f->vertex[0] = v1; f->vertex[1] = v1->next; f->vertex[2] = v1->next->next; /* Link edges to face. */ e1->adjface[0] = e2->adjface[0] = e3->adjface[0] = f; /* Find a fourth, non-coplanar point to form tetrahedron. */ v4 = v1->next->next->next; vol = Volume6( f, v4 ); while ( !vol ) { if ( ( v4 = v4->next ) == v1 ) { printf("All points are coplanar!\n"); exit(0); } vol = Volume6( f, v4 ); } v4->mark = PROCESSED; /* Store vertices in ccw order. */ if( vol < 0 ) { SWAP( t, f->vertex[1], f->vertex[2] ); SWAP( s, f->edge[1], f->edge[2] ); } /* Construct the faces and edges between the original triangle and the fourth point. */ e1->adjface[1] = MakeStructs( e1, v4 ); e2->adjface[1] = MakeStructs( e2, v4 ); e3->adjface[1] = MakeStructs( e3, v4 ); CleanUp(); } /*------------------------------------------------------------------------- ConstructHull adds the vertices to the hull one at a time. The hull vertices are those in the list marked as onhull. -------------------------------------------------------------------------*/ void ConstructHull( void ) { tVertex v; tFace f; int vol; bool changed; /* T if addition changes hull; not used. */ v = vertices; f = faces; do { if ( !v->mark ) { v->mark = PROCESSED; changed = AddOne( v ); CleanUp(); if ( debug ) { fprintf(stderr,"after cleanup:\n"); PrintOut( v ); Checks(); } } v = v->next; } while ( v != vertices ); } /*------------------------------------------------------------------------- AddOne is passed a vertex. It first determines all faces visible from that point. If none are visible then the point is marked as not onhull. Next is a loop over edges. If both faces adjacent to an edge are visible, then the edge is marked for deletion. If just one of the adjacent faces is visible then a new face is constructed. --------------------------------------------------------------------------*/ bool AddOne( tVertex p ) { tFace f; tEdge e; int vol; bool vis; /* Mark faces visible from p. */ f = faces; do { vol = Volume6( f, p ); /* if (debug) fprintf(stderr,"faddr: %6x paddr: %6x Vol = %d\n",f,p,vol); */ if ( vol < 0 ) { f->visible = VISIBLE; vis = TRUE; } f = f->next; } while ( f != faces ); /* If no faces are visible from p, then p is inside the hull. */ if ( !vis ) { p->onhull = !ONHULL; return FALSE; } /* Mark edges in interior of visible region for deletion. Erect a newface based on each border edge. */ e = edges; do { tEdge temp; temp = e->next; if ( e->adjface[0]->visible && e->adjface[1]->visible ) /* e interior: mark for deletion. */ e->delete = REMOVED; else if ( e->adjface[0]->visible || e->adjface[1]->visible ) /* e border: make a new face. */ e->newface = MakeStructs( e, p ); e = temp; } while ( e != edges ); return TRUE; } /*------------------------------------------------------------------------- Volume6 returns six times the volume of the tetrahedron determined by f and p. Volume6 is positive iff p is on the negative side of f, where the positive side is determined by the rh-rule. So the volume is positive if the ccw normal to f points outside the tetrahedron. --------------------------------------------------------------------------*/ int Volume6( tFace f, tVertex p ) { int vol; int ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz; double vold; int i; ax = f->vertex[0]->v[X]; ay = f->vertex[0]->v[Y]; az = f->vertex[0]->v[Z]; bx = f->vertex[1]->v[X]; by = f->vertex[1]->v[Y]; bz = f->vertex[1]->v[Z]; cx = f->vertex[2]->v[X]; cy = f->vertex[2]->v[Y]; cz = f->vertex[2]->v[Z]; dx = p->v[X]; dy = p->v[Y]; dz = p->v[Z]; vol = -az * by * cx + ay * bz * cx + az * bx * cy - ax * bz * cy - ay * bx * cz + ax * by * cz + az * by * dx - ay * bz * dx - az * cy * dx + bz * cy * dx + ay * cz * dx - by * cz * dx - az * bx * dy + ax * bz * dy + az * cx * dy - bz * cx * dy - ax * cz * dy + bx * cz * dy + ay * bx * dz - ax * by * dz - ay * cx * dz + by * cx * dz + ax * cy * dz - bx * cy * dz; /* Compare integer volume with double volume for saftey. */ vold = Volumed( f, p ); if (debug) { fprintf(stderr,"Face = %6x\tVertex = %d, vol(int) = %d\tvol(double) = %lf\n", f,p->vnum,vol,vold); } if ( fabs( vol - vold ) >= 1.0 ) { printf("Likely integer overflow in volume calculation\n"); printf("because coordinates are too large:\n"); printf("\tvol(int) = %d;\tvol(double) = %lf\n"); printf("\tfour points:\n"); for ( i=0; i < 3; i++ ) PrintPoint( f->vertex[i] ); PrintPoint( p ); printf("Basing decision based on double volume.\n"); /* Return based on double volume. */ if ( vold >= 1.0 ) return 1; else if ( vold <= -1.0 ) return -1; else return 0; } return vol; } /*------------------------------------------------------------------------- Volumed is the same as Volume6 but computed with doubles. For protection against overflow. --------------------------------------------------------------------------*/ double Volumed( tFace f, tVertex p ) { double vol; double ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz; ax = f->vertex[0]->v[X]; ay = f->vertex[0]->v[Y]; az = f->vertex[0]->v[Z]; bx = f->vertex[1]->v[X]; by = f->vertex[1]->v[Y]; bz = f->vertex[1]->v[Z]; cx = f->vertex[2]->v[X]; cy = f->vertex[2]->v[Y]; cz = f->vertex[2]->v[Z]; dx = p->v[X]; dy = p->v[Y]; dz = p->v[Z]; vol = -(az * by * cx) + ay * bz * cx + az * bx * cy - ax * bz * cy - ay * bx * cz + ax * by * cz + az * by * dx - ay * bz * dx - az * cy * dx + bz * cy * dx + ay * cz * dx - by * cz * dx - az * bx * dy + ax * bz * dy + az * cx * dy - bz * cx * dy - ax * cz * dy + bx * cz * dy + ay * bx * dz - ax * by * dz - ay * cx * dz + by * cx * dz + ax * cy * dz - bx * cy * dz; return vol; } /*-------------------------------------------------------------------------*/ void PrintPoint( tVertex p ) { int i; for ( i = 0; i < 3; i++ ) printf("\t%d", p->v[i]); putchar('\n'); } /*---------------------------------------------------------------------- MakeStructs makes a new face and two new edges between the edge and the point that are passed to it. It returns a pointer to the new face. ----------------------------------------------------------------------*/ tFace MakeStructs( tEdge e, tVertex p ) { tEdge new_edge[2]; tFace new_face; int i, j; /* Make two new edges (if don't already exist). */ for ( i=0; i < 2; ++i ) /* If the edge exists, copy it into new_edge. */ if ( !( new_edge[i] = e->endpts[i]->duplicate) ) { /* Otherwise (duplicate is NULL), MakeEdge. */ new_edge[i] = MakeEdge(); new_edge[i]->endpts[0] = e->endpts[i]; new_edge[i]->endpts[1] = p; e->endpts[i]->duplicate = new_edge[i]; } /* Make the new face. */ new_face = MakeFace(); new_face->edge[0] = e; new_face->edge[1] = new_edge[0]; new_face->edge[2] = new_edge[1]; MakeCcw( new_face, e, p ); /* Set the adjacent face pointers. */ for ( i=0; i < 2; ++i ) for ( j=0; j < 2; ++j ) /* Only one NULL link should be set to new_face. */ if ( !new_edge[i]->adjface[j] ) { new_edge[i]->adjface[j] = new_face; break; } return new_face; } /*------------------------------------------------------------------------ MakeCcw puts the vertices in the face structure in counterclockwise order. If there is no adjacent face[1] then we know that we are working with the first face of the initial tetrahedron. In this case we want to store the vertices in the opposite order from the initial face. Otherwise, we want to store the vertices in the same order as in the visible face. The third vertex is always p. ------------------------------------------------------------------------*/ void MakeCcw( tFace f, tEdge e, tVertex p ) { int i; /* Index */ tFace fi; /* The invisible face adjacent to e */ tEdge s; /* Temporary, for swapping */ /* If this is the initial tetrahedron, then e has only one adjacent face, and use that for fi. Otherwise, use the invisible face. */ if ( !e->adjface[1] ) fi = e->adjface[0]; else { if ( !e->adjface[0]->visible ) fi = e->adjface[0]; else fi = e->adjface[1]; } /* Set vertex[0] & [1] of f to have the opposite orientation as do the corresponding vertices of fi. */ /* Find the index i of e->endpoint[1] in fi. */ for ( i=0; fi->vertex[i] != e->endpts[1]; ++i ) ; /* Orient f opposite that of fi. */ if ( fi->vertex[ (i+1) % 3 ] != e->endpts[0] ) { f->vertex[0] = e->endpts[1]; f->vertex[1] = e->endpts[0]; } else { f->vertex[0] = e->endpts[0]; f->vertex[1] = e->endpts[1]; SWAP( s, f->edge[1], f->edge[2] ); } f->vertex[2] = p; } /*--------------------------------------------------------------------- MakeEdge creates a new cell and initializes all pointers to NULL and sets all flags to off. It returns a pointer to the empty cell. ---------------------------------------------------------------------*/ tEdge MakeEdge( void ) { tEdge e; NEW( e, tsEdge ); e->adjface[0] = e->adjface[1] = e->newface = NULL; e->endpts[0] = e->endpts[1] = NULL; e->delete = !REMOVED; ADD( edges, e ); return e; } /*--------------------------------------------------------------------- MakeFace creates a new face structure and initializes all of its flags to NULL and sets all the flags to off. It returns a pointer to the empty cell. ----------------------------------------------------------------------*/ tFace MakeFace( void ) { tFace f; int i; NEW( f, tsFace); for ( i=0; i < 3; ++i ) { f->edge[i] = NULL; f->vertex[i] = NULL; } f->visible = !VISIBLE; ADD( faces, f ); return f; } /*------------------------------------------------------------------------- CleanUp goes through each data structure list and clears all flags and NULLs out some pointers. The order of processing (edges, faces, vertices) is important. ------------------------------------------------------------------------*/ void CleanUp( void ) { CleanEdges(); CleanFaces(); CleanVertices(); } /*------------------------------------------------------------------------ CleanEdges runs through the edge list and cleans up the structure. If there is a newface then it will put that face in place of the visible face and NULL out newface. It also deletes so marked edges. -----------------------------------------------------------------------*/ void CleanEdges( void ) { tEdge e; /* Primary index into edge list. */ tEdge t; /* Temporary edge pointer. */ /* Integrate the newface's into the data structure. */ /* Check every edge. */ e = edges; do { if ( e->newface ) { if ( e->adjface[0]->visible ) e->adjface[0] = e->newface; else e->adjface[1] = e->newface; e->newface = NULL; } e = e->next; } while ( e != edges ); /* Delete any edges marked for deletion. */ while ( edges && edges->delete ) { e = edges; DELETE( edges, e ); } e = edges->next; do { if ( e->delete ) { t = e; e = e->next; DELETE( edges, t ); } else e = e->next; } while ( e != edges ); } /*------------------------------------------------------------------------ CleanFaces runs through the face list and deletes any face marked visible. -----------------------------------------------------------------------*/ void CleanFaces( void ) { tFace f; /* Primary pointer into face list. */ tFace t; /* Temporary pointer, for deleting. */ while ( faces && faces->visible ) { f = faces; DELETE( faces, f ); } f = faces->next; do { if ( f->visible ) { t = f; f = f->next; DELETE( faces, t ); } else f = f->next; } while ( f != faces ); } /*------------------------------------------------------------------------- CleanVertices runs through the vertex list and deletes the vertices that are marked as processed but are not incident to any undeleted edges. -------------------------------------------------------------------------*/ void CleanVertices( void ) { tEdge e; tVertex v, t; /* Mark all vertices incident to some undeleted edge as on the hull. */ e = edges; do { e->endpts[0]->onhull = e->endpts[1]->onhull = ONHULL; e = e->next; } while (e != edges); /* Delete all vertices that have been processed but are not on the hull. */ while ( vertices && vertices->mark && !vertices->onhull ) { v = vertices; DELETE( vertices, v ); } v = vertices->next; do { if ( v->mark && !v->onhull ) { t = v; v = v->next; DELETE( vertices, t ) } else v = v->next; } while ( v != vertices ); /* Reset flags. */ v = vertices; do { v->duplicate = NULL; v->onhull = !ONHULL; v = v->next; } while ( v != vertices ); } /*------------------------------------------------------------------ Collinear checks to see if the three points given are collinear, by checking to see if each element of the cross product is zero. ---------------------------------------------------------------------*/ bool Collinear( tVertex a, tVertex b, tVertex c ) { return ( c->v[Z] - a->v[Z] ) * ( b->v[Y] - a->v[Y] ) - ( b->v[Z] - a->v[Z] ) * ( c->v[Y] - a->v[Y] ) == 0 && ( b->v[Z] - a->v[Z] ) * ( c->v[X] - a->v[X] ) - ( b->v[X] - a->v[X] ) * ( c->v[Z] - a->v[Z] ) == 0 && ( b->v[X] - a->v[X] ) * ( c->v[Y] - a->v[Y] ) - ( b->v[Y] - a->v[Y] ) * ( c->v[X] - a->v[X] ) == 0 ; } /*------------------------------------------------------------------------ Consistency runs through the edge list and checks that all adjacent faces have their endpoints in opposite order. This verifies that the vertices are in counterclockwise order. -----------------------------------------------------------------------*/ void Consistency( void ) { register tEdge e; register int i, j; e = edges; do { /* find index of endpoint[0] in adjacent face[0] */ for ( i = 0; e->adjface[0]->vertex[i] != e->endpts[0]; ++i ) ; /* find index of endpoint[0] in adjacent face[1] */ for ( j = 0; e->adjface[1]->vertex[j] != e->endpts[0]; ++j ) ; /* check if the endpoints occur in opposite order */ if ( !( e->adjface[0]->vertex[ (i+1) % 3 ] == e->adjface[1]->vertex[ (j+2) % 3 ] || e->adjface[0]->vertex[ (i+2) % 3 ] == e->adjface[1]->vertex[ (j+1) % 3 ] ) ) break; e = e->next; } while ( e != edges ); if ( e != edges ) fprintf( stderr, "Checks: edges are NOT consistent.\n"); else fprintf( stderr, "Checks: edges consistent.\n"); } /*---------------------------------------------------------------------- Convexity checks that the volume between every face and every point is negative. This shows that each point is inside every face and therefore the hull is convex. ---------------------------------------------------------------------*/ void Convexity( void ) { register tFace f; register tVertex v; int vol; f = faces; do { v = vertices; do { if ( v->mark ) { vol = Volume6( f, v ); if ( vol < 0 ) break; } v = v->next; } while ( v != vertices ); f = f->next; } while ( f != faces ); if ( f != faces ) fprintf( stderr, "Checks: NOT convex.\n"); else if ( debug ) fprintf( stderr, "Checks: convex.\n"); } /*---------------------------------------------------------------------- CheckEuler checks Euler's relation, as well as its implications when all faces are known to be triangles. Only prints positive information when debug is true, but always prints negative information. ---------------------------------------------------------------------*/ void CheckEuler( int V, int E, int F ) { if ( debug ) fprintf( stderr, "Checks: V, E, F = %d %d %d:\t", V, E, F); if ( (V - E + F) != 2 ) fprintf( stderr, "Checks: V-E+F != 2\n"); else if ( debug ) fprintf( stderr, "V-E+F = 2\t"); if ( F != (2 * V - 4) ) fprintf( stderr, "Checks: F=%d != 2V-4=%d; V=%d\n", F, 2*V-4, V); else if ( debug ) fprintf( stderr, "F = 2V-4\t"); if ( (2 * E) != (3 * F) ) fprintf( stderr, "Checks: 2E=%d != 3F=%d; E=%d, F=%d\n", 2*E, 3*F, E, F ); else if ( debug ) fprintf( stderr, "2E = 3F\n"); } /*-----------------------------------------------------------------------*/ void Checks( void ) { tVertex v; tEdge e; tFace f; int V = 0, E = 0 , F = 0; Consistency(); Convexity(); if ( v = vertices ) do { if (v->mark) V++; v = v->next; } while ( v != vertices ); if ( e = edges ) do { E++; e = e->next; } while ( e != edges ); if ( f = faces ) do { F++; f = f ->next; } while ( f != faces ); CheckEuler( V, E, F ); } /*================================================================ These functions are used whenever the debug flag is set. They print out the entire contents of each data structure. Printing is to standard error. To grab the output in a file in the csh, use this: chull < i.file >&! o.file ==================================================================*/ /*-------------------------------------------------------------------*/ void PrintOut( tVertex v ) { fprintf( stderr, "\nAdding vertex %6x :\n", v ); PrintVertices(); PrintEdges(); PrintFaces(); } /*-------------------------------------------------------------------------*/ void PrintVertices( void ) { tVertex temp; temp = vertices; fprintf (stderr, "Vertex List\n"); if (vertices) do { fprintf(stderr," addr %6x\t", vertices ); fprintf(stderr," vnum %4d", vertices->vnum ); fprintf(stderr," (%6d,%6d,%6d)",vertices->v[X], vertices->v[Y], vertices->v[Z] ); fprintf(stderr," active:%3d", vertices->onhull ); fprintf(stderr," dup:%5x", vertices->duplicate ); fprintf(stderr," mark:%2d\n", vertices->mark ); vertices = vertices->next; } while ( vertices != temp ); } /*-------------------------------------------------------------------------*/ void PrintEdges( void ) { tEdge temp; int i; temp = edges; fprintf (stderr, "Edge List\n"); if (edges) do { fprintf( stderr, " addr: %6x\t", edges ); fprintf( stderr, "adj: "); for (i=0; i<2; ++i) fprintf( stderr, "%6x", edges->adjface[i] ); fprintf( stderr, " endpts:"); for (i=0; i<2; ++i) fprintf( stderr, "%4d", edges->endpts[i]->vnum); fprintf( stderr, " del:%3d\n", edges->delete ); edges = edges->next; } while (edges != temp ); } /*-------------------------------------------------------------------------*/ void PrintFaces( void ) { int i; tFace temp; temp = faces; fprintf (stderr, "Face List\n"); if (faces) do { fprintf(stderr, " addr: %6x\t", faces ); fprintf(stderr, " edges:"); for( i=0; i<3; ++i ) fprintf(stderr, "%6x", faces->edge[i] ); fprintf(stderr, " vert:"); for ( i=0; i<3; ++i) fprintf(stderr, "%4d", faces->vertex[i]->vnum ); fprintf(stderr, " vis: %d\n", faces->visible ); faces= faces->next; } while ( faces != temp ); } SHAR_EOF fi # end of overwriting check echo shar: extracting "'cube.c'" '(1065 characters)' if test -f 'cube.c' then echo shar: will not over-write existing file "'cube.c'" else cat << \SHAR_EOF > 'cube.c' /*----------------------------------------------------------------------- Cube Susan L. Weller ----------------------------------------------------------------------- This program will generate a given number of points that lie inside a cube. The number of points is given on command line. -----------------------------------------------------------------------*/ #include #include #define MAX_INT 2147483647 main( argc, argv ) int argc; char *argv[]; { int n; double x, y, z; srandom( (int) time( (long *) 0 ) ); if ( argc != 2 ) { printf("Usage: cube \n"); exit(1); } n = atoi( argv[1] ); while ( n-- ) { x = 2.0 * (double) random() / MAX_INT - 1.0; y = 2.0 * (double) random() / MAX_INT - 1.0; z = 2.0 * (double) random() / MAX_INT - 1.0; printf ("%6d %6d %6d\n", (int) (100*x ), (int) ( 100*y ), (int) (100 * z) ); } } SHAR_EOF fi # end of overwriting check echo shar: extracting "'dt.c'" '(1147 characters)' if test -f 'dt.c' then echo shar: will not over-write existing file "'dt.c'" else cat << \SHAR_EOF > 'dt.c' #include #define NMAX 100 main() { int x[NMAX],y[NMAX],z[NMAX];/* input points xy,z=x^2+y^2 */ int n; /* number of input points */ int i, j, k, m; /* indices of four points */ int xn, yn, zn; /* outward vector normal to (i,j,k) */ int flag; /* t if m above of (i,j,k) */ /* Input points and compute z = x^2 + y^2. */ scanf("%d", &n); for ( i = 0; i < n; i++ ) { scanf("%d %d", &x[i], &y[i]); z[i] = x[i] * x[i] + y[i] * y[i]; } /* For each triple (i,j,k) */ for ( i = 0; i < n - 2; i++ ) for ( j = i + 1; j < n; j++ ) for ( k = i + 1; k < n; k++ ) if ( j != k ) { /* Compute normal to triangle (i,j,k). */ xn = (y[j]-y[i])*(z[k]-z[i]) - (y[k]-y[i])*(z[j]-z[i]); yn = (x[k]-x[i])*(z[j]-z[i]) - (x[j]-x[i])*(z[k]-z[i]); zn = (x[j]-x[i])*(y[k]-y[i]) - (x[k]-x[i])*(y[j]-y[i]); /* Only examine faces on bottom of paraboloid: zn < 0. */ if ( flag = (zn < 0) ) /* For each other point m */ for (m = 0; m < n; m++) /* Check if m above (i,j,k). */ flag = flag && ((x[m]-x[i])*xn + (y[m]-y[i])*yn + (z[m]-z[i])*zn <= 0); if (flag) printf("%d\t%d\t%d\n", i, j, k); } } SHAR_EOF fi # end of overwriting check echo shar: extracting "'extreme.c'" '(5406 characters)' if test -f 'extreme.c' then echo shar: will not over-write existing file "'extreme.c'" else cat << \SHAR_EOF > 'extreme.c' /* extreme.c Written by Joseph O'Rourke orourke@cs.smith.edu */ #include #include #define ERROR 1 #define X 0 #define Y 1 typedef enum { FALSE, TRUE } bool; #define DIM 2 /* Dimension of points */ typedef int tPointi[DIM]; /* type integer point */ typedef double tPointd[DIM]; /* type double point */ #define PMAX 1000 /* Max # of pts in polygon */ typedef tPointi tPolygoni[PMAX];/* type integer polygon */ void SubVec( tPointi a, tPointi b, tPointi c ); int Dot( tPointi a, tPointi b ); void PrintPoly( int n, tPolygoni P, int labels[] ); void PrintPoint( tPointi p ); int Midway( int a, int b, int n ); int Extreme( tPointi u, tPolygoni P, int n ); main() { int n, c; tPolygoni P; tPointi u; n = ReadPoly( P ); u[0] = 0; u[1] = 1; c = Extreme( u, P, n ); printf("Extreme returns %d\n", c); u[0] = 0; u[1] = -1; c = Extreme( u, P, n ); printf("Extreme returns %d\n", c); u[0] = 1; u[1] = 0; c = Extreme( u, P, n ); printf("Extreme returns %d\n", c); u[0] = -1; u[1] = 0; c = Extreme( u, P, n ); printf("Extreme returns %d\n", c); } /* Returns index midway ccw from a to b, mod n. */ int Midway( int a, int b, int n ) { if (a < b) return ( a + b ) / 2; else return ( ( a + b + n ) / 2 ) % n; } /* Returns index of a point extreme in direction u. */ int Extreme( tPointi u, tPolygoni P, int n ) { int a,a1, b;/* [a,b] includes extreme; a1=a+1. */ int c, c1; /* index midway; c1 is c +- 1. */ tPointi A, C, B;/* edge vectors after a, after c, before c. */ int Adot, Cdot, Bdot; /* dots with u. */ int y; /* height difference: P[a] - P[c] in dir. u. */ tPointi ur; /* u rotated cw by pi/2 */ printf("\n==>Extreme: u = "); PrintPoint(u); putchar('\n'); ur[0] = u[1]; ur[1] = -u[0]; a = 0; b = 0; do { c = Midway( a, b, n ); printf("Extreme: a,b,c=%d\t%d\t%d\n", a, b, c); /* Compute basic vectors and dots. */ a1 = ( a + 1 ) % n; SubVec( P[a1], P[a], A ); c1 = ( c + 1 ) % n; SubVec( P[c1], P[c], C ); c1 = ( c + ( n-1 ) ) % n; SubVec( P[c], P[c1], B ); printf("Extreme: A,C,B: "); PrintPoint( A ); putchar('\t'); PrintPoint( C ); putchar('\t'); PrintPoint( B ); putchar('\n'); Adot = Dot(u,A); Cdot = Dot(u,C); Bdot = Dot(u,B); y = Dot(u,P[a]) - Dot(u,P[c]); printf("Extreme:\tAdot=%d\t", Adot); printf("Cdot=%d\t", Cdot); printf("Bdot=%d\t", Bdot); printf("y = %d\n", y); /* Termination conditions */ /* If either A or C points left of u, then at extreme. */ if ( (Adot == 0) && (Dot(ur,A) < 0) ) { printf("Extreme: A points left, so return a=%d\n",a); return a; } if ( (Cdot == 0) && (Dot(ur,C) < 0) ) { printf("Extreme: C points left, so return c=%d\n",c); return c; } /* From here on, can assume that zero dots put point on bot */ if ( (Cdot < 0) && (Bdot > 0) ) { printf("Extreme: unique extreme c=%d\n", c); return c; } if (a == c) { printf("Extreme: a=c=%d: choosing a or b=a+1\n", a); if (Adot > 0) return b; else return a; } /* Halving interval */ if ( (Adot >= 0) && (Cdot <= 0) ) b = c; /* new: (a,c) */ else if ( (Adot <= 0) && (Cdot >= 0) ) a = c; /* new: (c,b) */ else if ( (Adot > 0) && (Cdot > 0) ) { if ( y > 0 ) b = c; /* new: (a,c) */ else a = c; /* new: (c,b) */ } else if ( (Adot < 0) && (Cdot < 0) ) { if ( y < 0 ) b = c; /* new: (a,c) */ else a = c; /* new: (c,b) */ } else { printf("Error in Extreme: shouldn't reach here\n"); exit(1); } } while ( TRUE ); } /* Puts a - b into c. */ void SubVec( tPointi a, tPointi b, tPointi c ) { int i; for( i = 0; i < DIM; i++ ) c[i] = a[i] - b[i]; } /* Returns the dot product of the two input vectors. */ int Dot( tPointi a, tPointi b ) { int i; int sum; sum = 0; /* printf("Dot: a, b =\n"); PrintPoint(a); putchar('\n'); PrintPoint(b); putchar('\n'); */ for( i = 0; i < DIM; i++ ) { /*printf("before: i=%d, sum=%d, a[i]=%d, b[i]=%d\n", i, sum, a[i], b[i]);*/ sum += a[i] * b[i]; /*printf("after: i=%d, sum=%d, a[i]=%d, b[i]=%d\n", i, sum, a[i], b[i]);*/ } /*printf("Dot: returning %d\n", sum);*/ return sum; } /* Implements a = b, assignment of points/vectors. Assignment between arrays is not possible in C. */ void PointAssign( tPointi a, tPointi b ) { int i; for ( i = 0; i < DIM; i++ ) a[i] = b[i]; } void PrintPoint( tPointi p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%d", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } /* Reads in the coordinates of the vertices of a polygon from stdin, puts them into P, and returns n, the number of vertices. Formatting conventions: etc. */ int ReadPoly( tPolygoni P ) { int n = 0; printf("Polygon:\n"); printf(" i x y\n"); while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) { printf("%3d%4d%4d\n", n, P[n][0], P[n][1]); ++n; } if (n < PMAX) printf("n = %3d vertices read\n",n); else printf("Error in read_poly: too many points; max is %d\n", PMAX); putchar('\n'); return n; } void PrintPoly( int n, tPolygoni P, int labels[] ) { int i; printf("Polygon:\n"); printf(" i l x y\n"); for( i = 0; i < n; i++ ) printf("%3d%4d%4d%4d\n", i, labels[i], P[i][0], P[i][1]); } SHAR_EOF fi # end of overwriting check echo shar: extracting "'graham.c'" '(3937 characters)' if test -f 'graham.c' then echo shar: will not over-write existing file "'graham.c'" else cat << \SHAR_EOF > 'graham.c' /* Written by Joseph O'Rourke. Graham's algorithm for hull of 2-dim points. See "Computational Geometry in C." */ #include #include #include /* #include */ #define ERROR 1 #define X 0 #define Y 1 typedef enum { FALSE, TRUE } bool; #define DIM 2 /* Dimension of points */ typedef int tPointi[DIM]; /* type integer point */ typedef double tPointd[DIM]; /* type double point */ #define PMAX 1000 /* Max # of pts in polygon */ typedef tPointi tPolygoni[PMAX];/* type integer polygon */ typedef struct tCell *tStack; struct tCell { int i; tStack next; }; tStack Pop( tStack p ); tStack Push( tStack t, tStack p ); tStack GetCell( void ); void PrintStack( tStack t ); tStack GetPush( int i, tStack t ); tStack Graham( tPolygoni P ); /*----------from tri.c-------------*/ int Area2( tPointi a, tPointi b, tPointi c ); bool xor( bool x, bool y ); bool Left( tPointi a, tPointi b, tPointi c ); bool LeftOn( tPointi a, tPointi b, tPointi c ); void SubVec( tPointi a, tPointi b, tPointi c ); int Dot( tPointi a, tPointi b ); int Length2( tPointi p ); void PointAssign( tPointi a, tPointi b ); int ReadPoly( tPolygoni P ); void PrintPoly( int n, tPolygoni P ); void PrintPoint( tPointi p ); /*----------from tri.c-------------*/ int Compare( tPointi *p1, tPointi *p2 ); void FindLowest( int n, tPolygoni P ); static tPolygoni P; /* global so compare can access it. */ int n; main() { tStack top; n = ReadPoly( P ); FindLowest( n, P ); PrintPoly( n, P ); qsort( (void *) P[1], /* base: pointer to 1st elem */ n-1, /* count: # of elems */ sizeof( tPointi ), /* size of each elem */ Compare /* -1,0,+1 compare fnc */ ); PrintPoly( n, P ); top = NULL; top = Graham( P ); printf("Hull:\n"); PrintStack( top ); } /* FindLowest finds the rightmost lowest point and swaps with 0-th. The lowest point has the min y-coord, and amongst those, the max x-coord: so it is rightmost among the lowest. */ void FindLowest( int n, tPolygoni P ) { int i; int m; /* Index of lowest so far. */ tPointi low; /* To hold point when swapping. */ m = 0; for ( i = 1; i < n; i++ ) if ( (P[i][Y] < P[m][Y]) || ((P[i][Y] == P[m][Y]) && (P[i][X] > P[m][X])) ) m = i; /* swap */ PointAssign( low, P[m] ); PointAssign( P[m], P[0] ); PointAssign( P[0], low ); } /* Compare: returns -1,0,+1 if p1 < p2, =, or > respectively; here "<" means smaller angle. Follows the conventions of qsort. */ int Compare( tPointi *p1, tPointi *p2 ) { int a; /* area */ tPointi r1, r2; /* ri = pi - p0 */ int len1, len2; /* length of r1 & r2 */ a = Area2( P[0], *p1, *p2 ); if (a > 0) return -1; else if (a < 0) return 1; else { SubVec( *p1, P[0], r1 ); SubVec( *p2, P[0], r2 ); len1 = Length2( r1 ); len2 = Length2( r2 ); printf("compare:\n"); PrintPoint(r1);PrintPoint(r2); printf("\nlen1=%d, len2=%d\n", len1, len2); return (len1 < len2) ? -1 : (len1 > len2) ? 1 : 0; } } #include "stack.c" /* Get a cell, fill it with i, and push it onto the stack. Return pointer to stack top. */ tStack GetPush( int i, tStack top ) { tStack p; p = GetCell(); p->i = i; return Push( top, p ); } /* Performs the Graham scan on an array of angularly sorted points P. */ tStack Graham( tPolygoni P ) { int i; tStack top; /* Initialize stack to (P[n-1], P[0]). */ top = NULL; top = GetPush ( n-1, top ); top = GetPush ( 0, top ); /* Bottom two elements will never be removed. */ i = 1; while ( i < n ) { printf("Stack at top of while loop, i=%d:\n", i); PrintStack( top ); if ( Left( P[ (top->next)->i ], P[ top->i ], P[ i ] ) ) { top = GetPush ( i, top ); i++; } else top = Pop( top ); printf("Stack at bot of while loop, i=%d:\n", i); PrintStack( top ); putchar('\n'); } /* P[n-1] pushed twice, so pop it off. */ return Pop(top); } #include "vector.c" SHAR_EOF fi # end of overwriting check echo shar: extracting "'i.poly.18'" '(94 characters)' if test -f 'i.poly.18' then echo shar: will not over-write existing file "'i.poly.18'" else cat << \SHAR_EOF > 'i.poly.18' 0 0 10 7 12 3 20 8 13 17 10 12 12 14 13 11 7 11 6 14 10 15 6 18 -1 15 1 13 4 14 5 10 -2 9 5 5 SHAR_EOF fi # end of overwriting check echo shar: extracting "'inpoly.c'" '(2838 characters)' if test -f 'inpoly.c' then echo shar: will not over-write existing file "'inpoly.c'" else cat << \SHAR_EOF > 'inpoly.c' /* inpoly.c Written by Joseph O'Rourke orourke@sophia.smith.edu */ #include #include #define ERROR 1 #define X 0 #define Y 1 typedef enum { FALSE, TRUE } bool; #define DIM 2 /* Dimension of points */ typedef int tPointi[DIM]; /* type integer point */ typedef double tPointd[DIM]; /* type double point */ #define PMAX 1000 /* Max # of pts in polygon */ typedef tPointi tPolygoni[PMAX];/* type integer polygon */ void PrintPoly( int n, tPolygoni P ); void PrintPoint( tPointi p ); bool InPoly( tPointi q, tPolygoni P, int n ); main() { int n; tPolygoni P; tPointi q; n = ReadPoly( P ); /* Last point is the point q */ q[0] = P[n-1][0]; q[1] = P[n-1][1]; n = n - 1; PrintPoly( n, P ); printf("In = %d\n", InPoly( q, P, n )); } /* Returns true if q is inside polygon P. */ bool InPoly( tPointi q, tPolygoni P, int n ) { int i, i1; /* point index; i1 = i-1 mod n */ int d; /* dimension index */ double x; /* x intersection of e with ray */ int crossings = 0; /* number of edge/ray crossings */ printf("\n==>In: u = "); PrintPoint(q); putchar('\n'); /* Shift so that q is the origin. */ for( i = 0; i < n; i++ ) { for( d = 0; d < DIM; d++ ) P[i][d] = P[i][d] - q[d]; } /* For each edge e=(i-1,i), see if crosses ray. */ for( i = 0; i < n; i++ ) { i1 = ( i + n - 1 ) % n; printf("e=(%d,%d)\t", i1, i); /* if e straddles the x-axis... */ if( ( ( P[i] [Y] > 0 ) && ( P[i1][Y] <= 0 ) ) || ( ( P[i1][Y] > 0 ) && ( P[i] [Y] <= 0 ) ) ) { /* e straddles ray, so compute intersection with ray. */ x = (P[i][X] * P[i1][Y] - P[i1][X] * P[i][Y]) / (double)(P[i1][Y] - P[i][Y]); printf("straddles: x = %g\t", x); /* crosses ray if strictly positive intersection. */ if (x > 0) crossings++; } printf("crossings=%d\n", crossings); } /* q inside if an odd number of crossings. */ if( (crossings % 2) == 1 ) return TRUE; else return FALSE; } void PrintPoint( tPointi p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%d", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } /* Reads in the coordinates of the vertices of a polygon from stdin, puts them into P, and returns n, the number of vertices. Formatting conventions: etc. */ int ReadPoly( tPolygoni P ) { int n = 0; printf("Polygon:\n"); printf(" i x y\n"); while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) { printf("%3d%4d%4d\n", n, P[n][0], P[n][1]); ++n; } if (n < PMAX) printf("n = %3d vertices read\n",n); else printf("Error in read_poly: too many points; max is %d\n", PMAX); putchar('\n'); return n; } void PrintPoly( int n, tPolygoni P ) { int i; printf("Polygon:\n"); printf(" i x y\n"); for( i = 0; i < n; i++ ) printf("%3d%4d%4d\n", i, P[i][0], P[i][1]); } SHAR_EOF fi # end of overwriting check echo shar: extracting "'int.c'" '(7635 characters)' if test -f 'int.c' then echo shar: will not over-write existing file "'int.c'" else cat << \SHAR_EOF > 'int.c' #include #include #define ERROR 1 #define X 0 #define Y 1 typedef enum { FALSE, TRUE } bool; typedef enum { Pin, Qin, Unknown } tInFlag; #define DIM 2 /* Dimension of points */ typedef int tPointi[DIM]; /* type integer point */ typedef double tPointd[DIM]; /* type double point */ #define PMAX 1000 /* Max # of pts in polygon */ typedef tPointi tPolygoni[PMAX];/* type integer polygon */ int Area2( tPointi a, tPointi b, tPointi c ); void SubVec( tPointi a, tPointi b, tPointi c ); int Dot( tPointi a, tPointi b ); bool LeftOn( tPointi a, tPointi b, tPointi c ); void PrintPoly( int n, tPolygoni P ); void PrintPoint( tPointi p ); void PrintPointd( tPointd p ); void ConvexIntersect( tPolygoni P, tPolygoni Q, int n, int m ); tInFlag InOut( tPointd p, tInFlag inflag, bool aHB ); int Advance( int a, int *aadv, int n, bool inside, tPointi v ); bool Intersection( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p ); main() { int n, m; tPolygoni P, Q; n = ReadPoly( P ); m = ReadPoly( Q ); ConvexIntersect( P, Q, n, m ); } void ConvexIntersect( tPolygoni P, tPolygoni Q, int n, int m ) /* P has n vertices, Q has m vertices. */ { int a, b; /* indices on P and Q (resp.) */ int a1, b1; /* a-1, b-1 (resp.) */ tPointi A, B; /* directed edges on P and Q (resp.) */ int cross; /* A x B */ bool bHA, aHB; /* b in H(A); a in H(b). */ tPointi Origin = {0,0}; /* (0,0) */ tPointd p; /* double point of intersection */ int i; /* loop counter */ tInFlag inflag; /* {Pin, Qin, Unknown}: which polygon is known to be inside */ int aadv, badv; /* # advances on a & b indices (from first intersection) */ /* Initialize variables. */ a = 0; b = 0; i = 0; inflag = Unknown; do { /* Computations of key variables. */ a1 = (a + n - 1) % n; b1 = (b + m - 1) % m; SubVec( P[a], P[a1], A ); SubVec( Q[b], Q[b1], B ); cross = Area2( Origin, A, B ); bHA = LeftOn( P[a1], P[a], Q[b] ); aHB = LeftOn( Q[b1], Q[b], P[a] ); /* If A & B intersect, update inflag. */ if ( Intersection( P[a1], P[a], Q[b1], Q[b], p ) ) { if ( inflag == Unknown ) { aadv = 0; badv = 0; } inflag = InOut( p, inflag, aHB ); } /* Advance rules. */ if ( cross >= 0 ) if ( bHA ) a = Advance( a, &aadv, n, inflag == Pin, P[a] ); else { b = Advance( b, &badv, m, inflag == Qin, Q[b] ); } else { if ( aHB ) b = Advance( b, &badv, m, inflag == Qin, Q[b] ); else a = Advance( a, &aadv, n, inflag == Pin, P[a] ); } /* Quit when both indices have cycled. */ } while ( (aadv < n) || (badv < m) ); /* Deal with special cases: not implemented here. */ if ( inflag == Unknown) { printf("The boundaries of P and Q do not intersect:\n"); printf("\tEither P and Q do not intersect, or\n"); printf("\tone properly contains the other\n"); } } /* Prints out the double point of intersection, and toggles in/out flag. */ tInFlag InOut( tPointd p, tInFlag inflag, bool aHB ) { PrintPointd( p ); putchar('\n'); /* Update inflag. */ if ( aHB ) return Pin; else return Qin; } /* Advances and prints out an inside vertex if appropriate. */ int Advance( int a, int *aadv, int n, bool inside, tPointi v ) { if ( inside ) { PrintPoint( v ); putchar('\n'); } (*aadv)++; return (a+1) % n; } bool Intersection( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p ) { double s, t; /* The two parameters of the parametric eqns. */ double denom; /* Denoninator of solutions. */ denom = a[0] * ( d[1] - c[1] ) + b[0] * ( c[1] - d[1] ) + d[0] * ( b[1] - a[1] ) + c[0] * ( a[1] - b[1] ); /* If denom is zero, then the line segments are parallel. */ /* In this case, return false even though the segments might overlap. */ if (denom == 0.0) return FALSE; s = ( a[0] * ( d[1] - c[1] ) + c[0] * ( a[1] - d[1] ) + d[0] * ( c[1] - a[1] ) ) / denom; t = -( a[0] * ( c[1] - b[1] ) + b[0] * ( a[1] - c[1] ) + c[0] * ( b[1] - a[1] ) ) / denom; p[0] = a[0] + s * ( b[0] - a[0] ); p[1] = a[1] + s * ( b[1] - a[1] ); if ( (0.0 <= s) && (s <= 1.0) && (0.0 <= t) && (t <= 1.0) ) return TRUE; else return FALSE; } /* Implements a = b, assignment of points/vectors. Assignment between arrays is not possible in C. */ void PointAssign( tPointi a, tPointi b ) { int i; for ( i = 0; i < DIM; i++ ) a[i] = b[i]; } void PrintPoint( tPointi p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%d", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } void PrintPointd( tPointd p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%5.2f", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } /* Reads in the coordinates of the vertices of a polygon from stdin, puts them into P, and returns n, the number of vertices. Formatting conventions: etc. */ int ReadPoly( tPolygoni P ) { int n = 0; int nin; scanf("%d", &nin); printf("Polygon:\n"); printf(" i x y\n"); while ( (n < nin) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) { printf("%3d%4d%4d\n", n, P[n][0], P[n][1]); ++n; } if (n < PMAX) printf("n = %3d vertices read\n",n); else printf("Error in read_poly: too many points; max is %d\n", PMAX); putchar('\n'); return n; } void PrintPoly( int n, tPolygoni P ) { int i; printf("Polygon:\n"); printf(" i l x y\n"); for( i = 0; i < n; i++ ) printf("%3d%4d%4d%4d\n", i, P[i][0], P[i][1]); } /* Returns true iff c is strictly to the left of the directed line through a to b. */ bool Left( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) > 0; } bool LeftOn( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) >= 0; } bool Collinear( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) == 0; } /* Puts a - b into c. */ void SubVec( tPointi a, tPointi b, tPointi c ) { int i; for( i = 0; i < DIM; i++ ) c[i] = a[i] - b[i]; } /* Returns twice the signed area of the triangle determined by a,b,c, positive if a,b,c are oriented ccw, and negative if cw. */ int Area2( tPointi a, tPointi b, tPointi c ) { return a[0] * b[1] - a[1] * b[0] + a[1] * c[0] - a[0] * c[1] + b[0] * c[1] - c[0] * b[1]; } SHAR_EOF fi # end of overwriting check echo shar: extracting "'macros.h'" '(956 characters)' if test -f 'macros.h' then echo shar: will not over-write existing file "'macros.h'" else cat << \SHAR_EOF > 'macros.h' /*==================================================================== macros.h macros used to access data structures and perform quick tests. ====================================================================*/ /* general-purpose macros */ #define SWAP(t,x,y) { t = x; x = y; y = t; } char *malloc(); #define NEW(p,type) if ((p=(type *) malloc (sizeof(type))) == NULL) {\ printf ("Out of Memory!\n");\ exit(0);\ } #define FREE(p) if (p) { free ((char *) p); p = NULL; } #define ADD( head, p ) if ( head ) { \ p->next = head->next; \ p->prev = head; \ head->next = p; \ p->next->prev = p; \ } \ else { \ head = p; \ head->next = head->prev = p; \ } #define DELETE( head, p ) if ( head ) { \ if ( head == head->next ) \ head = NULL; \ else if ( p == head ) \ head = head->next; \ p->next->prev = p->prev; \ p->prev->next = p->next; \ FREE( p ); \ } SHAR_EOF fi # end of overwriting check echo shar: extracting "'quadratic.c'" '(10343 characters)' if test -f 'quadratic.c' then echo shar: will not over-write existing file "'quadratic.c'" else cat << \SHAR_EOF > 'quadratic.c' /* Written by Joseph O'Rourke. Last modified February 17, 1991. Code to triangulate a simple polygon. Quadratic algorithm. See "Computational Geometry in C." */ #include #include #define ERROR 1 #define X 0 #define Y 1 typedef enum { FALSE, TRUE } bool; #define DIM 2 /* Dimension of points */ typedef int tPointi[DIM]; /* type integer point */ typedef double tPointd[DIM]; /* type double point */ #define PMAX 1000 /* Max # of pts in polygon */ typedef tPointi tPolygoni[PMAX];/* type integer polygon */ int Area2( tPointi a, tPointi b, tPointi c ); int AreaPoly2( int n, tPolygoni P ); bool xor( bool x, bool y ); bool IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d ); bool Left( tPointi a, tPointi b, tPointi c ); bool LeftOn( tPointi a, tPointi b, tPointi c ); void SubVec( tPointi a, tPointi b, tPointi c ); int Dot( tPointi a, tPointi b ); int Length2( tPointi p ); bool Between( tPointi a, tPointi b, tPointi c ); bool Intersect( tPointi a, tPointi b, tPointi c, tPointi d ); bool Diagonalie( int i, int j, int n, tPolygoni P ); bool Diagonal( int i, int j, int n, tPolygoni P ); void PointAssign( tPointi a, tPointi b ); void Triangulate2( int n, tPolygoni P, bool Ear[], int labels[]); void Triangulate( int n, tPolygoni P ); void TriRecurse( int n, tPolygoni P, int labels[] ); bool InCone( int i, int j, int n, tPolygoni P ); void ClipEar1( int i, int n, tPolygoni P, int labels[] ); void ClipEar2( int i, int n, tPolygoni P, bool Ear[], int labels[] ); void ClipEar( int i, int n, tPolygoni P ); int ReadPoly( tPolygoni P ); void PrintPoly( int n, tPolygoni P, int labels[] ); void PrintPoint( tPointi p ); void UpdateArray( bool Ear[PMAX], int j, int n); main() { int n; tPolygoni P; n = ReadPoly( P ); Triangulate( n, P ); } /* Returns twice the signed area of the triangle determined by a,b,c. The area is positive if a,b,c are oriented ccw, negative if cw, and zero if the points are collinear. */ int Area2( tPointi a, tPointi b, tPointi c ) { return a[0] * b[1] - a[1] * b[0] + a[1] * c[0] - a[0] * c[1] + b[0] * c[1] - c[0] * b[1]; } /* Returns twice the area of polygon P. */ int AreaPoly2( int n, tPolygoni P ) { int i; int sum = 0; /* Partial area sum */ for (i = 1; i < n-1; i++) sum += Area2( P[0], P[i], P[i+1] ); return sum; } /* Exclusive or: true iff exactly one argument is true. The arguments are negated to ensure that they are 0/1 values. Then the bitwise xor operator may apply. (This idea is due to Michael Baldwin.) */ bool xor( bool x, bool y ) { return !x ^ !y; } /* Returns true iff ab properly intersects cd: they share a point interior to both segments. The properness of the intersection is ensured by using strict leftness. */ bool IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d ) { return xor( Left(a,b,c), Left(a,b,d) ) && xor( Left(c,d,a), Left(c,d,b) ); } /* Returns true iff c is strictly to the left of the directed line through a to b. */ bool Left( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) > 0; } bool LeftOn( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) >= 0; } bool Collinear( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) == 0; } /* Puts a - b into c. */ void SubVec( tPointi a, tPointi b, tPointi c ) { int i; for( i = 0; i < DIM; i++ ) c[i] = a[i] - b[i]; } /* Returns the dot product of the two input vectors. */ int Dot( tPointi a, tPointi b ) { int i; int sum = 0; for( i = 0; i < DIM; i++ ) sum += a[i] * b[i]; return sum; } /* Returns the square of the length of the vector p. */ int Length2( tPointi p ) { return Dot( p, p ); } /* Returns T iff point c lies on the closed segement ab. First checks that c is collinear with a and b. */ bool Between( tPointi a, tPointi b, tPointi c ) { tPointi ba, ca; if ( ! Collinear( a, b, c ) ) return FALSE; SubVec( b, a, ba ); SubVec( c, a, ca ); return Dot( ca, ba ) >= 0 && Length2( ca ) <= Length2( ba ); } /* Returns true iff segments ab and cd intersect, properly or improperly. */ bool Intersect( tPointi a, tPointi b, tPointi c, tPointi d ) { if ( IntersectProp( a, b, c, d ) ) return TRUE; else if ( Between( a, b, c ) || Between( a, b, d ) || Between( c, d, a ) || Between( c, d, b ) ) return TRUE; else return FALSE; } /* Returns T iff (v_i, v_j) is a proper internal *or* external diagonal of P, *ignoring edges incident to v_i and v_j*. */ bool Diagonalie( int i, int j, int n, tPolygoni P ) { int k; int k1; /* For each edge (k,k+1) of P */ for ( k = 0; k < n; k++ ) { k1 = (k+1) % n; /* Skip edges incident to i or j */ if ( ! ( ( k == i ) || ( k1 == i ) || ( k == j ) || ( k1 == j ) ) ) if ( Intersect( P[i], P[j], P[k], P[k1] ) ) return FALSE; } return TRUE; } /* Implements a = b, assignment of points/vectors. Assignment between arrays is not possible in C. */ void PointAssign( tPointi a, tPointi b ) { int i; for ( i = 0; i < DIM; i++ ) a[i] = b[i]; } void TriRecurse( int n, tPolygoni P, int labels[] ) { int i, i1, i2; printf("TriRecurse: n = %d\n", n); /*PrintPoly( n, P, labels );*/ if ( n > 3 ) for ( i = 0; i < n; i++ ) { i1 = (i+1) % n; i2 = (i+2) % n; if ( Diagonal( i, i2, n, P ) ) { printf("%d %d\n", labels[i], labels[i2]); ClipEar1( i1, n, P, labels ); TriRecurse( n-1, P, labels ); break; } } } void PrintPoint( tPointi p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%d", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } /* Prints out n-3 diagonals (as pairs of integer indices) which form a triangulation of P. This function initializes the data structures, and calls Triangulate2 to clip off the ears one by one. */ void Triangulate( int n, tPolygoni P ) { tPolygoni Pt; bool Ear[PMAX]; /* Ear[i] true iff (i,i+1,i+2) is an ear */ int labels[PMAX]; int i, i1, i2; /* Copy points into temporary array & initialize labels. */ for ( i = 0; i < n; i++ ){ PointAssign( Pt[i], P[i] ); labels[i] = i; } /* Ear[i] is true iff (i,i+1,i+2) is an ear. */ /* Initialize Ear[] for all i. */ if ( n > 3 ) for ( i = 0; i < n; i++ ) { i1 = (i+1) % n; i2 = (i+2) % n; Ear[i] = Diagonal( i, i2, n, P ); } Triangulate2( n, Pt, Ear, labels ); } /* Triangulate2 is an O(n^2) triangulation function (or it would be if the array squashing were replaced by pointer movements). Assumes all arrays (including Ear) have been initialized. */ void Triangulate2( int n, tPolygoni P , bool Ear[], int labels[] ) { int im1, i, i1, i2, i3; /* i-1,i,i+1,i+2,i+3 */ /* Each step of outer loop removes one ear. */ while ( n > 3 ) { /* Inner loop searches for an ear. */ for (i = 0; i < n; i++) if (Ear[i]) { /* Ear found. */ im1 = (i-1+n) % n; i1 = (i+1) % n; i2 = (i+2) % n; i3 = (i+3) % n; /* (i,i+2) is a diagonal */ printf("%d %d\n", labels[i], labels[i2]); /* Replace the two entries for segments incident to i+1: (i-1,i+1) ==> (i-1,i+2) (i+1,i+3) ==> (i,i+3) */ Ear[im1] = Diagonal( im1, i2, n, P); Ear[i] = Diagonal( i, i3, n, P); /* Squash out the i1 entry in all arrays */ ClipEar2( i1, n, P, Ear, labels ); n = n - 1; break; /* out of inner loop */ } } } void ClipEar2( int i, int n, tPolygoni P, bool Ear[], int labels[] ) { int k; for ( k = i; k < n-1; k++ ) { PointAssign( P[k], P[k+1] ); Ear[k] = Ear[k+1]; labels[k] = labels[k+1]; } } /* Returns true iff the diagonal (i,j) is striclty internal to the polygon P in the neighborhood of the i endpoint. */ bool InCone( int i, int j, int n, tPolygoni P ) { int i1; /* i + 1 */ int in1; /* i - 1 */ i1 = (i + 1) % n; in1 = (i + n - 1) % n; /* If P[i] is a convex vertex [ i+1 left of (i-1,i) ]. */ if ( Left( P[in1], P[i], P[i1] ) ) return Left( P[i], P[j], P[in1] ) && Left( P[j], P[i], P[i1] ); /* Assume (i-1,i,i+1) not collinear. */ /* else P[i] is reflex. */ else return !( LeftOn( P[i], P[j], P[i1] ) && LeftOn( P[j], P[i], P[in1] ) ); } /* Returns T iff (v_i, v_j) is a proper internal diagonal of P. */ bool Diagonal( int i, int j, int n, tPolygoni P ) { return InCone( i, j, n, P ) && Diagonalie( i, j, n, P ); } /* Removes P[i] by copying P[i+1]...P[n-1] left one index. */ void ClipEar1( int i, int n, tPolygoni P, int labels[] ) { int k; for ( k = i; k < n-1; k++ ) { PointAssign( P[k], P[k+1] ); labels[k] = labels[k+1]; } } void ClipEar( int i, int n, tPolygoni P ) { int k; for ( k = i; k < n-1; k++ ) PointAssign( P[k], P[k+1] ); } /* Reads in the coordinates of the vertices of a polygon from stdin, puts them into P, and returns n, the number of vertices. Formatting conventions: etc. */ int ReadPoly( tPolygoni P ) { int n = 0; printf("Polygon:\n"); printf(" i x y\n"); while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) { printf("%3d%4d%4d\n", n, P[n][0], P[n][1]); ++n; } if (n < PMAX) printf("n = %3d vertices read\n",n); else printf("Error in read_poly: too many points; max is %d\n", PMAX); putchar('\n'); return n; } void PrintPoly( int n, tPolygoni P, int labels[] ) { int i; printf("Polygon:\n"); printf(" i l x y\n"); for( i = 0; i < n; i++ ) printf("%3d%4d%4d%4d\n", i, labels[i], P[i][0], P[i][1]); } SHAR_EOF fi # end of overwriting check echo shar: extracting "'sphere.c'" '(1365 characters)' if test -f 'sphere.c' then echo shar: will not over-write existing file "'sphere.c'" else cat << \SHAR_EOF > 'sphere.c' /*---------------------------------------------------------------------- Sphere Susan L. Weller ---------------------------------------------------------------------- This program will generate a given number of points that lie on the unit sphere. The number of points is given on the command line. ---------------------------------------------------------------------- */ #include #include #define MAX_INT 2147483647 /* Max int is 2^31 - 1 */ #define FUZZ 0.000001 main( argc, argv ) int argc; char *argv[]; { int n; double x, y, z, r; srandom( (int) time((long *) 0 ) ); if ( argc != 2 ) printf( "Usage: sphere \n" ), exit(1); n = atoi( argv[1] ); /* n is number of points */ while (n--) { do { /* the do-while will generate points until one is inside the unit sphere. */ x = 2.0 * (double) random() / MAX_INT - 1.0; y = 2.0 * (double) random() / MAX_INT - 1.0; z = 2.0 * (double) random() / MAX_INT - 1.0; r = sqrt( x*x + y*y + z*z ); } while ( r <= FUZZ || r > 1.); printf ("%6d %6d %6d\n", (int) (100 * x/r), (int) (100 * y/r), (int) (100 * z/r) ); } } SHAR_EOF fi # end of overwriting check echo shar: extracting "'stack.c'" '(861 characters)' if test -f 'stack.c' then echo shar: will not over-write existing file "'stack.c'" else cat << \SHAR_EOF > 'stack.c' /*----------------------Stack routines------------------*/ /* Pushes cell p on top of stack t, and returns new top. */ tStack Push( tStack t, tStack p ) { p->next = t; return p; } /* Pops off top elment of stack p, frees up the cell, and returns new top. */ tStack Pop( tStack p ) { tStack top; top = p->next; free( (void *) p ); return top; } /* Prints the stack, both point index and point coordinates. */ void PrintStack( tStack t ) { if (t) { printf("i=%d:", t->i); PrintPoint( P[t->i] ); putchar('\n'); PrintStack( t->next ); } } /* GetCell returns a pointer to a newly allocated piece of storage of type tCell, or exits if no space is available. */ tStack GetCell( void ) { tStack p; p = (tStack) malloc( sizeof( struct tCell ) ); if (p == NULL) { printf("Error in GetCell: Out of memory!\n"); exit(1); } else return p; } SHAR_EOF fi # end of overwriting check echo shar: extracting "'tri.c'" '(9289 characters)' if test -f 'tri.c' then echo shar: will not over-write existing file "'tri.c'" else cat << \SHAR_EOF > 'tri.c' /* tri.c Written by Joseph O'Rourke orourke@cs.smith.edu Triangulation code for "Computational Geometry in C." Assumes polygon vertices are entered in ccw order. */ #include #include #define ERROR 1 #define X 0 #define Y 1 typedef enum { FALSE, TRUE } bool; #define DIM 2 /* Dimension of points */ typedef int tPointi[DIM]; /* type integer point */ typedef double tPointd[DIM]; /* type double point */ #define PMAX 1000 /* Max # of pts in polygon */ typedef tPointi tPolygoni[PMAX];/* type integer polygon */ int Area2( tPointi a, tPointi b, tPointi c ); int AreaPoly2( int n, tPolygoni P ); bool Xor( bool x, bool y ); bool IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d ); bool Left( tPointi a, tPointi b, tPointi c ); bool LeftOn( tPointi a, tPointi b, tPointi c ); bool Collinear( tPointi a, tPointi b, tPointi c ); void SubVec( tPointi a, tPointi b, tPointi c ); int Dot( tPointi a, tPointi b ); int Length2( tPointi p ); bool Between( tPointi a, tPointi b, tPointi c ); bool Intersect( tPointi a, tPointi b, tPointi c, tPointi d ); bool Diagonalie( int i, int j, int n, tPolygoni P ); bool Diagonal( int i, int j, int n, tPolygoni P ); void PointAssign( tPointi a, tPointi b ); void Triangulate1( int n, tPolygoni P ); void Triangulate( int n, tPolygoni P ); void TriRecurse( int n, tPolygoni P, int labels[] ); bool InCone( int i, int j, int n, tPolygoni P ); void ClipEar1( int i, int n, tPolygoni P, int labels[] ); void ClipEar( int i, int n, tPolygoni P ); int ReadPoints( tPolygoni P ); void PrintPoly( int n, tPolygoni P, int labels[] ); void PrintPoint( tPointi p ); main() { int n; tPolygoni P; n = ReadPoints( P ); Triangulate1( n, P ); Triangulate( n, P ); } /* Returns twice the signed area of the triangle determined by a,b,c, positive if a,b,c are oriented ccw, and negative if cw. */ int Area2( tPointi a, tPointi b, tPointi c ) { /* The text has this: return a[0] * b[1] - a[1] * b[0] + a[1] * c[0] - a[0] * c[1] + b[0] * c[1] - c[0] * b[1]; The following computation is algebraically equivalent but uses four fewer multiplications. It is obtained by shifting the coordinate system so that point a is the origin. */ return (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); } /* Returns twice the area of polygon P. */ int AreaPoly2( int n, tPolygoni P ) { int i; int sum = 0; /* Partial area sum */ for (i = 1; i < n-1; i++) sum += Area2( P[0], P[i], P[i+1] ); return sum; } /* Exclusive or: true iff exactly one argument is true. The arguments are negated to ensure that they are 0/1 values. Then the bitwise Xor operator may apply. (This idea is due to Michael Baldwin.) */ bool Xor( bool x, bool y ) { return !x ^ !y; } /* Returns true iff ab properly intersects cd: they share a point interior to both segments. The properness of the intersection is ensured by using strict leftness. */ bool IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d ) { /* Eliminate improper cases. */ if ( Collinear(a,b,c) || Collinear(a,b,d) || Collinear(c,d,a) || Collinear(c,d,b) ) return FALSE; return Xor( Left(a,b,c), Left(a,b,d) ) && Xor( Left(c,d,a), Left(c,d,b) ); } /* Returns true iff c is strictly to the left of the directed line through a to b. */ bool Left( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) > 0; } bool LeftOn( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) >= 0; } bool Collinear( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) == 0; } /* Returns T iff (a,b,c) are collinear and point c lies on the closed segement ab. */ bool Between( tPointi a, tPointi b, tPointi c ) { tPointi ba, ca; if ( ! Collinear( a, b, c ) ) return FALSE; /* If ab not vertical, check betweenness on x; else on y. */ if ( a[X] != b[X] ) return ((a[X] <= c[X]) && (c[X] <= b[X])) || ((a[X] >= c[X]) && (c[X] >= b[X])); else return ((a[Y] <= c[Y]) && (c[Y] <= b[Y])) || ((a[Y] >= c[Y]) && (c[Y] >= b[Y])); } /* Returns true iff segments ab and cd intersect, properly or improperly. */ bool Intersect( tPointi a, tPointi b, tPointi c, tPointi d ) { if ( IntersectProp( a, b, c, d ) ) return TRUE; else if ( Between( a, b, c ) || Between( a, b, d ) || Between( c, d, a ) || Between( c, d, b ) ) return TRUE; else return FALSE; } /* Returns T iff (v_i, v_j) is a proper internal *or* external diagonal of P, *ignoring edges incident to v_i and v_j*. */ bool Diagonalie( int i, int j, int n, tPolygoni P ) { int k; int k1; /* For each edge (k,k+1) of P */ for ( k = 0; k < n; k++ ) { k1 = (k+1) % n; /* Skip edges incident to i or j */ if ( ! ( ( k == i ) || ( k1 == i ) || ( k == j ) || ( k1 == j ) ) ) if ( Intersect( P[i], P[j], P[k], P[k1] ) ) return FALSE; } return TRUE; } /* Implements a = b, assignment of points/vectors. Assignment between arrays is not possible in C. */ void PointAssign( tPointi a, tPointi b ) { int i; for ( i = 0; i < DIM; i++ ) a[i] = b[i]; } /* Prints out n-3 diagonals (as pairs of integer indices) which form a triangulation of P. */ void Triangulate1( int n, tPolygoni P ) { tPolygoni Pt; int labels[PMAX]; int i; for ( i = 0; i < n; i++ ){ PointAssign( Pt[i], P[i] ); labels[i] = i; } TriRecurse( n, Pt, labels ); } void TriRecurse( int n, tPolygoni P, int labels[] ) { int i, i1, i2; printf("TriRecurse: n = %d\n", n); /*PrintPoly( n, P, labels );*/ if ( n > 3 ) for ( i = 0; i < n; i++ ) { i1 = (i+1) % n; i2 = (i+2) % n; if ( Diagonal( i, i2, n, P ) ) { printf("%d %d\n", labels[i], labels[i2]); ClipEar1( i1, n, P, labels ); TriRecurse( n-1, P, labels ); break; } } } void PrintPoint( tPointi p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%d", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } void Triangulate( int n, tPolygoni P ) { int i, i1, i2; if ( n > 3 ) for ( i = 0; i < n; i++ ) { i1 = (i+1) % n; i2 = (i+2) % n; if ( Diagonal( i, i2, n, P ) ) { PrintPoint( P[i] ); putchar('\t'); PrintPoint( P[i2] ); putchar('\n'); ClipEar( i1, n, P ); Triangulate( n-1, P ); break; } } } /* Returns true iff the diagonal (i,j) is striclty internal to the polygon P in the neighborhood of the i endpoint. */ bool InCone( int i, int j, int n, tPolygoni P ) { int i1; /* i + 1 */ int in1; /* i - 1 */ i1 = (i + 1) % n; in1 = (i + n - 1) % n; /* If P[i] is a convex vertex [ i+1 left or on (i-1,i) ]. */ if ( LeftOn( P[in1], P[i], P[i1] ) ) return Left( P[i], P[j], P[in1] ) && Left( P[j], P[i], P[i1] ); /* Assume (i-1,i,i+1) not collinear. */ /* else P[i] is reflex. */ else return !( LeftOn( P[i], P[j], P[i1] ) && LeftOn( P[j], P[i], P[in1] ) ); } /* Returns T iff (v_i, v_j) is a proper internal diagonal of P. */ bool Diagonal( int i, int j, int n, tPolygoni P ) { return InCone( i, j, n, P ) && Diagonalie( i, j, n, P ); } /* Removes P[i] by copying P[i+1]...P[n-1] left one index. */ void ClipEar1( int i, int n, tPolygoni P, int labels[] ) { int k; for ( k = i; k < n-1; k++ ) { PointAssign( P[k], P[k+1] ); labels[k] = labels[k+1]; } } void ClipEar( int i, int n, tPolygoni P ) { int k; for ( k = i; k < n-1; k++ ) PointAssign( P[k], P[k+1] ); } int ReadPoints( tPolygoni P ) /* Reads in the coordinates of the vertices of a polygon from stdin, puts them into P, and returns n, the number of vertices. The input is assumed to be pairs of whitespace-separated coordinates, one pair per line. The number of points is not part of the input. */ { int n = 0; printf("Polygon:\n"); printf(" i x y\n"); while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) { printf("%3d%4d%4d\n", n, P[n][0], P[n][1]); ++n; } if (n < PMAX) printf("n = %3d vertices read\n",n); else printf("Error in ReadPoints:\ too many points; max is %d\n", PMAX); putchar('\n'); return n; } void PrintPoly( int n, tPolygoni P, int labels[] ) { int i; printf("Polygon:\n"); printf(" i l x y\n"); for( i = 0; i < n; i++ ) printf("%3d%4d%4d%4d\n", i, labels[i], P[i][0], P[i][1]); } /* Puts a - b into c. */ void SubVec( tPointi a, tPointi b, tPointi c ) { int i; for( i = 0; i < DIM; i++ ) c[i] = a[i] - b[i]; } /* Returns the dot product of the two input vectors. */ int Dot( tPointi a, tPointi b ) { int i; int sum = 0; for( i = 0; i < DIM; i++ ) sum += a[i] * b[i]; return sum; } /* Returns the square of the length of the vector p. */ int Length2( tPointi p ) { return Dot( p, p ); } SHAR_EOF fi # end of overwriting check echo shar: extracting "'vector.c'" '(2570 characters)' if test -f 'vector.c' then echo shar: will not over-write existing file "'vector.c'" else cat << \SHAR_EOF > 'vector.c' /* Returns twice the signed area of the triangle determined by a,b,c. The area is positive if a,b,c are oriented ccw, negative if cw, and zero if the points are collinear. */ int Area2( tPointi a, tPointi b, tPointi c ) { return a[0] * b[1] - a[1] * b[0] + a[1] * c[0] - a[0] * c[1] + b[0] * c[1] - c[0] * b[1]; } /* Exclusive or: true iff exactly one argument is true. The arguments are negated to ensure that they are 0/1 values. Then the bitwise xor operator may apply. (This idea is due to Michael Baldwin.) */ bool xor( bool x, bool y ) { return !x ^ !y; } bool Left( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) > 0; } bool LeftOn( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) >= 0; } bool Collinear( tPointi a, tPointi b, tPointi c ) { return Area2( a, b, c ) == 0; } /* Puts a - b into c. */ void SubVec( tPointi a, tPointi b, tPointi c ) { int i; for( i = 0; i < DIM; i++ ) c[i] = a[i] - b[i]; } /* Returns the dot product of the two input vectors. */ int Dot( tPointi a, tPointi b ) { int i; int sum = 0; for( i = 0; i < DIM; i++ ) sum += sum + a[i] * b[i]; return sum; } /* Returns the square of the length of the vector p. */ int Length2( tPointi p ) { return Dot( p, p ); } /* Implements a = b, assignment of points/vectors. Assignment between arrays is not possible in C. */ void PointAssign( tPointi a, tPointi b ) { int i; for ( i = 0; i < DIM; i++ ) a[i] = b[i]; } void PrintPoint( tPointi p ) { int i; putchar('('); for ( i = 0; i < DIM; i++ ) { printf("%d", p[i]); if ( i != DIM-1 ) putchar(','); } putchar(')'); } /* Reads in the coordinates of the vertices of a polygon from stdin, puts them into P, and returns n, the number of vertices. Formatting conventions: etc. */ int ReadPoly( tPolygoni P ) { int n = 0; printf("Polygon:\n"); printf(" i x y\n"); while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) { printf("%3d%4d%4d\n", n, P[n][0], P[n][1]); ++n; } if (n < PMAX) printf("n = %3d vertices read\n",n); else printf("Error in read_poly: too many points; max is %d\n", PMAX); putchar('\n'); return n; } void PrintPoly( int n, tPolygoni P ) { int i; printf("Polygon:\n"); printf(" i x y\n"); for( i = 0; i < n; i++ ) printf("%3d%4d%4d\n", i, P[i][0], P[i][1]); } SHAR_EOF fi # end of overwriting check # End of shell archive exit 0