nuclear@0: #include <stdlib.h>
nuclear@0: #include "metasurf.h"
nuclear@0: #include "mcubes.h"
nuclear@0: 
nuclear@0: #undef USE_MTETRA
nuclear@0: #define USE_MCUBES
nuclear@0: 
nuclear@0: #if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))
nuclear@0: #error "pick either USE_MTETRA or USE_MCUBES, not both..."
nuclear@0: #endif
nuclear@0: 
nuclear@0: typedef float vec3[3];
nuclear@0: 
nuclear@0: struct metasurface {
nuclear@0: 	vec3 min, max;
nuclear@0: 	int res[3];
nuclear@0: 	float thres;
nuclear@0: 
nuclear@0: 	msurf_eval_func_t eval;
nuclear@0: 	msurf_vertex_func_t vertex;
nuclear@0: 	msurf_normal_func_t normal;
nuclear@0: 
nuclear@0: 	vec3 vbuf[3];
nuclear@0: 	int nverts;
nuclear@0: };
nuclear@0: 
nuclear@0: static int msurf_init(struct metasurface *ms);
nuclear@0: static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
nuclear@0: #ifdef USE_MTETRA
nuclear@0: static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
nuclear@0: #endif
nuclear@0: #ifdef USE_MCUBES
nuclear@0: static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
nuclear@0: #endif
nuclear@0: 
nuclear@0: 
nuclear@0: struct metasurface *msurf_create(void)
nuclear@0: {
nuclear@0: 	struct metasurface *ms;
nuclear@0: 
nuclear@0: 	if(!(ms = malloc(sizeof *ms))) {
nuclear@0: 		return 0;
nuclear@0: 	}
nuclear@0: 	if(msurf_init(ms) == -1) {
nuclear@0: 		free(ms);
nuclear@0: 	}
nuclear@0: 	return ms;
nuclear@0: }
nuclear@0: 
nuclear@0: void msurf_free(struct metasurface *ms)
nuclear@0: {
nuclear@0: 	free(ms);
nuclear@0: }
nuclear@0: 
nuclear@0: static int msurf_init(struct metasurface *ms)
nuclear@0: {
nuclear@0: 	ms->thres = 0.0;
nuclear@0: 	ms->eval = 0;
nuclear@0: 	ms->vertex = 0;
nuclear@0: 	ms->normal = 0;
nuclear@0: 	ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
nuclear@0: 	ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
nuclear@0: 	ms->res[0] = ms->res[1] = ms->res[2] = 32;
nuclear@0: 	ms->nverts = 0;
nuclear@0: 
nuclear@0: 	return 0;
nuclear@0: }
nuclear@0: 
nuclear@0: void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
nuclear@0: {
nuclear@0: 	ms->eval = func;
nuclear@0: }
nuclear@0: 
nuclear@0: void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
nuclear@0: {
nuclear@0: 	ms->vertex = func;
nuclear@0: }
nuclear@0: 
nuclear@0: void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
nuclear@0: {
nuclear@0: 	ms->normal = func;
nuclear@0: }
nuclear@0: 
nuclear@0: void msurf_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
nuclear@0: {
nuclear@0: 	ms->min[0] = xmin;
nuclear@0: 	ms->min[1] = ymin;
nuclear@0: 	ms->min[2] = zmin;
nuclear@0: 	ms->max[0] = xmax;
nuclear@0: 	ms->max[1] = ymax;
nuclear@0: 	ms->max[2] = zmax;
nuclear@0: }
nuclear@0: 
nuclear@0: void msurf_resolution(struct metasurface *ms, int xres, int yres, int zres)
nuclear@0: {
nuclear@0: 	ms->res[0] = xres;
nuclear@0: 	ms->res[1] = yres;
nuclear@0: 	ms->res[2] = zres;
nuclear@0: }
nuclear@0: 
nuclear@0: void msurf_threshold(struct metasurface *ms, float thres)
nuclear@0: {
nuclear@0: 	ms->thres = thres;
nuclear@0: }
nuclear@0: 
nuclear@0: 
nuclear@0: void msurf_polygonize(struct metasurface *ms)
nuclear@0: {
nuclear@0: 	int i, j, k;
nuclear@0: 	vec3 pos, delta;
nuclear@0: 
nuclear@0: 	for(i=0; i<3; i++) {
nuclear@0: 		delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	pos[0] = ms->min[0];
nuclear@0: 	for(i=0; i<ms->res[0] - 1; i++) {
nuclear@0: 
nuclear@0: 		pos[1] = ms->min[1];
nuclear@0: 		for(j=0; j<ms->res[1] - 1; j++) {
nuclear@0: 
nuclear@0: 			pos[2] = ms->min[2];
nuclear@0: 			for(k=0; k<ms->res[2] - 1; k++) {
nuclear@0: 
nuclear@0: 				process_cell(ms, pos, delta);
nuclear@0: 
nuclear@0: 				pos[2] += delta[2];
nuclear@0: 			}
nuclear@0: 			pos[1] += delta[1];
nuclear@0: 		}
nuclear@0: 		pos[0] += delta[0];
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: 
nuclear@0: static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
nuclear@0: {
nuclear@0: 	int i;
nuclear@0: 	vec3 p[8];
nuclear@0: 	float val[8];
nuclear@0: 
nuclear@0: #ifdef USE_MTETRA
nuclear@0: 	static int tetra[][4] = {
nuclear@0: 		{0, 2, 3, 7},
nuclear@0: 		{0, 2, 6, 7},
nuclear@0: 		{0, 4, 6, 7},
nuclear@0: 		{0, 6, 1, 2},
nuclear@0: 		{0, 6, 1, 4},
nuclear@0: 		{5, 6, 1, 4}
nuclear@0: 	};
nuclear@0: #endif
nuclear@0: 
nuclear@0: 	static const float offs[][3] = {
nuclear@0: 		{0.0f, 0.0f, 0.0f},
nuclear@0: 		{1.0f, 0.0f, 0.0f},
nuclear@0: 		{1.0f, 1.0f, 0.0f},
nuclear@0: 		{0.0f, 1.0f, 0.0f},
nuclear@0: 		{0.0f, 0.0f, 1.0f},
nuclear@0: 		{1.0f, 0.0f, 1.0f},
nuclear@0: 		{1.0f, 1.0f, 1.0f},
nuclear@0: 		{0.0f, 1.0f, 1.0f}
nuclear@0: 	};
nuclear@0: 
nuclear@0: 	for(i=0; i<8; i++) {
nuclear@0: 		p[i][0] = pos[0] + sz[0] * offs[i][2];
nuclear@0: 		p[i][1] = pos[1] + sz[1] * offs[i][1];
nuclear@0: 		p[i][2] = pos[2] + sz[2] * offs[i][0];
nuclear@0: 
nuclear@0: 		val[i] = ms->eval(p[i][0], p[i][1], p[i][2]);
nuclear@0: 	}
nuclear@0: 
nuclear@0: #ifdef USE_MTETRA
nuclear@0: 	for(i=0; i<6; i++) {
nuclear@0: 		process_tetra(ms, tetra[i], p, val);
nuclear@0: 	}
nuclear@0: #endif
nuclear@0: #ifdef USE_MCUBES
nuclear@0: 	process_cube(ms, p, val);
nuclear@0: #endif
nuclear@0: }
nuclear@0: 
nuclear@0: 
nuclear@0: /* ---- marching cubes implementation ---- */
nuclear@0: #ifdef USE_MCUBES
nuclear@0: 
nuclear@0: static unsigned int mc_bitcode(float *val, float thres);
nuclear@0: 
nuclear@0: static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
nuclear@0: {
nuclear@0: 	static const int pidx[12][2] = {
nuclear@0: 		{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
nuclear@0: 		{6, 7},	{7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
nuclear@0: 	};
nuclear@0: 	int i, j;
nuclear@0: 	vec3 vert[12];
nuclear@0: 	unsigned int code = mc_bitcode(val, ms->thres);
nuclear@0: 
nuclear@0: 	if(mc_edge_table[code] == 0) {
nuclear@0: 		return;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	for(i=0; i<12; i++) {
nuclear@0: 		if(mc_edge_table[code] & (1 << i)) {
nuclear@0: 			int p0 = pidx[i][0];
nuclear@0: 			int p1 = pidx[i][1];
nuclear@0: 
nuclear@0: 			float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
nuclear@0: 			vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
nuclear@0: 			vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
nuclear@0: 			vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
nuclear@0: 		}
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	for(i=0; mc_tri_table[code][i] != -1; i+=3) {
nuclear@0: 		for(j=0; j<3; j++) {
nuclear@0: 			float *v = vert[mc_tri_table[code][i + j]];
nuclear@0: 			ms->vertex(v[0], v[1], v[2]);
nuclear@0: 		}
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: static unsigned int mc_bitcode(float *val, float thres)
nuclear@0: {
nuclear@0: 	unsigned int i, res = 0;
nuclear@0: 
nuclear@0: 	for(i=0; i<8; i++) {
nuclear@0: 		if(val[i] > thres) {
nuclear@0: 			res |= 1 << i;
nuclear@0: 		}
nuclear@0: 	}
nuclear@0: 	return res;
nuclear@0: }
nuclear@0: #endif	/* USE_MCUBES */
nuclear@0: 
nuclear@0: 
nuclear@0: /* ---- marching tetrahedra implementation (incomplete) ---- */
nuclear@0: #ifdef USE_MTETRA
nuclear@0: 
nuclear@0: static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
nuclear@0: static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
nuclear@0: 
nuclear@0: 
nuclear@0: #define REVBIT(x)	((x) & 8)
nuclear@0: #define INV(x)		(~(x) & 0xf)
nuclear@0: #define EDGE(a, b)	emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
nuclear@0: static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
nuclear@0: {
nuclear@0: 	unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
nuclear@0: 
nuclear@0: 	switch(code) {
nuclear@0: 	/*case 1:
nuclear@0: 	case INV(1):*/
nuclear@0: 	case 0x0e:
nuclear@0: 	case 0x01:
nuclear@0: 		EDGE(0, 1);
nuclear@0: 		EDGE(0, 2);
nuclear@0: 		EDGE(0, 3);
nuclear@0: 		break;
nuclear@0: 
nuclear@0: 	/*case 2:
nuclear@0: 	case INV(2):*/
nuclear@0: 	case 0x0d:
nuclear@0: 	case 0x02:
nuclear@0: 		EDGE(1, 0);
nuclear@0: 		EDGE(1, 3);
nuclear@0: 		EDGE(1, 2);
nuclear@0: 		break;
nuclear@0: 
nuclear@0: 	/*case 3:
nuclear@0: 	case INV(3):*/
nuclear@0: 	case 0x0c:
nuclear@0: 	case 0x03:
nuclear@0: 		EDGE(0, 3);
nuclear@0: 		EDGE(0, 2);
nuclear@0: 		EDGE(1, 3);
nuclear@0: 
nuclear@0: 		EDGE(1, 3);
nuclear@0: 		EDGE(1, 2);
nuclear@0: 		EDGE(0, 2);
nuclear@0: 		break;
nuclear@0: 
nuclear@0: 	/*case 4:
nuclear@0: 	case INV(4):*/
nuclear@0: 	case 0x0b:
nuclear@0: 	case 0x04:
nuclear@0: 		EDGE(2, 0);
nuclear@0: 		EDGE(2, 1);
nuclear@0: 		EDGE(2, 3);
nuclear@0: 		break;
nuclear@0: 
nuclear@0: 	/*case 5:
nuclear@0: 	case INV(5):*/
nuclear@0: 	case 0x0a:
nuclear@0: 	case 0x05:
nuclear@0: 		EDGE(0, 1);
nuclear@0: 		EDGE(2, 3);
nuclear@0: 		EDGE(0, 3);
nuclear@0: 
nuclear@0: 		EDGE(0, 1);
nuclear@0: 		EDGE(1, 2);
nuclear@0: 		EDGE(2, 3);
nuclear@0: 		break;
nuclear@0: 
nuclear@0: 	/*case 6:
nuclear@0: 	case INV(6):*/
nuclear@0: 	case 0x09:
nuclear@0: 	case 0x06:
nuclear@0: 		EDGE(0, 1);
nuclear@0: 		EDGE(1, 3);
nuclear@0: 		EDGE(2, 3);
nuclear@0: 
nuclear@0: 		EDGE(0, 1);
nuclear@0: 		EDGE(0, 2);
nuclear@0: 		EDGE(2, 3);
nuclear@0: 		break;
nuclear@0: 
nuclear@0: 	/*case 7:
nuclear@0: 	case INV(7):*/
nuclear@0: 	case 0x07:
nuclear@0: 	case 0x08:
nuclear@0: 		EDGE(3, 0);
nuclear@0: 		EDGE(3, 2);
nuclear@0: 		EDGE(3, 1);
nuclear@0: 		break;
nuclear@0: 
nuclear@0: 	default:
nuclear@0: 		break;	/* cases 0 and 15 */
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: #define BIT(i)	((v##i > thres) ? (1 << i) : 0)
nuclear@0: static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
nuclear@0: {
nuclear@0: 	return BIT(0) | BIT(1) | BIT(2) | BIT(3);
nuclear@0: }
nuclear@0: 
nuclear@0: static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
nuclear@0: {
nuclear@0: 	int i;
nuclear@0: 	float t = (ms->thres - v0) / (v1 - v0);
nuclear@0: 
nuclear@0: 	vec3 p;
nuclear@0: 	for(i=0; i<3; i++) {
nuclear@0: 		p[i] = p0[i] + (p1[i] - p0[i]) * t;
nuclear@0: 	}
nuclear@0: 	ms->vertex(p[0], p[1], p[2]);
nuclear@0: 
nuclear@0: 	/*for(i=0; i<3; i++) {
nuclear@0: 		ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	if(++ms->nverts >= 3) {
nuclear@0: 		ms->nverts = 0;
nuclear@0: 
nuclear@0: 		for(i=0; i<3; i++) {
nuclear@0: 			int idx = rev ? (2 - i) : i;
nuclear@0: 			ms->vertex(ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
nuclear@0: 		}
nuclear@0: 	}*/
nuclear@0: }
nuclear@0: 
nuclear@0: #endif	/* USE_MTETRA */