/* * C code from the article * "An Implicit Surface Polygonizer" * by Jules Bloomenthal, jbloom@beauty.gmu.edu * in "Graphics Gems IV", Academic Press, 1994 */ /* Modified by Curtis Janssen: * 1. Eliminate memory leaks. * 2. Make main routine optional (-DMAIN to compile a main routine). */ /* implicit.c * an implicit surface polygonizer, translated from Mesa * applications should call polygonize() * * to compile a test program for ASCII output: * cc -DMAIN implicit.c -o implicit -lm * * to compile a test program for display on an SGI workstation: * cc -DMAIN -DSGIGFX implicit.c -o implicit -lgl_s -lm * * Authored by Jules Bloomenthal, Xerox PARC. * Copyright (c) Xerox Corporation, 1991. All rights reserved. * Permission is granted to reproduce, use and distribute this code for * any and all purposes, provided that this notice appears in all copies. */ #include #include #include #include #include #define TET 0 /* use tetrahedral decomposition */ #define NOTET 1 /* no tetrahedral decomposition */ #define RES 10 /* # converge iterations */ #define L 0 /* left direction: -x, -i */ #define R 1 /* right direction: +x, +i */ #define B 2 /* bottom direction: -y, -j */ #define T 3 /* top direction: +y, +j */ #define N 4 /* near direction: -z, -k */ #define F 5 /* far direction: +z, +k */ #define LBN 0 /* left bottom near corner */ #define LBF 1 /* left bottom far corner */ #define LTN 2 /* left top near corner */ #define LTF 3 /* left top far corner */ #define RBN 4 /* right bottom near corner */ #define RBF 5 /* right bottom far corner */ #define RTN 6 /* right top near corner */ #define RTF 7 /* right top far corner */ /* the LBN corner of cube (i, j, k), corresponds with location * (start.x+(i-.5)*size, start.y+(j-.5)*size, start.z+(k-.5)*size) */ #define RAND() ((rand()&32767)/32767.) /* random number between 0 and 1 */ #define HASHBIT (5) #define HASHSIZE (size_t)(1<<(3*HASHBIT)) /* hash table size (32768) */ #define MASK ((1<>(bit))&1) #define FLIP(i,bit) ((i)^1<<(bit)) /* flip the given bit of i */ typedef struct point { /* a three-dimensional point */ double x, y, z; /* its coordinates */ } POINT; typedef struct test { /* test the function for a signed value */ POINT p; /* location of test */ double value; /* function value at p */ int ok; /* if value is of correct sign */ } TEST; typedef struct vertex { /* surface vertex */ POINT position, normal; /* position and surface normal */ } VERTEX; typedef struct vertices { /* list of vertices in polygonization */ int count, max; /* # vertices, max # allowed */ VERTEX *ptr; /* dynamically allocated */ } VERTICES; typedef struct corner { /* corner of a cube */ int i, j, k; /* (i, j, k) is index within lattice */ double x, y, z, value; /* location and function value */ } CORNER; typedef struct cube { /* partitioning cell (cube) */ int i, j, k; /* lattice location of cube */ CORNER *corners[8]; /* eight corners */ } CUBE; typedef struct cubes { /* linked list of cubes acting as stack */ CUBE cube; /* a single cube */ struct cubes *next; /* remaining elements */ } CUBES; typedef struct centerlist { /* list of cube locations */ int i, j, k; /* cube location */ struct centerlist *next; /* remaining elements */ } CENTERLIST; typedef struct cornerlist { /* list of corners */ int i, j, k; /* corner id */ double value; /* corner value */ struct cornerlist *next; /* remaining elements */ } CORNERLIST; typedef struct edgelist { /* list of edges */ int i1, j1, k1, i2, j2, k2; /* edge corner ids */ int vid; /* vertex id */ struct edgelist *next; /* remaining elements */ } EDGELIST; typedef struct intlist { /* list of integers */ int i; /* an integer */ struct intlist *next; /* remaining elements */ } INTLIST; typedef struct intlists { /* list of list of integers */ INTLIST *list; /* a list of integers */ struct intlists *next; /* remaining elements */ } INTLISTS; typedef struct process { /* parameters, function, storage */ double (*function)(); /* implicit surface function */ int (*triproc)(); /* triangle output function */ double size, delta; /* cube size, normal delta */ int bounds; /* cube range within lattice */ POINT start; /* start point on surface */ CUBES *cubes; /* active cubes */ VERTICES vertices; /* surface vertices */ CENTERLIST **centers; /* cube center hash table */ CORNERLIST **corners; /* corner value hash table */ EDGELIST **edges; /* edge and vertex id hash table */ } PROCESS; void *calloc(); #define mycalloc(n,nbyte) _mycalloc(n,nbyte,__LINE__) #define myfree(ptr) _myfree(ptr,__LINE__) static void makecubetable (); static void free_cubetable(); static void converge(POINT*,POINT*,double,double(*f)(),POINT*); static CORNER *setcorner (PROCESS*, int, int, int); static int setcenter(CENTERLIST *table[], int, int, int); static int dotet (CUBE*, int, int, int, int, PROCESS*); static int docube(CUBE*,PROCESS*); static void testface (int,int,int,CUBE*,int,int,int,int,int,PROCESS*); static TEST find (int,PROCESS*,double,double,double); static void vnormal (POINT*,PROCESS*,POINT*); static void addtovertices (VERTICES*, VERTEX); static int vertid (CORNER*,CORNER*,PROCESS*); static void free_process_data(PROCESS *); static void clean_malloc(); static char *_mycalloc (int nitems, int nbytes, int line); static void _myfree(void*ptr, int lineno); #ifdef MAIN /**** A Test Program ****/ /* ffunction: a piece of an atomic f function */ double ffunction (x, y, z) double x, y, z; { return x*y*z; } /* torus: a torus with major, minor radii = 0.5, 0.1, try size = .05 */ double torus (x, y, z) double x, y, z; { double x2 = x*x, y2 = y*y, z2 = z*z; double a = x2+y2+z2+(0.5*0.5)-(0.1*0.1); return a*a-4.0*(0.5*0.5)*(y2+z2); } /* sphere: an inverse square function (always positive) */ double sphere (x, y, z) double x, y, z; { double rsq = x*x+y*y+z*z; return 1.0/(rsq < 0.00001? 0.00001 : rsq); } /* blob: a three-pole blend function, try size = .1 */ double blob (x, y, z) double x, y, z; { return 4.0-sphere(x+1.0,y,z)-sphere(x,y+1.0,z)-sphere(x,y,z+1.0); } #ifdef SGIGFX /**************************************************************/ #include /* triangle: called by polygonize() for each triangle; set SGI lines */ triangle (i1, i2, i3, vertices) int i1, i2, i3; VERTICES vertices; { float v[3]; int i, ids[3]; ids[0] = i1; ids[1] = i2; ids[2] = i3; bgnclosedline(); for (i = 0; i < 3; i++) { POINT *p = &vertices.ptr[ids[i]].position; v[0] = p->x; v[1] = p->y; v[2] = p->z; v3f(v); } endclosedline(); return 1; } /* main: call polygonize() with torus function * display lines on SGI */ main () { char *err, *polygonize(); keepaspect(1, 1); winopen("implicit"); doublebuffer(); gconfig(); perspective(450, 1.0/1.0, 0.1, 10.0); color(7); clear(); swapbuffers(); makeobj(1); if ((err = polygonize(torus, .1, 20, 0.,0.,0., triangle, TET)) != NULL) { fprintf(stderr, "%s\n", err); exit(1); } closeobj(); translate(0.0, 0.0, -2.0); pushmatrix(); while(1) { /* spin the object */ reshapeviewport(); color(7); clear(); color(0); callobj(1); rot(0.8, 'x'); rot(0.3, 'y'); rot(0.1, 'z'); swapbuffers(); } } #else /***********************************************************************/ int gntris; /* global needed by application */ VERTICES gvertices; /* global needed by application */ /* triangle: called by polygonize() for each triangle; write to stdout */ triangle (i1, i2, i3, vertices) int i1, i2, i3; VERTICES vertices; { gvertices = vertices; gntris++; fprintf(stdout, "%d %d %d\n", i1, i2, i3); return 1; } /* main: call polygonize() with torus function * write points-polygon formatted data to stdout */ main () { int i; char *err, *polygonize(); gntris = 0; fprintf(stdout, "triangles\n\n"); if ((err = polygonize(torus, .05, 20, 0.,0.,0., triangle, TET)) != NULL) { fprintf(stdout, "%s\n", err); exit(1); } fprintf(stdout, "\n%d triangles, %d vertices\n", gntris, gvertices.count); fprintf(stdout, "\nvertices\n\n"); for (i = 0; i < gvertices.count; i++) { VERTEX v; v = gvertices.ptr[i]; fprintf(stdout, "%f %f %f\t%f %f %f\n", v.position.x, v.position.y, v.position.z, v.normal.x, v.normal.y, v.normal.z); } fprintf(stderr, "%d triangles, %d vertices\n", gntris, gvertices.count); exit(0); } #endif /**********************************************************************/ #endif /* MAIN */ /**** An Implicit Surface Polygonizer ****/ /* polygonize: polygonize the implicit surface function * arguments are: * double function (x, y, z) * double x, y, z (an arbitrary 3D point) * the implicit surface function * return negative for inside, positive for outside * double size * width of the partitioning cube * int bounds * max. range of cubes (+/- on the three axes) from first cube * double x, y, z * coordinates of a starting point on or near the surface * may be defaulted to 0., 0., 0. * int triproc (i1, i2, i3, vertices) * int i1, i2, i3 (indices into the vertex array) * VERTICES vertices (the vertex array, indexed from 0) * called for each triangle * the triangle coordinates are (for i = i1, i2, i3): * vertices.ptr[i].position.x, .y, and .z * vertices are ccw when viewed from the out (positive) side * in a left-handed coordinate system * vertex normals point outwards * return 1 to continue, 0 to abort * int mode * TET: decompose cube and polygonize six tetrahedra * NOTET: polygonize cube directly * returns error or NULL */ char *polygonize (function, size, bounds, x, y, z, triproc, mode) double (*function)(), size, x, y, z; int bounds, (*triproc)(), mode; { PROCESS p; int n, noabort; CORNER *setcorner(); TEST in, out, find(); p.function = function; p.triproc = triproc; p.size = size; p.bounds = bounds; p.delta = size/(double)(RES*RES); /* allocate hash tables and build cube polygon table: */ p.centers = (CENTERLIST **) mycalloc(HASHSIZE,sizeof(CENTERLIST *)); p.corners = (CORNERLIST **) mycalloc(HASHSIZE,sizeof(CORNERLIST *)); p.edges = (EDGELIST **) mycalloc(2*HASHSIZE,sizeof(EDGELIST *)); makecubetable(); /* find point on surface, beginning search at (x, y, z): */ srand(1); in = find(1, &p, x, y, z); out = find(0, &p, x, y, z); if (!in.ok || !out.ok) { free_cubetable(); free_process_data(&p); clean_malloc(); return "can't find starting point"; } converge(&in.p, &out.p, in.value, p.function, &p.start); /* push initial cube on stack: */ p.cubes = (CUBES *) mycalloc(1, sizeof(CUBES)); /* list of 1 */ p.cubes->cube.i = p.cubes->cube.j = p.cubes->cube.k = 0; p.cubes->next = NULL; /* set corners of initial cube: */ for (n = 0; n < 8; n++) p.cubes->cube.corners[n] = setcorner(&p, BIT(n,2), BIT(n,1), BIT(n,0)); p.vertices.count = p.vertices.max = 0; /* no vertices yet */ p.vertices.ptr = NULL; setcenter(p.centers, 0, 0, 0); while (p.cubes != NULL) { /* process active cubes till none left */ int i; CUBE c; CUBES *temp = p.cubes; c = p.cubes->cube; noabort = mode == TET? /* either decompose into tetrahedra and polygonize: */ dotet(&c, LBN, LTN, RBN, LBF, &p) && dotet(&c, RTN, LTN, LBF, RBN, &p) && dotet(&c, RTN, LTN, LTF, LBF, &p) && dotet(&c, RTN, RBN, LBF, RBF, &p) && dotet(&c, RTN, LBF, LTF, RBF, &p) && dotet(&c, RTN, LTF, RTF, RBF, &p) : /* or polygonize the cube directly: */ docube(&c, &p); if (! noabort) { free_cubetable(); free_process_data(&p); clean_malloc(); return "aborted"; } /* pop current cube from stack */ p.cubes = p.cubes->next; /* test six face directions, maybe add to stack: */ testface(c.i-1, c.j, c.k, &c, L, LBN, LBF, LTN, LTF, &p); testface(c.i+1, c.j, c.k, &c, R, RBN, RBF, RTN, RTF, &p); testface(c.i, c.j-1, c.k, &c, B, LBN, LBF, RBN, RBF, &p); testface(c.i, c.j+1, c.k, &c, T, LTN, LTF, RTN, RTF, &p); testface(c.i, c.j, c.k-1, &c, N, LBN, LTN, RBN, RTN, &p); testface(c.i, c.j, c.k+1, &c, F, LBF, LTF, RBF, RTF, &p); /* get rid of the current cube */ for (i=0; i<8; i++) { myfree(temp->cube.corners[i]); temp->cube.corners[i]=0; } myfree(temp); } free_cubetable(); free_process_data(&p); clean_malloc(); return NULL; } static void free_process_data(p) PROCESS *p; { int i; CUBES *cubes,*nextcubes; if (p->vertices.ptr) myfree(p->vertices.ptr); for (i=0; icenters[i]; l; l=next) { next = l->next; myfree(l); } } for (i=0; icorners[i]; l; l=next) { next = l->next; myfree(l); } } for (i=0; i<2*HASHSIZE; i++) { EDGELIST *l,*next; for (l=p->edges[i]; l; l=next) { next = l->next; myfree(l); } } for (cubes=p->cubes; cubes; cubes=nextcubes) { nextcubes = cubes->next; for (i=0; i<8; i++) { myfree(cubes->cube.corners[i]); } myfree(cubes); } myfree(p->centers); myfree(p->corners); myfree(p->edges); } /* testface: given cube at lattice (i, j, k), and four corners of face, * if surface crosses face, compute other four corners of adjacent cube * and add new cube to cube stack */ static void testface (i, j, k, old, face, c1, c2, c3, c4, p) CUBE *old; PROCESS *p; int i, j, k, face, c1, c2, c3, c4; { CUBE new; CUBES *oldcubes = p->cubes; CORNER *setcorner(); int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0; /* static int facebit[6] = {2, 2, 1, 1, 0, 0}; */ /* int bit = facebit[face]; */ /* test if no surface crossing, cube out of bounds, or already visited: */ if ((old->corners[c2]->value > 0) == pos && (old->corners[c3]->value > 0) == pos && (old->corners[c4]->value > 0) == pos) return; if (abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) { static int have_been_warned = 0; if (!have_been_warned) { fprintf(stderr,"WARNING: testface: cube out of bounds\n"); have_been_warned = 1; } /* abort(); */ return; } if (setcenter(p->centers, i, j, k)) return; /* create new cube: */ new.i = i; new.j = j; new.k = k; /* CLJ: changed this to make memory management possible. */ /* for (n = 0; n < 8; n++) new.corners[n] = NULL; */ /* new.corners[FLIP(c1, bit)] = old->corners[c1]; */ /* new.corners[FLIP(c2, bit)] = old->corners[c2]; */ /* new.corners[FLIP(c3, bit)] = old->corners[c3]; */ /* new.corners[FLIP(c4, bit)] = old->corners[c4]; */ /* for (n = 0; n < 8; n++) */ /* if (new.corners[n] == NULL) */ /* new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0)); */ for (n = 0; n < 8; n++) new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0)); /*add cube to top of stack: */ p->cubes = (CUBES *) mycalloc(1, sizeof(CUBES)); p->cubes->cube = new; p->cubes->next = oldcubes; } /* setcorner: return corner with the given lattice location set (and cache) its function value */ static CORNER *setcorner (p, i, j, k) int i, j, k; PROCESS *p; { /* for speed, do corner value caching here */ CORNER *c = (CORNER *) mycalloc(1, sizeof(CORNER)); int index = HASH(i, j, k); CORNERLIST *l = p->corners[index]; c->i = i; c->x = p->start.x+((double)i-.5)*p->size; c->j = j; c->y = p->start.y+((double)j-.5)*p->size; c->k = k; c->z = p->start.z+((double)k-.5)*p->size; for (; l != NULL; l = l->next) if (l->i == i && l->j == j && l->k == k) { c->value = l->value; return c; } l = (CORNERLIST *) mycalloc(1, sizeof(CORNERLIST)); l->i = i; l->j = j; l->k = k; l->value = c->value = p->function(c->x, c->y, c->z); if (c->value > 100.0 || c->value < -100.0) { fprintf(stderr,"suspicious\n"); abort(); } l->next = p->corners[index]; p->corners[index] = l; return c; } /* find: search for point with value of given sign (0: neg, 1: pos) */ static TEST find (sign, p, x, y, z) int sign; PROCESS *p; double x, y, z; { int i; TEST test; double range = p->size; test.ok = 1; for (i = 0; i < 10000; i++) { test.p.x = x+range*(RAND()-0.5); test.p.y = y+range*(RAND()-0.5); test.p.z = z+range*(RAND()-0.5); test.value = p->function(test.p.x, test.p.y, test.p.z); if (sign == (test.value > 0.0)) return test; range = range*1.0005; /* slowly expand search outwards */ } test.ok = 0; return test; } /**** Tetrahedral Polygonization ****/ /* dotet: triangulate the tetrahedron * b, c, d should appear clockwise when viewed from a * return 0 if client aborts, 1 otherwise */ static int dotet (cube, c1, c2, c3, c4, p) CUBE *cube; int c1, c2, c3, c4; PROCESS *p; { CORNER *a = cube->corners[c1]; CORNER *b = cube->corners[c2]; CORNER *c = cube->corners[c3]; CORNER *d = cube->corners[c4]; int index = 0, apos, bpos, cpos, dpos, e1=0, e2=0, e3=0, e4=0, e5=0, e6=0; if ((apos = (a->value > 0.0))) index += 8; if ((bpos = (b->value > 0.0))) index += 4; if ((cpos = (c->value > 0.0))) index += 2; if ((dpos = (d->value > 0.0))) index += 1; /* index is now 4-bit number representing one of the 16 possible cases */ if (apos != bpos) e1 = vertid(a, b, p); if (apos != cpos) e2 = vertid(a, c, p); if (apos != dpos) e3 = vertid(a, d, p); if (bpos != cpos) e4 = vertid(b, c, p); if (bpos != dpos) e5 = vertid(b, d, p); if (cpos != dpos) e6 = vertid(c, d, p); /* 14 productive tetrahedral cases (0000 and 1111 do not yield polygons */ switch (index) { case 1: return p->triproc(e5, e6, e3, p->vertices); case 2: return p->triproc(e2, e6, e4, p->vertices); case 3: return p->triproc(e3, e5, e4, p->vertices) && p->triproc(e3, e4, e2, p->vertices); case 4: return p->triproc(e1, e4, e5, p->vertices); case 5: return p->triproc(e3, e1, e4, p->vertices) && p->triproc(e3, e4, e6, p->vertices); case 6: return p->triproc(e1, e2, e6, p->vertices) && p->triproc(e1, e6, e5, p->vertices); case 7: return p->triproc(e1, e2, e3, p->vertices); case 8: return p->triproc(e1, e3, e2, p->vertices); case 9: return p->triproc(e1, e5, e6, p->vertices) && p->triproc(e1, e6, e2, p->vertices); case 10: return p->triproc(e1, e3, e6, p->vertices) && p->triproc(e1, e6, e4, p->vertices); case 11: return p->triproc(e1, e5, e4, p->vertices); case 12: return p->triproc(e3, e2, e4, p->vertices) && p->triproc(e3, e4, e5, p->vertices); case 13: return p->triproc(e6, e2, e4, p->vertices); case 14: return p->triproc(e5, e3, e6, p->vertices); } return 1; } /**** Cubical Polygonization (optional) ****/ #define LB 0 /* left bottom edge */ #define LT 1 /* left top edge */ #define LN 2 /* left near edge */ #define LF 3 /* left far edge */ #define RB 4 /* right bottom edge */ #define RT 5 /* right top edge */ #define RN 6 /* right near edge */ #define RF 7 /* right far edge */ #define BN 8 /* bottom near edge */ #define BF 9 /* bottom far edge */ #define TN 10 /* top near edge */ #define TF 11 /* top far edge */ static INTLISTS *cubetable[256]; /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */ static int corner1[12] = {LBN,LTN,LBN,LBF,RBN,RTN,RBN,RBF,LBN,LBF,LTN,LTF}; static int corner2[12] = {LBF,LTF,LTN,LTF,RBF,RTF,RTN,RTF,RBN,RBF,RTN,RTF}; static int leftface[12] = {B, L, L, F, R, T, N, R, N, B, T, F}; /* face on left when going corner1 to corner2 */ static int rightface[12] = {L, T, N, L, B, R, R, F, B, F, N, T}; /* face on right when going corner1 to corner2 */ /* docube: triangulate the cube directly, without decomposition */ static int docube (cube, p) CUBE *cube; PROCESS *p; { INTLISTS *polys; int i, index = 0; for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0) index += (1<next) { INTLIST *edges; int a = -1, b = -1, count = 0; for (edges = polys->list; edges; edges = edges->next) { CORNER *c1 = cube->corners[corner1[edges->i]]; CORNER *c2 = cube->corners[corner2[edges->i]]; int c = vertid(c1, c2, p); if (++count > 2 && ! p->triproc(a, b, c, p->vertices)) return 0; if (count < 3) a = b; b = c; } } return 1; } /* nextcwedge: return next clockwise edge from given edge around given face */ static int nextcwedge (edge, face) int edge, face; { switch (edge) { case LB: return (face == L)? LF : BN; case LT: return (face == L)? LN : TF; case LN: return (face == L)? LB : TN; case LF: return (face == L)? LT : BF; case RB: return (face == R)? RN : BF; case RT: return (face == R)? RF : TN; case RN: return (face == R)? RT : BN; case RF: return (face == R)? RB : TF; case BN: return (face == B)? RB : LN; case BF: return (face == B)? LB : RF; case TN: return (face == T)? LT : RN; case TF: return (face == T)? RT : LF; } return -1; } /* otherface: return face adjoining edge that is not the given face */ static int otherface (edge, face) int edge, face; { int other = leftface[edge]; return face == other? rightface[edge] : other; } /* makecubetable: create the 256 entry table for cubical polygonization */ static void makecubetable () { int i, e, c, done[12], pos[8]; memset(cubetable, 0, sizeof(cubetable)); for (i = 0; i < 256; i++) { for (e = 0; e < 12; e++) done[e] = 0; for (c = 0; c < 8; c++) pos[c] = BIT(i, c); for (e = 0; e < 12; e++) if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) { INTLIST *ints = 0; INTLISTS *lists = (INTLISTS *) mycalloc(1, sizeof(INTLISTS)); int start = e, edge = e; /* get face that is to right of edge from pos to neg corner: */ int face = pos[corner1[e]]? rightface[e] : leftface[e]; while (1) { edge = nextcwedge(edge, face); done[edge] = 1; if (pos[corner1[edge]] != pos[corner2[edge]]) { INTLIST *tmp = ints; ints = (INTLIST *) mycalloc(1, sizeof(INTLIST)); ints->i = edge; ints->next = tmp; /* add edge to head of list */ if (edge == start) break; face = otherface(edge, face); } } lists->list = ints; /* add ints to head of table entry */ lists->next = cubetable[i]; cubetable[i] = lists; } } } static void free_cubetable() { int i; for (i=0; i<256; i++) { INTLISTS *l,*nextl; for (l=cubetable[i]; l; l=nextl) { INTLIST *m, *nextm; for (m=l->list; m; m=nextm) { nextm = m->next; myfree(m); } nextl = l->next; myfree(l); } } } /**** Storage ****/ #undef CHECK_MALLOC #ifdef CHECK_MALLOC static char allocwarn[10000]; static char delwarn[10000]; #endif /* mycalloc: return successful calloc or exit program */ typedef struct mallocdata { int lineno; char* ptr; size_t size; struct mallocdata* next; } MALLOCDATA; #ifdef CHECK_MALLOC static MALLOCDATA *malloc_list; static void add_mallocdata(char* ptr, int lineno, size_t size) { MALLOCDATA * old = malloc_list; malloc_list = (MALLOCDATA*) malloc(sizeof(MALLOCDATA)); malloc_list->next = old; malloc_list->ptr = ptr; malloc_list->size = size; malloc_list->lineno = lineno; } static size_t del_mallocdata(char* ptr,int lineno) { MALLOCDATA *i, *ilast = 0; int size; for (i=malloc_list; i; ilast=i,i=i->next) { if (i->ptr == ptr) { if (ilast) { MALLOCDATA * tmp = i->next; ilast->next = i->next; } else { malloc_list = i->next; } size = i->size; free(i); return size; } } if (!delwarn[lineno]) { fprintf(stderr,"tried to delete unknown data at line %d\n",lineno); delwarn[lineno] = 1; } return 0; } #endif static void clean_malloc() { #ifdef CHECK_MALLOC MALLOCDATA*i; int count=0; for (i=malloc_list; i; i=i->next) { if (!allocwarn[i->lineno]) { fprintf(stderr,"have memory allocated from line %d\n",i->lineno); allocwarn[i->lineno] = 1; } count++; } fprintf(stderr,"%d allocated pieces of memory remain\n",count); #endif } static char *_mycalloc (nitems, nbytes, line) int nitems, nbytes, line; { char *ptr = calloc(nitems, nbytes); #ifdef CHECK_MALLOC add_mallocdata(ptr,line,nitems*nbytes); #endif if (ptr != NULL) return ptr; fprintf(stderr, "can't calloc %d bytes\n", nitems*nbytes); abort(); return 0; } static void _myfree(ptr, lineno) void* ptr; int lineno; { #ifdef CHECK_MALLOC size_t size = del_mallocdata(ptr,lineno); char*tmp = ptr; for (int i=0; inext) if (l->i == i && l->j == j && l->k == k) return 1; new = (CENTERLIST *) mycalloc(1, sizeof(CENTERLIST)); new->i = i; new->j = j; new->k = k; new->next = q; table[index] = new; return 0; } /* setedge: set vertex id for edge */ static void setedge (table, i1, j1, k1, i2, j2, k2, vid) EDGELIST *table[]; int i1, j1, k1, i2, j2, k2, vid; { unsigned int index; EDGELIST *new; if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) { int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t; } index = HASH(i1, j1, k1) + HASH(i2, j2, k2); new = (EDGELIST *) mycalloc(1, sizeof(EDGELIST)); new->i1 = i1; new->j1 = j1; new->k1 = k1; new->i2 = i2; new->j2 = j2; new->k2 = k2; new->vid = vid; new->next = table[index]; table[index] = new; } /* getedge: return vertex id for edge; return -1 if not set */ static int getedge (table, i1, j1, k1, i2, j2, k2) EDGELIST *table[]; int i1, j1, k1, i2, j2, k2; { EDGELIST *q; if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) { int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t; }; q = table[HASH(i1, j1, k1)+HASH(i2, j2, k2)]; for (; q != NULL; q = q->next) if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 && q->i2 == i2 && q->j2 == j2 && q->k2 == k2) return q->vid; return -1; } /**** Vertices ****/ /* vertid: return index for vertex on edge: * c1->value and c2->value are presumed of different sign * return saved index if any; else compute vertex and save */ static int vertid (c1, c2, p) CORNER *c1, *c2; PROCESS *p; { VERTEX v; POINT a, b; int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k); if (vid != -1) return vid; /* previously computed */ a.x = c1->x; a.y = c1->y; a.z = c1->z; b.x = c2->x; b.y = c2->y; b.z = c2->z; converge(&a, &b, c1->value, p->function, &v.position); /* position */ vnormal(&v.position, p, &v.normal); /* normal */ addtovertices(&p->vertices, v); /* save vertex */ vid = p->vertices.count-1; setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid); return vid; } /* addtovertices: add v to sequence of vertices */ static void addtovertices (vertices, v) VERTICES *vertices; VERTEX v; { if (vertices->count == vertices->max) { int i; VERTEX *new; vertices->max = vertices->count == 0 ? 10 : 2*vertices->count; new = (VERTEX *) mycalloc(vertices->max, sizeof(VERTEX)); for (i = 0; i < vertices->count; i++) new[i] = vertices->ptr[i]; if (vertices->ptr != NULL) myfree(vertices->ptr); vertices->ptr = new; } vertices->ptr[vertices->count++] = v; } /* vnormal: compute unit length surface normal at point */ static void vnormal (point, p, v) POINT *point, *v; PROCESS *p; { double f = p->function(point->x, point->y, point->z); v->x = p->function(point->x+p->delta, point->y, point->z)-f; v->y = p->function(point->x, point->y+p->delta, point->z)-f; v->z = p->function(point->x, point->y, point->z+p->delta)-f; f = sqrt(v->x*v->x + v->y*v->y + v->z*v->z); if (f != 0.0) {v->x /= f; v->y /= f; v->z /= f;} } /* converge: from two points of differing sign, converge to zero crossing */ static void converge (p1, p2, v, function, p) double v; double (*function)(); POINT *p1, *p2, *p; { int i = 0; POINT pos, neg; if (v < 0) { pos.x = p2->x; pos.y = p2->y; pos.z = p2->z; neg.x = p1->x; neg.y = p1->y; neg.z = p1->z; } else { pos.x = p1->x; pos.y = p1->y; pos.z = p1->z; neg.x = p2->x; neg.y = p2->y; neg.z = p2->z; } while (1) { p->x = 0.5*(pos.x + neg.x); p->y = 0.5*(pos.y + neg.y); p->z = 0.5*(pos.z + neg.z); if (i++ == RES) return; if ((function(p->x, p->y, p->z)) > 0.0) {pos.x = p->x; pos.y = p->y; pos.z = p->z;} else {neg.x = p->x; neg.y = p->y; neg.z = p->z;} } }