metasurf
diff src/metasurf.c @ 0:7aa4627e492b
first commit
author | John Tsiombikas <nuclear@member.fsf.org> |
---|---|
date | Tue, 25 Oct 2011 07:34:31 +0300 |
parents | |
children | 9ab057fba0c5 |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/src/metasurf.c Tue Oct 25 07:34:31 2011 +0300 1.3 @@ -0,0 +1,362 @@ 1.4 +#include <stdlib.h> 1.5 +#include "metasurf.h" 1.6 +#include "mcubes.h" 1.7 + 1.8 +#undef USE_MTETRA 1.9 +#define USE_MCUBES 1.10 + 1.11 +#if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES)) 1.12 +#error "pick either USE_MTETRA or USE_MCUBES, not both..." 1.13 +#endif 1.14 + 1.15 +typedef float vec3[3]; 1.16 + 1.17 +struct metasurface { 1.18 + vec3 min, max; 1.19 + int res[3]; 1.20 + float thres; 1.21 + 1.22 + msurf_eval_func_t eval; 1.23 + msurf_vertex_func_t vertex; 1.24 + msurf_normal_func_t normal; 1.25 + 1.26 + vec3 vbuf[3]; 1.27 + int nverts; 1.28 +}; 1.29 + 1.30 +static int msurf_init(struct metasurface *ms); 1.31 +static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz); 1.32 +#ifdef USE_MTETRA 1.33 +static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val); 1.34 +#endif 1.35 +#ifdef USE_MCUBES 1.36 +static void process_cube(struct metasurface *ms, vec3 *pos, float *val); 1.37 +#endif 1.38 + 1.39 + 1.40 +struct metasurface *msurf_create(void) 1.41 +{ 1.42 + struct metasurface *ms; 1.43 + 1.44 + if(!(ms = malloc(sizeof *ms))) { 1.45 + return 0; 1.46 + } 1.47 + if(msurf_init(ms) == -1) { 1.48 + free(ms); 1.49 + } 1.50 + return ms; 1.51 +} 1.52 + 1.53 +void msurf_free(struct metasurface *ms) 1.54 +{ 1.55 + free(ms); 1.56 +} 1.57 + 1.58 +static int msurf_init(struct metasurface *ms) 1.59 +{ 1.60 + ms->thres = 0.0; 1.61 + ms->eval = 0; 1.62 + ms->vertex = 0; 1.63 + ms->normal = 0; 1.64 + ms->min[0] = ms->min[1] = ms->min[2] = -1.0; 1.65 + ms->max[0] = ms->max[1] = ms->max[2] = 1.0; 1.66 + ms->res[0] = ms->res[1] = ms->res[2] = 32; 1.67 + ms->nverts = 0; 1.68 + 1.69 + return 0; 1.70 +} 1.71 + 1.72 +void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func) 1.73 +{ 1.74 + ms->eval = func; 1.75 +} 1.76 + 1.77 +void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func) 1.78 +{ 1.79 + ms->vertex = func; 1.80 +} 1.81 + 1.82 +void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func) 1.83 +{ 1.84 + ms->normal = func; 1.85 +} 1.86 + 1.87 +void msurf_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax) 1.88 +{ 1.89 + ms->min[0] = xmin; 1.90 + ms->min[1] = ymin; 1.91 + ms->min[2] = zmin; 1.92 + ms->max[0] = xmax; 1.93 + ms->max[1] = ymax; 1.94 + ms->max[2] = zmax; 1.95 +} 1.96 + 1.97 +void msurf_resolution(struct metasurface *ms, int xres, int yres, int zres) 1.98 +{ 1.99 + ms->res[0] = xres; 1.100 + ms->res[1] = yres; 1.101 + ms->res[2] = zres; 1.102 +} 1.103 + 1.104 +void msurf_threshold(struct metasurface *ms, float thres) 1.105 +{ 1.106 + ms->thres = thres; 1.107 +} 1.108 + 1.109 + 1.110 +void msurf_polygonize(struct metasurface *ms) 1.111 +{ 1.112 + int i, j, k; 1.113 + vec3 pos, delta; 1.114 + 1.115 + for(i=0; i<3; i++) { 1.116 + delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i]; 1.117 + } 1.118 + 1.119 + pos[0] = ms->min[0]; 1.120 + for(i=0; i<ms->res[0] - 1; i++) { 1.121 + 1.122 + pos[1] = ms->min[1]; 1.123 + for(j=0; j<ms->res[1] - 1; j++) { 1.124 + 1.125 + pos[2] = ms->min[2]; 1.126 + for(k=0; k<ms->res[2] - 1; k++) { 1.127 + 1.128 + process_cell(ms, pos, delta); 1.129 + 1.130 + pos[2] += delta[2]; 1.131 + } 1.132 + pos[1] += delta[1]; 1.133 + } 1.134 + pos[0] += delta[0]; 1.135 + } 1.136 +} 1.137 + 1.138 + 1.139 +static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz) 1.140 +{ 1.141 + int i; 1.142 + vec3 p[8]; 1.143 + float val[8]; 1.144 + 1.145 +#ifdef USE_MTETRA 1.146 + static int tetra[][4] = { 1.147 + {0, 2, 3, 7}, 1.148 + {0, 2, 6, 7}, 1.149 + {0, 4, 6, 7}, 1.150 + {0, 6, 1, 2}, 1.151 + {0, 6, 1, 4}, 1.152 + {5, 6, 1, 4} 1.153 + }; 1.154 +#endif 1.155 + 1.156 + static const float offs[][3] = { 1.157 + {0.0f, 0.0f, 0.0f}, 1.158 + {1.0f, 0.0f, 0.0f}, 1.159 + {1.0f, 1.0f, 0.0f}, 1.160 + {0.0f, 1.0f, 0.0f}, 1.161 + {0.0f, 0.0f, 1.0f}, 1.162 + {1.0f, 0.0f, 1.0f}, 1.163 + {1.0f, 1.0f, 1.0f}, 1.164 + {0.0f, 1.0f, 1.0f} 1.165 + }; 1.166 + 1.167 + for(i=0; i<8; i++) { 1.168 + p[i][0] = pos[0] + sz[0] * offs[i][2]; 1.169 + p[i][1] = pos[1] + sz[1] * offs[i][1]; 1.170 + p[i][2] = pos[2] + sz[2] * offs[i][0]; 1.171 + 1.172 + val[i] = ms->eval(p[i][0], p[i][1], p[i][2]); 1.173 + } 1.174 + 1.175 +#ifdef USE_MTETRA 1.176 + for(i=0; i<6; i++) { 1.177 + process_tetra(ms, tetra[i], p, val); 1.178 + } 1.179 +#endif 1.180 +#ifdef USE_MCUBES 1.181 + process_cube(ms, p, val); 1.182 +#endif 1.183 +} 1.184 + 1.185 + 1.186 +/* ---- marching cubes implementation ---- */ 1.187 +#ifdef USE_MCUBES 1.188 + 1.189 +static unsigned int mc_bitcode(float *val, float thres); 1.190 + 1.191 +static void process_cube(struct metasurface *ms, vec3 *pos, float *val) 1.192 +{ 1.193 + static const int pidx[12][2] = { 1.194 + {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6}, 1.195 + {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7} 1.196 + }; 1.197 + int i, j; 1.198 + vec3 vert[12]; 1.199 + unsigned int code = mc_bitcode(val, ms->thres); 1.200 + 1.201 + if(mc_edge_table[code] == 0) { 1.202 + return; 1.203 + } 1.204 + 1.205 + for(i=0; i<12; i++) { 1.206 + if(mc_edge_table[code] & (1 << i)) { 1.207 + int p0 = pidx[i][0]; 1.208 + int p1 = pidx[i][1]; 1.209 + 1.210 + float t = (ms->thres - val[p0]) / (val[p1] - val[p0]); 1.211 + vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t; 1.212 + vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t; 1.213 + vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t; 1.214 + } 1.215 + } 1.216 + 1.217 + for(i=0; mc_tri_table[code][i] != -1; i+=3) { 1.218 + for(j=0; j<3; j++) { 1.219 + float *v = vert[mc_tri_table[code][i + j]]; 1.220 + ms->vertex(v[0], v[1], v[2]); 1.221 + } 1.222 + } 1.223 +} 1.224 + 1.225 +static unsigned int mc_bitcode(float *val, float thres) 1.226 +{ 1.227 + unsigned int i, res = 0; 1.228 + 1.229 + for(i=0; i<8; i++) { 1.230 + if(val[i] > thres) { 1.231 + res |= 1 << i; 1.232 + } 1.233 + } 1.234 + return res; 1.235 +} 1.236 +#endif /* USE_MCUBES */ 1.237 + 1.238 + 1.239 +/* ---- marching tetrahedra implementation (incomplete) ---- */ 1.240 +#ifdef USE_MTETRA 1.241 + 1.242 +static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres); 1.243 +static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev) 1.244 + 1.245 + 1.246 +#define REVBIT(x) ((x) & 8) 1.247 +#define INV(x) (~(x) & 0xf) 1.248 +#define EDGE(a, b) emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code)) 1.249 +static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val) 1.250 +{ 1.251 + unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres); 1.252 + 1.253 + switch(code) { 1.254 + /*case 1: 1.255 + case INV(1):*/ 1.256 + case 0x0e: 1.257 + case 0x01: 1.258 + EDGE(0, 1); 1.259 + EDGE(0, 2); 1.260 + EDGE(0, 3); 1.261 + break; 1.262 + 1.263 + /*case 2: 1.264 + case INV(2):*/ 1.265 + case 0x0d: 1.266 + case 0x02: 1.267 + EDGE(1, 0); 1.268 + EDGE(1, 3); 1.269 + EDGE(1, 2); 1.270 + break; 1.271 + 1.272 + /*case 3: 1.273 + case INV(3):*/ 1.274 + case 0x0c: 1.275 + case 0x03: 1.276 + EDGE(0, 3); 1.277 + EDGE(0, 2); 1.278 + EDGE(1, 3); 1.279 + 1.280 + EDGE(1, 3); 1.281 + EDGE(1, 2); 1.282 + EDGE(0, 2); 1.283 + break; 1.284 + 1.285 + /*case 4: 1.286 + case INV(4):*/ 1.287 + case 0x0b: 1.288 + case 0x04: 1.289 + EDGE(2, 0); 1.290 + EDGE(2, 1); 1.291 + EDGE(2, 3); 1.292 + break; 1.293 + 1.294 + /*case 5: 1.295 + case INV(5):*/ 1.296 + case 0x0a: 1.297 + case 0x05: 1.298 + EDGE(0, 1); 1.299 + EDGE(2, 3); 1.300 + EDGE(0, 3); 1.301 + 1.302 + EDGE(0, 1); 1.303 + EDGE(1, 2); 1.304 + EDGE(2, 3); 1.305 + break; 1.306 + 1.307 + /*case 6: 1.308 + case INV(6):*/ 1.309 + case 0x09: 1.310 + case 0x06: 1.311 + EDGE(0, 1); 1.312 + EDGE(1, 3); 1.313 + EDGE(2, 3); 1.314 + 1.315 + EDGE(0, 1); 1.316 + EDGE(0, 2); 1.317 + EDGE(2, 3); 1.318 + break; 1.319 + 1.320 + /*case 7: 1.321 + case INV(7):*/ 1.322 + case 0x07: 1.323 + case 0x08: 1.324 + EDGE(3, 0); 1.325 + EDGE(3, 2); 1.326 + EDGE(3, 1); 1.327 + break; 1.328 + 1.329 + default: 1.330 + break; /* cases 0 and 15 */ 1.331 + } 1.332 +} 1.333 + 1.334 +#define BIT(i) ((v##i > thres) ? (1 << i) : 0) 1.335 +static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres) 1.336 +{ 1.337 + return BIT(0) | BIT(1) | BIT(2) | BIT(3); 1.338 +} 1.339 + 1.340 +static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev) 1.341 +{ 1.342 + int i; 1.343 + float t = (ms->thres - v0) / (v1 - v0); 1.344 + 1.345 + vec3 p; 1.346 + for(i=0; i<3; i++) { 1.347 + p[i] = p0[i] + (p1[i] - p0[i]) * t; 1.348 + } 1.349 + ms->vertex(p[0], p[1], p[2]); 1.350 + 1.351 + /*for(i=0; i<3; i++) { 1.352 + ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t; 1.353 + } 1.354 + 1.355 + if(++ms->nverts >= 3) { 1.356 + ms->nverts = 0; 1.357 + 1.358 + for(i=0; i<3; i++) { 1.359 + int idx = rev ? (2 - i) : i; 1.360 + ms->vertex(ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]); 1.361 + } 1.362 + }*/ 1.363 +} 1.364 + 1.365 +#endif /* USE_MTETRA */