git.mutantstargoat.com Git - laserbrain_demo/blob - src/blobs/metasurf.c

   1 /*
   2 metasurf - a library for implicit surface polygonization
   3 Copyright (C) 2011-2016  John Tsiombikas <nuclear@member.fsf.org>
   4
   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.
   9
  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.
  14
  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 /* this is pulled from: https://github.com/jtsiomb/metasurf */
  19 #include <stdio.h>
  20 #include <stdlib.h>
  21 #include "metasurf.h"
  22 #include "mcubes.h"
  23
  24 #undef USE_MTETRA
  25 #define USE_MCUBES
  26
  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..."
  29 #endif
  30
  31 typedef float vec3[3];
  32
  33 struct metasurface {
  34         vec3 min, max;
  35         int res[3];
  36         float thres;
  37
  38         msurf_eval_func_t eval;
  39         msurf_vertex_func_t vertex;
  40         msurf_normal_func_t normal;
  41         void *udata;
  42
  43         float dx, dy, dz;
  44         int flip;
  45
  46         vec3 vbuf[3];
  47         int nverts;
  48 };
  49
  50 static int msurf_init(struct metasurface *ms);
  51 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
  52 #ifdef USE_MTETRA
  53 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
  54 #endif
  55 #ifdef USE_MCUBES
  56 static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
  57 #endif
  58
  59
  60 struct metasurface *msurf_create(void)
  61 {
  62         struct metasurface *ms;
  63
  64         if(!(ms = malloc(sizeof *ms))) {
  65                 return 0;
  66         }
  67         if(msurf_init(ms) == -1) {
  68                 free(ms);
  69         }
  70         return ms;
  71 }
  72
  73 void msurf_free(struct metasurface *ms)
  74 {
  75         free(ms);
  76 }
  77
  78 static int msurf_init(struct metasurface *ms)
  79 {
  80         ms->thres = 0.0;
  81         ms->eval = 0;
  82         ms->vertex = 0;
  83         ms->normal = 0;
  84         ms->udata = 0;
  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;
  88         ms->nverts = 0;
  89
  90         ms->dx = ms->dy = ms->dz = 0.001;
  91         ms->flip = 0;
  92
  93         return 0;
  94 }
  95
  96 void msurf_set_user_data(struct metasurface *ms, void *udata)
  97 {
  98         ms->udata = udata;
  99 }
 100
 101 void *msurf_get_user_data(struct metasurface *ms)
 102 {
 103         return ms->udata;
 104 }
 105
 106 void msurf_set_inside(struct metasurface *ms, int inside)
 107 {
 108         switch(inside) {
 109         case MSURF_GREATER:
 110                 ms->flip = 0;
 111                 break;
 112
 113         case MSURF_LESS:
 114                 ms->flip = 1;
 115                 break;
 116
 117         default:
 118                 fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
 119         }
 120 }
 121
 122 int msurf_get_inside(struct metasurface *ms)
 123 {
 124         return ms->flip ? MSURF_LESS : MSURF_GREATER;
 125 }
 126
 127 void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
 128 {
 129         ms->eval = func;
 130 }
 131
 132 void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
 133 {
 134         ms->vertex = func;
 135 }
 136
 137 void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
 138 {
 139         ms->normal = func;
 140 }
 141
 142 void msurf_set_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
 143 {
 144         ms->min[0] = xmin;
 145         ms->min[1] = ymin;
 146         ms->min[2] = zmin;
 147         ms->max[0] = xmax;
 148         ms->max[1] = ymax;
 149         ms->max[2] = zmax;
 150 }
 151
 152 void msurf_get_bounds(struct metasurface *ms, float *xmin, float *ymin, float *zmin, float *xmax, float *ymax, float *zmax)
 153 {
 154         *xmin = ms->min[0];
 155         *ymin = ms->min[1];
 156         *zmin = ms->min[2];
 157         *xmax = ms->max[0];
 158         *ymax = ms->max[1];
 159         *zmax = ms->max[2];
 160 }
 161
 162 void msurf_set_resolution(struct metasurface *ms, int xres, int yres, int zres)
 163 {
 164         ms->res[0] = xres;
 165         ms->res[1] = yres;
 166         ms->res[2] = zres;
 167 }
 168
 169 void msurf_get_resolution(struct metasurface *ms, int *xres, int *yres, int *zres)
 170 {
 171         *xres = ms->res[0];
 172         *yres = ms->res[1];
 173         *zres = ms->res[2];
 174 }
 175
 176 void msurf_set_threshold(struct metasurface *ms, float thres)
 177 {
 178         ms->thres = thres;
 179 }
 180
 181 float msurf_get_threshold(struct metasurface *ms)
 182 {
 183         return ms->thres;
 184 }
 185
 186
 187 int msurf_polygonize(struct metasurface *ms)
 188 {
 189         int i, j, k;
 190         vec3 pos, delta;
 191
 192         if(!ms->eval || !ms->vertex) {
 193                 fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
 194                 return -1;
 195         }
 196
 197         for(i=0; i<3; i++) {
 198                 delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
 199         }
 200
 201         pos[0] = ms->min[0];
 202         for(i=0; i<ms->res[0] - 1; i++) {
 203                 pos[1] = ms->min[1];
 204                 for(j=0; j<ms->res[1] - 1; j++) {
 205
 206                         pos[2] = ms->min[2];
 207                         for(k=0; k<ms->res[2] - 1; k++) {
 208
 209                                 process_cell(ms, pos, delta);
 210
 211                                 pos[2] += delta[2];
 212                         }
 213                         pos[1] += delta[1];
 214                 }
 215                 pos[0] += delta[0];
 216         }
 217         return 0;
 218 }
 219
 220
 221 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
 222 {
 223         int i;
 224         vec3 p[8];
 225         float val[8];
 226
 227 #ifdef USE_MTETRA
 228         static int tetra[][4] = {
 229                 {0, 2, 3, 7},
 230                 {0, 2, 6, 7},
 231                 {0, 4, 6, 7},
 232                 {0, 6, 1, 2},
 233                 {0, 6, 1, 4},
 234                 {5, 6, 1, 4}
 235         };
 236 #endif
 237
 238         static const float offs[][3] = {
 239                 {0.0f, 0.0f, 0.0f},
 240                 {1.0f, 0.0f, 0.0f},
 241                 {1.0f, 1.0f, 0.0f},
 242                 {0.0f, 1.0f, 0.0f},
 243                 {0.0f, 0.0f, 1.0f},
 244                 {1.0f, 0.0f, 1.0f},
 245                 {1.0f, 1.0f, 1.0f},
 246                 {0.0f, 1.0f, 1.0f}
 247         };
 248
 249         for(i=0; i<8; i++) {
 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];
 253
 254                 val[i] = ms->eval(ms, p[i][0], p[i][1], p[i][2]);
 255         }
 256
 257 #ifdef USE_MTETRA
 258         for(i=0; i<6; i++) {
 259                 process_tetra(ms, tetra[i], p, val);
 260         }
 261 #endif
 262 #ifdef USE_MCUBES
 263         process_cube(ms, p, val);
 264 #endif
 265 }
 266
 267
 268 /* ---- marching cubes implementation ---- */
 269 #ifdef USE_MCUBES
 270
 271 static unsigned int mc_bitcode(float *val, float thres);
 272
 273 static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
 274 {
 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}
 278         };
 279         int i, j;
 280         vec3 vert[12];
 281         unsigned int code = mc_bitcode(val, ms->thres);
 282
 283         if(ms->flip) {
 284                 code = ~code & 0xff;
 285         }
 286
 287         if(mc_edge_table[code] == 0) {
 288                 return;
 289         }
 290
 291         for(i=0; i<12; i++) {
 292                 if(mc_edge_table[code] & (1 << i)) {
 293                         int p0 = pidx[i][0];
 294                         int p1 = pidx[i][1];
 295
 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;
 300                 }
 301         }
 302
 303         for(i=0; mc_tri_table[code][i] != -1; i+=3) {
 304                 for(j=0; j<3; j++) {
 305                         float *v = vert[mc_tri_table[code][i + j]];
 306
 307                         if(ms->normal) {
 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);
 312
 313                                 if(ms->flip) {
 314                                         dfdx = -dfdx;
 315                                         dfdy = -dfdy;
 316                                         dfdz = -dfdz;
 317                                 }
 318                                 ms->normal(ms, dfdx, dfdy, dfdz);
 319                         }
 320
 321                         /* TODO multithreadied polygon emmit */
 322                         ms->vertex(ms, v[0], v[1], v[2]);
 323                 }
 324         }
 325 }
 326
 327 static unsigned int mc_bitcode(float *val, float thres)
 328 {
 329         unsigned int i, res = 0;
 330
 331         for(i=0; i<8; i++) {
 332                 if(val[i] > thres) {
 333                         res |= 1 << i;
 334                 }
 335         }
 336         return res;
 337 }
 338 #endif  /* USE_MCUBES */
 339
 340
 341 /* ---- marching tetrahedra implementation (incomplete) ---- */
 342 #ifdef USE_MTETRA
 343
 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)
 346
 347
 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)
 352 {
 353         unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
 354
 355         switch(code) {
 356         case 1:
 357         case INV(1):
 358                 EDGE(0, 1);
 359                 EDGE(0, 2);
 360                 EDGE(0, 3);
 361                 break;
 362
 363         case 2:
 364         case INV(2):
 365                 EDGE(1, 0);
 366                 EDGE(1, 3);
 367                 EDGE(1, 2);
 368                 break;
 369
 370         case 3:
 371         case INV(3):
 372                 EDGE(0, 3);
 373                 EDGE(0, 2);
 374                 EDGE(1, 3);
 375
 376                 EDGE(1, 3);
 377                 EDGE(1, 2);
 378                 EDGE(0, 2);
 379                 break;
 380
 381         case 4:
 382         case INV(4):
 383                 EDGE(2, 0);
 384                 EDGE(2, 1);
 385                 EDGE(2, 3);
 386                 break;
 387
 388         case 5:
 389         case INV(5):
 390                 EDGE(0, 1);
 391                 EDGE(2, 3);
 392                 EDGE(0, 3);
 393
 394                 EDGE(0, 1);
 395                 EDGE(1, 2);
 396                 EDGE(2, 3);
 397                 break;
 398
 399         case 6:
 400         case INV(6):
 401                 EDGE(0, 1);
 402                 EDGE(1, 3);
 403                 EDGE(2, 3);
 404
 405                 EDGE(0, 1);
 406                 EDGE(0, 2);
 407                 EDGE(2, 3);
 408                 break;
 409
 410         case 7:
 411         case INV(7):
 412                 EDGE(3, 0);
 413                 EDGE(3, 2);
 414                 EDGE(3, 1);
 415                 break;
 416
 417         default:
 418                 break;  /* cases 0 and 15 */
 419         }
 420 }
 421
 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)
 424 {
 425         return BIT(0) | BIT(1) | BIT(2) | BIT(3);
 426 }
 427
 428 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
 429 {
 430         int i;
 431         float t = (ms->thres - v0) / (v1 - v0);
 432
 433         vec3 p;
 434         for(i=0; i<3; i++) {
 435                 p[i] = p0[i] + (p1[i] - p0[i]) * t;
 436         }
 437         ms->vertex(ms, p[0], p[1], p[2]);
 438
 439         /*for(i=0; i<3; i++) {
 440                 ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
 441         }
 442
 443         if(++ms->nverts >= 3) {
 444                 ms->nverts = 0;
 445
 446                 for(i=0; i<3; i++) {
 447                         int idx = rev ? (2 - i) : i;
 448                         ms->vertex(ms, ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
 449                 }
 450         }*/
 451 }
 452 #endif  /* USE_MTETRA */