Revision: 71570
http://sourceforge.net/p/brlcad/code/71570
Author: starseeker
Date: 2018-08-23 19:56:43 +0000 (Thu, 23 Aug 2018)
Log Message:
-----------
Start cleaning up Bloomenthal's polygonizer for potential use in turning a
Winding Number field into triangles. Untested.
Modified Paths:
--------------
brlcad/trunk/src/libanalyze/CMakeLists.txt
Added Paths:
-----------
brlcad/trunk/src/libanalyze/polygonizer.c
brlcad/trunk/src/libanalyze/polygonizer.h
Modified: brlcad/trunk/src/libanalyze/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/libanalyze/CMakeLists.txt 2018-08-23 17:54:46 UTC (rev
71569)
+++ brlcad/trunk/src/libanalyze/CMakeLists.txt 2018-08-23 19:56:43 UTC (rev
71570)
@@ -52,6 +52,8 @@
analyze_private.h
find_subtracted_shapes.cpp
wnsurface.cxx
+ polygonizer.c
+ polygonizer.h
)
# Local Variables:
Added: brlcad/trunk/src/libanalyze/polygonizer.c
===================================================================
--- brlcad/trunk/src/libanalyze/polygonizer.c (rev 0)
+++ brlcad/trunk/src/libanalyze/polygonizer.c 2018-08-23 19:56:43 UTC (rev
71570)
@@ -0,0 +1,678 @@
+/*
+ * C code from the article
+ * "An Implicit Surface Polygonizer"
+ * http::www.unchainedgeometry.com/jbloom/papers/polygonizer.pdf
+ * by Jules Bloomenthal, ju...@bloomenthal.com
+ * in "Graphics Gems IV", Academic Press, 1994
+
+ * 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. */
+
+/* A Brief Explanation
+ *
+ * The main data structures in the polygonizer represent a hexahedral lattice,
+ * ie, a collection of semi-adjacent cubes, represented as cube centers,
corners,
+ * and edges. The centers and corners are three-dimensional indices rerpesented
+ * by integer i,j,k. The edges are two three-dimensional indices, represented
+ * by integer i1,j1,k1,i2,j2,k2. These indices and associated data are stored
+ * in hash tables.
+ *
+ * The polygonize() routine first allocates memory for the hash tables for the
+ * cube centers, corners, and edges that define the polygonizing lattice. It
+ * then finds a start point, ie, the center of the first lattice cell. It
+ * pushes this cell onto an initially empty stack of cells awaiting processing.
+ * It creates the first cell by computing its eight corners and assigning them
+ * an implicit value.
+ *
+ * polygonize() then enters a loop in which a cell is popped from the stack,
+ * becoming the 'active' cell c. c is (optionally) decomposed (ie, subdivided)
+ * into six tetrahedra; within each transverse tetrahedron (ie, those that
+ * intersect the surface), one or two triangles are produced.
+ *
+ * The six faces of c are tested for intersection with the implicit surface;
for
+ * a transverse face, a new cube is generated and placed on the stack.
+ *
+ * Some of the more important routines include:
+ *
+ * testface (called by polygonize): test given face for surface intersection;
+ * if transverse, create new cube by creating four new corners.
+ * setcorner (called by polygonize, testface): create new cell corner at given
+ * (i,j,k), compute its implicit value, and add to corners hash table.
+ * find (called by polygonize): search for point with given polarity
+ * dotet (called by polygonize) set edge vertices, output triangle by
+ * invoking callback
+ *
+ * The section Cubical Polygonization contains routines to polygonize directly
+ * from the lattice cell rather than first decompose it into tetrahedra;
+ * dotet, however, is recommended over docube.
+ *
+ * The section Storage provides routines to handle the linked lists
+ * in the hash tables.
+ *
+ * The section Vertices contains the following routines.
+ * vertid (called by dotet): given two corner indices defining a cell edge,
+ * test whether the edge has been stored in the hash table; if so, return
its
+ * associated vertex index. If not, compute intersection of edge and
implicit
+ * surface, compute associated surface normal, add vertex to mesh array, and
+ * update hash tables
+ * converge (called by polygonize, vertid): find surface crossing on edge */
+
+#include "common.h"
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <sys/types.h>
+#include "bu/log.h"
+#include "bu/malloc.h"
+#include "./polygonizer.h"
+
+#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<<HASHBIT)-1)
+#define HASH(i,j,k) ((((((i)&MASK)<<HASHBIT)|((j)&MASK))<<HASHBIT)|((k)&MASK))
+#define BIT(i, bit) (((i)>>(bit))&1)
+#define FLIP(i,bit) ((i)^1<<(bit)) /* flip the given bit of i */
+
+typedef struct test { /* test the function for a signed value */
+ point_t p; /* location of test */
+ double value; /* function value at p */
+ int ok; /* if value is of correct sign */
+} TEST;
+
+typedef struct polygonizer_vertex VERTEX; /* surface vertex */
+
+typedef struct polygonizer_vertices VERTICES;
+
+typedef struct triangle {
+ int i1, i2, i3;
+} TRIANGLE;
+
+typedef struct triangles {
+ int count, max;
+ TRIANGLE *ptr;
+} TRIANGLES;
+
+typedef struct corner { /* corner of a cube */
+ int i, j, k; /* (i, j, k) is index within lattice */
+ point_t p; /* location */
+ double value; /* 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 */
+ polygonize_func_t function; /* implicit surface function */
+ void *d; /* implicit surface function data */
+ polygonize_triproc_t triproc; /* triangle output function */
+ double size, delta; /* cube size, normal delta */
+ int bounds; /* cube range within lattice */
+ point_t 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;
+
+
+/* setcorner: return corner with the given lattice location
+ set (and cache) its function value */
+CORNER *
+setcorner (PROCESS *p, int i, int j, int k)
+{
+ /* for speed, do corner value caching here */
+ CORNER *c = (CORNER *) bu_calloc(1, sizeof(CORNER), "corner");
+ int index = HASH(i, j, k);
+ CORNERLIST *l = p->corners[index];
+ c->i = i; c->p[X] = p->start[X]+((double)i-.5)*p->size;
+ c->j = j; c->p[Y] = p->start[Y]+((double)j-.5)*p->size;
+ c->k = k; c->p[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 *) bu_calloc(1, sizeof(CORNERLIST), "cornerlist");
+ l->i = i; l->j = j; l->k = k;
+ l->value = c->value = p->function(c->p, p->d);
+ l->next = p->corners[index];
+ p->corners[index] = l;
+ return c;
+}
+
+/* setcenter: set (i,j,k) entry of table[]
+ * return 1 if already set; otherwise, set and return 0 */
+
+int
+setcenter(CENTERLIST *table[], int i, int j, int k)
+{
+ int index = HASH(i, j, k);
+ CENTERLIST *snew, *l, *q = table[index];
+ for (l = q; l != NULL; l = l->next)
+ if (l->i == i && l->j == j && l->k == k) return 1;
+ snew = (CENTERLIST *) bu_calloc(1, sizeof(CENTERLIST), "centerlist");
+ snew->i = i; snew->j = j; snew->k = k; snew->next = q;
+ table[index] = snew;
+ return 0;
+}
+
+
+/* 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 */
+void
+testface(int i, int j, int k, CUBE *old, int face, int c1, int c2, int c3, int
c4, PROCESS *p)
+{
+ CUBE cnew;
+ CUBES *oldcubes = p->cubes;
+ static int facebit[6] = {2, 2, 1, 1, 0, 0};
+ int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0, 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) return;
+ if (setcenter(p->centers, i, j, k)) return;
+
+ /* create new cube: */
+ cnew.i = i;
+ cnew.j = j;
+ cnew.k = k;
+ for (n = 0; n < 8; n++) cnew.corners[n] = NULL;
+ cnew.corners[FLIP(c1, bit)] = old->corners[c1];
+ cnew.corners[FLIP(c2, bit)] = old->corners[c2];
+ cnew.corners[FLIP(c3, bit)] = old->corners[c3];
+ cnew.corners[FLIP(c4, bit)] = old->corners[c4];
+ for (n = 0; n < 8; n++)
+ if (cnew.corners[n] == NULL)
+ cnew.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 *) bu_calloc(1, sizeof(CUBES), "cubes");
+ p->cubes->cube = cnew;
+ p->cubes->next = oldcubes;
+}
+
+/* find: search for point with value of given sign (0: neg, 1: pos) */
+
+TEST
+find(int sign, PROCESS *pr, point_t p)
+{
+ int i;
+ TEST test;
+ double range = pr->size;
+ test.ok = 1;
+ for (i = 0; i < 10000; i++) {
+ test.p[X] = p[X]+range*(RAND()-0.5);
+ test.p[Y] = p[Y]+range*(RAND()-0.5);
+ test.p[Z] = p[Z]+range*(RAND()-0.5);
+ test.value = pr->function(test.p, pr->d);
+ if (sign == (test.value > 0.0)) return test;
+ range = range*1.0005; /* slowly expand search outwards */
+ }
+ test.ok = 0;
+ return test;
+}
+
+
+/* getedge: return vertex id for edge; return -1 if not set */
+
+int
+getedge(EDGELIST *table[], int i1, int j1, int k1, int i2, int j2, int 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;
+}
+
+
+/* converge: from two points of differing sign, converge to zero crossing */
+
+void
+converge(point_t *p1, point_t *p2, double v, polygonize_func_t function,
point_t *p, void *d)
+{
+ int i = 0;
+ point_t pos, neg;
+ if (v < 0) {
+ VMOVE(pos, *p2);
+ VMOVE(neg, *p1);
+ }
+ else {
+ VMOVE(pos, *p1);
+ VMOVE(neg, *p2);
+ }
+ 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, d)) > 0.0) {
+ VMOVE(pos, *p);
+ } else {
+ VMOVE(neg, *p);
+ }
+ }
+}
+
+
+/* vnormal: compute unit length surface normal at point */
+
+void
+vnormal(point_t *point, PROCESS *p, point_t *v)
+{
+ point_t temp;
+ double f = p->function(*point, p->d);
+ VMOVE(temp, *point);
+ temp[X] = temp[X] + p->delta;
+ (*v)[X] = p->function(temp, p->d)-f;
+ VMOVE(temp, *point);
+ temp[Y] = temp[Y] + p->delta;
+ (*v)[Y] = p->function(temp, p->d)-f;
+ VMOVE(temp, *point);
+ temp[Z] = temp[Z] + p->delta;
+ (*v)[Z] = p->function(temp, p->d)-f;
+ f = MAGNITUDE(*v);
+ if (!NEAR_ZERO(f, VUNITIZE_TOL)) {
+ VSCALE(*v, *v, 1.0/f);
+ }
+}
+
+
+/* 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 */
+
+
+/* addtovertices: add v to sequence of vertices */
+
+void
+addtovertices(VERTICES *vertices, VERTEX v)
+{
+ if (vertices->count == vertices->max) {
+ int i;
+ VERTEX *vnew;
+ vertices->max = vertices->count == 0 ? 10 : 2*vertices->count;
+ vnew = (VERTEX *) bu_calloc(vertices->max, sizeof(VERTEX), "vertex");
+ for (i = 0; i < vertices->count; i++) vnew[i] = vertices->ptr[i];
+ if (vertices->ptr != NULL) free((char *) vertices->ptr);
+ vertices->ptr = vnew;
+ }
+ vertices->ptr[vertices->count++] = v;
+}
+
+
+/* setedge: set vertex id for edge */
+
+void
+setedge(EDGELIST *table[], int i1, int j1, int k1, int i2, int j2, int k2, int
vid)
+{
+ unsigned int index;
+ EDGELIST *enew;
+ 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);
+ enew = (EDGELIST *) bu_calloc(1, sizeof(EDGELIST), "edgelist");
+ enew->i1 = i1; enew->j1 = j1; enew->k1 = k1;
+ enew->i2 = i2; enew->j2 = j2; enew->k2 = k2;
+ enew->vid = vid;
+ enew->next = table[index];
+ table[index] = enew;
+}
+
+int
+vertid(CORNER *c1, CORNER *c2, PROCESS *p)
+{
+ VERTEX v;
+ point_t 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 */
+ VMOVE(a, c1->p);
+ VMOVE(b, c2->p);
+ converge(&a, &b, c1->value, p->function, &v.position, p->d); /* 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;
+}
+
+
+/**** 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 *cube, int c1, int c2, int c3, int 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;
+ apos = (a->value > 0.0);
+ bpos = (b->value > 0.0);
+ cpos = (c->value > 0.0);
+ dpos = (d->value > 0.0);
+ if (apos) index += 8;
+ if (bpos) index += 4;
+ if (cpos) index += 2;
+ if (dpos) 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 (int edge, int 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 */
+
+int
+otherface (int edge, int face)
+{
+ int other = leftface[edge];
+ return face == other? rightface[edge] : other;
+}
+
+
+/* makecubetable: create the 256 entry table for cubical polygonization */
+
+void
+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 *) bu_calloc(1, sizeof(INTLISTS),
"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 *) bu_calloc(1, sizeof(INTLIST),
"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;
+ }
+ }
+}
+
+/**** An Implicit Surface Polygonizer ****/
+
+int
+polygonalize(
+ polygonize_func_t pf,
+ void *pf_d,
+ fastf_t size,
+ int bounds,
+ point_t p_s,
+ polygonize_triproc_t triproc,
+ int mode)
+{
+ PROCESS p;
+ int n, noabort;
+ TEST in, out;
+
+ p.function = pf;
+ p.triproc = triproc;
+ p.size = size;
+ p.bounds = bounds;
+ p.delta = size/(double)(RES*RES);
+ p.d = pf_d;
+
+ /* allocate hash tables and build cube polygon table: */
+ p.centers = (CENTERLIST **) bu_calloc(HASHSIZE, sizeof(CENTERLIST *),
"hashsize centerlist");
+ p.corners = (CORNERLIST **) bu_calloc(HASHSIZE,sizeof(CORNERLIST *),
"hashsize, cornerlist");
+ p.edges = (EDGELIST **) bu_calloc(2*HASHSIZE,sizeof(EDGELIST *),
"2*hashsize, edgelist");
+ makecubetable();
+
+ /* find point on surface, beginning search at point p_s: */
+ srand(1);
+ in = find(1, &p, p_s);
+ out = find(0, &p, p_s);
+ if (!in.ok || !out.ok) {
+ bu_log ("polygonizer: Error, can't find starting point");
+ return -1;
+ }
+ converge(&in.p, &out.p, in.value, p.function, &p.start, p.d);
+
+ /* push initial cube on stack: */
+ p.cubes = (CUBES *) bu_calloc(1, sizeof(CUBES), "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 */
+ CUBE c;
+ CUBES *temp = p.cubes;
+ c = p.cubes->cube;
+
+ noabort = mode == POLYGONIZE_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) return -1;
+
+ /* pop current cube from stack */
+ p.cubes = p.cubes->next;
+ free((char *) temp);
+ /* 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);
+ }
+ return 0;
+}
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
+
Property changes on: brlcad/trunk/src/libanalyze/polygonizer.c
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Added: brlcad/trunk/src/libanalyze/polygonizer.h
===================================================================
--- brlcad/trunk/src/libanalyze/polygonizer.h (rev 0)
+++ brlcad/trunk/src/libanalyze/polygonizer.h 2018-08-23 19:56:43 UTC (rev
71570)
@@ -0,0 +1,96 @@
+/*
+ * C code from the article
+ * "An Implicit Surface Polygonizer"
+ * http::www.unchainedgeometry.com/jbloom/papers/polygonizer.pdf
+ * by Jules Bloomenthal, ju...@bloomenthal.com
+ * in "Graphics Gems IV", Academic Press, 1994
+
+ * 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 "common.h"
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <sys/types.h>
+
+#include "vmath.h"
+
+
+#define POLYGONIZE_TET 0 /** decompose cube and polygonize six tetrahedra */
+#define POLYGONIZE_NOTET 1 /** polygonize cube directly */
+
+/**
+ * Callback function signature for the function used to decide if the query
+ * point q is inside or outside the surface. The value in d is intended to
+ * hold any information needed when evaluating the function.
+ */
+typedef int (*polygonize_func_t)(point_t q, void *d);
+
+
+
+struct polygonizer_vertex { /* surface vertex */
+ point_t position;
+ point_t normal; /* surface normal */
+};
+
+struct polygonizer_vertices { /* list of vertices in polygonization
*/
+ int count, max; /* # vertices, max # allowed */
+ struct polygonizer_vertex *ptr; /* dynamically allocated */
+};
+
+/**
+ * Callback function signature for the function called when a triangle is
+ * generated - this is how the results of polygonalization are communicated
+ * to the calling application.
+ *
+ * i1, i2, i3 (indices into the vertex array defining the triangle)
+ * vertices (the vertex array, indexed from 0)
+ *
+ * vertices are ccw when viewed from the out (positive) side in a left-handed
+ * coordinate system
+ *
+ * vertex normals point outwards
+ *
+ * Function should return 1 to continue the polygonalization, 0 to abort
+ */
+typedef int (*polygonize_triproc_t)(int i1, int i2, int i3, struct
polygonizer_vertices *vertices);
+
+/**
+ * @brief
+ * An Implicit Surface Polygonizer
+ *
+ * Given a function to determine inside/outside for a surface, construct a
+ * triangle approximation of that surface.
+ *
+ * @param[in] pf The implicit surface function - must return negative
for inside, positive for outside
+ * @param[in] pf_d Caller supplied data for pf evaluation - may be NULL
+ * @param[in] size Width of the partitioning cube
+ * @param[in] bounds Maximum range of cubes (+/- on the three axes) from
first cube
+ * @param[in] p_s Coordinates of a starting point on or near the surface
+ * @param[in] triproc Callback called for each triangle generated by the
polygonalizer
+ * @param[in] mode Either POLYGONALIZE_TET (default) or POLYGONALIZE_NOTET
+ */
+int
+polygonalize(
+ polygonize_func_t pf,
+ void *pf_d,
+ fastf_t size,
+ int bounds,
+ point_t p_s,
+ polygonize_triproc_t triproc,
+ int mode);
+
+
+/*
+ * Local Variables:
+ * mode: C
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Property changes on: brlcad/trunk/src/libanalyze/polygonizer.h
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
BRL-CAD Source Commits mailing list
brlcad-commits@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
