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