take mass into account
[dosdemo] / src / metasurf.c
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include "metasurf.h"
5 #include "mcubes.h"
6
7 typedef float vec3[3];
8
9 struct metasurface {
10         vec3 min, max;
11         int res[3], newres[3];
12         float thres;
13
14         float dx, dy, dz;
15         unsigned int flags;
16
17         float *voxels;
18
19         int varr_size, varr_alloc_size;
20         struct g3d_vertex *varr;
21 };
22
23 static int msurf_init(struct metasurface *ms);
24 static void process_cell(struct metasurface *ms, int xcell, int ycell, int zcell, vec3 pos, vec3 sz);
25
26
27 struct metasurface *msurf_create(void)
28 {
29         struct metasurface *ms;
30
31         if(!(ms = malloc(sizeof *ms))) {
32                 return 0;
33         }
34         if(msurf_init(ms) == -1) {
35                 free(ms);
36         }
37         return ms;
38 }
39
40 void msurf_free(struct metasurface *ms)
41 {
42         if(ms) {
43                 free(ms->voxels);
44                 free(ms->varr);
45                 free(ms);
46         }
47 }
48
49 static int msurf_init(struct metasurface *ms)
50 {
51         ms->voxels = 0;
52         ms->thres = 0.0;
53         ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
54         ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
55         ms->res[0] = ms->res[1] = ms->res[2] = 0;
56         ms->newres[0] = ms->newres[1] = ms->newres[2] = 40;
57
58         ms->varr_alloc_size = ms->varr_size = 0;
59         ms->varr = 0;
60
61         ms->dx = ms->dy = ms->dz = 0.001;
62         ms->flags = 0;
63
64         return 0;
65 }
66
67 void msurf_enable(struct metasurface *ms, unsigned int opt)
68 {
69         ms->flags |= opt;
70 }
71
72 void msurf_disable(struct metasurface *ms, unsigned int opt)
73 {
74         ms->flags &= ~opt;
75 }
76
77 int msurf_is_enabled(struct metasurface *ms, unsigned int opt)
78 {
79         return ms->flags & opt;
80 }
81
82 void msurf_set_inside(struct metasurface *ms, int inside)
83 {
84         switch(inside) {
85         case MSURF_GREATER:
86                 msurf_enable(ms, MSURF_FLIP);
87                 break;
88
89         case MSURF_LESS:
90                 msurf_disable(ms, MSURF_FLIP);
91                 break;
92
93         default:
94                 fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
95         }
96 }
97
98 int msurf_get_inside(struct metasurface *ms)
99 {
100         return msurf_is_enabled(ms, MSURF_FLIP) ? MSURF_LESS : MSURF_GREATER;
101 }
102
103 void msurf_set_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
104 {
105         ms->min[0] = xmin;
106         ms->min[1] = ymin;
107         ms->min[2] = zmin;
108         ms->max[0] = xmax;
109         ms->max[1] = ymax;
110         ms->max[2] = zmax;
111 }
112
113 void msurf_get_bounds(struct metasurface *ms, float *xmin, float *ymin, float *zmin, float *xmax, float *ymax, float *zmax)
114 {
115         *xmin = ms->min[0];
116         *ymin = ms->min[1];
117         *zmin = ms->min[2];
118         *xmax = ms->max[0];
119         *ymax = ms->max[1];
120         *zmax = ms->max[2];
121 }
122
123 void msurf_set_resolution(struct metasurface *ms, int xres, int yres, int zres)
124 {
125         ms->newres[0] = xres;
126         ms->newres[1] = yres;
127         ms->newres[2] = zres;
128 }
129
130 void msurf_get_resolution(struct metasurface *ms, int *xres, int *yres, int *zres)
131 {
132         *xres = ms->res[0];
133         *yres = ms->res[1];
134         *zres = ms->res[2];
135 }
136
137 void msurf_set_threshold(struct metasurface *ms, float thres)
138 {
139         ms->thres = thres;
140 }
141
142 float msurf_get_threshold(struct metasurface *ms)
143 {
144         return ms->thres;
145 }
146
147
148 float *msurf_voxels(struct metasurface *ms)
149 {
150         if(ms->res[0] != ms->newres[0] || ms->res[1] != ms->newres[1] || ms->res[2] != ms->newres[2]) {
151                 int count;
152                 ms->res[0] = ms->newres[0];
153                 ms->res[1] = ms->newres[1];
154                 ms->res[2] = ms->newres[2];
155                 count = ms->res[0] * ms->res[1] * ms->res[2];
156                 free(ms->voxels);
157                 if(!(ms->voxels = malloc(count * sizeof *ms->voxels))) {
158                         return 0;
159                 }
160         }
161         return ms->voxels;
162 }
163
164 float *msurf_slice(struct metasurface *ms, int idx)
165 {
166         float *vox = msurf_voxels(ms);
167         if(!vox) return 0;
168
169         return vox + ms->res[0] * ms->res[1] * idx;
170 }
171
172 int msurf_polygonize(struct metasurface *ms)
173 {
174         int i, j, k;
175         vec3 pos, delta;
176
177         if(!ms->voxels) return -1;
178
179         ms->varr_size = 0;
180
181         for(i=0; i<3; i++) {
182                 delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
183         }
184
185         for(i=0; i<ms->res[0] - 1; i++) {
186                 for(j=0; j<ms->res[1] - 1; j++) {
187                         for(k=0; k<ms->res[2] - 1; k++) {
188
189                                 pos[0] = ms->min[0] + i * delta[0];
190                                 pos[1] = ms->min[1] + j * delta[1];
191                                 pos[2] = ms->min[2] + k * delta[2];
192
193                                 process_cell(ms, i, j, k, pos, delta);
194                         }
195                 }
196         }
197         return 0;
198 }
199
200 int msurf_vertex_count(struct metasurface *ms)
201 {
202         return ms->varr_size;
203 }
204
205 struct g3d_vertex *msurf_vertices(struct metasurface *ms)
206 {
207         return ms->varr;
208 }
209
210 static unsigned int mc_bitcode(float *val, float thres);
211
212 static void process_cell(struct metasurface *ms, int xcell, int ycell, int zcell, vec3 cellpos, vec3 cellsz)
213 {
214         int i, j, slice_size;
215         vec3 pos[8];
216         float dfdx[8], dfdy[8], dfdz[8];
217         vec3 vert[12], norm[12];
218         float val[8];
219         float *cellptr;
220         unsigned int code;
221
222         static const int offs[][3] = {
223                 {0, 0, 0},
224                 {1, 0, 0},
225                 {1, 1, 0},
226                 {0, 1, 0},
227                 {0, 0, 1},
228                 {1, 0, 1},
229                 {1, 1, 1},
230                 {0, 1, 1}
231         };
232
233         static const int pidx[12][2] = {
234                 {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
235                 {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
236         };
237
238         slice_size = ms->res[0] * ms->res[1];
239         cellptr = ms->voxels + slice_size * zcell + ms->res[0] * ycell + xcell;
240
241 #define GRIDOFFS(x, y, z)       ((z) * slice_size + (y) * ms->res[0] + (x))
242
243         for(i=0; i<8; i++) {
244                 val[i] = cellptr[GRIDOFFS(offs[i][0], offs[i][1], offs[i][2])];
245         }
246
247         code = mc_bitcode(val, ms->thres);
248         if(ms->flags & MSURF_FLIP) {
249                 code = ~code & 0xff;
250         }
251         if(mc_edge_table[code] == 0) {
252                 return;
253         }
254
255         /* calculate normals at the 8 corners */
256         for(i=0; i<8; i++) {
257                 float *ptr = cellptr + GRIDOFFS(offs[i][0], offs[i][1], offs[i][2]);
258
259                 if(xcell < ms->res[0] - 1) {
260                         dfdx[i] = ptr[GRIDOFFS(1, 0, 0)] - *ptr;
261                 } else {
262                         dfdx[i] = *ptr - ptr[GRIDOFFS(-1, 0, 0)];
263                 }
264                 if(ycell < ms->res[1] - 1) {
265                         dfdy[i] = ptr[GRIDOFFS(0, 1, 0)] - *ptr;
266                 } else {
267                         dfdy[i] = *ptr - ptr[GRIDOFFS(0, -1, 0)];
268                 }
269                 if(zcell < ms->res[2] - 1) {
270                         dfdz[i] = ptr[GRIDOFFS(0, 0, 1)] - *ptr;
271                 } else {
272                         dfdz[i] = *ptr - ptr[GRIDOFFS(0, 0, -1)];
273                 }
274         }
275
276         /* calculate the world-space position of each corner */
277         for(i=0; i<8; i++) {
278                 pos[i][0] = cellpos[0] + cellsz[0] * offs[i][0];
279                 pos[i][1] = cellpos[1] + cellsz[1] * offs[i][1];
280                 pos[i][2] = cellpos[2] + cellsz[2] * offs[i][2];
281         }
282
283         /* generate up to a max of 12 vertices per cube. interpolate positions and normals for each one */
284         for(i=0; i<12; i++) {
285                 if(mc_edge_table[code] & (1 << i)) {
286                         float nx, ny, nz;
287                         int p0 = pidx[i][0];
288                         int p1 = pidx[i][1];
289
290                         float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
291                         vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
292                         vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
293                         vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
294
295                         nx = dfdx[p0] + (dfdx[p1] - dfdx[p0]) * t;
296                         ny = dfdy[p0] + (dfdy[p1] - dfdy[p0]) * t;
297                         nz = dfdz[p0] + (dfdz[p1] - dfdz[p0]) * t;
298
299                         if(ms->flags & MSURF_FLIP) {
300                                 nx = -nx;
301                                 ny = -ny;
302                                 nz = -nz;
303                         }
304
305                         if(ms->flags & MSURF_NORMALIZE) {
306                                 float len = sqrt(nx * nx + ny * ny + nz * nz);
307                                 if(len != 0.0) {
308                                         float s = 1.0f / len;
309                                         nx *= s;
310                                         ny *= s;
311                                         nz *= s;
312                                 }
313                         }
314
315                         norm[i][0] = nx;
316                         norm[i][1] = ny;
317                         norm[i][2] = nz;
318                 }
319         }
320
321         /* for each triangle of the cube add the appropriate vertices to the vertex buffer */
322         for(i=0; mc_tri_table[code][i] != -1; i+=3) {
323                 for(j=0; j<3; j++) {
324                         int idx = mc_tri_table[code][i + j];
325                         float *v = vert[idx];
326                         float *n = norm[idx];
327                         struct g3d_vertex *vdest;
328
329                         if(ms->varr_size >= ms->varr_alloc_size) {
330                                 int newsz = ms->varr_alloc_size ? ms->varr_alloc_size * 2 : 1024;
331                                 struct g3d_vertex *new_varr;
332
333                                 if(!(new_varr = realloc(ms->varr, newsz * sizeof *new_varr))) {
334                                         fprintf(stderr, "msurf_polygonize: failed to grow vertex buffers to %d elements\n", newsz);
335                                         return;
336                                 }
337                                 ms->varr = new_varr;
338                                 ms->varr_alloc_size = newsz;
339                         }
340
341                         vdest = ms->varr + ms->varr_size++;
342                         vdest->x = v[0];
343                         vdest->y = v[1];
344                         vdest->z = v[2];
345                         vdest->w = 1.0f;
346                         vdest->nx = n[0];
347                         vdest->ny = n[1];
348                         vdest->nz = n[2];
349                 }
350         }
351 }
352
353 static unsigned int mc_bitcode(float *val, float thres)
354 {
355         unsigned int i, res = 0;
356
357         for(i=0; i<8; i++) {
358                 if(val[i] > thres) {
359                         res |= 1 << i;
360                 }
361         }
362         return res;
363 }