git.mutantstargoat.com Git - metatoy/blob - src/metasurf.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include "metasurf.h"
   5 #include "mcubes.h"
   6
   7 typedef float vec3[3];
   8
   9 struct metasurface {
  10         vec3 min, max;
  11         int res[3], newres[3];
  12         float thres;
  13
  14         float dx, dy, dz;
  15         unsigned int flags;
  16
  17         float *voxels;
  18
  19         int varr_size, varr_alloc_size;
  20         float *varr, *narr;
  21 };
  22
  23 static int msurf_init(struct metasurface *ms);
  24 static void process_cell(struct metasurface *ms, int xcell, int ycell, int zcell, vec3 pos, vec3 sz);
  25
  26 static int **mc_tri_table;
  27
  28 static int decompress_tables(void)
  29 {
  30         int i, j, run, nib, *data;
  31         unsigned char *indata;
  32
  33         if(!(mc_tri_table = malloc(256 * 16 * sizeof(int) + 256 * sizeof(void*)))) {
  34                 fprintf(stderr, "decompress_tables: failed to allocate memory\n");
  35                 return -1;
  36         }
  37
  38         indata = tritab_data;
  39         nib = 0;
  40         data = (int*)(mc_tri_table + 256);
  41         for(i=0; i<256; i++) {
  42                 mc_tri_table[i] = data + i * 16;
  43
  44                 run = ((i & 1) == 0) ? tritab_runlen[i/2] & 0xf : tritab_runlen[i/2] >> 4;
  45                 for(j=0; j<16; j++) {
  46                         if(j < run) {
  47                                 mc_tri_table[i][j] = ((nib++ & 1) == 0) ? *indata & 0xf : *indata++ >> 4;
  48                         } else {
  49                                 mc_tri_table[i][j] = -1;
  50                         }
  51                 }
  52         }
  53
  54         return 0;
  55 }
  56
  57 struct metasurface *msurf_create(void)
  58 {
  59         struct metasurface *ms;
  60
  61         if(!mc_tri_table) {
  62                 if(decompress_tables() == -1) {
  63                         return 0;
  64                 }
  65         }
  66
  67         if(!(ms = malloc(sizeof *ms))) {
  68                 return 0;
  69         }
  70         if(msurf_init(ms) == -1) {
  71                 free(ms);
  72         }
  73         return ms;
  74 }
  75
  76 void msurf_free(struct metasurface *ms)
  77 {
  78         if(ms) {
  79                 free(ms->voxels);
  80                 free(ms->varr);
  81                 free(ms->narr);
  82                 free(ms);
  83         }
  84 }
  85
  86 static int msurf_init(struct metasurface *ms)
  87 {
  88         ms->voxels = 0;
  89         ms->thres = 0.0;
  90         ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
  91         ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
  92         ms->res[0] = ms->res[1] = ms->res[2] = 0;
  93         ms->newres[0] = ms->newres[1] = ms->newres[2] = 40;
  94
  95         ms->varr_alloc_size = ms->varr_size = 0;
  96         ms->varr = ms->narr = 0;
  97
  98         ms->dx = ms->dy = ms->dz = 0.001;
  99         ms->flags = 0;
 100
 101         return 0;
 102 }
 103
 104 void msurf_enable(struct metasurface *ms, unsigned int opt)
 105 {
 106         ms->flags |= opt;
 107 }
 108
 109 void msurf_disable(struct metasurface *ms, unsigned int opt)
 110 {
 111         ms->flags &= ~opt;
 112 }
 113
 114 int msurf_is_enabled(struct metasurface *ms, unsigned int opt)
 115 {
 116         return ms->flags & opt;
 117 }
 118
 119 void msurf_set_inside(struct metasurface *ms, int inside)
 120 {
 121         switch(inside) {
 122         case MSURF_GREATER:
 123                 msurf_enable(ms, MSURF_FLIP);
 124                 break;
 125
 126         case MSURF_LESS:
 127                 msurf_disable(ms, MSURF_FLIP);
 128                 break;
 129
 130         default:
 131                 fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
 132         }
 133 }
 134
 135 int msurf_get_inside(struct metasurface *ms)
 136 {
 137         return msurf_is_enabled(ms, MSURF_FLIP) ? MSURF_LESS : MSURF_GREATER;
 138 }
 139
 140 void msurf_set_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
 141 {
 142         ms->min[0] = xmin;
 143         ms->min[1] = ymin;
 144         ms->min[2] = zmin;
 145         ms->max[0] = xmax;
 146         ms->max[1] = ymax;
 147         ms->max[2] = zmax;
 148 }
 149
 150 void msurf_get_bounds(struct metasurface *ms, float *xmin, float *ymin, float *zmin, float *xmax, float *ymax, float *zmax)
 151 {
 152         *xmin = ms->min[0];
 153         *ymin = ms->min[1];
 154         *zmin = ms->min[2];
 155         *xmax = ms->max[0];
 156         *ymax = ms->max[1];
 157         *zmax = ms->max[2];
 158 }
 159
 160 void msurf_set_resolution(struct metasurface *ms, int xres, int yres, int zres)
 161 {
 162         ms->newres[0] = xres;
 163         ms->newres[1] = yres;
 164         ms->newres[2] = zres;
 165 }
 166
 167 void msurf_get_resolution(struct metasurface *ms, int *xres, int *yres, int *zres)
 168 {
 169         *xres = ms->res[0];
 170         *yres = ms->res[1];
 171         *zres = ms->res[2];
 172 }
 173
 174 void msurf_set_threshold(struct metasurface *ms, float thres)
 175 {
 176         ms->thres = thres;
 177 }
 178
 179 float msurf_get_threshold(struct metasurface *ms)
 180 {
 181         return ms->thres;
 182 }
 183
 184
 185 float *msurf_voxels(struct metasurface *ms)
 186 {
 187         if(ms->res[0] != ms->newres[0] || ms->res[1] != ms->newres[1] || ms->res[2] != ms->newres[2]) {
 188                 int count;
 189                 ms->res[0] = ms->newres[0];
 190                 ms->res[1] = ms->newres[1];
 191                 ms->res[2] = ms->newres[2];
 192                 count = ms->res[0] * ms->res[1] * ms->res[2];
 193                 free(ms->voxels);
 194                 if(!(ms->voxels = malloc(count * sizeof *ms->voxels))) {
 195                         return 0;
 196                 }
 197         }
 198         return ms->voxels;
 199 }
 200
 201 float *msurf_slice(struct metasurface *ms, int idx)
 202 {
 203         float *vox = msurf_voxels(ms);
 204         if(!vox) return 0;
 205
 206         return vox + ms->res[0] * ms->res[1] * idx;
 207 }
 208
 209 int msurf_polygonize(struct metasurface *ms)
 210 {
 211         int i, j, k;
 212         vec3 pos, delta;
 213
 214         if(!ms->voxels) return -1;
 215
 216         ms->varr_size = 0;
 217
 218         for(i=0; i<3; i++) {
 219                 delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
 220         }
 221
 222         for(i=0; i<ms->res[0] - 1; i++) {
 223                 for(j=0; j<ms->res[1] - 1; j++) {
 224                         for(k=0; k<ms->res[2] - 1; k++) {
 225
 226                                 pos[0] = ms->min[0] + i * delta[0];
 227                                 pos[1] = ms->min[1] + j * delta[1];
 228                                 pos[2] = ms->min[2] + k * delta[2];
 229
 230                                 process_cell(ms, i, j, k, pos, delta);
 231                         }
 232                 }
 233         }
 234         return 0;
 235 }
 236
 237 int msurf_vertex_count(struct metasurface *ms)
 238 {
 239         return ms->varr_size / 3;
 240 }
 241
 242 float *msurf_vertices(struct metasurface *ms)
 243 {
 244         return ms->varr;
 245 }
 246
 247 float *msurf_normals(struct metasurface *ms)
 248 {
 249         return ms->narr;
 250 }
 251
 252 static unsigned int mc_bitcode(float *val, float thres);
 253
 254 static void process_cell(struct metasurface *ms, int xcell, int ycell, int zcell, vec3 cellpos, vec3 cellsz)
 255 {
 256         int i, j, k, slice_size;
 257         vec3 pos[8];
 258         float dfdx[8], dfdy[8], dfdz[8];
 259         vec3 vert[12], norm[12];
 260         float val[8];
 261         float *cellptr;
 262         unsigned int code;
 263
 264         static const int offs[][3] = {
 265                 {0, 0, 0},
 266                 {1, 0, 0},
 267                 {1, 1, 0},
 268                 {0, 1, 0},
 269                 {0, 0, 1},
 270                 {1, 0, 1},
 271                 {1, 1, 1},
 272                 {0, 1, 1}
 273         };
 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
 280         slice_size = ms->res[0] * ms->res[1];
 281         cellptr = ms->voxels + slice_size * zcell + ms->res[0] * ycell + xcell;
 282
 283 #define GRIDOFFS(x, y, z)       ((z) * slice_size + (y) * ms->res[0] + (x))
 284
 285         for(i=0; i<8; i++) {
 286                 val[i] = cellptr[GRIDOFFS(offs[i][0], offs[i][1], offs[i][2])];
 287         }
 288
 289         code = mc_bitcode(val, ms->thres);
 290         if(ms->flags & MSURF_FLIP) {
 291                 code = ~code & 0xff;
 292         }
 293         if(mc_edge_table[code] == 0) {
 294                 return;
 295         }
 296
 297         /* calculate normals at the 8 corners */
 298         for(i=0; i<8; i++) {
 299                 float *ptr = cellptr + GRIDOFFS(offs[i][0], offs[i][1], offs[i][2]);
 300
 301                 if(xcell < ms->res[0] - 1) {
 302                         dfdx[i] = ptr[GRIDOFFS(1, 0, 0)] - *ptr;
 303                 } else {
 304                         dfdx[i] = *ptr - ptr[GRIDOFFS(-1, 0, 0)];
 305                 }
 306                 if(ycell < ms->res[1] - 1) {
 307                         dfdy[i] = ptr[GRIDOFFS(0, 1, 0)] - *ptr;
 308                 } else {
 309                         dfdy[i] = *ptr - ptr[GRIDOFFS(0, -1, 0)];
 310                 }
 311                 if(zcell < ms->res[2] - 1) {
 312                         dfdz[i] = ptr[GRIDOFFS(0, 0, 1)] - *ptr;
 313                 } else {
 314                         dfdz[i] = *ptr - ptr[GRIDOFFS(0, 0, -1)];
 315                 }
 316         }
 317
 318         /* calculate the world-space position of each corner */
 319         for(i=0; i<8; i++) {
 320                 pos[i][0] = cellpos[0] + cellsz[0] * offs[i][0];
 321                 pos[i][1] = cellpos[1] + cellsz[1] * offs[i][1];
 322                 pos[i][2] = cellpos[2] + cellsz[2] * offs[i][2];
 323         }
 324
 325         /* generate up to a max of 12 vertices per cube. interpolate positions and normals for each one */
 326         for(i=0; i<12; i++) {
 327                 if(mc_edge_table[code] & (1 << i)) {
 328                         float nx, ny, nz;
 329                         int p0 = pidx[i][0];
 330                         int p1 = pidx[i][1];
 331
 332                         float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
 333                         vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
 334                         vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
 335                         vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
 336
 337                         nx = dfdx[p0] + (dfdx[p1] - dfdx[p0]) * t;
 338                         ny = dfdy[p0] + (dfdy[p1] - dfdy[p0]) * t;
 339                         nz = dfdz[p0] + (dfdz[p1] - dfdz[p0]) * t;
 340
 341                         if(ms->flags & MSURF_FLIP) {
 342                                 nx = -nx;
 343                                 ny = -ny;
 344                                 nz = -nz;
 345                         }
 346
 347                         if(ms->flags & MSURF_NORMALIZE) {
 348                                 float len = sqrt(nx * nx + ny * ny + nz * nz);
 349                                 if(len != 0.0) {
 350                                         float s = 1.0f / len;
 351                                         nx *= s;
 352                                         ny *= s;
 353                                         nz *= s;
 354                                 }
 355                         }
 356
 357                         norm[i][0] = nx;
 358                         norm[i][1] = ny;
 359                         norm[i][2] = nz;
 360                 }
 361         }
 362
 363         /* for each triangle of the cube add the appropriate vertices to the vertex buffer */
 364         for(i=0; mc_tri_table[code][i] != -1; i+=3) {
 365                 for(j=0; j<3; j++) {
 366                         int idx = mc_tri_table[code][i + j];
 367                         float *v = vert[idx];
 368                         float *n = norm[idx];
 369
 370                         /* TODO multithreadied polygon emit */
 371                         if(ms->varr_size + 3 > ms->varr_alloc_size) {
 372                                 int newsz = ms->varr_alloc_size ? ms->varr_alloc_size * 2 : 1024;
 373                                 float *new_varr, *new_narr;
 374
 375                                 if(!(new_varr = realloc(ms->varr, newsz * sizeof *new_varr)) ||
 376                                                 !(new_narr = realloc(ms->narr, newsz * sizeof *new_narr))) {
 377                                         free(new_varr);
 378                                         fprintf(stderr, "msurf_polygonize: failed to grow vertex buffers to %d elements\n", newsz);
 379                                         return;
 380                                 }
 381                                 ms->varr = new_varr;
 382                                 ms->narr = new_narr;
 383                                 ms->varr_alloc_size = newsz;
 384                         }
 385
 386                         for(k=0; k<3; k++) {
 387                                 ms->varr[ms->varr_size] = v[k];
 388                                 ms->narr[ms->varr_size] = n[k];
 389                                 ++ms->varr_size;
 390                         }
 391                 }
 392         }
 393 }
 394
 395 static unsigned int mc_bitcode(float *val, float thres)
 396 {
 397         unsigned int i, res = 0;
 398
 399         for(i=0; i<8; i++) {
 400                 if(val[i] > thres) {
 401                         res |= 1 << i;
 402                 }
 403         }
 404         return res;
 405 }