1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/libs/metasurf/metasurf.c	Mon Feb 18 03:46:52 2013 +0200
     1.3 @@ -0,0 +1,413 @@
     1.4 +/*
     1.5 +metasurf - a library for implicit surface polygonization
     1.6 +Copyright (C) 2011  John Tsiombikas <nuclear@member.fsf.org>
     1.7 +
     1.8 +This program is free software: you can redistribute it and/or modify
     1.9 +it under the terms of the GNU Lesser General Public License as published by
    1.10 +the Free Software Foundation, either version 3 of the License, or
    1.11 +(at your option) any later version.
    1.12 +
    1.13 +This program is distributed in the hope that it will be useful,
    1.14 +but WITHOUT ANY WARRANTY; without even the implied warranty of
    1.15 +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    1.16 +GNU Lesser General Public License for more details.
    1.17 +
    1.18 +You should have received a copy of the GNU Lesser General Public License
    1.19 +along with this program.  If not, see <http://www.gnu.org/licenses/>.
    1.20 +*/
    1.21 +#include <stdio.h>
    1.22 +#include <stdlib.h>
    1.23 +#include "metasurf.h"
    1.24 +#include "mcubes.h"
    1.25 +
    1.26 +#undef USE_MTETRA
    1.27 +#define USE_MCUBES
    1.28 +
    1.29 +#if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))
    1.30 +#error "pick either USE_MTETRA or USE_MCUBES, not both..."
    1.31 +#endif
    1.32 +
    1.33 +typedef float vec3[3];
    1.34 +
    1.35 +struct metasurface {
    1.36 +	vec3 min, max;
    1.37 +	int res[3];
    1.38 +	float thres;
    1.39 +
    1.40 +	msurf_eval_func_t eval;
    1.41 +	msurf_vertex_func_t vertex;
    1.42 +	msurf_normal_func_t normal;
    1.43 +
    1.44 +	float dx, dy, dz;
    1.45 +	int flip;
    1.46 +
    1.47 +	vec3 vbuf[3];
    1.48 +	int nverts;
    1.49 +};
    1.50 +
    1.51 +static int msurf_init(struct metasurface *ms);
    1.52 +static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
    1.53 +#ifdef USE_MTETRA
    1.54 +static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
    1.55 +#endif
    1.56 +#ifdef USE_MCUBES
    1.57 +static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
    1.58 +#endif
    1.59 +
    1.60 +
    1.61 +struct metasurface *msurf_create(void)
    1.62 +{
    1.63 +	struct metasurface *ms;
    1.64 +
    1.65 +	if(!(ms = malloc(sizeof *ms))) {
    1.66 +		return 0;
    1.67 +	}
    1.68 +	if(msurf_init(ms) == -1) {
    1.69 +		free(ms);
    1.70 +	}
    1.71 +	return ms;
    1.72 +}
    1.73 +
    1.74 +void msurf_free(struct metasurface *ms)
    1.75 +{
    1.76 +	free(ms);
    1.77 +}
    1.78 +
    1.79 +static int msurf_init(struct metasurface *ms)
    1.80 +{
    1.81 +	ms->thres = 0.0;
    1.82 +	ms->eval = 0;
    1.83 +	ms->vertex = 0;
    1.84 +	ms->normal = 0;
    1.85 +	ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
    1.86 +	ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
    1.87 +	ms->res[0] = ms->res[1] = ms->res[2] = 40;
    1.88 +	ms->nverts = 0;
    1.89 +
    1.90 +	ms->dx = ms->dy = ms->dz = 0.001;
    1.91 +	ms->flip = 0;
    1.92 +
    1.93 +	return 0;
    1.94 +}
    1.95 +
    1.96 +void msurf_inside(struct metasurface *ms, int inside)
    1.97 +{
    1.98 +	switch(inside) {
    1.99 +	case MSURF_GREATER:
   1.100 +		ms->flip = 0;
   1.101 +		break;
   1.102 +
   1.103 +	case MSURF_LESS:
   1.104 +		ms->flip = 1;
   1.105 +		break;
   1.106 +
   1.107 +	default:
   1.108 +		fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
   1.109 +	}
   1.110 +}
   1.111 +
   1.112 +void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
   1.113 +{
   1.114 +	ms->eval = func;
   1.115 +}
   1.116 +
   1.117 +void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
   1.118 +{
   1.119 +	ms->vertex = func;
   1.120 +}
   1.121 +
   1.122 +void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
   1.123 +{
   1.124 +	ms->normal = func;
   1.125 +}
   1.126 +
   1.127 +void msurf_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
   1.128 +{
   1.129 +	ms->min[0] = xmin;
   1.130 +	ms->min[1] = ymin;
   1.131 +	ms->min[2] = zmin;
   1.132 +	ms->max[0] = xmax;
   1.133 +	ms->max[1] = ymax;
   1.134 +	ms->max[2] = zmax;
   1.135 +}
   1.136 +
   1.137 +void msurf_resolution(struct metasurface *ms, int xres, int yres, int zres)
   1.138 +{
   1.139 +	ms->res[0] = xres;
   1.140 +	ms->res[1] = yres;
   1.141 +	ms->res[2] = zres;
   1.142 +}
   1.143 +
   1.144 +void msurf_threshold(struct metasurface *ms, float thres)
   1.145 +{
   1.146 +	ms->thres = thres;
   1.147 +}
   1.148 +
   1.149 +
   1.150 +int msurf_polygonize(struct metasurface *ms)
   1.151 +{
   1.152 +	int i, j, k;
   1.153 +	vec3 pos, delta;
   1.154 +
   1.155 +	if(!ms->eval || !ms->vertex) {
   1.156 +		fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
   1.157 +		return -1;
   1.158 +	}
   1.159 +
   1.160 +	for(i=0; i<3; i++) {
   1.161 +		delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
   1.162 +	}
   1.163 +
   1.164 +	pos[0] = ms->min[0];
   1.165 +	for(i=0; i<ms->res[0] - 1; i++) {
   1.166 +
   1.167 +		pos[1] = ms->min[1];
   1.168 +		for(j=0; j<ms->res[1] - 1; j++) {
   1.169 +
   1.170 +			pos[2] = ms->min[2];
   1.171 +			for(k=0; k<ms->res[2] - 1; k++) {
   1.172 +
   1.173 +				process_cell(ms, pos, delta);
   1.174 +
   1.175 +				pos[2] += delta[2];
   1.176 +			}
   1.177 +			pos[1] += delta[1];
   1.178 +		}
   1.179 +		pos[0] += delta[0];
   1.180 +	}
   1.181 +	return 0;
   1.182 +}
   1.183 +
   1.184 +
   1.185 +static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
   1.186 +{
   1.187 +	int i;
   1.188 +	vec3 p[8];
   1.189 +	float val[8];
   1.190 +
   1.191 +#ifdef USE_MTETRA
   1.192 +	static int tetra[][4] = {
   1.193 +		{0, 2, 3, 7},
   1.194 +		{0, 2, 6, 7},
   1.195 +		{0, 4, 6, 7},
   1.196 +		{0, 6, 1, 2},
   1.197 +		{0, 6, 1, 4},
   1.198 +		{5, 6, 1, 4}
   1.199 +	};
   1.200 +#endif
   1.201 +
   1.202 +	static const float offs[][3] = {
   1.203 +		{0.0f, 0.0f, 0.0f},
   1.204 +		{1.0f, 0.0f, 0.0f},
   1.205 +		{1.0f, 1.0f, 0.0f},
   1.206 +		{0.0f, 1.0f, 0.0f},
   1.207 +		{0.0f, 0.0f, 1.0f},
   1.208 +		{1.0f, 0.0f, 1.0f},
   1.209 +		{1.0f, 1.0f, 1.0f},
   1.210 +		{0.0f, 1.0f, 1.0f}
   1.211 +	};
   1.212 +
   1.213 +	for(i=0; i<8; i++) {
   1.214 +		p[i][0] = pos[0] + sz[0] * offs[i][2];
   1.215 +		p[i][1] = pos[1] + sz[1] * offs[i][1];
   1.216 +		p[i][2] = pos[2] + sz[2] * offs[i][0];
   1.217 +
   1.218 +		val[i] = ms->eval(p[i][0], p[i][1], p[i][2]);
   1.219 +	}
   1.220 +
   1.221 +#ifdef USE_MTETRA
   1.222 +	for(i=0; i<6; i++) {
   1.223 +		process_tetra(ms, tetra[i], p, val);
   1.224 +	}
   1.225 +#endif
   1.226 +#ifdef USE_MCUBES
   1.227 +	process_cube(ms, p, val);
   1.228 +#endif
   1.229 +}
   1.230 +
   1.231 +
   1.232 +/* ---- marching cubes implementation ---- */
   1.233 +#ifdef USE_MCUBES
   1.234 +
   1.235 +static unsigned int mc_bitcode(float *val, float thres);
   1.236 +
   1.237 +static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
   1.238 +{
   1.239 +	static const int pidx[12][2] = {
   1.240 +		{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
   1.241 +		{6, 7},	{7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
   1.242 +	};
   1.243 +	int i, j;
   1.244 +	vec3 vert[12];
   1.245 +	unsigned int code = mc_bitcode(val, ms->thres);
   1.246 +
   1.247 +	if(ms->flip) {
   1.248 +		code = ~code & 0xff;
   1.249 +	}
   1.250 +
   1.251 +	if(mc_edge_table[code] == 0) {
   1.252 +		return;
   1.253 +	}
   1.254 +
   1.255 +	for(i=0; i<12; i++) {
   1.256 +		if(mc_edge_table[code] & (1 << i)) {
   1.257 +			int p0 = pidx[i][0];
   1.258 +			int p1 = pidx[i][1];
   1.259 +
   1.260 +			float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
   1.261 +			vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
   1.262 +			vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
   1.263 +			vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
   1.264 +		}
   1.265 +	}
   1.266 +
   1.267 +	for(i=0; mc_tri_table[code][i] != -1; i+=3) {
   1.268 +		for(j=0; j<3; j++) {
   1.269 +			float *v = vert[mc_tri_table[code][i + j]];
   1.270 +
   1.271 +			if(ms->normal) {
   1.272 +				float dfdx, dfdy, dfdz;
   1.273 +				dfdx = ms->eval(v[0] - ms->dx, v[1], v[2]) - ms->eval(v[0] + ms->dx, v[1], v[2]);
   1.274 +				dfdy = ms->eval(v[0], v[1] - ms->dy, v[2]) - ms->eval(v[0], v[1] + ms->dy, v[2]);
   1.275 +				dfdz = ms->eval(v[0], v[1], v[2] - ms->dz) - ms->eval(v[0], v[1], v[2] + ms->dz);
   1.276 +
   1.277 +				if(ms->flip) {
   1.278 +					dfdx = -dfdx;
   1.279 +					dfdy = -dfdy;
   1.280 +					dfdz = -dfdz;
   1.281 +				}
   1.282 +				ms->normal(dfdx, dfdy, dfdz);
   1.283 +			}
   1.284 +
   1.285 +			ms->vertex(v[0], v[1], v[2]);
   1.286 +		}
   1.287 +	}
   1.288 +}
   1.289 +
   1.290 +static unsigned int mc_bitcode(float *val, float thres)
   1.291 +{
   1.292 +	unsigned int i, res = 0;
   1.293 +
   1.294 +	for(i=0; i<8; i++) {
   1.295 +		if(val[i] > thres) {
   1.296 +			res |= 1 << i;
   1.297 +		}
   1.298 +	}
   1.299 +	return res;
   1.300 +}
   1.301 +#endif	/* USE_MCUBES */
   1.302 +
   1.303 +
   1.304 +/* ---- marching tetrahedra implementation (incomplete) ---- */
   1.305 +#ifdef USE_MTETRA
   1.306 +
   1.307 +static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
   1.308 +static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
   1.309 +
   1.310 +
   1.311 +#define REVBIT(x)	((x) & 8)
   1.312 +#define INV(x)		(~(x) & 0xf)
   1.313 +#define EDGE(a, b)	emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
   1.314 +static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
   1.315 +{
   1.316 +	unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
   1.317 +
   1.318 +	switch(code) {
   1.319 +	case 1:
   1.320 +	case INV(1):
   1.321 +		EDGE(0, 1);
   1.322 +		EDGE(0, 2);
   1.323 +		EDGE(0, 3);
   1.324 +		break;
   1.325 +
   1.326 +	case 2:
   1.327 +	case INV(2):
   1.328 +		EDGE(1, 0);
   1.329 +		EDGE(1, 3);
   1.330 +		EDGE(1, 2);
   1.331 +		break;
   1.332 +
   1.333 +	case 3:
   1.334 +	case INV(3):
   1.335 +		EDGE(0, 3);
   1.336 +		EDGE(0, 2);
   1.337 +		EDGE(1, 3);
   1.338 +
   1.339 +		EDGE(1, 3);
   1.340 +		EDGE(1, 2);
   1.341 +		EDGE(0, 2);
   1.342 +		break;
   1.343 +
   1.344 +	case 4:
   1.345 +	case INV(4):
   1.346 +		EDGE(2, 0);
   1.347 +		EDGE(2, 1);
   1.348 +		EDGE(2, 3);
   1.349 +		break;
   1.350 +
   1.351 +	case 5:
   1.352 +	case INV(5):
   1.353 +		EDGE(0, 1);
   1.354 +		EDGE(2, 3);
   1.355 +		EDGE(0, 3);
   1.356 +
   1.357 +		EDGE(0, 1);
   1.358 +		EDGE(1, 2);
   1.359 +		EDGE(2, 3);
   1.360 +		break;
   1.361 +
   1.362 +	case 6:
   1.363 +	case INV(6):
   1.364 +		EDGE(0, 1);
   1.365 +		EDGE(1, 3);
   1.366 +		EDGE(2, 3);
   1.367 +
   1.368 +		EDGE(0, 1);
   1.369 +		EDGE(0, 2);
   1.370 +		EDGE(2, 3);
   1.371 +		break;
   1.372 +
   1.373 +	case 7:
   1.374 +	case INV(7):
   1.375 +		EDGE(3, 0);
   1.376 +		EDGE(3, 2);
   1.377 +		EDGE(3, 1);
   1.378 +		break;
   1.379 +
   1.380 +	default:
   1.381 +		break;	/* cases 0 and 15 */
   1.382 +	}
   1.383 +}
   1.384 +
   1.385 +#define BIT(i)	((v##i > thres) ? (1 << i) : 0)
   1.386 +static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
   1.387 +{
   1.388 +	return BIT(0) | BIT(1) | BIT(2) | BIT(3);
   1.389 +}
   1.390 +
   1.391 +static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
   1.392 +{
   1.393 +	int i;
   1.394 +	float t = (ms->thres - v0) / (v1 - v0);
   1.395 +
   1.396 +	vec3 p;
   1.397 +	for(i=0; i<3; i++) {
   1.398 +		p[i] = p0[i] + (p1[i] - p0[i]) * t;
   1.399 +	}
   1.400 +	ms->vertex(p[0], p[1], p[2]);
   1.401 +
   1.402 +	/*for(i=0; i<3; i++) {
   1.403 +		ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
   1.404 +	}
   1.405 +
   1.406 +	if(++ms->nverts >= 3) {
   1.407 +		ms->nverts = 0;
   1.408 +
   1.409 +		for(i=0; i<3; i++) {
   1.410 +			int idx = rev ? (2 - i) : i;
   1.411 +			ms->vertex(ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
   1.412 +		}
   1.413 +	}*/
   1.414 +}
   1.415 +
   1.416 +#endif	/* USE_MTETRA */