Logo Search packages:      
Sourcecode: viewmol version File versions  Download package

marcub.c

/*******************************************************************************
*                                                                              *
*                                   Viewmol                                    *
*                                                                              *
*                               M A R C U B . C                                *
*                                                                              *
*                 Copyright (c) Joerg-R. Hill, October 2003                    *
*                                                                              *
********************************************************************************
*
* $Id: marcub.c,v 1.6 2003/11/07 11:06:36 jrh Exp $
* $Log: marcub.c,v $
* Revision 1.6  2003/11/07 11:06:36  jrh
* Release 2.4
*
* Revision 1.5  2000/12/10 15:10:52  jrh
* Release 2.3
*
* Revision 1.4  1999/05/24 01:26:15  jrh
* Release 2.2.1
*
* Revision 1.3  1999/02/07 21:51:14  jrh
* Release 2.2
*
* Revision 1.2  1998/01/26 00:48:24  jrh
* Release 2.1
*
* Revision 1.1  1996/12/10  18:41:59  jrh
* Initial revision
*
*
* This function implements the drawing of the basic patterns of the improved
* marching cube algorithm as described in W. Heiden, T. Goetze, and J. Brick-
* mann: J. Comp. Chem. 14 (1993), 246-250. The calculation of the grid has to be
* done in a separate function. The parameters passed to this function are:
* xstart, ystart, and zstart - the cartesian coordinates of the lower left front
* corner of the volume which is used to draw the isosurface (the origin of the
* grid); grid - a pointer to an array which contains the values of the function
* at each grid point, three dimensional, stored row by row, step - the distance
* between two grid points; xres, yres, and zres - the number of times step has
* to be added to xstart, ystart, or zstart to walk through the whole grid, is
* equal to the number of grid points minus one in each direction; level - a
* floating point number which determines which isosurface to select;
* begin, end, vert, and normal - pointers to the drawing functions.
* The numbering of the cube's corners is as follows:
*
*                              7+------+6
*          y                   /|     /|
*          ^  z              3+------+2|
*          | /                |4+----|-+5
*          |/                 |/     |/
*          +----> x          0+------+1 
*                                      
*
*/
#include<math.h>
#include<stdio.h>
#include<GL/gl.h>
#include "marcub.h"
#include "viewmol.h"

#define FALSE 0
#define TRUE  1
#define index(i, j, k) ((grid->ngridpoints[1]+1)*(grid->ngridpoints[2]+1)*(k)+(grid->ngridpoints[2]+1)*(j)+(i))

/* static double norms[27][3], vertices[27][3];
static int vcount=0, recordit=FALSE; */
/* Prototypes */
void set_vertex(double, double, double, int, int, struct GRIDOBJECT *, int, int, int,
                int, void (*)(), void (*)());

/* Global variables */
static int corner[] = {
0, 1, 0, 3, 0, 4,                                    /* 1 */
1, 5, 1, 2, 1, 0, 
0, 3, 0, 4, 1, 2, 1, 5,
3, 7, 3, 0, 3, 2, 
0, 4, 0, 1, 3, 7, 3, 2,                              /* 5 */
1, 2, 1, 0, 1, 5, -1, 3, 2, 3, 7, 3, 0, 
0, 4, 1, 5, 3, 7, 1, 2, 3, 2, 
2, 3, 2, 1, 2, 6, 
0, 1, 0, 3, 0, 4, -1, 2, 3, 2, 1, 2, 6, 
1, 5, 2, 6, 1, 0, 2, 3,                              /* 10 */
1, 5, 2, 6, 0, 4, 2, 3, 0, 3, 
2, 6, 3, 7, 2, 1, 3, 0, 
3, 7, 0, 4, 2, 6, 0, 1, 2, 1, 
2, 6, 3, 7, 1, 5, 3, 0, 1, 0, 
1, 5, 2, 6, 0, 4, 3, 7,                              /* 15 */
4, 5, 4, 0, 4, 7, 
0, 3, 4, 7, 0, 1, 4, 5, 
1, 0, 1, 5, 1, 2, -1, 4, 5, 4, 0, 4, 7, 
0, 3, 4, 7, 1, 2, 4, 5, 1, 5, 
3, 2, 3, 7, 3, 0, -1, 4, 5, 4, 0, 4, 7,              /* 20 */
0, 1, 3, 2, 4, 5, 3, 7, 4, 7, 
1, 0, 1, 5, 1, 2, -1, 3, 2, 3, 7, 3, 0, -1, 4, 5, 4, 0, 4, 7, 
1, 5, 1, 2, 4, 5, 3, 2, 4, 7, 3, 7, 
2, 3, 2, 1, 2, 6, -1, 4, 5, 4, 0, 4, 7, 
0, 3, 4, 7, 0, 1, 4, 5, -1, 2, 3, 2, 1, 2, 6,        /* 25 */
1, 5, 2, 6, 1, 0, 2, 3, -1, 4, 5, 4, 0, 4, 7, 
2, 6, 2, 3, 1, 5, 0, 3, 4, 7, -1, 1, 5, 4, 7, 4, 5, 
2, 6, 3, 7, 2, 1, 3, 0, -1, 4, 5, 4, 0, 4, 7, 
2, 6, 3, 7, 2, 1, 4, 5, 0, 1, -1, 3, 7, 4, 7, 4, 5, 
4, 5, 4, 0, 4, 7, -1, 2, 6, 3, 7, 1, 5, 3, 0, 1, 0,  /* 30 */
6, 2, 7, 3, 5, 1, 7, 4, 5, 4, 
5, 4, 5, 6, 5, 1, 
0, 1, 0, 3, 0, 4, -1, 5, 4, 5, 6, 5, 1, 
1, 2, 1, 0, 5, 6, 5, 4, 
1, 2, 0, 3, 5, 6, 0, 4, 5, 4,                        /* 35 */
3, 2, 3, 7, 3, 0, -1, 5, 4, 5, 6, 5, 1, 
0, 4, 0, 1, 3, 7, 3, 2, -1, 5, 4, 5, 6, 5, 1, 
1, 2, 1, 0, 5, 6, 5, 4, -1, 3, 2, 3, 7, 3, 0, 
3, 7, 0, 4, 3, 2, 5, 6, 1, 2, -1, 0, 4, 5, 4, 5, 6, 
2, 3, 2, 1, 2, 6, -1, 5, 4, 5, 6, 5, 1,              /* 40 */
0, 1, 0, 3, 0, 4, -1, 2, 3, 2, 1, 2, 6, -1, 5, 4, 5, 6, 5, 1, 
1, 0, 5, 4, 2, 3, 5, 6, 2, 6, 
0, 4, 5, 4, 0, 3, 5, 6, 2, 3, 2, 6, 
2, 6, 3, 7, 2, 1, 3, 0, -1, 5, 4, 5, 6, 5, 1, 
5, 4, 5, 6, 5, 1, -1, 3, 7, 0, 4, 2, 6, 0, 1, 2, 1,  /* 45 */
3, 7, 3, 0, 2, 6, 1, 0, 5, 4, -1, 2, 6, 5, 4, 5, 6, 
7, 3, 4, 0, 6, 2, 4, 5, 6, 5, 
4, 0, 4, 7, 5, 1, 5, 6, 
4, 7, 5, 6, 0, 3, 5, 1, 0, 1, 
5, 6, 1, 2, 4, 7, 1, 0, 4, 0,                        /* 50 */
1, 2, 0, 3, 5, 6, 4, 7, 
4, 0, 4, 7, 5, 1, 5, 6, -1, 3, 2, 3, 7, 3, 0, 
3, 2, 3, 7, 0, 1, 4, 7, 5, 6, -1, 0, 1, 5, 6, 5, 1, 
3, 2, 3, 7, 3, 0, -1, 5, 6, 1, 2, 4, 7, 1, 0, 4, 0, 
6, 5, 2, 1, 7, 4, 2, 3, 7, 3,                        /* 55 */
4, 0, 4, 7, 5, 1, 5, 6, -1, 2, 3, 2, 1, 2, 6, 
2, 3, 2, 1, 2, 6, -1, 4, 7, 5, 6, 0, 3, 5, 1, 0, 1, 
2, 3, 1, 0, 2, 6, 4, 7, 5, 6, -1, 1, 0, 4, 0, 4, 7, 
7, 4, 6, 5, 3, 0, 6, 2, 3, 2, 
2, 6, 3, 7, 2, 1, 3, 0, -1, 4, 0, 4, 7, 5, 1, 5, 6,  /* 60 */
5, 6, 1, 5, 4, 7, 0, 1, 3, 7, 1, 2, 2, 6,
3, 7, 0, 3, 2, 6, 0, 1, 5, 6, 0, 4, 4, 7,
6, 2, 7, 3, 6, 5, 7, 4, 
7, 6, 7, 4, 7, 3, 
0, 1, 0, 3, 0, 4, -1, 7, 6, 7, 4, 7, 3,              /* 65 */
1, 0, 1, 5, 1, 2, -1, 7, 6, 7, 4, 7, 3, 
0, 4, 1, 5, 0, 3, 1, 2, -1, 7, 6, 7, 4, 7, 3, 
3, 0, 3, 2, 7, 4, 7, 6, 
3, 2, 7, 6, 0, 1, 7, 4, 0, 4, 
3, 0, 3, 2, 7, 4, 7, 6, -1, 1, 0, 1, 5, 1, 2,        /* 70 */
1, 5, 1, 2, 0, 4, 3, 2, 7, 6, -1, 0, 4, 7, 6, 7, 4, 
2, 3, 2, 1, 2, 6, -1, 7, 6, 7, 4, 7, 3, 
0, 1, 0, 3, 0, 4, -1, 2, 3, 2, 1, 2, 6, -1, 7, 6, 7, 4, 7, 3, 
1, 5, 2, 6, 1, 0, 2, 3, -1, 7, 6, 7, 4, 7, 3, 
7, 6, 7, 4, 7, 3, -1, 1, 5, 2, 6, 0, 4, 2, 3, 0, 3,  /* 75 */
3, 0, 2, 1, 7, 4, 2, 6, 7, 6, 
0, 4, 0, 1, 7, 4, 2, 1, 7, 6, 2, 6, 
1, 5, 2, 6, 1, 0, 7, 4, 3, 0, -1, 2, 6, 7, 6, 7, 4, 
5, 1, 6, 2, 4, 0, 6, 7, 4, 7, 
4, 0, 7, 3, 4, 5, 7, 6,                              /* 80 */
4, 5, 0, 1, 7, 6, 0, 3, 7, 3, 
4, 0, 7, 3, 4, 5, 7, 6, -1, 1, 0, 1, 5, 1, 2, 
1, 2, 0, 3, 1, 5, 7, 6, 4, 5, -1, 0, 3, 7, 3, 7, 6, 
7, 6, 4, 5, 3, 2, 4, 0, 3, 0, 
3, 2, 7, 6, 0, 1, 4, 5,                              /* 85 */
1, 0, 1, 5, 1, 2, -1, 7, 6, 4, 5, 3, 2, 4, 0, 3, 0, 
6, 7, 5, 4, 2, 3, 5, 1, 2, 1, 
4, 0, 7, 3, 4, 5, 7, 6, -1, 2, 3, 2, 1, 2, 6, 
2, 3, 2, 1, 2, 6, -1, 4, 5, 0, 1, 7, 6, 0, 3, 7, 3, 
1, 5, 2, 6, 1, 0, 2, 3, -1, 4, 0, 7, 3, 4, 5, 7, 6,  /* 90 */
2, 6, 2, 3, 1, 5, 0, 3, 4, 5, 3, 7, 6, 7,
2, 1, 2, 6, 3, 0, 7, 6, 4, 5, -1, 3, 0, 4, 5, 4, 0, 
5, 4, 1, 0, 6, 7, 1, 2, 6, 2, 
4, 5, 0, 4, 6, 7, 0, 3, 2, 6, 0, 1, 1, 5,
5, 1, 6, 2, 5, 4, 6, 7,                              /* 95 */
5, 4, 5, 6, 5, 1, -1, 7, 6, 7, 4, 7, 3, 
0, 1, 0, 3, 0, 4, -1, 5, 4, 5, 6, 5, 1, -1, 7, 6, 7, 4, 7, 3, 
1, 2, 1, 0, 5, 6, 5, 4, -1, 7, 6, 7, 4, 7, 3, 
7, 6, 7, 4, 7, 3, -1, 1, 2, 0, 3, 5, 6, 0, 4, 5, 4, 
3, 0, 3, 2, 7, 4, 7, 6, -1, 5, 4, 5, 6, 5, 1,        /* 100 */
5, 4, 5, 6, 5, 1, -1, 3, 2, 7, 6, 0, 1, 7, 4, 0, 4, 
1, 2, 1, 0, 5, 6, 5, 4, -1, 3, 0, 3, 2, 7, 4, 7, 6, 
6, 7, 4, 7, 2, 3, 0, 4, 1, 2, 4, 5, 5, 6,
2, 3, 2, 1, 2, 6, -1, 5, 4, 5, 6, 5, 1, -1, 7, 6, 7, 4, 7, 3, 
0, 1, 0, 3, 0, 4, -1, 2, 3, 2, 1, 2, 6, -1, 5, 4, 5, 6, 5, 1, -1, 7, 6, 7, 4, 7, 3,  /* 105 */
7, 6, 7, 4, 7, 3, -1, 1, 0, 5, 4, 2, 3, 5, 6, 2, 6, 
3, 0, 4, 0, 3, 2, 4, 5, 3, 7, 4, 7, -1, 6, 7, 6, 5, 6, 2, 
5, 4, 5, 6, 5, 1, -1, 3, 0, 2, 1, 7, 4, 2, 6, 7, 6, 
1, 5, 4, 5, 1, 2, 4, 7, 0, 1, 0, 4, -1, 6, 7, 6, 5, 6, 2,
4, 5, 6, 5, 0, 1, 6, 2, 0, 3, 6, 7, 4, 7,            /* 110 */
4, 5, 6, 5, 4, 0, 6, 2, 4, 7, 6, 7, 
4, 0, 7, 3, 5, 1, 7, 6, 5, 6, 
0, 3, 7, 3, 0, 1, 7, 6, 5, 1, 5, 6, 
1, 2, 1, 0, 5, 6, 4, 0, 7, 3, -1, 5, 6, 7, 3, 7, 6, 
2, 1, 3, 0, 6, 5, 3, 7, 6, 7,                        /* 115 */
3, 2, 7, 6, 3, 0, 5, 1, 4, 0, -1, 7, 6, 5, 6, 5, 1, 
2, 3, 6, 7, 1, 0, 6, 5, 1, 5, 
1, 2, 0, 1, 5, 6, 0, 4, 6, 7, 0, 3, 2, 3,
2, 1, 2, 3, 6, 5, 6, 7, 
2, 3, 2, 1, 2, 6, -1, 4, 0, 7, 3, 5, 1, 7, 6, 5, 6,  /* 120 */
1, 0, 3, 0, 1, 5, 3, 7, 1, 2, 3, 2, -1, 6, 7, 6, 5, 6, 2, 
3, 7, 6, 7, 0, 4, 5, 6, 0, 1, 2, 6, 2, 3,
2, 6, 2, 3, 5, 6, 0, 3, 6, 7, 3, 7,
1, 2, 2, 6, 0, 3, 6, 7, 0, 4, 5, 6, 1, 5,
1, 2, 6, 2, 1, 0, 6, 7, 1, 5, 6, 5,                  /* 125 */
0, 1, 0, 4, 0, 3, -1, 6, 7, 6, 5, 6, 2, 
6, 7, 6, 5, 6, 2, 
6, 7, 6, 2, 6, 5, 
0, 1, 0, 3, 0, 4, -1, 6, 7, 6, 2, 6, 5, 
1, 0, 1, 5, 1, 2, -1, 6, 7, 6, 2, 6, 5,              /* 130 */
0, 4, 1, 5, 0, 3, 1, 2, -1, 6, 7, 6, 2, 6, 5, 
3, 2, 3, 7, 3, 0, -1, 6, 7, 6, 2, 6, 5, 
0, 4, 0, 1, 3, 7, 3, 2, -1, 6, 7, 6, 2, 6, 5, 
1, 0, 1, 5, 1, 2, -1, 3, 2, 3, 7, 3, 0, -1, 6, 7, 6, 2, 6, 5, 
6, 7, 6, 2, 6, 5, -1, 0, 4, 1, 5, 3, 7, 1, 2, 3, 2,  /* 135 */
2, 1, 6, 5, 2, 3, 6, 7, 
2, 1, 6, 5, 2, 3, 6, 7, -1, 0, 1, 0, 3, 0, 4, 
2, 3, 1, 0, 6, 7, 1, 5, 6, 5, 
0, 4, 1, 5, 0, 3, 6, 7, 2, 3, -1, 1, 5, 6, 5, 6, 7, 
2, 1, 6, 5, 3, 0, 6, 7, 3, 7,                        /* 140 */
0, 4, 0, 1, 3, 7, 2, 1, 6, 5, -1, 3, 7, 6, 5, 6, 7, 
1, 5, 6, 5, 1, 0, 6, 7, 3, 0, 3, 7, 
4, 0, 5, 1, 7, 3, 5, 6, 7, 6, 
4, 5, 4, 0, 4, 7, -1, 6, 7, 6, 2, 6, 5, 
0, 3, 4, 7, 0, 1, 4, 5, -1, 6, 7, 6, 2, 6, 5,        /* 145 */
1, 0, 1, 5, 1, 2, -1, 4, 5, 4, 0, 4, 7, -1, 6, 7, 6, 2, 6, 5, 
6, 7, 6, 2, 6, 5, -1, 0, 3, 4, 7, 1, 2, 4, 5, 1, 5, 
3, 2, 3, 7, 3, 0, -1, 4, 5, 4, 0, 4, 7, -1, 6, 7, 6, 2, 6, 5, 
6, 7, 6, 2, 6, 5, -1, 0, 1, 3, 2, 4, 5, 3, 7, 4, 7, 
1, 0, 1, 5, 1, 2, -1, 3, 2, 3, 7, 3, 0, -1, 4, 5, 4, 0, 4, 7, -1, 6, 7, 6, 2, 6, 5,  /* 150 */
1, 5, 1, 2, 4, 5, 2, 3, 5, 6, 2, 6, -1, 7, 6, 7, 3, 7, 4, 
2, 1, 6, 5, 2, 3, 6, 7, -1, 4, 5, 4, 0, 4, 7, 
0, 3, 4, 7, 0, 1, 4, 5, -1, 2, 1, 6, 5, 2, 3, 6, 7, 
4, 5, 4, 0, 4, 7, -1, 2, 3, 1, 0, 6, 7, 1, 5, 6, 5, 
4, 7, 4, 5, 0, 3, 1, 5, 2, 3, 5, 6, 6, 7,            /* 155 */
4, 5, 4, 0, 4, 7, -1, 2, 1, 6, 5, 3, 0, 6, 7, 3, 7, 
5, 6, 6, 7, 1, 2, 3, 7, 0, 1, 4, 7, 4, 5,
0, 1, 5, 1, 0, 3, 5, 6, 0, 4, 5, 4, -1, 7, 6, 7, 3, 7, 4, 
5, 6, 6, 7, 1, 5, 3, 7, 4, 5, 4, 7,
5, 1, 5, 4, 6, 2, 6, 7,                              /* 160 */
5, 1, 5, 4, 6, 2, 6, 7, -1, 0, 1, 0, 3, 0, 4, 
5, 4, 6, 7, 1, 0, 6, 2, 1, 2, 
0, 3, 0, 4, 1, 2, 5, 4, 6, 7, -1, 1, 2, 6, 7, 6, 2, 
5, 1, 5, 4, 6, 2, 6, 7, -1, 3, 2, 3, 7, 3, 0, 
0, 4, 0, 1, 3, 7, 3, 2, -1, 5, 1, 5, 4, 6, 2, 6, 7,  /* 165 */
3, 2, 3, 7, 3, 0, -1, 5, 4, 6, 7, 1, 0, 6, 2, 1, 2, 
6, 7, 2, 6, 4, 5, 1, 2, 0, 4, 2, 3, 3, 7,
6, 7, 2, 3, 5, 4, 2, 1, 5, 1, 
0, 1, 0, 3, 0, 4, -1, 6, 7, 2, 3, 5, 4, 2, 1, 5, 1, 
2, 3, 1, 0, 6, 7, 5, 4,                              /* 170 */
7, 6, 3, 2, 4, 5, 3, 0, 4, 0, 
3, 0, 2, 1, 3, 7, 5, 4, 6, 7, -1, 2, 1, 5, 1, 5, 4, 
0, 4, 0, 1, 3, 7, 1, 2, 6, 7, 1, 5, 4, 5,
4, 5, 7, 6, 0, 1, 7, 3, 0, 3, 
4, 0, 4, 5, 7, 3, 7, 6,                              /* 175 */
5, 1, 4, 0, 6, 2, 4, 7, 6, 7, 
0, 3, 4, 7, 0, 1, 6, 2, 5, 1, -1, 4, 7, 6, 7, 6, 2, 
1, 2, 1, 0, 6, 2, 4, 0, 6, 7, 4, 7, 
3, 0, 7, 4, 2, 1, 7, 6, 2, 6, 
3, 2, 3, 7, 3, 0, -1, 5, 1, 4, 0, 6, 2, 4, 7, 6, 7,  /* 180 */
2, 3, 3, 7, 0, 1, 4, 7, 1, 5, 6, 7, 2, 6,
0, 3, 2, 3, 0, 4, 2, 6, 0, 1, 1, 2, -1, 7, 6, 7, 3, 7, 4,
2, 3, 7, 3, 2, 1, 7, 4, 2, 6, 7, 6, 
2, 3, 2, 1, 6, 7, 5, 1, 4, 0, -1, 6, 7, 4, 0, 4, 7, 
2, 3, 1, 2, 6, 7, 1, 5, 4, 7, 0, 1, 0, 3,            /* 185 */
3, 2, 0, 1, 7, 6, 0, 4, 7, 4, 
3, 0, 7, 4, 3, 2, 7, 6, 
0, 4, 4, 7, 1, 5, 6, 7, 1, 2, 3, 7, 0, 3,
1, 0, 1, 2, 1, 5, -1, 7, 6, 7, 3, 7, 4, 
3, 7, 0, 3, 6, 7, 0, 1, 4, 7, 0, 4,                  /* 190 */
7, 6, 7, 3, 7, 4, 
6, 2, 6, 5, 7, 3, 7, 4, 
6, 2, 6, 5, 7, 3, 7, 4, -1, 0, 1, 0, 3, 0, 4, 
6, 2, 6, 5, 7, 3, 7, 4, -1, 1, 0, 1, 5, 1, 2, 
0, 4, 1, 5, 0, 3, 1, 2, -1, 6, 2, 6, 5, 7, 3, 7, 4,  /* 195 */
7, 4, 3, 0, 6, 5, 3, 2, 6, 2, 
0, 1, 3, 2, 0, 4, 6, 5, 7, 4, -1, 3, 2, 6, 2, 6, 5, 
1, 0, 1, 5, 1, 2, -1, 7, 4, 3, 0, 6, 5, 3, 2, 6, 2, 
1, 5, 1, 2, 0, 4, 2, 3, 4, 7, 2, 6, 5, 6,
6, 5, 7, 4, 2, 1, 7, 3, 2, 3,                        /* 200 */
0, 1, 0, 3, 0, 4, -1, 6, 5, 7, 4, 2, 1, 7, 3, 2, 3, 
1, 0, 1, 5, 2, 3, 6, 5, 7, 4, -1, 2, 3, 7, 4, 7, 3, 
4, 7, 3, 7, 5, 6, 2, 3, 1, 5, 0, 3, 0, 4,
3, 0, 2, 1, 7, 4, 6, 5, 
5, 6, 4, 7, 1, 2, 4, 0, 1, 0,                        /* 205 */
4, 7, 0, 3, 5, 6, 0, 1, 5, 1, 
4, 0, 5, 1, 4, 7, 5, 6, 
7, 3, 6, 2, 4, 0, 6, 5, 4, 5, 
0, 1, 0, 3, 4, 5, 7, 3, 6, 2, -1, 4, 5, 6, 2, 6, 5, 
1, 0, 1, 5, 1, 2, -1, 7, 3, 6, 2, 4, 0, 6, 5, 4, 5,  /* 210 */
2, 6, 5, 6, 3, 7, 4, 5, 0, 3, 1, 5, 1, 2,
3, 2, 6, 2, 3, 0, 6, 5, 4, 0, 4, 5, 
1, 0, 2, 3, 5, 4, 2, 6, 5, 6, 
1, 2, 0, 1, 2, 6, 0, 4, 2, 3, 0, 3, -1, 5, 4, 5, 1, 5, 6,
2, 6, 5, 6, 2, 3, 4, 5, 1, 2, 1, 5,                  /* 215 */
2, 1, 6, 5, 2, 3, 4, 0, 7, 3,-1, 6, 5, 4, 5, 4, 0,
0, 1, 0, 3, 4, 5, 3, 7, 5, 6, 2, 3, 1, 2,
0, 1, 1, 5, 2, 3, 5, 6, 3, 7, 4, 5, 0, 4,
3, 2, 3, 0, 3, 7, -1, 5, 4, 5, 1, 5, 6, 
1, 2, 5, 6, 0, 3, 5, 4, 0, 4,                        /* 220 */
1, 2, 5, 6, 1, 0, 5, 4, 
0, 1, 5, 1, 0, 3, 5, 6, 0, 4, 5, 4, 
5, 4, 5, 1, 5, 6, 
6, 2, 5, 1, 7, 3, 5, 4, 7, 4, 
0, 1, 0, 3, 0, 4, -1, 6, 2, 5, 1, 7, 3, 5, 4, 7, 4,  /* 225 */
1, 0, 5, 4, 1, 2, 7, 3, 6, 2, -1, 5, 4, 7, 4, 7, 3, 
0, 3, 0, 4, 1, 2, 4, 5, 2, 6, 4, 7, 3, 7,
3, 0, 3, 2, 7, 4, 6, 2, 5, 1, -1, 7, 4, 5, 1, 5, 4, 
1, 5, 4, 5, 2, 6, 4, 7, 2, 3, 0, 4, 0, 1,
0, 3, 2, 3, 4, 7, 2, 6, 4, 5, 1, 2, 0, 1,            /* 230 */
2, 3, 2, 6, 2, 1, -1, 4, 5, 4, 7, 4, 0, 
2, 3, 2, 1, 7, 3, 5, 1, 7, 4, 5, 4, 
1, 0, 3, 0, 1, 5, 3, 7, 1, 2, 3, 2, -1, 4, 5, 4, 7, 4, 0, 
0, 1, 4, 5, 3, 2, 4, 7, 3, 7, 
3, 0, 4, 0, 3, 2, 4, 5, 3, 7, 4, 7,                  /* 235 */
0, 3, 1, 2, 4, 7, 1, 5, 4, 5, 
0, 4, 0, 1, 4, 7, 1, 2, 4, 5, 1, 5,
0, 3, 0, 1, 4, 7, 4, 5, 
4, 5, 4, 7, 4, 0, 
5, 1, 4, 0, 6, 2, 7, 3,                              /* 240 */
2, 6, 1, 5, 3, 7, 1, 0, 3, 0, 
3, 7, 2, 6, 0, 4, 2, 1, 0, 1, 
2, 6, 2, 1, 3, 7, 3, 0, 
1, 5, 0, 4, 2, 6, 0, 3, 2, 3, 
1, 5, 1, 0, 2, 6, 2, 3,                              /* 245 */
1, 2, 0, 1, 2, 6, 0, 4, 2, 3, 0, 3,
2, 3, 2, 6, 2, 1, 
0, 4, 3, 7, 1, 5, 3, 2, 1, 2, 
1, 0, 3, 0, 1, 5, 3, 7, 1, 2, 3, 2, 
0, 4, 3, 7, 0, 1, 3, 2,                              /* 250 */
3, 2, 3, 0, 3, 7, 
0, 4, 0, 3, 1, 5, 1, 2, 
1, 0, 1, 2, 1, 5, 
0, 1, 0, 4, 0, 3};
static int count[] = {0,
    0,    6,   12,   20,   26,   34,   47,   57,   63,   76,   84,   94,  102,
  112,  122,  130,  136,  144,  157,  167,  180,  190,  210,  222,  235,  250,
  265,  282,  297,  314,  331,  341,  347,  360,  368,  378,  391,  406,  421,
  438,  451,  471,  481,  493,  508,  525,  542,  552,  560,  570,  580,  588,
  603,  620,  637,  647,  662,  679,  696,  706,  723,  737,  751,  759,  765,
  778,  791,  806,  814,  824,  839,  856,  869,  889,  904,  921,  931,  943,
  960,  970,  978,  988, 1003, 1020, 1030, 1038, 1055, 1065, 1080, 1097, 1114,
 1128, 1145, 1155, 1169, 1177, 1190, 1210, 1225, 1242, 1257, 1274, 1291, 1305,
 1325, 1352, 1369, 1388, 1405, 1424, 1438, 1450, 1460, 1472, 1489, 1499, 1516,
 1526, 1540, 1548, 1565, 1584, 1598, 1610, 1624, 1636, 1649, 1655, 1661, 1674,
 1687, 1702, 1715, 1730, 1750, 1767, 1775, 1790, 1800, 1817, 1827, 1844, 1856,
 1866, 1879, 1894, 1914, 1931, 1951, 1968, 1995, 2014, 2029, 2046, 2063, 2077,
 2094, 2108, 2127, 2139, 2147, 2162, 2172, 2189, 2204, 2221, 2238, 2252, 2262,
 2279, 2287, 2297, 2314, 2328, 2338, 2346, 2356, 2373, 2385, 2395, 2412, 2426,
 2445, 2457, 2474, 2488, 2498, 2506, 2520, 2533, 2545, 2551, 2559, 2574, 2589,
 2606, 2616, 2633, 2650, 2664, 2674, 2691, 2708, 2722, 2730, 2740, 2750, 2758,
 2768, 2785, 2802, 2816, 2828, 2838, 2857, 2869, 2886, 2900, 2914, 2927, 2937,
 2945, 2957, 2963, 2973, 2990, 3007, 3021, 3038, 3052, 3066, 3079, 3091, 3110,
 3120, 3132, 3142, 3154, 3162, 3168, 3176, 3186, 3196, 3204, 3214, 3222, 3234,
 3240, 3250, 3262, 3270, 3276, 3284, 3290, 3296};

static double cube[8], level;
static int interpol;

void marchingCube(struct GRIDOBJECT *grid, double surface, int primitive,
                  void (*begin)(int), void (*end)(void),
                  void (*vert)(double, double, double),
                  void (*norm)(double, double, double), int ip, int debug)
{
  register double x, y, z;
  int i, j, k, l, set, cc;
  register int ii, jj, kk;

  level=surface;
  interpol=ip;

/*printf("%f %f %f\n", grid->start[0], grid->start[1], grid->start[2]);
  printf("%f %f %f\n", grid->step[0], grid->step[1], grid->step[2]);*/
  for (k=0; k<grid->npoints[2]; k++)
  {
    kk=k % grid->ngridpoints[2];
    for (j=0; j<grid->npoints[1]; j++)
    {
      jj=j % grid->ngridpoints[1];
      for (i=0; i<grid->npoints[0]; i++)
      {
        ii=i % grid->ngridpoints[0];
        set=0;
        if (grid->grid[index(ii,jj,kk)]       < level) set|=0x01;
        if (grid->grid[index(ii+1,jj,kk)]     < level) set|=0x02;
        if (grid->grid[index(ii,jj+1,kk)]     < level) set|=0x04;
        if (grid->grid[index(ii+1,jj+1,kk)]   < level) set|=0x08;
        if (grid->grid[index(ii,jj,kk+1)]     < level) set|=0x10;
        if (grid->grid[index(ii+1,jj,kk+1)]   < level) set|=0x20;
        if (grid->grid[index(ii,jj+1,kk+1)]   < level) set|=0x40;
        if (grid->grid[index(ii+1,jj+1,kk+1)] < level) set|=0x80;

        if (set == 0 || set == 255) continue;

        cube[0]=grid->grid[index(ii,jj,kk)];
        cube[1]=grid->grid[index(ii+1,jj,kk)];
        cube[2]=grid->grid[index(ii+1,jj+1,kk)];
        cube[3]=grid->grid[index(ii,jj+1,kk)];
        cube[4]=grid->grid[index(ii,jj,kk+1)];
        cube[5]=grid->grid[index(ii+1,jj,kk+1)];
        cube[6]=grid->grid[index(ii+1,jj+1,kk+1)];
        cube[7]=grid->grid[index(ii,jj+1,kk+1)];
  
        x=grid->vector1[0]*(double)i*grid->step[0]
         +grid->vector2[0]*(double)j*grid->step[1]
         +grid->vector3[0]*(double)k*grid->step[2]+grid->start[0];
        y=grid->vector1[1]*(double)i*grid->step[0]
         +grid->vector2[1]*(double)j*grid->step[1]
         +grid->vector3[1]*(double)k*grid->step[2]+grid->start[1];
        z=grid->vector1[2]*(double)i*grid->step[0]
         +grid->vector2[2]*(double)j*grid->step[1]
         +grid->vector3[2]*(double)k*grid->step[2]+grid->start[2];
/*
        glColor4f(0.0, 0.0, 0.0, 0.0);
        glBegin(GL_POINTS);
        glVertex3d(x, y, z);
        glEnd();
*/
        if (primitive == GL_TRIANGLE_STRIP)
        {
          (*begin)(primitive);
          for (l=count[set]; l<count[set+1]; l+=2)
          {
            if (corner[l] == -1)
            {
              (*end)();
              (*begin)(primitive);
              l--;
            }
            else
              set_vertex(x, y, z, corner[l], corner[l+1], grid, ii, jj, kk,
                         primitive, vert, norm);
          }
          (*end)();
        }
        else
        {
          cc=0;
          for (l=count[set]; l<count[set+1]; l+=2)
          {
            if (corner[l] == -1)
            {
                cc=0;
              l++;
            }
            if (cc == 0 || cc > 2) (*begin)(primitive);
            set_vertex(x, y, z, corner[l], corner[l+1], grid, ii, jj, kk,
                       primitive, vert, norm);
            cc++;
            if (cc > 3)
            {
              set_vertex(x, y, z, corner[l-2], corner[l-1], grid, ii, jj, kk,
                         primitive, vert, norm);
              set_vertex(x, y, z, corner[l-4], corner[l-3], grid, ii, jj, kk,
                         primitive, vert, norm);
              cc+=2;
            }
            if (cc % 3 == 0) (*end)();
          }
        }
      }
    }
  }
}

void set_vertex(double x, double y, double z, int c1, int c2, struct GRIDOBJECT *grid,
                int ig, int jg, int kg, int primitive, void (*vert)(double, double, double),
                void (*norm)(double, double, double))
{
  static double vertex[8][8][3]={{{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.0, 0.5, 0.0},
                                  {0.0, 0.0, 0.5}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}},
                                 {{-0.5, 0.0, 0.0}, {0.0, 0.0, 0.0},
                                  {1.0, 0.5, 0.0}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {1.0, 0.0, 0.5},
                                  {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}},
                                 {{0.0, 0.0, 0.0}, {1.0, -0.5, 0.0},
                                  {0.0, 0.0, 0.0}, {-0.5, 1.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0},
                                  {1.0, 1.0, 0.5}, {0.0, 0.0, 0.0}},
                                 {{0.0, -0.5, 0.0}, {0.0, 0.0, 0.0},
                                  {0.5, 1.0, 0.0}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.0, 1.0, 0.5}},
                                 {{0.0, 0.0, -0.5}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.5, 0.0, 1.0},
                                  {0.0, 0.0, 0.0}, {0.0, 0.5, 1.0}},
                                 {{0.0, 0.0, 0.0}, {1.0, 0.0, -0.5},
                                  {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0},
                                  {-0.5, 0.0, 1.0}, {0.0, 0.0, 0.0},
                                  {1.0, 0.5, 1.0}, {0.0, 0.0, 0.0}},
                                 {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0},
                                  {1.0, 1.0, -0.5}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {1.0, -0.5, 1.0},
                                  {0.0, 0.0, 0.0}, {-0.5, 1.0, 1.0}},
                                 {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0},
                                  {0.0, 0.0, 0.0}, {0.0, 1.0, -0.5},
                                  {0.0, -0.5, 1.0}, {0.0, 0.0, 0.0},
                                  {0.5, 1.0, 1.0}, {0.0, 0.0, 0.0}},
                                  };
  double v[3], x1=cube[c1], x2=cube[c2], l=level, n1[3], n2[3];
  int ii, jj, kk;
  register int i;
 
  switch (interpol)
  {
    case IP_NONE:   v[0]=fabs(vertex[c1][c2][0]);
                    v[1]=fabs(vertex[c1][c2][1]);
                    v[2]=fabs(vertex[c1][c2][2]);
                    break;
    case IP_LOG:    if (x1 > 0.0) x1=log(x1);
                    if (x2 > 0.0) x2=log(x2);
                    if (l  > 0.0) l=log(l);
    case IP_LINEAR: x1-=l;
                    x2-=l;
                    for (i=0; i<3; i++)
                    {
                      v[i]=vertex[c1][c2][i];
                      if (v[i] ==  0.5) v[i]=fabs(x1)/(fabs(x1)+fabs(x2));
                      if (v[i] == -0.5) v[i]=1.0-fabs(x1)/(fabs(x1)+fabs(x2));
                    }
                    break;
  }
  if (primitive == GL_TRIANGLE_STRIP)
  {
    ii=ig;
    jj=jg;
    kk=kg;
    switch (c1)
    {
      case 1: ii++;
              break;
      case 2: ii++;
              jj++;
              break;
      case 3: jj++;
              break;
      case 4: kk++;
              break;
      case 5: ii++;
              kk++;
              break;
      case 6: ii++;
              jj++;
              kk++;
              break;
      case 7: jj++;
              kk++;
              break;
    }
    if (ii >= grid->ngridpoints[0]-1)
      n1[0]=(grid->grid[index(ii-1, jj, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[0];
    else if (ii == 0)
      n1[0]=(grid->grid[index(ii+1, jj, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[0];
    else
      n1[0]=(grid->grid[index(ii-1, jj, kk)]-grid->grid[index(ii+1, jj, kk)])/(2.0*grid->step[0]);
    if (jj >= grid->ngridpoints[1]-1)
      n1[1]=(grid->grid[index(ii, jj-1, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[1];
    else if (jj == 0)
      n1[1]=(grid->grid[index(ii, jj+1, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[1];
    else
      n1[1]=(grid->grid[index(ii, jj-1, kk)]-grid->grid[index(ii, jj+1, kk)])/(2.0*grid->step[1]);
    if (kk >= grid->ngridpoints[2]-1)
      n1[2]=(grid->grid[index(ii, jj, kk-1)]-grid->grid[index(ii, jj, kk)])/grid->step[2];
    else if (kk == 0)
      n1[2]=(grid->grid[index(ii, jj, kk+1)]-grid->grid[index(ii, jj, kk)])/grid->step[2];
    else
      n1[2]=(grid->grid[index(ii, jj, kk-1)]-grid->grid[index(ii, jj, kk+1)])/(2.0*grid->step[2]);
    ii=ig;
    jj=jg;
    kk=kg;
    switch (c2)
    {
      case 1: ii++;
              break;
      case 2: ii++;
              jj++;
              break;
      case 3: jj++;
              break;
      case 4: kk++;
              break;
      case 5: ii++;
              kk++;
              break;
      case 6: ii++;
              jj++;
              kk++;
              break;
      case 7: jj++;
              kk++;
              break;
    }
    if (ii >= grid->ngridpoints[0]-1)
      n2[0]=(grid->grid[index(ii-1, jj, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[0];
    else if (ii == 0)
      n2[0]=(grid->grid[index(ii+1, jj, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[0];
    else
      n2[0]=(grid->grid[index(ii-1, jj, kk)]-grid->grid[index(ii+1, jj, kk)])/(2.0*grid->step[0]);
    if (jj >= grid->ngridpoints[1]-1)
      n2[1]=(grid->grid[index(ii, jj-1, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[1];
    else if (jj == 0)
      n2[1]=(grid->grid[index(ii, jj+1, kk)]-grid->grid[index(ii, jj, kk)])/grid->step[1];
    else
      n2[1]=(grid->grid[index(ii, jj-1, kk)]-grid->grid[index(ii, jj+1, kk)])/(2.0*grid->step[1]);
    if (kk >= grid->ngridpoints[2]-1)
      n2[2]=(grid->grid[index(ii, jj, kk-1)]-grid->grid[index(ii, jj, kk)])/grid->step[2];
    else if (kk == 0)
      n2[2]=(grid->grid[index(ii, jj, kk+1)]-grid->grid[index(ii, jj, kk)])/grid->step[2];
    else
      n2[2]=(grid->grid[index(ii, jj, kk-1)]-grid->grid[index(ii, jj, kk+1)])/(2.0*grid->step[2]);
    n1[0]=0.5*(n1[0]+n2[0]);
    n1[1]=0.5*(n1[1]+n2[1]);
    n1[2]=0.5*(n1[2]+n2[2]);
    if (level > 0.0)
      l=sqrt(n1[0]*n1[0]+n1[1]*n1[1]+n1[2]*n1[2]);
    else
      l=-sqrt(n1[0]*n1[0]+n1[1]*n1[1]+n1[2]*n1[2]);
    n1[0]/=l;
    n1[1]/=l;
    n1[2]/=l;
    (*norm)(n1[0], n1[1], n1[2]);
  }
/*recordit=TRUE; */
  v[0]=x+v[0]*grid->vector1[0]*grid->step[0]
      +v[1]*grid->vector2[0]*grid->step[1]
      +v[2]*grid->vector3[0]*grid->step[2];
  v[1]=y+v[0]*grid->vector1[1]*grid->step[0]
      +v[1]*grid->vector2[1]*grid->step[1]
      +v[2]*grid->vector3[1]*grid->step[2];
  v[2]=z+v[0]*grid->vector1[2]*grid->step[0]
      +v[1]*grid->vector2[2]*grid->step[1]
      +v[2]*grid->vector3[2]*grid->step[2];
  (*vert)(v[0], v[1], v[2]);
/*printf("%f %f %f %f %f %f\n", v[0], v[1], v[2], n1[0], n1[1], n1[2]); */
}
/*
void testNorm(double x, double y, double z)
{
  norms[vcount][0]=0.2*x;
  norms[vcount][1]=0.2*y;
  norms[vcount][2]=0.2*z;
  glNormal3d(x, y, z);
}

void testVertex(double x, double y, double z)
{
  if (recordit)
  {
    vertices[vcount][0]=x;
    vertices[vcount][1]=y;
    vertices[vcount][2]=z;
    vcount++;
  }
  glVertex3d(x, y, z);
}

void testBegin(GLenum what)
{
  vcount=0;
  glBegin(what);
}

void testEnd(void)
{
  int i;

  glEnd();
  if (recordit)
  {
    for (i=0; i<vcount; i++)
    {
      glColor4f(0.0, 0.0, 0.0, 0.5);
      glBegin(GL_LINES);
      glVertex3d(vertices[i][0], vertices[i][1], vertices[i][2]);
      glVertex3d(vertices[i][0]+norms[i][0], vertices[i][1]+norms[i][1], vertices[i][2]+norms[i][2]);
      glEnd();
      if (level > 0.0)
        glColor4f(1.0, 0.0, 0.0, 0.0);
      else
        glColor4f(0.0, 0.0, 1.0, 0.0);
    }
    recordit=FALSE;
  }
} */

Generated by  Doxygen 1.6.0   Back to index