source: ThirdParty/mpqc_open/src/lib/math/isosurf/implicit.c@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 29.3 KB
RevLine 
[0b990d]1/*
2 * C code from the article
3 * "An Implicit Surface Polygonizer"
4 * by Jules Bloomenthal, jbloom@beauty.gmu.edu
5 * in "Graphics Gems IV", Academic Press, 1994
6 */
7
8/* Modified by Curtis Janssen:
9 * 1. Eliminate memory leaks.
10 * 2. Make main routine optional (-DMAIN to compile a main routine).
11 */
12
13/* implicit.c
14 * an implicit surface polygonizer, translated from Mesa
15 * applications should call polygonize()
16 *
17 * to compile a test program for ASCII output:
18 * cc -DMAIN implicit.c -o implicit -lm
19 *
20 * to compile a test program for display on an SGI workstation:
21 * cc -DMAIN -DSGIGFX implicit.c -o implicit -lgl_s -lm
22 *
23 * Authored by Jules Bloomenthal, Xerox PARC.
24 * Copyright (c) Xerox Corporation, 1991. All rights reserved.
25 * Permission is granted to reproduce, use and distribute this code for
26 * any and all purposes, provided that this notice appears in all copies. */
27
28#include <stdlib.h>
29#include <string.h>
30#include <math.h>
31#include <stdio.h>
32#include <sys/types.h>
33
34#define TET 0 /* use tetrahedral decomposition */
35#define NOTET 1 /* no tetrahedral decomposition */
36
37#define RES 10 /* # converge iterations */
38
39#define L 0 /* left direction: -x, -i */
40#define R 1 /* right direction: +x, +i */
41#define B 2 /* bottom direction: -y, -j */
42#define T 3 /* top direction: +y, +j */
43#define N 4 /* near direction: -z, -k */
44#define F 5 /* far direction: +z, +k */
45#define LBN 0 /* left bottom near corner */
46#define LBF 1 /* left bottom far corner */
47#define LTN 2 /* left top near corner */
48#define LTF 3 /* left top far corner */
49#define RBN 4 /* right bottom near corner */
50#define RBF 5 /* right bottom far corner */
51#define RTN 6 /* right top near corner */
52#define RTF 7 /* right top far corner */
53
54/* the LBN corner of cube (i, j, k), corresponds with location
55 * (start.x+(i-.5)*size, start.y+(j-.5)*size, start.z+(k-.5)*size) */
56
57#define RAND() ((rand()&32767)/32767.) /* random number between 0 and 1 */
58#define HASHBIT (5)
59#define HASHSIZE (size_t)(1<<(3*HASHBIT)) /* hash table size (32768) */
60#define MASK ((1<<HASHBIT)-1)
61#define HASH(i,j,k) ((((((i)&MASK)<<HASHBIT)|((j)&MASK))<<HASHBIT)|((k)&MASK))
62#define BIT(i, bit) (((i)>>(bit))&1)
63#define FLIP(i,bit) ((i)^1<<(bit)) /* flip the given bit of i */
64
65typedef struct point { /* a three-dimensional point */
66 double x, y, z; /* its coordinates */
67} POINT;
68
69typedef struct test { /* test the function for a signed value */
70 POINT p; /* location of test */
71 double value; /* function value at p */
72 int ok; /* if value is of correct sign */
73} TEST;
74
75typedef struct vertex { /* surface vertex */
76 POINT position, normal; /* position and surface normal */
77} VERTEX;
78
79typedef struct vertices { /* list of vertices in polygonization */
80 int count, max; /* # vertices, max # allowed */
81 VERTEX *ptr; /* dynamically allocated */
82} VERTICES;
83
84typedef struct corner { /* corner of a cube */
85 int i, j, k; /* (i, j, k) is index within lattice */
86 double x, y, z, value; /* location and function value */
87} CORNER;
88
89typedef struct cube { /* partitioning cell (cube) */
90 int i, j, k; /* lattice location of cube */
91 CORNER *corners[8]; /* eight corners */
92} CUBE;
93
94typedef struct cubes { /* linked list of cubes acting as stack */
95 CUBE cube; /* a single cube */
96 struct cubes *next; /* remaining elements */
97} CUBES;
98
99typedef struct centerlist { /* list of cube locations */
100 int i, j, k; /* cube location */
101 struct centerlist *next; /* remaining elements */
102} CENTERLIST;
103
104typedef struct cornerlist { /* list of corners */
105 int i, j, k; /* corner id */
106 double value; /* corner value */
107 struct cornerlist *next; /* remaining elements */
108} CORNERLIST;
109
110typedef struct edgelist { /* list of edges */
111 int i1, j1, k1, i2, j2, k2; /* edge corner ids */
112 int vid; /* vertex id */
113 struct edgelist *next; /* remaining elements */
114} EDGELIST;
115
116typedef struct intlist { /* list of integers */
117 int i; /* an integer */
118 struct intlist *next; /* remaining elements */
119} INTLIST;
120
121typedef struct intlists { /* list of list of integers */
122 INTLIST *list; /* a list of integers */
123 struct intlists *next; /* remaining elements */
124} INTLISTS;
125
126typedef struct process { /* parameters, function, storage */
127 double (*function)(); /* implicit surface function */
128 int (*triproc)(); /* triangle output function */
129 double size, delta; /* cube size, normal delta */
130 int bounds; /* cube range within lattice */
131 POINT start; /* start point on surface */
132 CUBES *cubes; /* active cubes */
133 VERTICES vertices; /* surface vertices */
134 CENTERLIST **centers; /* cube center hash table */
135 CORNERLIST **corners; /* corner value hash table */
136 EDGELIST **edges; /* edge and vertex id hash table */
137} PROCESS;
138
139void *calloc();
140
141#define mycalloc(n,nbyte) _mycalloc(n,nbyte,__LINE__)
142#define myfree(ptr) _myfree(ptr,__LINE__)
143
144static void makecubetable ();
145static void free_cubetable();
146static void converge(POINT*,POINT*,double,double(*f)(),POINT*);
147static CORNER *setcorner (PROCESS*, int, int, int);
148static int setcenter(CENTERLIST *table[], int, int, int);
149static int dotet (CUBE*, int, int, int, int, PROCESS*);
150static int docube(CUBE*,PROCESS*);
151static void testface (int,int,int,CUBE*,int,int,int,int,int,PROCESS*);
152static TEST find (int,PROCESS*,double,double,double);
153static void vnormal (POINT*,PROCESS*,POINT*);
154static void addtovertices (VERTICES*, VERTEX);
155static int vertid (CORNER*,CORNER*,PROCESS*);
156static void free_process_data(PROCESS *);
157static void clean_malloc();
158static char *_mycalloc (int nitems, int nbytes, int line);
159static void _myfree(void*ptr, int lineno);
160
161#ifdef MAIN
162
163/**** A Test Program ****/
164
165
166/* ffunction: a piece of an atomic f function */
167
168double ffunction (x, y, z)
169double x, y, z;
170{
171 return x*y*z;
172}
173
174/* torus: a torus with major, minor radii = 0.5, 0.1, try size = .05 */
175
176double torus (x, y, z)
177double x, y, z;
178{
179 double x2 = x*x, y2 = y*y, z2 = z*z;
180 double a = x2+y2+z2+(0.5*0.5)-(0.1*0.1);
181 return a*a-4.0*(0.5*0.5)*(y2+z2);
182}
183
184
185/* sphere: an inverse square function (always positive) */
186
187double sphere (x, y, z)
188double x, y, z;
189{
190 double rsq = x*x+y*y+z*z;
191 return 1.0/(rsq < 0.00001? 0.00001 : rsq);
192}
193
194
195/* blob: a three-pole blend function, try size = .1 */
196
197double blob (x, y, z)
198double x, y, z;
199{
200 return 4.0-sphere(x+1.0,y,z)-sphere(x,y+1.0,z)-sphere(x,y,z+1.0);
201}
202
203#ifdef SGIGFX /**************************************************************/
204
205#include <math/isosurf/gl.h>
206
207/* triangle: called by polygonize() for each triangle; set SGI lines */
208
209triangle (i1, i2, i3, vertices)
210int i1, i2, i3;
211VERTICES vertices;
212{
213 float v[3];
214 int i, ids[3];
215 ids[0] = i1;
216 ids[1] = i2;
217 ids[2] = i3;
218 bgnclosedline();
219 for (i = 0; i < 3; i++) {
220 POINT *p = &vertices.ptr[ids[i]].position;
221 v[0] = p->x; v[1] = p->y; v[2] = p->z;
222 v3f(v);
223 }
224 endclosedline();
225 return 1;
226}
227
228
229/* main: call polygonize() with torus function
230 * display lines on SGI */
231
232main ()
233{
234 char *err, *polygonize();
235
236 keepaspect(1, 1);
237 winopen("implicit");
238 doublebuffer();
239 gconfig();
240 perspective(450, 1.0/1.0, 0.1, 10.0);
241 color(7);
242 clear();
243 swapbuffers();
244 makeobj(1);
245 if ((err = polygonize(torus, .1, 20, 0.,0.,0., triangle, TET)) != NULL) {
246 fprintf(stderr, "%s\n", err);
247 exit(1);
248 }
249 closeobj();
250 translate(0.0, 0.0, -2.0);
251 pushmatrix();
252 while(1) { /* spin the object */
253 reshapeviewport();
254 color(7);
255 clear();
256 color(0);
257 callobj(1);
258 rot(0.8, 'x');
259 rot(0.3, 'y');
260 rot(0.1, 'z');
261 swapbuffers();
262
263 }
264}
265
266#else /***********************************************************************/
267
268int gntris; /* global needed by application */
269VERTICES gvertices; /* global needed by application */
270
271
272/* triangle: called by polygonize() for each triangle; write to stdout */
273
274triangle (i1, i2, i3, vertices)
275int i1, i2, i3;
276VERTICES vertices;
277{
278 gvertices = vertices;
279 gntris++;
280 fprintf(stdout, "%d %d %d\n", i1, i2, i3);
281 return 1;
282}
283
284
285/* main: call polygonize() with torus function
286 * write points-polygon formatted data to stdout */
287
288main ()
289 {
290 int i;
291 char *err, *polygonize();
292 gntris = 0;
293 fprintf(stdout, "triangles\n\n");
294 if ((err = polygonize(torus, .05, 20, 0.,0.,0., triangle, TET)) != NULL) {
295 fprintf(stdout, "%s\n", err);
296 exit(1);
297 }
298 fprintf(stdout, "\n%d triangles, %d vertices\n", gntris, gvertices.count);
299 fprintf(stdout, "\nvertices\n\n");
300 for (i = 0; i < gvertices.count; i++) {
301 VERTEX v;
302 v = gvertices.ptr[i];
303 fprintf(stdout, "%f %f %f\t%f %f %f\n",
304 v.position.x, v.position.y, v.position.z,
305 v.normal.x, v.normal.y, v.normal.z);
306 }
307 fprintf(stderr, "%d triangles, %d vertices\n", gntris, gvertices.count);
308 exit(0);
309}
310
311#endif /**********************************************************************/
312
313#endif /* MAIN */
314
315
316/**** An Implicit Surface Polygonizer ****/
317
318
319/* polygonize: polygonize the implicit surface function
320 * arguments are:
321 * double function (x, y, z)
322 * double x, y, z (an arbitrary 3D point)
323 * the implicit surface function
324 * return negative for inside, positive for outside
325 * double size
326 * width of the partitioning cube
327 * int bounds
328 * max. range of cubes (+/- on the three axes) from first cube
329 * double x, y, z
330 * coordinates of a starting point on or near the surface
331 * may be defaulted to 0., 0., 0.
332 * int triproc (i1, i2, i3, vertices)
333 * int i1, i2, i3 (indices into the vertex array)
334 * VERTICES vertices (the vertex array, indexed from 0)
335 * called for each triangle
336 * the triangle coordinates are (for i = i1, i2, i3):
337 * vertices.ptr[i].position.x, .y, and .z
338 * vertices are ccw when viewed from the out (positive) side
339 * in a left-handed coordinate system
340 * vertex normals point outwards
341 * return 1 to continue, 0 to abort
342 * int mode
343 * TET: decompose cube and polygonize six tetrahedra
344 * NOTET: polygonize cube directly
345 * returns error or NULL
346 */
347
348char *polygonize (function, size, bounds, x, y, z, triproc, mode)
349double (*function)(), size, x, y, z;
350int bounds, (*triproc)(), mode;
351{
352 PROCESS p;
353 int n, noabort;
354 CORNER *setcorner();
355 TEST in, out, find();
356
357 p.function = function;
358 p.triproc = triproc;
359 p.size = size;
360 p.bounds = bounds;
361 p.delta = size/(double)(RES*RES);
362
363 /* allocate hash tables and build cube polygon table: */
364 p.centers = (CENTERLIST **) mycalloc(HASHSIZE,sizeof(CENTERLIST *));
365 p.corners = (CORNERLIST **) mycalloc(HASHSIZE,sizeof(CORNERLIST *));
366 p.edges = (EDGELIST **) mycalloc(2*HASHSIZE,sizeof(EDGELIST *));
367 makecubetable();
368
369 /* find point on surface, beginning search at (x, y, z): */
370 srand(1);
371 in = find(1, &p, x, y, z);
372 out = find(0, &p, x, y, z);
373 if (!in.ok || !out.ok) {
374 free_cubetable();
375 free_process_data(&p);
376 clean_malloc();
377 return "can't find starting point";
378 }
379 converge(&in.p, &out.p, in.value, p.function, &p.start);
380
381 /* push initial cube on stack: */
382 p.cubes = (CUBES *) mycalloc(1, sizeof(CUBES)); /* list of 1 */
383 p.cubes->cube.i = p.cubes->cube.j = p.cubes->cube.k = 0;
384 p.cubes->next = NULL;
385
386 /* set corners of initial cube: */
387 for (n = 0; n < 8; n++)
388 p.cubes->cube.corners[n] = setcorner(&p, BIT(n,2), BIT(n,1), BIT(n,0));
389
390 p.vertices.count = p.vertices.max = 0; /* no vertices yet */
391 p.vertices.ptr = NULL;
392
393 setcenter(p.centers, 0, 0, 0);
394
395 while (p.cubes != NULL) { /* process active cubes till none left */
396 int i;
397 CUBE c;
398 CUBES *temp = p.cubes;
399 c = p.cubes->cube;
400
401 noabort = mode == TET?
402 /* either decompose into tetrahedra and polygonize: */
403 dotet(&c, LBN, LTN, RBN, LBF, &p) &&
404 dotet(&c, RTN, LTN, LBF, RBN, &p) &&
405 dotet(&c, RTN, LTN, LTF, LBF, &p) &&
406 dotet(&c, RTN, RBN, LBF, RBF, &p) &&
407 dotet(&c, RTN, LBF, LTF, RBF, &p) &&
408 dotet(&c, RTN, LTF, RTF, RBF, &p)
409 :
410 /* or polygonize the cube directly: */
411 docube(&c, &p);
412 if (! noabort) {
413 free_cubetable();
414 free_process_data(&p);
415 clean_malloc();
416 return "aborted";
417 }
418
419 /* pop current cube from stack */
420 p.cubes = p.cubes->next;
421
422 /* test six face directions, maybe add to stack: */
423 testface(c.i-1, c.j, c.k, &c, L, LBN, LBF, LTN, LTF, &p);
424 testface(c.i+1, c.j, c.k, &c, R, RBN, RBF, RTN, RTF, &p);
425 testface(c.i, c.j-1, c.k, &c, B, LBN, LBF, RBN, RBF, &p);
426 testface(c.i, c.j+1, c.k, &c, T, LTN, LTF, RTN, RTF, &p);
427 testface(c.i, c.j, c.k-1, &c, N, LBN, LTN, RBN, RTN, &p);
428 testface(c.i, c.j, c.k+1, &c, F, LBF, LTF, RBF, RTF, &p);
429
430 /* get rid of the current cube */
431 for (i=0; i<8; i++) {
432 myfree(temp->cube.corners[i]);
433 temp->cube.corners[i]=0;
434 }
435 myfree(temp);
436 }
437 free_cubetable();
438 free_process_data(&p);
439 clean_malloc();
440 return NULL;
441}
442
443static void
444free_process_data(p)
445 PROCESS *p;
446{
447 int i;
448 CUBES *cubes,*nextcubes;
449
450 if (p->vertices.ptr) myfree(p->vertices.ptr);
451
452 for (i=0; i<HASHSIZE; i++) {
453 CENTERLIST *l,*next;
454 for (l=p->centers[i]; l; l=next) {
455 next = l->next;
456 myfree(l);
457 }
458 }
459
460 for (i=0; i<HASHSIZE; i++) {
461 CORNERLIST *l,*next;
462 for (l=p->corners[i]; l; l=next) {
463 next = l->next;
464 myfree(l);
465 }
466 }
467
468 for (i=0; i<2*HASHSIZE; i++) {
469 EDGELIST *l,*next;
470 for (l=p->edges[i]; l; l=next) {
471 next = l->next;
472 myfree(l);
473 }
474 }
475
476 for (cubes=p->cubes; cubes; cubes=nextcubes) {
477 nextcubes = cubes->next;
478 for (i=0; i<8; i++) {
479 myfree(cubes->cube.corners[i]);
480 }
481 myfree(cubes);
482 }
483
484 myfree(p->centers);
485 myfree(p->corners);
486 myfree(p->edges);
487}
488
489
490/* testface: given cube at lattice (i, j, k), and four corners of face,
491 * if surface crosses face, compute other four corners of adjacent cube
492 * and add new cube to cube stack */
493
494static void
495testface (i, j, k, old, face, c1, c2, c3, c4, p)
496CUBE *old;
497PROCESS *p;
498int i, j, k, face, c1, c2, c3, c4;
499{
500 CUBE new;
501 CUBES *oldcubes = p->cubes;
502 CORNER *setcorner();
503 int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0;
504 /* static int facebit[6] = {2, 2, 1, 1, 0, 0}; */
505 /* int bit = facebit[face]; */
506
507 /* test if no surface crossing, cube out of bounds, or already visited: */
508 if ((old->corners[c2]->value > 0) == pos &&
509 (old->corners[c3]->value > 0) == pos &&
510 (old->corners[c4]->value > 0) == pos) return;
511 if (abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) {
512 static int have_been_warned = 0;
513 if (!have_been_warned) {
514 fprintf(stderr,"WARNING: testface: cube out of bounds\n");
515 have_been_warned = 1;
516 }
517 /* abort(); */
518 return;
519 }
520 if (setcenter(p->centers, i, j, k)) return;
521
522 /* create new cube: */
523 new.i = i;
524 new.j = j;
525 new.k = k;
526 /* CLJ: changed this to make memory management possible. */
527/* for (n = 0; n < 8; n++) new.corners[n] = NULL; */
528/* new.corners[FLIP(c1, bit)] = old->corners[c1]; */
529/* new.corners[FLIP(c2, bit)] = old->corners[c2]; */
530/* new.corners[FLIP(c3, bit)] = old->corners[c3]; */
531/* new.corners[FLIP(c4, bit)] = old->corners[c4]; */
532/* for (n = 0; n < 8; n++) */
533/* if (new.corners[n] == NULL) */
534/* new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0)); */
535 for (n = 0; n < 8; n++)
536 new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0));
537
538 /*add cube to top of stack: */
539 p->cubes = (CUBES *) mycalloc(1, sizeof(CUBES));
540 p->cubes->cube = new;
541 p->cubes->next = oldcubes;
542}
543
544
545/* setcorner: return corner with the given lattice location
546 set (and cache) its function value */
547
548static CORNER *setcorner (p, i, j, k)
549int i, j, k;
550PROCESS *p;
551{
552 /* for speed, do corner value caching here */
553 CORNER *c = (CORNER *) mycalloc(1, sizeof(CORNER));
554 int index = HASH(i, j, k);
555 CORNERLIST *l = p->corners[index];
556 c->i = i; c->x = p->start.x+((double)i-.5)*p->size;
557 c->j = j; c->y = p->start.y+((double)j-.5)*p->size;
558 c->k = k; c->z = p->start.z+((double)k-.5)*p->size;
559 for (; l != NULL; l = l->next)
560 if (l->i == i && l->j == j && l->k == k) {
561 c->value = l->value;
562 return c;
563 }
564 l = (CORNERLIST *) mycalloc(1, sizeof(CORNERLIST));
565 l->i = i; l->j = j; l->k = k;
566 l->value = c->value = p->function(c->x, c->y, c->z);
567 if (c->value > 100.0 || c->value < -100.0) {
568 fprintf(stderr,"suspicious\n");
569 abort();
570 }
571 l->next = p->corners[index];
572 p->corners[index] = l;
573 return c;
574}
575
576
577/* find: search for point with value of given sign (0: neg, 1: pos) */
578
579static TEST find (sign, p, x, y, z)
580int sign;
581PROCESS *p;
582double x, y, z;
583{
584 int i;
585 TEST test;
586 double range = p->size;
587 test.ok = 1;
588 for (i = 0; i < 10000; i++) {
589 test.p.x = x+range*(RAND()-0.5);
590 test.p.y = y+range*(RAND()-0.5);
591 test.p.z = z+range*(RAND()-0.5);
592 test.value = p->function(test.p.x, test.p.y, test.p.z);
593 if (sign == (test.value > 0.0)) return test;
594 range = range*1.0005; /* slowly expand search outwards */
595 }
596 test.ok = 0;
597 return test;
598}
599
600
601/**** Tetrahedral Polygonization ****/
602
603
604/* dotet: triangulate the tetrahedron
605 * b, c, d should appear clockwise when viewed from a
606 * return 0 if client aborts, 1 otherwise */
607
608static int dotet (cube, c1, c2, c3, c4, p)
609CUBE *cube;
610int c1, c2, c3, c4;
611PROCESS *p;
612{
613 CORNER *a = cube->corners[c1];
614 CORNER *b = cube->corners[c2];
615 CORNER *c = cube->corners[c3];
616 CORNER *d = cube->corners[c4];
617 int index = 0, apos, bpos, cpos, dpos, e1=0, e2=0, e3=0, e4=0, e5=0, e6=0;
618 if ((apos = (a->value > 0.0))) index += 8;
619 if ((bpos = (b->value > 0.0))) index += 4;
620 if ((cpos = (c->value > 0.0))) index += 2;
621 if ((dpos = (d->value > 0.0))) index += 1;
622 /* index is now 4-bit number representing one of the 16 possible cases */
623 if (apos != bpos) e1 = vertid(a, b, p);
624 if (apos != cpos) e2 = vertid(a, c, p);
625 if (apos != dpos) e3 = vertid(a, d, p);
626 if (bpos != cpos) e4 = vertid(b, c, p);
627 if (bpos != dpos) e5 = vertid(b, d, p);
628 if (cpos != dpos) e6 = vertid(c, d, p);
629 /* 14 productive tetrahedral cases (0000 and 1111 do not yield polygons */
630 switch (index) {
631 case 1: return p->triproc(e5, e6, e3, p->vertices);
632 case 2: return p->triproc(e2, e6, e4, p->vertices);
633 case 3: return p->triproc(e3, e5, e4, p->vertices) &&
634 p->triproc(e3, e4, e2, p->vertices);
635 case 4: return p->triproc(e1, e4, e5, p->vertices);
636 case 5: return p->triproc(e3, e1, e4, p->vertices) &&
637 p->triproc(e3, e4, e6, p->vertices);
638 case 6: return p->triproc(e1, e2, e6, p->vertices) &&
639 p->triproc(e1, e6, e5, p->vertices);
640 case 7: return p->triproc(e1, e2, e3, p->vertices);
641 case 8: return p->triproc(e1, e3, e2, p->vertices);
642 case 9: return p->triproc(e1, e5, e6, p->vertices) &&
643 p->triproc(e1, e6, e2, p->vertices);
644 case 10: return p->triproc(e1, e3, e6, p->vertices) &&
645 p->triproc(e1, e6, e4, p->vertices);
646 case 11: return p->triproc(e1, e5, e4, p->vertices);
647 case 12: return p->triproc(e3, e2, e4, p->vertices) &&
648 p->triproc(e3, e4, e5, p->vertices);
649 case 13: return p->triproc(e6, e2, e4, p->vertices);
650 case 14: return p->triproc(e5, e3, e6, p->vertices);
651 }
652 return 1;
653}
654
655
656/**** Cubical Polygonization (optional) ****/
657
658
659#define LB 0 /* left bottom edge */
660#define LT 1 /* left top edge */
661#define LN 2 /* left near edge */
662#define LF 3 /* left far edge */
663#define RB 4 /* right bottom edge */
664#define RT 5 /* right top edge */
665#define RN 6 /* right near edge */
666#define RF 7 /* right far edge */
667#define BN 8 /* bottom near edge */
668#define BF 9 /* bottom far edge */
669#define TN 10 /* top near edge */
670#define TF 11 /* top far edge */
671
672static INTLISTS *cubetable[256];
673
674/* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
675static int corner1[12] = {LBN,LTN,LBN,LBF,RBN,RTN,RBN,RBF,LBN,LBF,LTN,LTF};
676static int corner2[12] = {LBF,LTF,LTN,LTF,RBF,RTF,RTN,RTF,RBN,RBF,RTN,RTF};
677static int leftface[12] = {B, L, L, F, R, T, N, R, N, B, T, F};
678 /* face on left when going corner1 to corner2 */
679static int rightface[12] = {L, T, N, L, B, R, R, F, B, F, N, T};
680 /* face on right when going corner1 to corner2 */
681
682
683/* docube: triangulate the cube directly, without decomposition */
684
685static int docube (cube, p)
686CUBE *cube;
687PROCESS *p;
688{
689 INTLISTS *polys;
690 int i, index = 0;
691 for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0) index += (1<<i);
692 for (polys = cubetable[index]; polys; polys = polys->next) {
693 INTLIST *edges;
694 int a = -1, b = -1, count = 0;
695 for (edges = polys->list; edges; edges = edges->next) {
696 CORNER *c1 = cube->corners[corner1[edges->i]];
697 CORNER *c2 = cube->corners[corner2[edges->i]];
698 int c = vertid(c1, c2, p);
699 if (++count > 2 && ! p->triproc(a, b, c, p->vertices)) return 0;
700 if (count < 3) a = b;
701 b = c;
702 }
703 }
704 return 1;
705}
706
707
708/* nextcwedge: return next clockwise edge from given edge around given face */
709
710static int nextcwedge (edge, face)
711int edge, face;
712{
713 switch (edge) {
714 case LB: return (face == L)? LF : BN;
715 case LT: return (face == L)? LN : TF;
716 case LN: return (face == L)? LB : TN;
717 case LF: return (face == L)? LT : BF;
718 case RB: return (face == R)? RN : BF;
719 case RT: return (face == R)? RF : TN;
720 case RN: return (face == R)? RT : BN;
721 case RF: return (face == R)? RB : TF;
722 case BN: return (face == B)? RB : LN;
723 case BF: return (face == B)? LB : RF;
724 case TN: return (face == T)? LT : RN;
725 case TF: return (face == T)? RT : LF;
726 }
727
728 return -1;
729}
730
731
732/* otherface: return face adjoining edge that is not the given face */
733
734static int otherface (edge, face)
735int edge, face;
736{
737 int other = leftface[edge];
738 return face == other? rightface[edge] : other;
739}
740
741
742/* makecubetable: create the 256 entry table for cubical polygonization */
743
744static void makecubetable ()
745{
746 int i, e, c, done[12], pos[8];
747 memset(cubetable, 0, sizeof(cubetable));
748 for (i = 0; i < 256; i++) {
749 for (e = 0; e < 12; e++) done[e] = 0;
750 for (c = 0; c < 8; c++) pos[c] = BIT(i, c);
751 for (e = 0; e < 12; e++)
752 if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
753 INTLIST *ints = 0;
754 INTLISTS *lists = (INTLISTS *) mycalloc(1, sizeof(INTLISTS));
755 int start = e, edge = e;
756 /* get face that is to right of edge from pos to neg corner: */
757 int face = pos[corner1[e]]? rightface[e] : leftface[e];
758 while (1) {
759 edge = nextcwedge(edge, face);
760 done[edge] = 1;
761 if (pos[corner1[edge]] != pos[corner2[edge]]) {
762 INTLIST *tmp = ints;
763 ints = (INTLIST *) mycalloc(1, sizeof(INTLIST));
764 ints->i = edge;
765 ints->next = tmp; /* add edge to head of list */
766 if (edge == start) break;
767 face = otherface(edge, face);
768 }
769 }
770 lists->list = ints; /* add ints to head of table entry */
771 lists->next = cubetable[i];
772 cubetable[i] = lists;
773 }
774 }
775}
776
777static void
778free_cubetable()
779{
780 int i;
781 for (i=0; i<256; i++) {
782 INTLISTS *l,*nextl;
783 for (l=cubetable[i]; l; l=nextl) {
784 INTLIST *m, *nextm;
785 for (m=l->list; m; m=nextm) {
786 nextm = m->next;
787 myfree(m);
788 }
789 nextl = l->next;
790 myfree(l);
791 }
792 }
793}
794
795/**** Storage ****/
796
797#undef CHECK_MALLOC
798
799#ifdef CHECK_MALLOC
800static char allocwarn[10000];
801static char delwarn[10000];
802#endif
803
804/* mycalloc: return successful calloc or exit program */
805
806typedef struct mallocdata {
807 int lineno;
808 char* ptr;
809 size_t size;
810 struct mallocdata* next;
811} MALLOCDATA;
812
813#ifdef CHECK_MALLOC
814static MALLOCDATA *malloc_list;
815static void add_mallocdata(char* ptr, int lineno, size_t size)
816{
817 MALLOCDATA * old = malloc_list;
818 malloc_list = (MALLOCDATA*) malloc(sizeof(MALLOCDATA));
819 malloc_list->next = old;
820 malloc_list->ptr = ptr;
821 malloc_list->size = size;
822 malloc_list->lineno = lineno;
823}
824
825static size_t del_mallocdata(char* ptr,int lineno)
826{
827 MALLOCDATA *i, *ilast = 0;
828 int size;
829 for (i=malloc_list; i; ilast=i,i=i->next) {
830 if (i->ptr == ptr) {
831 if (ilast) {
832 MALLOCDATA * tmp = i->next;
833 ilast->next = i->next;
834 }
835 else {
836 malloc_list = i->next;
837 }
838 size = i->size;
839 free(i);
840 return size;
841 }
842 }
843 if (!delwarn[lineno]) {
844 fprintf(stderr,"tried to delete unknown data at line %d\n",lineno);
845 delwarn[lineno] = 1;
846 }
847 return 0;
848}
849#endif
850
851static void clean_malloc()
852{
853#ifdef CHECK_MALLOC
854 MALLOCDATA*i;
855 int count=0;
856 for (i=malloc_list; i; i=i->next) {
857 if (!allocwarn[i->lineno]) {
858 fprintf(stderr,"have memory allocated from line %d\n",i->lineno);
859 allocwarn[i->lineno] = 1;
860 }
861 count++;
862 }
863 fprintf(stderr,"%d allocated pieces of memory remain\n",count);
864#endif
865}
866
867static char *_mycalloc (nitems, nbytes, line)
868int nitems, nbytes, line;
869{
870 char *ptr = calloc(nitems, nbytes);
871#ifdef CHECK_MALLOC
872 add_mallocdata(ptr,line,nitems*nbytes);
873#endif
874 if (ptr != NULL) return ptr;
875 fprintf(stderr, "can't calloc %d bytes\n", nitems*nbytes);
876 abort();
877 return 0;
878}
879
880static void _myfree(ptr, lineno)
881 void* ptr;
882 int lineno;
883{
884#ifdef CHECK_MALLOC
885 size_t size = del_mallocdata(ptr,lineno);
886 char*tmp = ptr;
887 for (int i=0; i<size; i++) {
888 *tmp++ = 0x00;
889 }
890#endif
891
892 free(ptr);
893}
894
895
896/* setcenter: set (i,j,k) entry of table[]
897 * return 1 if already set; otherwise, set and return 0 */
898
899static int setcenter(table, i, j, k)
900CENTERLIST *table[];
901int i, j, k;
902{
903 int index = HASH(i, j, k);
904 CENTERLIST *new, *l, *q = table[index];
905 for (l = q; l != NULL; l = l->next)
906 if (l->i == i && l->j == j && l->k == k) return 1;
907 new = (CENTERLIST *) mycalloc(1, sizeof(CENTERLIST));
908 new->i = i; new->j = j; new->k = k; new->next = q;
909 table[index] = new;
910 return 0;
911}
912
913
914/* setedge: set vertex id for edge */
915
916static void setedge (table, i1, j1, k1, i2, j2, k2, vid)
917EDGELIST *table[];
918int i1, j1, k1, i2, j2, k2, vid;
919{
920 unsigned int index;
921 EDGELIST *new;
922 if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
923 int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
924 }
925 index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
926 new = (EDGELIST *) mycalloc(1, sizeof(EDGELIST));
927 new->i1 = i1; new->j1 = j1; new->k1 = k1;
928 new->i2 = i2; new->j2 = j2; new->k2 = k2;
929 new->vid = vid;
930 new->next = table[index];
931 table[index] = new;
932}
933
934
935/* getedge: return vertex id for edge; return -1 if not set */
936
937static int getedge (table, i1, j1, k1, i2, j2, k2)
938EDGELIST *table[];
939int i1, j1, k1, i2, j2, k2;
940{
941 EDGELIST *q;
942 if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
943 int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
944 };
945 q = table[HASH(i1, j1, k1)+HASH(i2, j2, k2)];
946 for (; q != NULL; q = q->next)
947 if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
948 q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
949 return q->vid;
950 return -1;
951}
952
953
954/**** Vertices ****/
955
956
957/* vertid: return index for vertex on edge:
958 * c1->value and c2->value are presumed of different sign
959 * return saved index if any; else compute vertex and save */
960
961static int vertid (c1, c2, p)
962CORNER *c1, *c2;
963PROCESS *p;
964{
965 VERTEX v;
966 POINT a, b;
967 int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
968 if (vid != -1) return vid; /* previously computed */
969 a.x = c1->x; a.y = c1->y; a.z = c1->z;
970 b.x = c2->x; b.y = c2->y; b.z = c2->z;
971 converge(&a, &b, c1->value, p->function, &v.position); /* position */
972 vnormal(&v.position, p, &v.normal); /* normal */
973 addtovertices(&p->vertices, v); /* save vertex */
974 vid = p->vertices.count-1;
975 setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
976 return vid;
977}
978
979
980/* addtovertices: add v to sequence of vertices */
981
982static void addtovertices (vertices, v)
983VERTICES *vertices;
984VERTEX v;
985{
986 if (vertices->count == vertices->max) {
987 int i;
988 VERTEX *new;
989 vertices->max = vertices->count == 0 ? 10 : 2*vertices->count;
990 new = (VERTEX *) mycalloc(vertices->max, sizeof(VERTEX));
991 for (i = 0; i < vertices->count; i++) new[i] = vertices->ptr[i];
992 if (vertices->ptr != NULL) myfree(vertices->ptr);
993 vertices->ptr = new;
994 }
995 vertices->ptr[vertices->count++] = v;
996}
997
998
999/* vnormal: compute unit length surface normal at point */
1000
1001static void vnormal (point, p, v)
1002POINT *point, *v;
1003PROCESS *p;
1004{
1005 double f = p->function(point->x, point->y, point->z);
1006 v->x = p->function(point->x+p->delta, point->y, point->z)-f;
1007 v->y = p->function(point->x, point->y+p->delta, point->z)-f;
1008 v->z = p->function(point->x, point->y, point->z+p->delta)-f;
1009 f = sqrt(v->x*v->x + v->y*v->y + v->z*v->z);
1010 if (f != 0.0) {v->x /= f; v->y /= f; v->z /= f;}
1011}
1012
1013
1014/* converge: from two points of differing sign, converge to zero crossing */
1015
1016static void converge (p1, p2, v, function, p)
1017double v;
1018double (*function)();
1019POINT *p1, *p2, *p;
1020{
1021 int i = 0;
1022 POINT pos, neg;
1023 if (v < 0) {
1024 pos.x = p2->x; pos.y = p2->y; pos.z = p2->z;
1025 neg.x = p1->x; neg.y = p1->y; neg.z = p1->z;
1026 }
1027 else {
1028 pos.x = p1->x; pos.y = p1->y; pos.z = p1->z;
1029 neg.x = p2->x; neg.y = p2->y; neg.z = p2->z;
1030 }
1031 while (1) {
1032 p->x = 0.5*(pos.x + neg.x);
1033 p->y = 0.5*(pos.y + neg.y);
1034 p->z = 0.5*(pos.z + neg.z);
1035 if (i++ == RES) return;
1036 if ((function(p->x, p->y, p->z)) > 0.0)
1037 {pos.x = p->x; pos.y = p->y; pos.z = p->z;}
1038 else {neg.x = p->x; neg.y = p->y; neg.z = p->z;}
1039 }
1040}
Note: See TracBrowser for help on using the repository browser.