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