implicit.c

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

/*
 * C code from the article
 * "An Implicit Surface Polygonizer"
 * by Jules Bloomenthal, [e-mail]
 * 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 <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <sys/types.h>

#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<<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 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 <math/isosurf/gl.h>

/* 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; i<HASHSIZE; i++) {
      CENTERLIST *l,*next;
      for (l=p->centers[i]; l; l=next) {
          next = l->next;
          myfree(l);
        }
    }

  for (i=0; i<HASHSIZE; i++) {
      CORNERLIST *l,*next;
      for (l=p->corners[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;