metasurf

annotate src/metasurf.c @ 3:52664d3451ad

find changesets by author, revision, files, or words in the commit message
added volume rendering example
author John Tsiombikas <nuclear@member.fsf.org>
date Tue, 25 Oct 2011 13:30:03 +0300
parents 7aa4627e492b
children 2c575855f707
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@0 41 vec3 vbuf[3];
nuclear@0 42 int nverts;
nuclear@0 43 };
nuclear@0 44
nuclear@0 45 static int msurf_init(struct metasurface *ms);
nuclear@0 46 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
nuclear@0 47 #ifdef USE_MTETRA
nuclear@0 48 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
nuclear@0 49 #endif
nuclear@0 50 #ifdef USE_MCUBES
nuclear@0 51 static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
nuclear@0 52 #endif
nuclear@0 53
nuclear@0 54
nuclear@0 55 struct metasurface *msurf_create(void)
nuclear@0 56 {
nuclear@0 57 struct metasurface *ms;
nuclear@0 58
nuclear@0 59 if(!(ms = malloc(sizeof *ms))) {
nuclear@0 60 return 0;
nuclear@0 61 }
nuclear@0 62 if(msurf_init(ms) == -1) {
nuclear@0 63 free(ms);
nuclear@0 64 }
nuclear@0 65 return ms;
nuclear@0 66 }
nuclear@0 67
nuclear@0 68 void msurf_free(struct metasurface *ms)
nuclear@0 69 {
nuclear@0 70 free(ms);
nuclear@0 71 }
nuclear@0 72
nuclear@0 73 static int msurf_init(struct metasurface *ms)
nuclear@0 74 {
nuclear@0 75 ms->thres = 0.0;
nuclear@0 76 ms->eval = 0;
nuclear@0 77 ms->vertex = 0;
nuclear@0 78 ms->normal = 0;
nuclear@0 79 ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
nuclear@0 80 ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
nuclear@2 81 ms->res[0] = ms->res[1] = ms->res[2] = 40;
nuclear@0 82 ms->nverts = 0;
nuclear@0 83
nuclear@0 84 return 0;
nuclear@0 85 }
nuclear@0 86
nuclear@0 87 void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
nuclear@0 88 {
nuclear@0 89 ms->eval = func;
nuclear@0 90 }
nuclear@0 91
nuclear@0 92 void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
nuclear@0 93 {
nuclear@0 94 ms->vertex = func;
nuclear@0 95 }
nuclear@0 96
nuclear@0 97 void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
nuclear@0 98 {
nuclear@0 99 ms->normal = func;
nuclear@0 100 }
nuclear@0 101
nuclear@0 102 void msurf_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
nuclear@0 103 {
nuclear@0 104 ms->min[0] = xmin;
nuclear@0 105 ms->min[1] = ymin;
nuclear@0 106 ms->min[2] = zmin;
nuclear@0 107 ms->max[0] = xmax;
nuclear@0 108 ms->max[1] = ymax;
nuclear@0 109 ms->max[2] = zmax;
nuclear@0 110 }
nuclear@0 111
nuclear@0 112 void msurf_resolution(struct metasurface *ms, int xres, int yres, int zres)
nuclear@0 113 {
nuclear@0 114 ms->res[0] = xres;
nuclear@0 115 ms->res[1] = yres;
nuclear@0 116 ms->res[2] = zres;
nuclear@0 117 }
nuclear@0 118
nuclear@0 119 void msurf_threshold(struct metasurface *ms, float thres)
nuclear@0 120 {
nuclear@0 121 ms->thres = thres;
nuclear@0 122 }
nuclear@0 123
nuclear@0 124
nuclear@2 125 int msurf_polygonize(struct metasurface *ms)
nuclear@0 126 {
nuclear@0 127 int i, j, k;
nuclear@0 128 vec3 pos, delta;
nuclear@0 129
nuclear@2 130 if(!ms->eval || !ms->vertex) {
nuclear@2 131 fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
nuclear@2 132 return -1;
nuclear@2 133 }
nuclear@2 134
nuclear@0 135 for(i=0; i<3; i++) {
nuclear@0 136 delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
nuclear@0 137 }
nuclear@0 138
nuclear@0 139 pos[0] = ms->min[0];
nuclear@0 140 for(i=0; i<ms->res[0] - 1; i++) {
nuclear@0 141
nuclear@0 142 pos[1] = ms->min[1];
nuclear@0 143 for(j=0; j<ms->res[1] - 1; j++) {
nuclear@0 144
nuclear@0 145 pos[2] = ms->min[2];
nuclear@0 146 for(k=0; k<ms->res[2] - 1; k++) {
nuclear@0 147
nuclear@0 148 process_cell(ms, pos, delta);
nuclear@0 149
nuclear@0 150 pos[2] += delta[2];
nuclear@0 151 }
nuclear@0 152 pos[1] += delta[1];
nuclear@0 153 }
nuclear@0 154 pos[0] += delta[0];
nuclear@0 155 }
nuclear@2 156 return 0;
nuclear@0 157 }
nuclear@0 158
nuclear@0 159
nuclear@0 160 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
nuclear@0 161 {
nuclear@0 162 int i;
nuclear@0 163 vec3 p[8];
nuclear@0 164 float val[8];
nuclear@0 165
nuclear@0 166 #ifdef USE_MTETRA
nuclear@0 167 static int tetra[][4] = {
nuclear@0 168 {0, 2, 3, 7},
nuclear@0 169 {0, 2, 6, 7},
nuclear@0 170 {0, 4, 6, 7},
nuclear@0 171 {0, 6, 1, 2},
nuclear@0 172 {0, 6, 1, 4},
nuclear@0 173 {5, 6, 1, 4}
nuclear@0 174 };
nuclear@0 175 #endif
nuclear@0 176
nuclear@0 177 static const float offs[][3] = {
nuclear@0 178 {0.0f, 0.0f, 0.0f},
nuclear@0 179 {1.0f, 0.0f, 0.0f},
nuclear@0 180 {1.0f, 1.0f, 0.0f},
nuclear@0 181 {0.0f, 1.0f, 0.0f},
nuclear@0 182 {0.0f, 0.0f, 1.0f},
nuclear@0 183 {1.0f, 0.0f, 1.0f},
nuclear@0 184 {1.0f, 1.0f, 1.0f},
nuclear@0 185 {0.0f, 1.0f, 1.0f}
nuclear@0 186 };
nuclear@0 187
nuclear@0 188 for(i=0; i<8; i++) {
nuclear@0 189 p[i][0] = pos[0] + sz[0] * offs[i][2];
nuclear@0 190 p[i][1] = pos[1] + sz[1] * offs[i][1];
nuclear@0 191 p[i][2] = pos[2] + sz[2] * offs[i][0];
nuclear@0 192
nuclear@0 193 val[i] = ms->eval(p[i][0], p[i][1], p[i][2]);
nuclear@0 194 }
nuclear@0 195
nuclear@0 196 #ifdef USE_MTETRA
nuclear@0 197 for(i=0; i<6; i++) {
nuclear@0 198 process_tetra(ms, tetra[i], p, val);
nuclear@0 199 }
nuclear@0 200 #endif
nuclear@0 201 #ifdef USE_MCUBES
nuclear@0 202 process_cube(ms, p, val);
nuclear@0 203 #endif
nuclear@0 204 }
nuclear@0 205
nuclear@0 206
nuclear@0 207 /* ---- marching cubes implementation ---- */
nuclear@0 208 #ifdef USE_MCUBES
nuclear@0 209
nuclear@0 210 static unsigned int mc_bitcode(float *val, float thres);
nuclear@0 211
nuclear@0 212 static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
nuclear@0 213 {
nuclear@0 214 static const int pidx[12][2] = {
nuclear@0 215 {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
nuclear@0 216 {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
nuclear@0 217 };
nuclear@0 218 int i, j;
nuclear@0 219 vec3 vert[12];
nuclear@0 220 unsigned int code = mc_bitcode(val, ms->thres);
nuclear@0 221
nuclear@0 222 if(mc_edge_table[code] == 0) {
nuclear@0 223 return;
nuclear@0 224 }
nuclear@0 225
nuclear@0 226 for(i=0; i<12; i++) {
nuclear@0 227 if(mc_edge_table[code] & (1 << i)) {
nuclear@0 228 int p0 = pidx[i][0];
nuclear@0 229 int p1 = pidx[i][1];
nuclear@0 230
nuclear@0 231 float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
nuclear@0 232 vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
nuclear@0 233 vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
nuclear@0 234 vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
nuclear@0 235 }
nuclear@0 236 }
nuclear@0 237
nuclear@0 238 for(i=0; mc_tri_table[code][i] != -1; i+=3) {
nuclear@0 239 for(j=0; j<3; j++) {
nuclear@0 240 float *v = vert[mc_tri_table[code][i + j]];
nuclear@0 241 ms->vertex(v[0], v[1], v[2]);
nuclear@0 242 }
nuclear@0 243 }
nuclear@0 244 }
nuclear@0 245
nuclear@0 246 static unsigned int mc_bitcode(float *val, float thres)
nuclear@0 247 {
nuclear@0 248 unsigned int i, res = 0;
nuclear@0 249
nuclear@0 250 for(i=0; i<8; i++) {
nuclear@0 251 if(val[i] > thres) {
nuclear@0 252 res |= 1 << i;
nuclear@0 253 }
nuclear@0 254 }
nuclear@0 255 return res;
nuclear@0 256 }
nuclear@0 257 #endif /* USE_MCUBES */
nuclear@0 258
nuclear@0 259
nuclear@0 260 /* ---- marching tetrahedra implementation (incomplete) ---- */
nuclear@0 261 #ifdef USE_MTETRA
nuclear@0 262
nuclear@0 263 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
nuclear@0 264 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
nuclear@0 265
nuclear@0 266
nuclear@0 267 #define REVBIT(x) ((x) & 8)
nuclear@0 268 #define INV(x) (~(x) & 0xf)
nuclear@0 269 #define EDGE(a, b) emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
nuclear@0 270 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
nuclear@0 271 {
nuclear@0 272 unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
nuclear@0 273
nuclear@0 274 switch(code) {
nuclear@2 275 case 1:
nuclear@2 276 case INV(1):
nuclear@0 277 EDGE(0, 1);
nuclear@0 278 EDGE(0, 2);
nuclear@0 279 EDGE(0, 3);
nuclear@0 280 break;
nuclear@0 281
nuclear@2 282 case 2:
nuclear@2 283 case INV(2):
nuclear@0 284 EDGE(1, 0);
nuclear@0 285 EDGE(1, 3);
nuclear@0 286 EDGE(1, 2);
nuclear@0 287 break;
nuclear@0 288
nuclear@2 289 case 3:
nuclear@2 290 case INV(3):
nuclear@0 291 EDGE(0, 3);
nuclear@0 292 EDGE(0, 2);
nuclear@0 293 EDGE(1, 3);
nuclear@0 294
nuclear@0 295 EDGE(1, 3);
nuclear@0 296 EDGE(1, 2);
nuclear@0 297 EDGE(0, 2);
nuclear@0 298 break;
nuclear@0 299
nuclear@2 300 case 4:
nuclear@2 301 case INV(4):
nuclear@0 302 EDGE(2, 0);
nuclear@0 303 EDGE(2, 1);
nuclear@0 304 EDGE(2, 3);
nuclear@0 305 break;
nuclear@0 306
nuclear@2 307 case 5:
nuclear@2 308 case INV(5):
nuclear@0 309 EDGE(0, 1);
nuclear@0 310 EDGE(2, 3);
nuclear@0 311 EDGE(0, 3);
nuclear@0 312
nuclear@0 313 EDGE(0, 1);
nuclear@0 314 EDGE(1, 2);
nuclear@0 315 EDGE(2, 3);
nuclear@0 316 break;
nuclear@0 317
nuclear@2 318 case 6:
nuclear@2 319 case INV(6):
nuclear@0 320 EDGE(0, 1);
nuclear@0 321 EDGE(1, 3);
nuclear@0 322 EDGE(2, 3);
nuclear@0 323
nuclear@0 324 EDGE(0, 1);
nuclear@0 325 EDGE(0, 2);
nuclear@0 326 EDGE(2, 3);
nuclear@0 327 break;
nuclear@0 328
nuclear@2 329 case 7:
nuclear@2 330 case INV(7):
nuclear@0 331 EDGE(3, 0);
nuclear@0 332 EDGE(3, 2);
nuclear@0 333 EDGE(3, 1);
nuclear@0 334 break;
nuclear@0 335
nuclear@0 336 default:
nuclear@0 337 break; /* cases 0 and 15 */
nuclear@0 338 }
nuclear@0 339 }
nuclear@0 340
nuclear@0 341 #define BIT(i) ((v##i > thres) ? (1 << i) : 0)
nuclear@0 342 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
nuclear@0 343 {
nuclear@0 344 return BIT(0) | BIT(1) | BIT(2) | BIT(3);
nuclear@0 345 }
nuclear@0 346
nuclear@0 347 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
nuclear@0 348 {
nuclear@0 349 int i;
nuclear@0 350 float t = (ms->thres - v0) / (v1 - v0);
nuclear@0 351
nuclear@0 352 vec3 p;
nuclear@0 353 for(i=0; i<3; i++) {
nuclear@0 354 p[i] = p0[i] + (p1[i] - p0[i]) * t;
nuclear@0 355 }
nuclear@0 356 ms->vertex(p[0], p[1], p[2]);
nuclear@0 357
nuclear@0 358 /*for(i=0; i<3; i++) {
nuclear@0 359 ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
nuclear@0 360 }
nuclear@0 361
nuclear@0 362 if(++ms->nverts >= 3) {
nuclear@0 363 ms->nverts = 0;
nuclear@0 364
nuclear@0 365 for(i=0; i<3; i++) {
nuclear@0 366 int idx = rev ? (2 - i) : i;
nuclear@0 367 ms->vertex(ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
nuclear@0 368 }
nuclear@0 369 }*/
nuclear@0 370 }
nuclear@0 371
nuclear@0 372 #endif /* USE_MTETRA */