nuclear@5: /*
nuclear@5: metasurf - a library for implicit surface polygonization
nuclear@5: Copyright (C) 2011  John Tsiombikas <nuclear@member.fsf.org>
nuclear@5: 
nuclear@5: This program is free software: you can redistribute it and/or modify
nuclear@5: it under the terms of the GNU Lesser General Public License as published by
nuclear@5: the Free Software Foundation, either version 3 of the License, or
nuclear@5: (at your option) any later version.
nuclear@5: 
nuclear@5: This program is distributed in the hope that it will be useful,
nuclear@5: but WITHOUT ANY WARRANTY; without even the implied warranty of
nuclear@5: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
nuclear@5: GNU Lesser General Public License for more details.
nuclear@5: 
nuclear@5: You should have received a copy of the GNU Lesser General Public License
nuclear@5: along with this program.  If not, see <http://www.gnu.org/licenses/>.
nuclear@5: */
nuclear@5: #include <stdio.h>
nuclear@5: #include <stdlib.h>
nuclear@5: #include "metasurf.h"
nuclear@5: #include "mcubes.h"
nuclear@5: 
nuclear@5: #undef USE_MTETRA
nuclear@5: #define USE_MCUBES
nuclear@5: 
nuclear@5: #if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))
nuclear@5: #error "pick either USE_MTETRA or USE_MCUBES, not both..."
nuclear@5: #endif
nuclear@5: 
nuclear@5: typedef float vec3[3];
nuclear@5: 
nuclear@5: struct metasurface {
nuclear@5: 	vec3 min, max;
nuclear@5: 	int res[3];
nuclear@5: 	float thres;
nuclear@5: 
nuclear@5: 	msurf_eval_func_t eval;
nuclear@5: 	msurf_vertex_func_t vertex;
nuclear@5: 	msurf_normal_func_t normal;
nuclear@5: 
nuclear@5: 	float dx, dy, dz;
nuclear@5: 	int flip;
nuclear@5: 
nuclear@5: 	vec3 vbuf[3];
nuclear@5: 	int nverts;
nuclear@5: };
nuclear@5: 
nuclear@5: static int msurf_init(struct metasurface *ms);
nuclear@5: static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
nuclear@5: #ifdef USE_MTETRA
nuclear@5: static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
nuclear@5: #endif
nuclear@5: #ifdef USE_MCUBES
nuclear@5: static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
nuclear@5: #endif
nuclear@5: 
nuclear@5: 
nuclear@5: struct metasurface *msurf_create(void)
nuclear@5: {
nuclear@5: 	struct metasurface *ms;
nuclear@5: 
nuclear@5: 	if(!(ms = malloc(sizeof *ms))) {
nuclear@5: 		return 0;
nuclear@5: 	}
nuclear@5: 	if(msurf_init(ms) == -1) {
nuclear@5: 		free(ms);
nuclear@5: 	}
nuclear@5: 	return ms;
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_free(struct metasurface *ms)
nuclear@5: {
nuclear@5: 	free(ms);
nuclear@5: }
nuclear@5: 
nuclear@5: static int msurf_init(struct metasurface *ms)
nuclear@5: {
nuclear@5: 	ms->thres = 0.0;
nuclear@5: 	ms->eval = 0;
nuclear@5: 	ms->vertex = 0;
nuclear@5: 	ms->normal = 0;
nuclear@5: 	ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
nuclear@5: 	ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
nuclear@5: 	ms->res[0] = ms->res[1] = ms->res[2] = 40;
nuclear@5: 	ms->nverts = 0;
nuclear@5: 
nuclear@5: 	ms->dx = ms->dy = ms->dz = 0.001;
nuclear@5: 	ms->flip = 0;
nuclear@5: 
nuclear@5: 	return 0;
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_inside(struct metasurface *ms, int inside)
nuclear@5: {
nuclear@5: 	switch(inside) {
nuclear@5: 	case MSURF_GREATER:
nuclear@5: 		ms->flip = 0;
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	case MSURF_LESS:
nuclear@5: 		ms->flip = 1;
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	default:
nuclear@5: 		fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");
nuclear@5: 	}
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
nuclear@5: {
nuclear@5: 	ms->eval = func;
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
nuclear@5: {
nuclear@5: 	ms->vertex = func;
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
nuclear@5: {
nuclear@5: 	ms->normal = func;
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
nuclear@5: {
nuclear@5: 	ms->min[0] = xmin;
nuclear@5: 	ms->min[1] = ymin;
nuclear@5: 	ms->min[2] = zmin;
nuclear@5: 	ms->max[0] = xmax;
nuclear@5: 	ms->max[1] = ymax;
nuclear@5: 	ms->max[2] = zmax;
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_resolution(struct metasurface *ms, int xres, int yres, int zres)
nuclear@5: {
nuclear@5: 	ms->res[0] = xres;
nuclear@5: 	ms->res[1] = yres;
nuclear@5: 	ms->res[2] = zres;
nuclear@5: }
nuclear@5: 
nuclear@5: void msurf_threshold(struct metasurface *ms, float thres)
nuclear@5: {
nuclear@5: 	ms->thres = thres;
nuclear@5: }
nuclear@5: 
nuclear@5: 
nuclear@5: int msurf_polygonize(struct metasurface *ms)
nuclear@5: {
nuclear@5: 	int i, j, k;
nuclear@5: 	vec3 pos, delta;
nuclear@5: 
nuclear@5: 	if(!ms->eval || !ms->vertex) {
nuclear@5: 		fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
nuclear@5: 		return -1;
nuclear@5: 	}
nuclear@5: 
nuclear@5: 	for(i=0; i<3; i++) {
nuclear@5: 		delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
nuclear@5: 	}
nuclear@5: 
nuclear@5: 	pos[0] = ms->min[0];
nuclear@5: 	for(i=0; i<ms->res[0] - 1; i++) {
nuclear@5: 
nuclear@5: 		pos[1] = ms->min[1];
nuclear@5: 		for(j=0; j<ms->res[1] - 1; j++) {
nuclear@5: 
nuclear@5: 			pos[2] = ms->min[2];
nuclear@5: 			for(k=0; k<ms->res[2] - 1; k++) {
nuclear@5: 
nuclear@5: 				process_cell(ms, pos, delta);
nuclear@5: 
nuclear@5: 				pos[2] += delta[2];
nuclear@5: 			}
nuclear@5: 			pos[1] += delta[1];
nuclear@5: 		}
nuclear@5: 		pos[0] += delta[0];
nuclear@5: 	}
nuclear@5: 	return 0;
nuclear@5: }
nuclear@5: 
nuclear@5: 
nuclear@5: static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
nuclear@5: {
nuclear@5: 	int i;
nuclear@5: 	vec3 p[8];
nuclear@5: 	float val[8];
nuclear@5: 
nuclear@5: #ifdef USE_MTETRA
nuclear@5: 	static int tetra[][4] = {
nuclear@5: 		{0, 2, 3, 7},
nuclear@5: 		{0, 2, 6, 7},
nuclear@5: 		{0, 4, 6, 7},
nuclear@5: 		{0, 6, 1, 2},
nuclear@5: 		{0, 6, 1, 4},
nuclear@5: 		{5, 6, 1, 4}
nuclear@5: 	};
nuclear@5: #endif
nuclear@5: 
nuclear@5: 	static const float offs[][3] = {
nuclear@5: 		{0.0f, 0.0f, 0.0f},
nuclear@5: 		{1.0f, 0.0f, 0.0f},
nuclear@5: 		{1.0f, 1.0f, 0.0f},
nuclear@5: 		{0.0f, 1.0f, 0.0f},
nuclear@5: 		{0.0f, 0.0f, 1.0f},
nuclear@5: 		{1.0f, 0.0f, 1.0f},
nuclear@5: 		{1.0f, 1.0f, 1.0f},
nuclear@5: 		{0.0f, 1.0f, 1.0f}
nuclear@5: 	};
nuclear@5: 
nuclear@5: 	for(i=0; i<8; i++) {
nuclear@5: 		p[i][0] = pos[0] + sz[0] * offs[i][2];
nuclear@5: 		p[i][1] = pos[1] + sz[1] * offs[i][1];
nuclear@5: 		p[i][2] = pos[2] + sz[2] * offs[i][0];
nuclear@5: 
nuclear@5: 		val[i] = ms->eval(p[i][0], p[i][1], p[i][2]);
nuclear@5: 	}
nuclear@5: 
nuclear@5: #ifdef USE_MTETRA
nuclear@5: 	for(i=0; i<6; i++) {
nuclear@5: 		process_tetra(ms, tetra[i], p, val);
nuclear@5: 	}
nuclear@5: #endif
nuclear@5: #ifdef USE_MCUBES
nuclear@5: 	process_cube(ms, p, val);
nuclear@5: #endif
nuclear@5: }
nuclear@5: 
nuclear@5: 
nuclear@5: /* ---- marching cubes implementation ---- */
nuclear@5: #ifdef USE_MCUBES
nuclear@5: 
nuclear@5: static unsigned int mc_bitcode(float *val, float thres);
nuclear@5: 
nuclear@5: static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
nuclear@5: {
nuclear@5: 	static const int pidx[12][2] = {
nuclear@5: 		{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
nuclear@5: 		{6, 7},	{7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
nuclear@5: 	};
nuclear@5: 	int i, j;
nuclear@5: 	vec3 vert[12];
nuclear@5: 	unsigned int code = mc_bitcode(val, ms->thres);
nuclear@5: 
nuclear@5: 	if(ms->flip) {
nuclear@5: 		code = ~code & 0xff;
nuclear@5: 	}
nuclear@5: 
nuclear@5: 	if(mc_edge_table[code] == 0) {
nuclear@5: 		return;
nuclear@5: 	}
nuclear@5: 
nuclear@5: 	for(i=0; i<12; i++) {
nuclear@5: 		if(mc_edge_table[code] & (1 << i)) {
nuclear@5: 			int p0 = pidx[i][0];
nuclear@5: 			int p1 = pidx[i][1];
nuclear@5: 
nuclear@5: 			float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);
nuclear@5: 			vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
nuclear@5: 			vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
nuclear@5: 			vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
nuclear@5: 		}
nuclear@5: 	}
nuclear@5: 
nuclear@5: 	for(i=0; mc_tri_table[code][i] != -1; i+=3) {
nuclear@5: 		for(j=0; j<3; j++) {
nuclear@5: 			float *v = vert[mc_tri_table[code][i + j]];
nuclear@5: 
nuclear@5: 			if(ms->normal) {
nuclear@5: 				float dfdx, dfdy, dfdz;
nuclear@5: 				dfdx = ms->eval(v[0] - ms->dx, v[1], v[2]) - ms->eval(v[0] + ms->dx, v[1], v[2]);
nuclear@5: 				dfdy = ms->eval(v[0], v[1] - ms->dy, v[2]) - ms->eval(v[0], v[1] + ms->dy, v[2]);
nuclear@5: 				dfdz = ms->eval(v[0], v[1], v[2] - ms->dz) - ms->eval(v[0], v[1], v[2] + ms->dz);
nuclear@5: 
nuclear@5: 				if(ms->flip) {
nuclear@5: 					dfdx = -dfdx;
nuclear@5: 					dfdy = -dfdy;
nuclear@5: 					dfdz = -dfdz;
nuclear@5: 				}
nuclear@5: 				ms->normal(dfdx, dfdy, dfdz);
nuclear@5: 			}
nuclear@5: 
nuclear@5: 			ms->vertex(v[0], v[1], v[2]);
nuclear@5: 		}
nuclear@5: 	}
nuclear@5: }
nuclear@5: 
nuclear@5: static unsigned int mc_bitcode(float *val, float thres)
nuclear@5: {
nuclear@5: 	unsigned int i, res = 0;
nuclear@5: 
nuclear@5: 	for(i=0; i<8; i++) {
nuclear@5: 		if(val[i] > thres) {
nuclear@5: 			res |= 1 << i;
nuclear@5: 		}
nuclear@5: 	}
nuclear@5: 	return res;
nuclear@5: }
nuclear@5: #endif	/* USE_MCUBES */
nuclear@5: 
nuclear@5: 
nuclear@5: /* ---- marching tetrahedra implementation (incomplete) ---- */
nuclear@5: #ifdef USE_MTETRA
nuclear@5: 
nuclear@5: static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
nuclear@5: static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
nuclear@5: 
nuclear@5: 
nuclear@5: #define REVBIT(x)	((x) & 8)
nuclear@5: #define INV(x)		(~(x) & 0xf)
nuclear@5: #define EDGE(a, b)	emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
nuclear@5: static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
nuclear@5: {
nuclear@5: 	unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
nuclear@5: 
nuclear@5: 	switch(code) {
nuclear@5: 	case 1:
nuclear@5: 	case INV(1):
nuclear@5: 		EDGE(0, 1);
nuclear@5: 		EDGE(0, 2);
nuclear@5: 		EDGE(0, 3);
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	case 2:
nuclear@5: 	case INV(2):
nuclear@5: 		EDGE(1, 0);
nuclear@5: 		EDGE(1, 3);
nuclear@5: 		EDGE(1, 2);
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	case 3:
nuclear@5: 	case INV(3):
nuclear@5: 		EDGE(0, 3);
nuclear@5: 		EDGE(0, 2);
nuclear@5: 		EDGE(1, 3);
nuclear@5: 
nuclear@5: 		EDGE(1, 3);
nuclear@5: 		EDGE(1, 2);
nuclear@5: 		EDGE(0, 2);
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	case 4:
nuclear@5: 	case INV(4):
nuclear@5: 		EDGE(2, 0);
nuclear@5: 		EDGE(2, 1);
nuclear@5: 		EDGE(2, 3);
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	case 5:
nuclear@5: 	case INV(5):
nuclear@5: 		EDGE(0, 1);
nuclear@5: 		EDGE(2, 3);
nuclear@5: 		EDGE(0, 3);
nuclear@5: 
nuclear@5: 		EDGE(0, 1);
nuclear@5: 		EDGE(1, 2);
nuclear@5: 		EDGE(2, 3);
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	case 6:
nuclear@5: 	case INV(6):
nuclear@5: 		EDGE(0, 1);
nuclear@5: 		EDGE(1, 3);
nuclear@5: 		EDGE(2, 3);
nuclear@5: 
nuclear@5: 		EDGE(0, 1);
nuclear@5: 		EDGE(0, 2);
nuclear@5: 		EDGE(2, 3);
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	case 7:
nuclear@5: 	case INV(7):
nuclear@5: 		EDGE(3, 0);
nuclear@5: 		EDGE(3, 2);
nuclear@5: 		EDGE(3, 1);
nuclear@5: 		break;
nuclear@5: 
nuclear@5: 	default:
nuclear@5: 		break;	/* cases 0 and 15 */
nuclear@5: 	}
nuclear@5: }
nuclear@5: 
nuclear@5: #define BIT(i)	((v##i > thres) ? (1 << i) : 0)
nuclear@5: static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
nuclear@5: {
nuclear@5: 	return BIT(0) | BIT(1) | BIT(2) | BIT(3);
nuclear@5: }
nuclear@5: 
nuclear@5: static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
nuclear@5: {
nuclear@5: 	int i;
nuclear@5: 	float t = (ms->thres - v0) / (v1 - v0);
nuclear@5: 
nuclear@5: 	vec3 p;
nuclear@5: 	for(i=0; i<3; i++) {
nuclear@5: 		p[i] = p0[i] + (p1[i] - p0[i]) * t;
nuclear@5: 	}
nuclear@5: 	ms->vertex(p[0], p[1], p[2]);
nuclear@5: 
nuclear@5: 	/*for(i=0; i<3; i++) {
nuclear@5: 		ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
nuclear@5: 	}
nuclear@5: 
nuclear@5: 	if(++ms->nverts >= 3) {
nuclear@5: 		ms->nverts = 0;
nuclear@5: 
nuclear@5: 		for(i=0; i<3; i++) {
nuclear@5: 			int idx = rev ? (2 - i) : i;
nuclear@5: 			ms->vertex(ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
nuclear@5: 		}
nuclear@5: 	}*/
nuclear@5: }
nuclear@5: 
nuclear@5: #endif	/* USE_MTETRA */