2 metasurf - a library for implicit surface polygonization
3 Copyright (C) 2011-2016 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/>.
18 /* this is pulled from: https://github.com/jtsiomb/metasurf */
27 #if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))
28 #error "pick either USE_MTETRA or USE_MCUBES, not both..."
31 typedef float vec3[3];
38 msurf_eval_func_t eval;
39 msurf_vertex_func_t vertex;
40 msurf_normal_func_t normal;
50 static int msurf_init(struct metasurface *ms);
51 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
53 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
56 static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
60 struct metasurface *msurf_create(void)
62 struct metasurface *ms;
64 if(!(ms = malloc(sizeof *ms))) {
67 if(msurf_init(ms) == -1) {
73 void msurf_free(struct metasurface *ms)
78 static int msurf_init(struct metasurface *ms)
85 ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
86 ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
87 ms->res[0] = ms->res[1] = ms->res[2] = 40;
90 ms->dx = ms->dy = ms->dz = 0.001;
96 void msurf_set_user_data(struct metasurface *ms, void *udata)
101 void *msurf_get_user_data(struct metasurface *ms)
106 void msurf_set_inside(struct metasurface *ms, int inside)
118 fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
122 int msurf_get_inside(struct metasurface *ms)
124 return ms->flip ? MSURF_LESS : MSURF_GREATER;
127 void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
132 void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
137 void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
142 void msurf_set_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
152 void msurf_get_bounds(struct metasurface *ms, float *xmin, float *ymin, float *zmin, float *xmax, float *ymax, float *zmax)
162 void msurf_set_resolution(struct metasurface *ms, int xres, int yres, int zres)
169 void msurf_get_resolution(struct metasurface *ms, int *xres, int *yres, int *zres)
176 void msurf_set_threshold(struct metasurface *ms, float thres)
181 float msurf_get_threshold(struct metasurface *ms)
187 int msurf_polygonize(struct metasurface *ms)
192 if(!ms->eval || !ms->vertex) {
193 fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
198 delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
202 for(i=0; i<ms->res[0] - 1; i++) {
204 for(j=0; j<ms->res[1] - 1; j++) {
207 for(k=0; k<ms->res[2] - 1; k++) {
209 process_cell(ms, pos, delta);
221 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
228 static int tetra[][4] = {
238 static const float offs[][3] = {
250 p[i][0] = pos[0] + sz[0] * offs[i][2];
251 p[i][1] = pos[1] + sz[1] * offs[i][1];
252 p[i][2] = pos[2] + sz[2] * offs[i][0];
254 val[i] = ms->eval(ms, p[i][0], p[i][1], p[i][2]);
259 process_tetra(ms, tetra[i], p, val);
263 process_cube(ms, p, val);
268 /* ---- marching cubes implementation ---- */
271 static unsigned int mc_bitcode(float *val, float thres);
273 static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
275 static const int pidx[12][2] = {
276 {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
277 {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
281 unsigned int code = mc_bitcode(val, ms->thres);
287 if(mc_edge_table[code] == 0) {
291 for(i=0; i<12; i++) {
292 if(mc_edge_table[code] & (1 << i)) {
296 float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
297 vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
298 vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
299 vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
303 for(i=0; mc_tri_table[code][i] != -1; i+=3) {
305 float *v = vert[mc_tri_table[code][i + j]];
308 float dfdx, dfdy, dfdz;
309 dfdx = ms->eval(ms, v[0] - ms->dx, v[1], v[2]) - ms->eval(ms, v[0] + ms->dx, v[1], v[2]);
310 dfdy = ms->eval(ms, v[0], v[1] - ms->dy, v[2]) - ms->eval(ms, v[0], v[1] + ms->dy, v[2]);
311 dfdz = ms->eval(ms, v[0], v[1], v[2] - ms->dz) - ms->eval(ms, v[0], v[1], v[2] + ms->dz);
318 ms->normal(ms, dfdx, dfdy, dfdz);
321 /* TODO multithreadied polygon emmit */
322 ms->vertex(ms, v[0], v[1], v[2]);
327 static unsigned int mc_bitcode(float *val, float thres)
329 unsigned int i, res = 0;
338 #endif /* USE_MCUBES */
341 /* ---- marching tetrahedra implementation (incomplete) ---- */
344 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
345 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
348 #define REVBIT(x) ((x) & 8)
349 #define INV(x) (~(x) & 0xf)
350 #define EDGE(a, b) emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
351 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
353 unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
418 break; /* cases 0 and 15 */
422 #define BIT(i) ((v##i > thres) ? (1 << i) : 0)
423 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
425 return BIT(0) | BIT(1) | BIT(2) | BIT(3);
428 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
431 float t = (ms->thres - v0) / (v1 - v0);
435 p[i] = p0[i] + (p1[i] - p0[i]) * t;
437 ms->vertex(ms, p[0], p[1], p[2]);
439 /*for(i=0; i<3; i++) {
440 ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
443 if(++ms->nverts >= 3) {
447 int idx = rev ? (2 - i) : i;
448 ms->vertex(ms, ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
452 #endif /* USE_MTETRA */