git.hungrycats.org Git - xscreensaver/blob - hacks/glx/marching.c

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