implicit.c

[Game.Programming].Academic - Graphics Gems (6 books source code)

    /* create new cube: */
    new.i = i;
    new.j = j;
    new.k = k;
    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));

    /*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 */

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);
    l->next = p->corners[index];
    p->corners[index] = l;
    return c;
}


/* find: search for point with value of given sign (0: neg, 1: pos) */

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 */

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, e2, e3, e4, e5, e6;
    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 */

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<<i);
    for (polys = cubetable[index]; polys; polys = polys->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 */

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;
    }
}


/* otherface: return face adjoining edge that is not the given face */

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 */

makecubetable ()
{
    int i, e, c, done[12], pos[8];
    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;
	    }
    }
}


/**** Storage ****/


/* mycalloc: return successful calloc or exit program */

char *mycalloc (nitems, nbytes)
int nitems, nbytes;
{
   char *ptr = calloc(nitems, nbytes);
   if (ptr != NULL) return ptr;
   fprintf(stderr, "can't calloc %d bytes\n", nitems*nbytes);
   exit(1);
}


/* setcenter: set (i,j,k) entry of table[]
 * return 1 if already set; otherwise, set and return 0 */

int setcenter(table, i, j, k)
CENTERLIST *table[];
int i, j, k;
{
    int index = HASH(i, j, k);
    CENTERLIST *new, *l, *q = table[index];
    for (l = q; l != NULL; l = l->next)
	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 */

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 */

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 */

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 */

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) free((char *) vertices->ptr);
	vertices->ptr = new;
    }
    vertices->ptr[vertices->count++] = v;
}


/* vnormal: compute unit length surface normal at point */

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 */

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;}
    }
}