metasurf

annotate src/metasurf.c @ 4:2c575855f707

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