bloomenthal.cpp

Go to the documentation of this file.

/***** bloomenthal.cpp */

/*
 * C++ code adapted from the
 * ANSI C code from the article
 * "An Implicit Surface Polygonizer"
 * by Jules Bloomenthal
 * 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
 *
 * Last modified 11 Jan 95 by Jules Bloomenthal and Paul Heckbert
 * 25 Sep 2000 by John Hart to C++
 */

#include <math.h>
#include <malloc.h>
#include <stdio.h>
#include <sys/types.h>
#include <stdlib.h> // for (int)abs(int)

#include "bloomenthal.h"

/**** An Implicit Surface Polygonizer ****/


/* polygonize: polygonize the implicit surface function
 *   arguments are:
 *       double function (double x, double y, double z)
 *           the implicit surface function given an arbitrary point
 *           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)
 *               jbVertices 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
 *   returns error or NULL
 */

char * jbProcess::polygonize (
//    double (*in_function)(gmVector3),
    Implicit *in_f,
    double in_size,
    int in_bounds,
    gmVector3 x,
    int (*in_triproc)(int i1, int i2, int i3, jbVertices vertices))
{
    int n;
    jbTest in, out;
    
    f = in_f;
    triproc = in_triproc;
    size = in_size;
    bounds = in_bounds;
    delta = size/(double)(RES*RES);

    /* allocate hash tables: */
    centers = (jbCenterlist **) mycalloc(HSIZE, sizeof(jbCenterlist *));
    corners = (jbCornerlist **) mycalloc(HSIZE,   sizeof(jbCornerlist *));
    edges   = (jbEdgelist   **) mycalloc(2*HSIZE, sizeof(jbEdgelist *));

    vertices.count = vertices.max = 0; /* no vertices yet */
    vertices.ptr = NULL;
    
    /* find point on surface, beginning search at (x, y, z):  */
    srand(1);
    in = find(1 == 1, x);
    out = find(0 == 1, x);
    if (!in.ok || !out.ok) {
        freeprocess();
        if (!in.ok) printf ("in not ok\n");
        if (!out.ok) printf ("out not ok\n");
        return "can't find starting point";
    }
    start = converge(in.pt, out.pt, in.value);

    /* push initial cube on stack: */
    cubes = (jbCubes *) mycalloc(1, sizeof(jbCubes)); /* list of 1 */
    cubes->cube.i = cubes->cube.j = cubes->cube.k = 0;
    cubes->next = NULL;

    /* set corners of initial cube: */
    for (n = 0; n < 8; n++)
        cubes->cube.corners[n] = \
            setcorner(BIT(n,2), BIT(n,1), BIT(n,0));

    setcenter(0, 0, 0);

    while (cubes != NULL) { /* process active cubes till none left */
        jbCube c;
        jbCubes *temp = cubes;
        c = cubes->cube;

        /* decompose into tetrahedra and polygonize: */
        if (!(dotet(&c, JB_LBN, JB_LTN, JB_RBN, JB_LBF) &&
              dotet(&c, JB_RTN, JB_LTN, JB_LBF, JB_RBN) &&
              dotet(&c, JB_RTN, JB_LTN, JB_LTF, JB_LBF) &&
              dotet(&c, JB_RTN, JB_RBN, JB_LBF, JB_RBF) &&
              dotet(&c, JB_RTN, JB_LBF, JB_LTF, JB_RBF) &&
              dotet(&c, JB_RTN, JB_LTF, JB_RTF, JB_RBF)))
         {
              freeprocess();
              return "aborted";
         }

        /* pop current cube from stack */
        cubes = cubes->next;
        free((char *) temp);
        /* test six face directions, maybe add to stack: */
        testface(c.i-1, c.j, c.k, &c, JB_L, JB_LBN, JB_LBF, JB_LTN, JB_LTF);
        testface(c.i+1, c.j, c.k, &c, JB_R, JB_RBN, JB_RBF, JB_RTN, JB_RTF);
        testface(c.i, c.j-1, c.k, &c, JB_B, JB_LBN, JB_LBF, JB_RBN, JB_RBF);
        testface(c.i, c.j+1, c.k, &c, JB_T, JB_LTN, JB_LTF, JB_RTN, JB_RTF);
        testface(c.i, c.j, c.k-1, &c, JB_N, JB_LBN, JB_LTN, JB_RBN, JB_RTN);
        testface(c.i, c.j, c.k+1, &c, JB_F, JB_LBF, JB_LTF, JB_RBF, JB_RTF);
    }
    freeprocess();

    return NULL;
}

char * jbProcess::marchingcubes (
    Implicit *in_f,
    double in_size,
    gmVector3 x0, gmVector3 x1,
    int (*in_triproc)(int i1, int i2, int i3, jbVertices vertices))
{
    int i,j,k,n,i1,j1,k1;
        jbCube cube;
        int pos;
    
    f = in_f;
    triproc = in_triproc;
    size = in_size;
    bounds = 0; // not needed this time
    delta = size/(double)(RES*RES);

    /* allocate hash tables: */
    centers = (jbCenterlist **) mycalloc(HSIZE, sizeof(jbCenterlist *));
    corners = (jbCornerlist **) mycalloc(HSIZE,   sizeof(jbCornerlist *));
    edges   = (jbEdgelist   **) mycalloc(2*HSIZE, sizeof(jbEdgelist *));

    vertices.count = vertices.max = 0; /* no vertices yet */
    vertices.ptr = NULL;


        // add noise to hopefully nudge precise surfaces off lattice

        srand(1);
        start = x1 - size*gmVector3(.5,.5,.5) - 0.001*size*gmVector3(RAND(),RAND(),RAND());

        i1 = (int)((x1[0] - x0[0])/size);
        j1 = (int)((x1[1] - x0[1])/size);
        k1 = (int)((x1[2] - x0[2])/size);

    for (k = 0; k <= k1; k++) {
                for (j = 0; j <= j1; j++) {
                        for (i = 0; i <= i1; i++) {

                                cube.i = i;
                                cube.j = j;
                                cube.k = k;

                                pos = 0;

                            /* set corners of cube: */
                                for (n = 0; n < 8; n++) {
                                        cube.corners[n] = setcorner(i+BIT(n,2), j+BIT(n,1), k+BIT(n,0));
                                        if (cube.corners[n]->value > 0.0) pos++;
                                }

                                if (pos == 0 || pos == 8)
                                        continue;       /* cube probably empty */

                                dotet(&cube, JB_LBN, JB_LTN, JB_RBN, JB_LBF);
                                dotet(&cube, JB_RTN, JB_LTN, JB_LBF, JB_RBN);
                                dotet(&cube, JB_RTN, JB_LTN, JB_LTF, JB_LBF);
                                dotet(&cube, JB_RTN, JB_RBN, JB_LBF, JB_RBF);
                                dotet(&cube, JB_RTN, JB_LBF, JB_LTF, JB_RBF);
                                dotet(&cube, JB_RTN, JB_LTF, JB_RTF, JB_RBF);
                        }
                }
        }

    freeprocess();

    return NULL;
}


/* freeprocess: free all allocated memory */

void jbProcess::freeprocess() {
    int index;
    jbCornerlist *l, *lnext;
    jbCenterlist *cl, *clnext;
    jbEdgelist *edge, *edgenext;

    for (index = 0; index < HSIZE; index++)
        for (l = corners[index]; l; l = lnext) {
            lnext = l->next;
            free((char *) l);           /* free CORNERLIST */
        }
    for (index = 0; index < HSIZE; index++)
        for (cl = centers[index]; cl; cl = clnext) {
            clnext = cl->next;
            free((char *) cl);          /* free CENTERLIST */
        }
    for (index = 0; index < 2*HSIZE; index++)
        for (edge = edges[index]; edge; edge = edgenext) {
            edgenext = edge->next;
            free((char *) edge);        /* free EDGELIST */
        }
    free((char *) edges);            /* free array EDGELIST ptrs */
    free((char *) corners);          /* free array CORNERLIST ptrs */
    free((char *) centers);          /* free array CENTERLIST ptrs */
    if (vertices.ptr)
        free((char *) vertices.ptr); /* free VERTEX array */
}

/* 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 jbProcess::testface (
    int i, int j, int k,
    jbCube *old,
    int face, int c1, int c2, int c3, int c4)
{
    jbCube newcube;
    jbCubes *oldcubes = cubes;
    static int facebit[6] = {2, 2, 1, 1, 0, 0};
    int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0;
    int bit = facebit[face];

    /* test if  no surface crossing, cube out of bounds, or prev. visited? */
    if ((old->corners[c2]->value > 0) == pos &&
        (old->corners[c3]->value > 0) == pos &&
        (old->corners[c4]->value > 0) == pos) return;
    if (abs(i) > bounds || abs(j) > bounds || abs(k) > bounds)
        return;
    if (setcenter(i, j, k)) return;

    /* create new cube: */
    newcube.i = i;
    newcube.j = j;
    newcube.k = k;
    for (n = 0; n < 8; n++) newcube.corners[n] = NULL;
    newcube.corners[FLIP(c1, bit)] = old->corners[c1];
    newcube.corners[FLIP(c2, bit)] = old->corners[c2];
    newcube.corners[FLIP(c3, bit)] = old->corners[c3];
    newcube.corners[FLIP(c4, bit)] = old->corners[c4];
    for (n = 0; n < 8; n++)
        if (newcube.corners[n] == NULL) newcube.corners[n] =
            setcorner(i+BIT(n,2), j+BIT(n,1), k+BIT(n,0));

    /*add cube to top of stack: */
    cubes = (jbCubes *) mycalloc(1, sizeof(jbCubes));
    cubes->cube = newcube;
    cubes->next = oldcubes;
}


/* setpoint: set point location given lattice location */

gmVector3 jbProcess::setpoint (int i, int j, int k)
{
        return start - size*gmVector3(i-0.5,j-0.5,k-0.5);
}


/* setcorner: return corner with the given lattice location
   set (and cache) its function value */

jbCornerlist * jbProcess::setcorner (int i, int j, int k)
{
    /* for speed, do corner value caching here */
    int index = HASH(i, j, k);
    jbCornerlist *l = corners[index];
    gmVector3 pt;
    
    for (; l != NULL; l = l->next)
        if (l->i == i && l->j == j && l->k == k) return l;
            
    pt = setpoint (i, j, k);
    l = (jbCornerlist *) mycalloc(1, sizeof(jbCornerlist));
    l->i = i; l->j = j; l->k = k;
    l->value = f->proc(pt);
    l->next = corners[index];
    corners[index] = l;
    return l;
}


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

jbTest jbProcess::find(int sign, gmVector3 x)
{
    int i;
        jbTest test;
    double range = size;

    test.ok = 1;
    for (i = 0; i < 10000; i++) {
                test.pt = x + range*gmVector3(RAND()-0.5,RAND()-0.5,RAND()-0.5);
        test.value = f->proc(test.pt);
        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 jbProcess::dotet (jbCube *cube, int c1, int c2, int c3, int c4) {
    jbCornerlist *a = cube->corners[c1],
                                 *b = cube->corners[c2],
                                 *c = cube->corners[c3],
                                 *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 now 4-bit number equal to one of the 16 possible cases */
    if (apos != bpos) e1 = vertid(a, b);
    if (apos != cpos) e2 = vertid(a, c);
    if (apos != dpos) e3 = vertid(a, d);
    if (bpos != cpos) e4 = vertid(b, c);
    if (bpos != dpos) e5 = vertid(b, d);
    if (cpos != dpos) e6 = vertid(c, d);
    /* 14 productive tet. cases (0000 and 1111 do not yield polygons */
    switch (index) {
        case 1:  return triproc(e5, e6, e3, vertices);
        case 2:  return triproc(e2, e6, e4, vertices);
        case 3:  return triproc(e3, e5, e4, vertices) &&
                        triproc(e3, e4, e2, vertices);
        case 4:  return triproc(e1, e4, e5, vertices);
        case 5:  return triproc(e3, e1, e4, vertices) &&
                        triproc(e3, e4, e6, vertices);
        case 6:  return triproc(e1, e2, e6, vertices) &&
                        triproc(e1, e6, e5, vertices);
        case 7:  return triproc(e1, e2, e3, vertices);
        case 8:  return triproc(e1, e3, e2, vertices);
        case 9:  return triproc(e1, e5, e6, vertices) &&
                        triproc(e1, e6, e2, vertices);
        case 10: return triproc(e1, e3, e6, vertices) &&
                        triproc(e1, e6, e4, vertices);
        case 11: return triproc(e1, e5, e4, vertices);
        case 12: return triproc(e3, e2, e4, vertices) &&
                        triproc(e3, e4, e5, vertices);
        case 13: return triproc(e6, e2, e4, vertices);
        case 14: return triproc(e5, e3, e6, vertices);
    }
    return 1;
}


/**** Storage ****/


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

char *mycalloc (int nitems, int nbytes) {
   char *ptr = (char *) 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 jbProcess::setcenter(int i, int j, int k) {
    int index = HASH(i, j, k);
    jbCenterlist *newlist, *l, *q = centers[index];

    for (l = q; l != NULL; l = l->next)
        if (l->i == i && l->j == j && l->k == k) return 1;
    newlist = (jbCenterlist *) mycalloc(1, sizeof(jbCenterlist));
    newlist->i = i; newlist->j = j; newlist->k = k; newlist->next = q;
    centers[index] = newlist;
    return 0;
}


/* setedge: set vertex id for edge */

void jbProcess::setedge (int i1, int j1, int k1, int i2, int j2, int k2, int vid) {
    unsigned int index;
    jbEdgelist *newlist;

    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);
    newlist = (jbEdgelist *) mycalloc(1, sizeof(jbEdgelist));
    newlist->i1 = i1; newlist->j1 = j1; newlist->k1 = k1;
    newlist->i2 = i2; newlist->j2 = j2; newlist->k2 = k2;
    newlist->vid = vid;
    newlist->next = edges[index];
    edges[index] = newlist;
}


/* getedge: return vertex id for edge; return -1 if not set */

int jbProcess::getedge (int i1, int j1, int k1, int i2, int j2, int k2)
{
    jbEdgelist *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 = edges[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 jbProcess::vertid (jbCornerlist *c1, jbCornerlist *c2) {
    jbVertex v;
    gmVector3 a, b;
    int vid =
        getedge(c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
    if (vid != -1) return vid;                /* previously computed */
    a = setpoint (c1->i, c1->j, c1->k);
    b = setpoint (c2->i, c2->j, c2->k);
    v.position = converge (a, b, c1->value); /* posn.  */
//    vid = addtovertices(&vertices, v);                   /* save   */
        vid = vertices.add(v);
    setedge(c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
    return vid;
}


/* addtovertices: add v to sequence of vertices and return its id */

int jbVertices::add(jbVertex v) {
        int i;
        jbVertex *newv;
        
        if (count == max) {
      max = count == 0 ? 10 : 2*count;
      newv = (jbVertex *) mycalloc((unsigned)max,sizeof(jbVertex));
      for (i = 0; i < count; i++) newv[i] = ptr[i];
      if (ptr != NULL) free((char *) ptr);
      ptr = newv;
   }
   ptr[count++] = v;
   return (count-1);
}

/* converge: from two points of differing sign, converge to surface */

gmVector3 jbProcess::converge (gmVector3 p1, gmVector3 p2, double v)
{
    int i = 0;
    gmVector3 pos, neg, p;

    if (v < 0) {
                pos = p2;
                neg = p1;
    }
    else {
                pos = p1;
                neg = p2;
    }
    for (;;) {
                p = 0.5*(pos+neg);
        if (i++ == RES) return p;
        if ((f->proc(p)) > 0.0) {
                        pos = p;
        } else {
                        neg = p;
                }
    }
}

---