dbf-udg

view libs/metasurf/metasurf.c @ 5:e09cbb2e9d4f

metaballs
author John Tsiombikas <nuclear@member.fsf.org>
date Mon, 18 Feb 2013 03:46:52 +0200
parents
children
line source
1 /*
2 metasurf - a library for implicit surface polygonization
3 Copyright (C) 2011 John Tsiombikas <nuclear@member.fsf.org>
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
15 You should have received a copy of the GNU Lesser General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include "metasurf.h"
21 #include "mcubes.h"
23 #undef USE_MTETRA
24 #define USE_MCUBES
26 #if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))
27 #error "pick either USE_MTETRA or USE_MCUBES, not both..."
28 #endif
30 typedef float vec3[3];
32 struct metasurface {
33 vec3 min, max;
34 int res[3];
35 float thres;
37 msurf_eval_func_t eval;
38 msurf_vertex_func_t vertex;
39 msurf_normal_func_t normal;
41 float dx, dy, dz;
42 int flip;
44 vec3 vbuf[3];
45 int nverts;
46 };
48 static int msurf_init(struct metasurface *ms);
49 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
50 #ifdef USE_MTETRA
51 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
52 #endif
53 #ifdef USE_MCUBES
54 static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
55 #endif
58 struct metasurface *msurf_create(void)
59 {
60 struct metasurface *ms;
62 if(!(ms = malloc(sizeof *ms))) {
63 return 0;
64 }
65 if(msurf_init(ms) == -1) {
66 free(ms);
67 }
68 return ms;
69 }
71 void msurf_free(struct metasurface *ms)
72 {
73 free(ms);
74 }
76 static int msurf_init(struct metasurface *ms)
77 {
78 ms->thres = 0.0;
79 ms->eval = 0;
80 ms->vertex = 0;
81 ms->normal = 0;
82 ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
83 ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
84 ms->res[0] = ms->res[1] = ms->res[2] = 40;
85 ms->nverts = 0;
87 ms->dx = ms->dy = ms->dz = 0.001;
88 ms->flip = 0;
90 return 0;
91 }
93 void msurf_inside(struct metasurface *ms, int inside)
94 {
95 switch(inside) {
96 case MSURF_GREATER:
97 ms->flip = 0;
98 break;
100 case MSURF_LESS:
101 ms->flip = 1;
102 break;
104 default:
105 fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
106 }
107 }
109 void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
110 {
111 ms->eval = func;
112 }
114 void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
115 {
116 ms->vertex = func;
117 }
119 void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
120 {
121 ms->normal = func;
122 }
124 void msurf_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
125 {
126 ms->min[0] = xmin;
127 ms->min[1] = ymin;
128 ms->min[2] = zmin;
129 ms->max[0] = xmax;
130 ms->max[1] = ymax;
131 ms->max[2] = zmax;
132 }
134 void msurf_resolution(struct metasurface *ms, int xres, int yres, int zres)
135 {
136 ms->res[0] = xres;
137 ms->res[1] = yres;
138 ms->res[2] = zres;
139 }
141 void msurf_threshold(struct metasurface *ms, float thres)
142 {
143 ms->thres = thres;
144 }
147 int msurf_polygonize(struct metasurface *ms)
148 {
149 int i, j, k;
150 vec3 pos, delta;
152 if(!ms->eval || !ms->vertex) {
153 fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
154 return -1;
155 }
157 for(i=0; i<3; i++) {
158 delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
159 }
161 pos[0] = ms->min[0];
162 for(i=0; i<ms->res[0] - 1; i++) {
164 pos[1] = ms->min[1];
165 for(j=0; j<ms->res[1] - 1; j++) {
167 pos[2] = ms->min[2];
168 for(k=0; k<ms->res[2] - 1; k++) {
170 process_cell(ms, pos, delta);
172 pos[2] += delta[2];
173 }
174 pos[1] += delta[1];
175 }
176 pos[0] += delta[0];
177 }
178 return 0;
179 }
182 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
183 {
184 int i;
185 vec3 p[8];
186 float val[8];
188 #ifdef USE_MTETRA
189 static int tetra[][4] = {
190 {0, 2, 3, 7},
191 {0, 2, 6, 7},
192 {0, 4, 6, 7},
193 {0, 6, 1, 2},
194 {0, 6, 1, 4},
195 {5, 6, 1, 4}
196 };
197 #endif
199 static const float offs[][3] = {
200 {0.0f, 0.0f, 0.0f},
201 {1.0f, 0.0f, 0.0f},
202 {1.0f, 1.0f, 0.0f},
203 {0.0f, 1.0f, 0.0f},
204 {0.0f, 0.0f, 1.0f},
205 {1.0f, 0.0f, 1.0f},
206 {1.0f, 1.0f, 1.0f},
207 {0.0f, 1.0f, 1.0f}
208 };
210 for(i=0; i<8; i++) {
211 p[i][0] = pos[0] + sz[0] * offs[i][2];
212 p[i][1] = pos[1] + sz[1] * offs[i][1];
213 p[i][2] = pos[2] + sz[2] * offs[i][0];
215 val[i] = ms->eval(p[i][0], p[i][1], p[i][2]);
216 }
218 #ifdef USE_MTETRA
219 for(i=0; i<6; i++) {
220 process_tetra(ms, tetra[i], p, val);
221 }
222 #endif
223 #ifdef USE_MCUBES
224 process_cube(ms, p, val);
225 #endif
226 }
229 /* ---- marching cubes implementation ---- */
230 #ifdef USE_MCUBES
232 static unsigned int mc_bitcode(float *val, float thres);
234 static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
235 {
236 static const int pidx[12][2] = {
237 {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
238 {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
239 };
240 int i, j;
241 vec3 vert[12];
242 unsigned int code = mc_bitcode(val, ms->thres);
244 if(ms->flip) {
245 code = ~code & 0xff;
246 }
248 if(mc_edge_table[code] == 0) {
249 return;
250 }
252 for(i=0; i<12; i++) {
253 if(mc_edge_table[code] & (1 << i)) {
254 int p0 = pidx[i][0];
255 int p1 = pidx[i][1];
257 float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
258 vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
259 vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
260 vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
261 }
262 }
264 for(i=0; mc_tri_table[code][i] != -1; i+=3) {
265 for(j=0; j<3; j++) {
266 float *v = vert[mc_tri_table[code][i + j]];
268 if(ms->normal) {
269 float dfdx, dfdy, dfdz;
270 dfdx = ms->eval(v[0] - ms->dx, v[1], v[2]) - ms->eval(v[0] + ms->dx, v[1], v[2]);
271 dfdy = ms->eval(v[0], v[1] - ms->dy, v[2]) - ms->eval(v[0], v[1] + ms->dy, v[2]);
272 dfdz = ms->eval(v[0], v[1], v[2] - ms->dz) - ms->eval(v[0], v[1], v[2] + ms->dz);
274 if(ms->flip) {
275 dfdx = -dfdx;
276 dfdy = -dfdy;
277 dfdz = -dfdz;
278 }
279 ms->normal(dfdx, dfdy, dfdz);
280 }
282 ms->vertex(v[0], v[1], v[2]);
283 }
284 }
285 }
287 static unsigned int mc_bitcode(float *val, float thres)
288 {
289 unsigned int i, res = 0;
291 for(i=0; i<8; i++) {
292 if(val[i] > thres) {
293 res |= 1 << i;
294 }
295 }
296 return res;
297 }
298 #endif /* USE_MCUBES */
301 /* ---- marching tetrahedra implementation (incomplete) ---- */
302 #ifdef USE_MTETRA
304 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
305 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
308 #define REVBIT(x) ((x) & 8)
309 #define INV(x) (~(x) & 0xf)
310 #define EDGE(a, b) emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
311 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
312 {
313 unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
315 switch(code) {
316 case 1:
317 case INV(1):
318 EDGE(0, 1);
319 EDGE(0, 2);
320 EDGE(0, 3);
321 break;
323 case 2:
324 case INV(2):
325 EDGE(1, 0);
326 EDGE(1, 3);
327 EDGE(1, 2);
328 break;
330 case 3:
331 case INV(3):
332 EDGE(0, 3);
333 EDGE(0, 2);
334 EDGE(1, 3);
336 EDGE(1, 3);
337 EDGE(1, 2);
338 EDGE(0, 2);
339 break;
341 case 4:
342 case INV(4):
343 EDGE(2, 0);
344 EDGE(2, 1);
345 EDGE(2, 3);
346 break;
348 case 5:
349 case INV(5):
350 EDGE(0, 1);
351 EDGE(2, 3);
352 EDGE(0, 3);
354 EDGE(0, 1);
355 EDGE(1, 2);
356 EDGE(2, 3);
357 break;
359 case 6:
360 case INV(6):
361 EDGE(0, 1);
362 EDGE(1, 3);
363 EDGE(2, 3);
365 EDGE(0, 1);
366 EDGE(0, 2);
367 EDGE(2, 3);
368 break;
370 case 7:
371 case INV(7):
372 EDGE(3, 0);
373 EDGE(3, 2);
374 EDGE(3, 1);
375 break;
377 default:
378 break; /* cases 0 and 15 */
379 }
380 }
382 #define BIT(i) ((v##i > thres) ? (1 << i) : 0)
383 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
384 {
385 return BIT(0) | BIT(1) | BIT(2) | BIT(3);
386 }
388 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
389 {
390 int i;
391 float t = (ms->thres - v0) / (v1 - v0);
393 vec3 p;
394 for(i=0; i<3; i++) {
395 p[i] = p0[i] + (p1[i] - p0[i]) * t;
396 }
397 ms->vertex(p[0], p[1], p[2]);
399 /*for(i=0; i<3; i++) {
400 ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
401 }
403 if(++ms->nverts >= 3) {
404 ms->nverts = 0;
406 for(i=0; i<3; i++) {
407 int idx = rev ? (2 - i) : i;
408 ms->vertex(ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
409 }
410 }*/
411 }
413 #endif /* USE_MTETRA */