Pileus Git - grits/blob - src/marching.c

   1 /**
   2  * Copyright (c) 2002 Jamie Zawinski <jwz@jwz.org>
   3  * Copyright (c) 2009 Andy Spencer <andy753421@gmail.com>
   4  *
   5  * Utility functions to create "marching cubes" meshes from 3d fields.
   6  *
   7  * Permission to use, copy, modify, distribute, and sell this software and its
   8  * documentation for any purpose is hereby granted without fee, provided that
   9  * the above copyright notice appear in all copies and that both that
  10  * copyright notice and this permission notice appear in supporting
  11  * documentation.  No representations are made about the suitability of this
  12  * software for any purpose.  It is provided "as is" without express or
  13  * implied warranty.
  14  *
  15  * Marching cubes implementation by Paul Bourke <pbourke@swin.edu.au>
  16  * http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/
  17  */
  18
  19 #include <math.h>
  20 #include <glib.h>
  21
  22 #include <GL/gl.h>
  23
  24 #include "marching.h"
  25
  26 /* Calculate the unit normal at p given two other points p1,p2 on the
  27  * surface. The normal points in the direction of p1 crossproduct p2
  28  */
  29 XYZ calc_normal (XYZ p, XYZ p1, XYZ p2)
  30 {
  31         XYZ n, pa, pb;
  32         pa.x = p1.x - p.x;
  33         pa.y = p1.y - p.y;
  34         pa.z = p1.z - p.z;
  35         pb.x = p2.x - p.x;
  36         pb.y = p2.y - p.y;
  37         pb.z = p2.z - p.z;
  38         n.x = pa.y * pb.z - pa.z * pb.y;
  39         n.y = pa.z * pb.x - pa.x * pb.z;
  40         n.z = pa.x * pb.y - pa.y * pb.x;
  41         return (n);
  42 }
  43
  44 /* Call glNormal3f() with a normal of the indicated triangle.
  45  */
  46 void do_normal(GLfloat x1, GLfloat y1, GLfloat z1,
  47                 GLfloat x2, GLfloat y2, GLfloat z2,
  48                 GLfloat x3, GLfloat y3, GLfloat z3)
  49 {
  50         XYZ p1, p2, p3, p;
  51         p1.x = x1; p1.y = y1; p1.z = z1;
  52         p2.x = x2; p2.y = y2; p2.z = z2;
  53         p3.x = x3; p3.y = y3; p3.z = z3;
  54         p = calc_normal (p1, p2, p3);
  55         glNormal3f (p.x, p.y, p.z);
  56 }
  57
  58 /* Indexing convention:
  59  *
  60  *           Vertices:                    Edges:
  61  *
  62  *        4  ______________ 5           ______________
  63  *         /|            /|           /|     4      /|
  64  *        / |         6 / |       7  / |8        5 / |
  65  *    7  /_____________/  |        /______________/  | 9
  66  *      |   |         |   |        |   |   6     |   |
  67  *      | 0 |_________|___| 1      |   |_________|10_|
  68  *      |  /          |  /      11 | 3/     0    |  /
  69  *      | /           | /          | /           | / 1
  70  *    3 |/____________|/ 2         |/____________|/
  71  *                                        2
  72  */
  73
  74 static const int edgeTable[256] = {
  75         0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
  76         0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
  77         0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
  78         0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
  79         0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
  80         0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
  81         0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
  82         0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
  83         0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
  84         0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
  85         0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
  86         0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
  87         0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
  88         0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
  89         0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
  90         0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
  91         0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
  92         0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
  93         0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
  94         0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
  95         0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
  96         0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
  97         0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
  98         0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
  99         0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
 100         0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
 101         0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
 102         0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
 103         0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
 104         0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
 105         0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
 106         0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0
 107 };
 108
 109 static const int triTable[256][16] = {
 110         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 111         { 0,  8,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 112         { 0,  1,  9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 113         { 1,  8,  3,  9,  8,  1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 114         { 1,  2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 115         { 0,  8,  3,  1,  2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 116         { 9,  2, 10,  0,  2,  9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 117         { 2,  8,  3,  2, 10,  8, 10,  9,  8, -1, -1, -1, -1, -1, -1, -1},
 118         { 3, 11,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 119         { 0, 11,  2,  8, 11,  0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 120         { 1,  9,  0,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 121         { 1, 11,  2,  1,  9, 11,  9,  8, 11, -1, -1, -1, -1, -1, -1, -1},
 122         { 3, 10,  1, 11, 10,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 123         { 0, 10,  1,  0,  8, 10,  8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
 124         { 3,  9,  0,  3, 11,  9, 11, 10,  9, -1, -1, -1, -1, -1, -1, -1},
 125         { 9,  8, 10, 10,  8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 126         { 4,  7,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 127         { 4,  3,  0,  7,  3,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 128         { 0,  1,  9,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 129         { 4,  1,  9,  4,  7,  1,  7,  3,  1, -1, -1, -1, -1, -1, -1, -1},
 130         { 1,  2, 10,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 131         { 3,  4,  7,  3,  0,  4,  1,  2, 10, -1, -1, -1, -1, -1, -1, -1},
 132         { 9,  2, 10,  9,  0,  2,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1},
 133         { 2, 10,  9,  2,  9,  7,  2,  7,  3,  7,  9,  4, -1, -1, -1, -1},
 134         { 8,  4,  7,  3, 11,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 135         {11,  4,  7, 11,  2,  4,  2,  0,  4, -1, -1, -1, -1, -1, -1, -1},
 136         { 9,  0,  1,  8,  4,  7,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1},
 137         { 4,  7, 11,  9,  4, 11,  9, 11,  2,  9,  2,  1, -1, -1, -1, -1},
 138         { 3, 10,  1,  3, 11, 10,  7,  8,  4, -1, -1, -1, -1, -1, -1, -1},
 139         { 1, 11, 10,  1,  4, 11,  1,  0,  4,  7, 11,  4, -1, -1, -1, -1},
 140         { 4,  7,  8,  9,  0, 11,  9, 11, 10, 11,  0,  3, -1, -1, -1, -1},
 141         { 4,  7, 11,  4, 11,  9,  9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
 142         { 9,  5,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 143         { 9,  5,  4,  0,  8,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 144         { 0,  5,  4,  1,  5,  0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 145         { 8,  5,  4,  8,  3,  5,  3,  1,  5, -1, -1, -1, -1, -1, -1, -1},
 146         { 1,  2, 10,  9,  5,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 147         { 3,  0,  8,  1,  2, 10,  4,  9,  5, -1, -1, -1, -1, -1, -1, -1},
 148         { 5,  2, 10,  5,  4,  2,  4,  0,  2, -1, -1, -1, -1, -1, -1, -1},
 149         { 2, 10,  5,  3,  2,  5,  3,  5,  4,  3,  4,  8, -1, -1, -1, -1},
 150         { 9,  5,  4,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 151         { 0, 11,  2,  0,  8, 11,  4,  9,  5, -1, -1, -1, -1, -1, -1, -1},
 152         { 0,  5,  4,  0,  1,  5,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1},
 153         { 2,  1,  5,  2,  5,  8,  2,  8, 11,  4,  8,  5, -1, -1, -1, -1},
 154         {10,  3, 11, 10,  1,  3,  9,  5,  4, -1, -1, -1, -1, -1, -1, -1},
 155         { 4,  9,  5,  0,  8,  1,  8, 10,  1,  8, 11, 10, -1, -1, -1, -1},
 156         { 5,  4,  0,  5,  0, 11,  5, 11, 10, 11,  0,  3, -1, -1, -1, -1},
 157         { 5,  4,  8,  5,  8, 10, 10,  8, 11, -1, -1, -1, -1, -1, -1, -1},
 158         { 9,  7,  8,  5,  7,  9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 159         { 9,  3,  0,  9,  5,  3,  5,  7,  3, -1, -1, -1, -1, -1, -1, -1},
 160         { 0,  7,  8,  0,  1,  7,  1,  5,  7, -1, -1, -1, -1, -1, -1, -1},
 161         { 1,  5,  3,  3,  5,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 162         { 9,  7,  8,  9,  5,  7, 10,  1,  2, -1, -1, -1, -1, -1, -1, -1},
 163         {10,  1,  2,  9,  5,  0,  5,  3,  0,  5,  7,  3, -1, -1, -1, -1},
 164         { 8,  0,  2,  8,  2,  5,  8,  5,  7, 10,  5,  2, -1, -1, -1, -1},
 165         { 2, 10,  5,  2,  5,  3,  3,  5,  7, -1, -1, -1, -1, -1, -1, -1},
 166         { 7,  9,  5,  7,  8,  9,  3, 11,  2, -1, -1, -1, -1, -1, -1, -1},
 167         { 9,  5,  7,  9,  7,  2,  9,  2,  0,  2,  7, 11, -1, -1, -1, -1},
 168         { 2,  3, 11,  0,  1,  8,  1,  7,  8,  1,  5,  7, -1, -1, -1, -1},
 169         {11,  2,  1, 11,  1,  7,  7,  1,  5, -1, -1, -1, -1, -1, -1, -1},
 170         { 9,  5,  8,  8,  5,  7, 10,  1,  3, 10,  3, 11, -1, -1, -1, -1},
 171         { 5,  7,  0,  5,  0,  9,  7, 11,  0,  1,  0, 10, 11, 10,  0, -1},
 172         {11, 10,  0, 11,  0,  3, 10,  5,  0,  8,  0,  7,  5,  7,  0, -1},
 173         {11, 10,  5,  7, 11,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 174         {10,  6,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 175         { 0,  8,  3,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 176         { 9,  0,  1,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 177         { 1,  8,  3,  1,  9,  8,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1},
 178         { 1,  6,  5,  2,  6,  1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 179         { 1,  6,  5,  1,  2,  6,  3,  0,  8, -1, -1, -1, -1, -1, -1, -1},
 180         { 9,  6,  5,  9,  0,  6,  0,  2,  6, -1, -1, -1, -1, -1, -1, -1},
 181         { 5,  9,  8,  5,  8,  2,  5,  2,  6,  3,  2,  8, -1, -1, -1, -1},
 182         { 2,  3, 11, 10,  6,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 183         {11,  0,  8, 11,  2,  0, 10,  6,  5, -1, -1, -1, -1, -1, -1, -1},
 184         { 0,  1,  9,  2,  3, 11,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1},
 185         { 5, 10,  6,  1,  9,  2,  9, 11,  2,  9,  8, 11, -1, -1, -1, -1},
 186         { 6,  3, 11,  6,  5,  3,  5,  1,  3, -1, -1, -1, -1, -1, -1, -1},
 187         { 0,  8, 11,  0, 11,  5,  0,  5,  1,  5, 11,  6, -1, -1, -1, -1},
 188         { 3, 11,  6,  0,  3,  6,  0,  6,  5,  0,  5,  9, -1, -1, -1, -1},
 189         { 6,  5,  9,  6,  9, 11, 11,  9,  8, -1, -1, -1, -1, -1, -1, -1},
 190         { 5, 10,  6,  4,  7,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 191         { 4,  3,  0,  4,  7,  3,  6,  5, 10, -1, -1, -1, -1, -1, -1, -1},
 192         { 1,  9,  0,  5, 10,  6,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1},
 193         {10,  6,  5,  1,  9,  7,  1,  7,  3,  7,  9,  4, -1, -1, -1, -1},
 194         { 6,  1,  2,  6,  5,  1,  4,  7,  8, -1, -1, -1, -1, -1, -1, -1},
 195         { 1,  2,  5,  5,  2,  6,  3,  0,  4,  3,  4,  7, -1, -1, -1, -1},
 196         { 8,  4,  7,  9,  0,  5,  0,  6,  5,  0,  2,  6, -1, -1, -1, -1},
 197         { 7,  3,  9,  7,  9,  4,  3,  2,  9,  5,  9,  6,  2,  6,  9, -1},
 198         { 3, 11,  2,  7,  8,  4, 10,  6,  5, -1, -1, -1, -1, -1, -1, -1},
 199         { 5, 10,  6,  4,  7,  2,  4,  2,  0,  2,  7, 11, -1, -1, -1, -1},
 200         { 0,  1,  9,  4,  7,  8,  2,  3, 11,  5, 10,  6, -1, -1, -1, -1},
 201         { 9,  2,  1,  9, 11,  2,  9,  4, 11,  7, 11,  4,  5, 10,  6, -1},
 202         { 8,  4,  7,  3, 11,  5,  3,  5,  1,  5, 11,  6, -1, -1, -1, -1},
 203         { 5,  1, 11,  5, 11,  6,  1,  0, 11,  7, 11,  4,  0,  4, 11, -1},
 204         { 0,  5,  9,  0,  6,  5,  0,  3,  6, 11,  6,  3,  8,  4,  7, -1},
 205         { 6,  5,  9,  6,  9, 11,  4,  7,  9,  7, 11,  9, -1, -1, -1, -1},
 206         {10,  4,  9,  6,  4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 207         { 4, 10,  6,  4,  9, 10,  0,  8,  3, -1, -1, -1, -1, -1, -1, -1},
 208         {10,  0,  1, 10,  6,  0,  6,  4,  0, -1, -1, -1, -1, -1, -1, -1},
 209         { 8,  3,  1,  8,  1,  6,  8,  6,  4,  6,  1, 10, -1, -1, -1, -1},
 210         { 1,  4,  9,  1,  2,  4,  2,  6,  4, -1, -1, -1, -1, -1, -1, -1},
 211         { 3,  0,  8,  1,  2,  9,  2,  4,  9,  2,  6,  4, -1, -1, -1, -1},
 212         { 0,  2,  4,  4,  2,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 213         { 8,  3,  2,  8,  2,  4,  4,  2,  6, -1, -1, -1, -1, -1, -1, -1},
 214         {10,  4,  9, 10,  6,  4, 11,  2,  3, -1, -1, -1, -1, -1, -1, -1},
 215         { 0,  8,  2,  2,  8, 11,  4,  9, 10,  4, 10,  6, -1, -1, -1, -1},
 216         { 3, 11,  2,  0,  1,  6,  0,  6,  4,  6,  1, 10, -1, -1, -1, -1},
 217         { 6,  4,  1,  6,  1, 10,  4,  8,  1,  2,  1, 11,  8, 11,  1, -1},
 218         { 9,  6,  4,  9,  3,  6,  9,  1,  3, 11,  6,  3, -1, -1, -1, -1},
 219         { 8, 11,  1,  8,  1,  0, 11,  6,  1,  9,  1,  4,  6,  4,  1, -1},
 220         { 3, 11,  6,  3,  6,  0,  0,  6,  4, -1, -1, -1, -1, -1, -1, -1},
 221         { 6,  4,  8, 11,  6,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 222         { 7, 10,  6,  7,  8, 10,  8,  9, 10, -1, -1, -1, -1, -1, -1, -1},
 223         { 0,  7,  3,  0, 10,  7,  0,  9, 10,  6,  7, 10, -1, -1, -1, -1},
 224         {10,  6,  7,  1, 10,  7,  1,  7,  8,  1,  8,  0, -1, -1, -1, -1},
 225         {10,  6,  7, 10,  7,  1,  1,  7,  3, -1, -1, -1, -1, -1, -1, -1},
 226         { 1,  2,  6,  1,  6,  8,  1,  8,  9,  8,  6,  7, -1, -1, -1, -1},
 227         { 2,  6,  9,  2,  9,  1,  6,  7,  9,  0,  9,  3,  7,  3,  9, -1},
 228         { 7,  8,  0,  7,  0,  6,  6,  0,  2, -1, -1, -1, -1, -1, -1, -1},
 229         { 7,  3,  2,  6,  7,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 230         { 2,  3, 11, 10,  6,  8, 10,  8,  9,  8,  6,  7, -1, -1, -1, -1},
 231         { 2,  0,  7,  2,  7, 11,  0,  9,  7,  6,  7, 10,  9, 10,  7, -1},
 232         { 1,  8,  0,  1,  7,  8,  1, 10,  7,  6,  7, 10,  2,  3, 11, -1},
 233         {11,  2,  1, 11,  1,  7, 10,  6,  1,  6,  7,  1, -1, -1, -1, -1},
 234         { 8,  9,  6,  8,  6,  7,  9,  1,  6, 11,  6,  3,  1,  3,  6, -1},
 235         { 0,  9,  1, 11,  6,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 236         { 7,  8,  0,  7,  0,  6,  3, 11,  0, 11,  6,  0, -1, -1, -1, -1},
 237         { 7, 11,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 238         { 7,  6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 239         { 3,  0,  8, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 240         { 0,  1,  9, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 241         { 8,  1,  9,  8,  3,  1, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1},
 242         {10,  1,  2,  6, 11,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 243         { 1,  2, 10,  3,  0,  8,  6, 11,  7, -1, -1, -1, -1, -1, -1, -1},
 244         { 2,  9,  0,  2, 10,  9,  6, 11,  7, -1, -1, -1, -1, -1, -1, -1},
 245         { 6, 11,  7,  2, 10,  3, 10,  8,  3, 10,  9,  8, -1, -1, -1, -1},
 246         { 7,  2,  3,  6,  2,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 247         { 7,  0,  8,  7,  6,  0,  6,  2,  0, -1, -1, -1, -1, -1, -1, -1},
 248         { 2,  7,  6,  2,  3,  7,  0,  1,  9, -1, -1, -1, -1, -1, -1, -1},
 249         { 1,  6,  2,  1,  8,  6,  1,  9,  8,  8,  7,  6, -1, -1, -1, -1},
 250         {10,  7,  6, 10,  1,  7,  1,  3,  7, -1, -1, -1, -1, -1, -1, -1},
 251         {10,  7,  6,  1,  7, 10,  1,  8,  7,  1,  0,  8, -1, -1, -1, -1},
 252         { 0,  3,  7,  0,  7, 10,  0, 10,  9,  6, 10,  7, -1, -1, -1, -1},
 253         { 7,  6, 10,  7, 10,  8,  8, 10,  9, -1, -1, -1, -1, -1, -1, -1},
 254         { 6,  8,  4, 11,  8,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 255         { 3,  6, 11,  3,  0,  6,  0,  4,  6, -1, -1, -1, -1, -1, -1, -1},
 256         { 8,  6, 11,  8,  4,  6,  9,  0,  1, -1, -1, -1, -1, -1, -1, -1},
 257         { 9,  4,  6,  9,  6,  3,  9,  3,  1, 11,  3,  6, -1, -1, -1, -1},
 258         { 6,  8,  4,  6, 11,  8,  2, 10,  1, -1, -1, -1, -1, -1, -1, -1},
 259         { 1,  2, 10,  3,  0, 11,  0,  6, 11,  0,  4,  6, -1, -1, -1, -1},
 260         { 4, 11,  8,  4,  6, 11,  0,  2,  9,  2, 10,  9, -1, -1, -1, -1},
 261         {10,  9,  3, 10,  3,  2,  9,  4,  3, 11,  3,  6,  4,  6,  3, -1},
 262         { 8,  2,  3,  8,  4,  2,  4,  6,  2, -1, -1, -1, -1, -1, -1, -1},
 263         { 0,  4,  2,  4,  6,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 264         { 1,  9,  0,  2,  3,  4,  2,  4,  6,  4,  3,  8, -1, -1, -1, -1},
 265         { 1,  9,  4,  1,  4,  2,  2,  4,  6, -1, -1, -1, -1, -1, -1, -1},
 266         { 8,  1,  3,  8,  6,  1,  8,  4,  6,  6, 10,  1, -1, -1, -1, -1},
 267         {10,  1,  0, 10,  0,  6,  6,  0,  4, -1, -1, -1, -1, -1, -1, -1},
 268         { 4,  6,  3,  4,  3,  8,  6, 10,  3,  0,  3,  9, 10,  9,  3, -1},
 269         {10,  9,  4,  6, 10,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 270         { 4,  9,  5,  7,  6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 271         { 0,  8,  3,  4,  9,  5, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1},
 272         { 5,  0,  1,  5,  4,  0,  7,  6, 11, -1, -1, -1, -1, -1, -1, -1},
 273         {11,  7,  6,  8,  3,  4,  3,  5,  4,  3,  1,  5, -1, -1, -1, -1},
 274         { 9,  5,  4, 10,  1,  2,  7,  6, 11, -1, -1, -1, -1, -1, -1, -1},
 275         { 6, 11,  7,  1,  2, 10,  0,  8,  3,  4,  9,  5, -1, -1, -1, -1},
 276         { 7,  6, 11,  5,  4, 10,  4,  2, 10,  4,  0,  2, -1, -1, -1, -1},
 277         { 3,  4,  8,  3,  5,  4,  3,  2,  5, 10,  5,  2, 11,  7,  6, -1},
 278         { 7,  2,  3,  7,  6,  2,  5,  4,  9, -1, -1, -1, -1, -1, -1, -1},
 279         { 9,  5,  4,  0,  8,  6,  0,  6,  2,  6,  8,  7, -1, -1, -1, -1},
 280         { 3,  6,  2,  3,  7,  6,  1,  5,  0,  5,  4,  0, -1, -1, -1, -1},
 281         { 6,  2,  8,  6,  8,  7,  2,  1,  8,  4,  8,  5,  1,  5,  8, -1},
 282         { 9,  5,  4, 10,  1,  6,  1,  7,  6,  1,  3,  7, -1, -1, -1, -1},
 283         { 1,  6, 10,  1,  7,  6,  1,  0,  7,  8,  7,  0,  9,  5,  4, -1},
 284         { 4,  0, 10,  4, 10,  5,  0,  3, 10,  6, 10,  7,  3,  7, 10, -1},
 285         { 7,  6, 10,  7, 10,  8,  5,  4, 10,  4,  8, 10, -1, -1, -1, -1},
 286         { 6,  9,  5,  6, 11,  9, 11,  8,  9, -1, -1, -1, -1, -1, -1, -1},
 287         { 3,  6, 11,  0,  6,  3,  0,  5,  6,  0,  9,  5, -1, -1, -1, -1},
 288         { 0, 11,  8,  0,  5, 11,  0,  1,  5,  5,  6, 11, -1, -1, -1, -1},
 289         { 6, 11,  3,  6,  3,  5,  5,  3,  1, -1, -1, -1, -1, -1, -1, -1},
 290         { 1,  2, 10,  9,  5, 11,  9, 11,  8, 11,  5,  6, -1, -1, -1, -1},
 291         { 0, 11,  3,  0,  6, 11,  0,  9,  6,  5,  6,  9,  1,  2, 10, -1},
 292         {11,  8,  5, 11,  5,  6,  8,  0,  5, 10,  5,  2,  0,  2,  5, -1},
 293         { 6, 11,  3,  6,  3,  5,  2, 10,  3, 10,  5,  3, -1, -1, -1, -1},
 294         { 5,  8,  9,  5,  2,  8,  5,  6,  2,  3,  8,  2, -1, -1, -1, -1},
 295         { 9,  5,  6,  9,  6,  0,  0,  6,  2, -1, -1, -1, -1, -1, -1, -1},
 296         { 1,  5,  8,  1,  8,  0,  5,  6,  8,  3,  8,  2,  6,  2,  8, -1},
 297         { 1,  5,  6,  2,  1,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 298         { 1,  3,  6,  1,  6, 10,  3,  8,  6,  5,  6,  9,  8,  9,  6, -1},
 299         {10,  1,  0, 10,  0,  6,  9,  5,  0,  5,  6,  0, -1, -1, -1, -1},
 300         { 0,  3,  8,  5,  6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 301         {10,  5,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 302         {11,  5, 10,  7,  5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 303         {11,  5, 10, 11,  7,  5,  8,  3,  0, -1, -1, -1, -1, -1, -1, -1},
 304         { 5, 11,  7,  5, 10, 11,  1,  9,  0, -1, -1, -1, -1, -1, -1, -1},
 305         {10,  7,  5, 10, 11,  7,  9,  8,  1,  8,  3,  1, -1, -1, -1, -1},
 306         {11,  1,  2, 11,  7,  1,  7,  5,  1, -1, -1, -1, -1, -1, -1, -1},
 307         { 0,  8,  3,  1,  2,  7,  1,  7,  5,  7,  2, 11, -1, -1, -1, -1},
 308         { 9,  7,  5,  9,  2,  7,  9,  0,  2,  2, 11,  7, -1, -1, -1, -1},
 309         { 7,  5,  2,  7,  2, 11,  5,  9,  2,  3,  2,  8,  9,  8,  2, -1},
 310         { 2,  5, 10,  2,  3,  5,  3,  7,  5, -1, -1, -1, -1, -1, -1, -1},
 311         { 8,  2,  0,  8,  5,  2,  8,  7,  5, 10,  2,  5, -1, -1, -1, -1},
 312         { 9,  0,  1,  5, 10,  3,  5,  3,  7,  3, 10,  2, -1, -1, -1, -1},
 313         { 9,  8,  2,  9,  2,  1,  8,  7,  2, 10,  2,  5,  7,  5,  2, -1},
 314         { 1,  3,  5,  3,  7,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 315         { 0,  8,  7,  0,  7,  1,  1,  7,  5, -1, -1, -1, -1, -1, -1, -1},
 316         { 9,  0,  3,  9,  3,  5,  5,  3,  7, -1, -1, -1, -1, -1, -1, -1},
 317         { 9,  8,  7,  5,  9,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 318         { 5,  8,  4,  5, 10,  8, 10, 11,  8, -1, -1, -1, -1, -1, -1, -1},
 319         { 5,  0,  4,  5, 11,  0,  5, 10, 11, 11,  3,  0, -1, -1, -1, -1},
 320         { 0,  1,  9,  8,  4, 10,  8, 10, 11, 10,  4,  5, -1, -1, -1, -1},
 321         {10, 11,  4, 10,  4,  5, 11,  3,  4,  9,  4,  1,  3,  1,  4, -1},
 322         { 2,  5,  1,  2,  8,  5,  2, 11,  8,  4,  5,  8, -1, -1, -1, -1},
 323         { 0,  4, 11,  0, 11,  3,  4,  5, 11,  2, 11,  1,  5,  1, 11, -1},
 324         { 0,  2,  5,  0,  5,  9,  2, 11,  5,  4,  5,  8, 11,  8,  5, -1},
 325         { 9,  4,  5,  2, 11,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 326         { 2,  5, 10,  3,  5,  2,  3,  4,  5,  3,  8,  4, -1, -1, -1, -1},
 327         { 5, 10,  2,  5,  2,  4,  4,  2,  0, -1, -1, -1, -1, -1, -1, -1},
 328         { 3, 10,  2,  3,  5, 10,  3,  8,  5,  4,  5,  8,  0,  1,  9, -1},
 329         { 5, 10,  2,  5,  2,  4,  1,  9,  2,  9,  4,  2, -1, -1, -1, -1},
 330         { 8,  4,  5,  8,  5,  3,  3,  5,  1, -1, -1, -1, -1, -1, -1, -1},
 331         { 0,  4,  5,  1,  0,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 332         { 8,  4,  5,  8,  5,  3,  9,  0,  5,  0,  3,  5, -1, -1, -1, -1},
 333         { 9,  4,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 334         { 4, 11,  7,  4,  9, 11,  9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
 335         { 0,  8,  3,  4,  9,  7,  9, 11,  7,  9, 10, 11, -1, -1, -1, -1},
 336         { 1, 10, 11,  1, 11,  4,  1,  4,  0,  7,  4, 11, -1, -1, -1, -1},
 337         { 3,  1,  4,  3,  4,  8,  1, 10,  4,  7,  4, 11, 10, 11,  4, -1},
 338         { 4, 11,  7,  9, 11,  4,  9,  2, 11,  9,  1,  2, -1, -1, -1, -1},
 339         { 9,  7,  4,  9, 11,  7,  9,  1, 11,  2, 11,  1,  0,  8,  3, -1},
 340         {11,  7,  4, 11,  4,  2,  2,  4,  0, -1, -1, -1, -1, -1, -1, -1},
 341         {11,  7,  4, 11,  4,  2,  8,  3,  4,  3,  2,  4, -1, -1, -1, -1},
 342         { 2,  9, 10,  2,  7,  9,  2,  3,  7,  7,  4,  9, -1, -1, -1, -1},
 343         { 9, 10,  7,  9,  7,  4, 10,  2,  7,  8,  7,  0,  2,  0,  7, -1},
 344         { 3,  7, 10,  3, 10,  2,  7,  4, 10,  1, 10,  0,  4,  0, 10, -1},
 345         { 1, 10,  2,  8,  7,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 346         { 4,  9,  1,  4,  1,  7,  7,  1,  3, -1, -1, -1, -1, -1, -1, -1},
 347         { 4,  9,  1,  4,  1,  7,  0,  8,  1,  8,  7,  1, -1, -1, -1, -1},
 348         { 4,  0,  3,  7,  4,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 349         { 4,  8,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 350         { 9, 10,  8, 10, 11,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 351         { 3,  0,  9,  3,  9, 11, 11,  9, 10, -1, -1, -1, -1, -1, -1, -1},
 352         { 0,  1, 10,  0, 10,  8,  8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
 353         { 3,  1, 10, 11,  3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 354         { 1,  2, 11,  1, 11,  9,  9, 11,  8, -1, -1, -1, -1, -1, -1, -1},
 355         { 3,  0,  9,  3,  9, 11,  1,  2,  9,  2, 11,  9, -1, -1, -1, -1},
 356         { 0,  2, 11,  8,  0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 357         { 3,  2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 358         { 2,  3,  8,  2,  8, 10, 10,  8,  9, -1, -1, -1, -1, -1, -1, -1},
 359         { 9, 10,  2,  0,  9,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 360         { 2,  3,  8,  2,  8, 10,  0,  1,  8,  1, 10,  8, -1, -1, -1, -1},
 361         { 1, 10,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 362         { 1,  3,  8,  9,  1,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 363         { 0,  9,  1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 364         { 0,  3,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 365         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
 366 };
 367
 368
 369
 370 /* Linearly interpolate the position where an isosurface cuts
 371  * an edge between two vertices, each with their own scalar value
 372  */
 373 static XYZ interp_vertex (double isolevel, XYZ p1, XYZ p2, double valp1, double valp2)
 374 {
 375         double mu;
 376         XYZ p;
 377
 378         if (ABS(isolevel-valp1) < 0.00001)
 379                 return(p1);
 380         if (ABS(isolevel-valp2) < 0.00001)
 381                 return(p2);
 382         if (ABS(valp1-valp2) < 0.00001)
 383                 return(p1);
 384         mu = (isolevel - valp1) / (valp2 - valp1);
 385         p.x = p1.x + mu * (p2.x - p1.x);
 386         p.y = p1.y + mu * (p2.y - p1.y);
 387         p.z = p1.z + mu * (p2.z - p1.z);
 388
 389         return(p);
 390 }
 391
 392
 393 /* Given a grid cell and an isolevel, calculate the triangular
 394  * facets required to represent the isosurface through the cell.
 395  * Return the number of triangular facets.
 396  * `triangles' will be loaded up with the vertices at most 5 triangular facets.
 397  * 0 will be returned if the grid cell is either totally above
 398  * of totally below the isolevel.
 399  *
 400  * By Paul Bourke <pbourke@swin.edu.au>
 401  */
 402 int march_one_cube (GRIDCELL grid, double isolevel, TRIANGLE *triangles)
 403 {
 404         int i, ntriang;
 405         int cubeindex;
 406         XYZ vertlist[12];
 407
 408         /*
 409          * Determine the index into the edge table which
 410          * tells us which vertices are inside of the surface
 411          */
 412         cubeindex = 0;
 413         if (grid.val[0] < isolevel) cubeindex |= 1;
 414         if (grid.val[1] < isolevel) cubeindex |= 2;
 415         if (grid.val[2] < isolevel) cubeindex |= 4;
 416         if (grid.val[3] < isolevel) cubeindex |= 8;
 417         if (grid.val[4] < isolevel) cubeindex |= 16;
 418         if (grid.val[5] < isolevel) cubeindex |= 32;
 419         if (grid.val[6] < isolevel) cubeindex |= 64;
 420         if (grid.val[7] < isolevel) cubeindex |= 128;
 421
 422         /* Cube is entirely in/out of the surface */
 423         if (edgeTable[cubeindex] == 0)
 424                 return(0);
 425
 426         /* Find the vertices where the surface intersects the cube */
 427         if (edgeTable[cubeindex] & 1)
 428                 vertlist[0] =
 429                         interp_vertex (isolevel,grid.p[0],grid.p[1],grid.val[0],grid.val[1]);
 430         if (edgeTable[cubeindex] & 2)
 431                 vertlist[1] =
 432                         interp_vertex (isolevel,grid.p[1],grid.p[2],grid.val[1],grid.val[2]);
 433         if (edgeTable[cubeindex] & 4)
 434                 vertlist[2] =
 435                         interp_vertex (isolevel,grid.p[2],grid.p[3],grid.val[2],grid.val[3]);
 436         if (edgeTable[cubeindex] & 8)
 437                 vertlist[3] =
 438                         interp_vertex (isolevel,grid.p[3],grid.p[0],grid.val[3],grid.val[0]);
 439         if (edgeTable[cubeindex] & 16)
 440                 vertlist[4] =
 441                         interp_vertex (isolevel,grid.p[4],grid.p[5],grid.val[4],grid.val[5]);
 442         if (edgeTable[cubeindex] & 32)
 443                 vertlist[5] =
 444                         interp_vertex (isolevel,grid.p[5],grid.p[6],grid.val[5],grid.val[6]);
 445         if (edgeTable[cubeindex] & 64)
 446                 vertlist[6] =
 447                         interp_vertex (isolevel,grid.p[6],grid.p[7],grid.val[6],grid.val[7]);
 448         if (edgeTable[cubeindex] & 128)
 449                 vertlist[7] =
 450                         interp_vertex (isolevel,grid.p[7],grid.p[4],grid.val[7],grid.val[4]);
 451         if (edgeTable[cubeindex] & 256)
 452                 vertlist[8] =
 453                         interp_vertex (isolevel,grid.p[0],grid.p[4],grid.val[0],grid.val[4]);
 454         if (edgeTable[cubeindex] & 512)
 455                 vertlist[9] =
 456                         interp_vertex (isolevel,grid.p[1],grid.p[5],grid.val[1],grid.val[5]);
 457         if (edgeTable[cubeindex] & 1024)
 458                 vertlist[10] =
 459                         interp_vertex (isolevel,grid.p[2],grid.p[6],grid.val[2],grid.val[6]);
 460         if (edgeTable[cubeindex] & 2048)
 461                 vertlist[11] =
 462                         interp_vertex (isolevel,grid.p[3],grid.p[7],grid.val[3],grid.val[7]);
 463
 464         /* Create the triangle */
 465         ntriang = 0;
 466         for (i=0; triTable[cubeindex][i] != -1; i+=3)
 467         {
 468                 triangles[ntriang].p[0] = vertlist[triTable[cubeindex][i  ]];
 469                 triangles[ntriang].p[1] = vertlist[triTable[cubeindex][i+1]];
 470                 triangles[ntriang].p[2] = vertlist[triTable[cubeindex][i+2]];
 471                 ntriang++;
 472         }
 473
 474         return(ntriang);
 475 }
 476
 477
 478 /* Walking the grid.  By jwz.
 479  */
 480
 481
 482 /* Computes the normal of the scalar field at the given point,
 483  * for vertex normals (as opposed to face normals.)
 484  */
 485 void do_function_normal (double x, double y, double z,
 486                 double (*compute_fn) (double x, double y, double z,
 487                         void *closure),
 488                 void *c)
 489 {
 490         XYZ n;
 491         double off = 0.5;
 492         n.x = compute_fn (x-off, y, z, c) - compute_fn (x+off, y, z, c);
 493         n.y = compute_fn (x, y-off, z, c) - compute_fn (x, y+off, z, c);
 494         n.z = compute_fn (x, y, z-off, c) - compute_fn (x, y, z+off, c);
 495         /* normalize (&n); */
 496         glNormal3f (n.x, n.y, n.z);
 497 }
 498
 499
 500 /* Given a function capable of generating a value at any XYZ position,
 501  * creates OpenGL faces for the solids defined.
 502  *
 503  * init_fn is called at the beginning for initial, and returns an object.
 504  * free_fn is called at the end.
 505  *
 506  * compute_fn is called for each XYZ in the specified grid, and returns
 507  * the double value of that coordinate.  If smoothing is on, then
 508  * compute_fn will also be called twice more for each emitted vertex,
 509  * in order to calculate vertex normals (so don't assume it will only
 510  * be called with values falling on the grid boundaries.)
 511  *
 512  * Points are inside an object if the are less than `isolevel', and
 513  * outside otherwise.
 514  */
 515 void
 516 marching_cubes (int grid_size,     /* density of the mesh */
 517                 double isolevel,   /* cutoff point for "in" versus "out" */
 518                 int wireframe_p,   /* wireframe, or solid */
 519                 int smooth_p,      /* smooth, or faceted */
 520
 521                 void * (*init_fn)    (double grid_size, void *closure1),
 522                 double (*compute_fn) (double x, double y, double z,
 523                         void *closure2),
 524                 void   (*free_fn)    (void *closure2),
 525                 void *closure1,
 526
 527                 unsigned long *polygon_count)
 528 {
 529         int planesize = grid_size * grid_size;
 530         int x, y, z;
 531         void *closure2 = 0;
 532         unsigned long polys = 0;
 533         double *layers;
 534
 535         layers = (double *) g_malloc0(sizeof (*layers) * planesize * 2);
 536         if (!layers)
 537         {
 538                 g_error("out of memory for %dx%dx%d grid\n",
 539                                 grid_size, grid_size, 2);
 540         }
 541
 542         if (init_fn)
 543                 closure2 = init_fn (grid_size, closure1);
 544
 545         glFrontFace(GL_CCW);
 546         if (!wireframe_p)
 547                 glBegin (GL_TRIANGLES);
 548
 549         for (z = 0; z < grid_size; z++)
 550         {
 551                 double *layer0 = (z & 1 ? layers+planesize : layers);
 552                 double *layer1 = (z & 1 ? layers           : layers+planesize);
 553                 double *row;
 554
 555                 /* Fill in the XY grid on the currently-bottommost layer. */
 556                 row = layer1;
 557                 for (y = 0; y < grid_size; y++, row += grid_size)
 558                 {
 559                         double *cell = row;
 560                         for (x = 0; x < grid_size; x++, cell++)
 561                                 *cell = compute_fn (x, y, z, closure2);
 562                 }
 563
 564                 /* Now we've completed one layer (an XY slice of Z.)  Now we can
 565                    generate the polygons that fill the space between this layer
 566                    and the previous one (unless this is the first layer.)
 567                  */
 568                 if (z == 0) continue;
 569
 570                 for (y = 1; y < grid_size; y += 1)
 571                         for (x = 1; x < grid_size; x += 1)
 572                         {
 573                                 TRIANGLE tri[6];
 574                                 int i, ntri;
 575                                 GRIDCELL cell;
 576
 577                                 /* This is kinda hokey, there ought to be a more efficient
 578                                    way to do this... */
 579                                 cell.p[0].x = x-1; cell.p[0].y = y-1; cell.p[0].z = z-1;
 580                                 cell.p[1].x = x  ; cell.p[1].y = y-1; cell.p[1].z = z-1;
 581                                 cell.p[2].x = x  ; cell.p[2].y = y  ; cell.p[2].z = z-1;
 582                                 cell.p[3].x = x-1; cell.p[3].y = y  ; cell.p[3].z = z-1;
 583                                 cell.p[4].x = x-1; cell.p[4].y = y-1; cell.p[4].z = z  ;
 584                                 cell.p[5].x = x  ; cell.p[5].y = y-1; cell.p[5].z = z  ;
 585                                 cell.p[6].x = x  ; cell.p[6].y = y  ; cell.p[6].z = z  ;
 586                                 cell.p[7].x = x-1; cell.p[7].y = y  ; cell.p[7].z = z  ;
 587
 588 # define GRID(X,Y,WHICH) ((WHICH) \
 589                 ? layer1[((Y)*grid_size) + ((X))] \
 590                 : layer0[((Y)*grid_size) + ((X))])
 591
 592                                 cell.val[0] = GRID (x-1, y-1, 0);
 593                                 cell.val[1] = GRID (x  , y-1, 0);
 594                                 cell.val[2] = GRID (x  , y  , 0);
 595                                 cell.val[3] = GRID (x-1, y  , 0);
 596                                 cell.val[4] = GRID (x-1, y-1, 1);
 597                                 cell.val[5] = GRID (x  , y-1, 1);
 598                                 cell.val[6] = GRID (x  , y  , 1);
 599                                 cell.val[7] = GRID (x-1, y  , 1);
 600 # undef GRID
 601
 602                                 /* Now generate the triangles for this cubic segment,
 603                                    and emit the GL faces.
 604                                  */
 605                                 ntri = march_one_cube (cell, isolevel, tri);
 606                                 polys += ntri;
 607                                 for (i = 0; i < ntri; i++)
 608                                 {
 609                                         if (wireframe_p) glBegin (GL_LINE_LOOP);
 610
 611                                         /* If we're smoothing, we need to call the field function
 612                                            again for each vertex (via function_normal().)  If we're
 613                                            not smoothing, then we can just compute the normal from
 614                                            this triangle.
 615                                          */
 616                                         if (!smooth_p)
 617                                                 do_normal (tri[i].p[0].x, tri[i].p[0].y, tri[i].p[0].z,
 618                                                                 tri[i].p[1].x, tri[i].p[1].y, tri[i].p[1].z,
 619                                                                 tri[i].p[2].x, tri[i].p[2].y, tri[i].p[2].z);
 620
 621 # define VERT(X,Y,Z) \
 622                                         if (smooth_p) \
 623                                         do_function_normal ((X), (Y), (Z), compute_fn, closure2); \
 624                                         glVertex3f ((X), (Y), (Z))
 625
 626                                         VERT (tri[i].p[0].x, tri[i].p[0].y, tri[i].p[0].z);
 627                                         VERT (tri[i].p[1].x, tri[i].p[1].y, tri[i].p[1].z);
 628                                         VERT (tri[i].p[2].x, tri[i].p[2].y, tri[i].p[2].z);
 629 # undef VERT
 630                                         if (wireframe_p) glEnd ();
 631                                 }
 632                         }
 633         }
 634
 635         if (!wireframe_p)
 636                 glEnd ();
 637
 638         g_free(layers);
 639
 640         if (free_fn)
 641                 free_fn (closure2);
 642
 643         if (polygon_count)
 644                 *polygon_count = polys;
 645 }