a4cafb4860c43745865865beee252205adc10476
[xscreensaver] / 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 }