dbf-udg

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

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