X-Git-Url: http://git.mutantstargoat.com/user/nuclear/?p=laserbrain_demo;a=blobdiff_plain;f=src%2Fblobs%2Fmetasurf.c;h=0298e0561fc7f9b20015e9f5631b460e19b25f57;hp=2207b3df50f6d5253100d28c29166a61b553eab6;hb=b75b6703809abf0dc02f2e557ba73e9efbefa4d1;hpb=82a1d96d9a18e94ccdc13b4bda0c470e81e70768

diff --git a/src/blobs/metasurf.c b/src/blobs/metasurf.c
index 2207b3d..0298e05 100644
--- a/src/blobs/metasurf.c
+++ b/src/blobs/metasurf.c
@@ -1,60 +1,28 @@
-/*
-metasurf - a library for implicit surface polygonization
-Copyright (C) 2011-2016  John Tsiombikas <nuclear@member.fsf.org>
-
-This program is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public License
-along with this program.  If not, see <http://www.gnu.org/licenses/>.
-*/
-/* this is pulled from: https://github.com/jtsiomb/metasurf */
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
 #include "metasurf.h"
 #include "mcubes.h"
-
-#undef USE_MTETRA
-#define USE_MCUBES
-
-#if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))
-#error "pick either USE_MTETRA or USE_MCUBES, not both..."
-#endif
+#include "logger.h"
 
 typedef float vec3[3];
 
 struct metasurface {
 	vec3 min, max;
-	int res[3];
+	int res[3], newres[3];
 	float thres;
 
-	msurf_eval_func_t eval;
-	msurf_vertex_func_t vertex;
-	msurf_normal_func_t normal;
-	void *udata;
-
 	float dx, dy, dz;
-	int flip;
+	unsigned int flags;
+
+	float *voxels;
 
-	vec3 vbuf[3];
-	int nverts;
+	int varr_size, varr_alloc_size;
+	float *varr, *narr;
 };
 
 static int msurf_init(struct metasurface *ms);
-static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);
-#ifdef USE_MTETRA
-static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);
-#endif
-#ifdef USE_MCUBES
-static void process_cube(struct metasurface *ms, vec3 *pos, float *val);
-#endif
+static void process_cell(struct metasurface *ms, int xcell, int ycell, int zcell, vec3 pos, vec3 sz);
 
 
 struct metasurface *msurf_create(void)
@@ -72,46 +40,56 @@ struct metasurface *msurf_create(void)
 
 void msurf_free(struct metasurface *ms)
 {
-	free(ms);
+	if(ms) {
+		free(ms->voxels);
+		free(ms->varr);
+		free(ms->narr);
+		free(ms);
+	}
 }
 
 static int msurf_init(struct metasurface *ms)
 {
+	ms->voxels = 0;
 	ms->thres = 0.0;
-	ms->eval = 0;
-	ms->vertex = 0;
-	ms->normal = 0;
-	ms->udata = 0;
 	ms->min[0] = ms->min[1] = ms->min[2] = -1.0;
 	ms->max[0] = ms->max[1] = ms->max[2] = 1.0;
-	ms->res[0] = ms->res[1] = ms->res[2] = 40;
-	ms->nverts = 0;
+	ms->res[0] = ms->res[1] = ms->res[2] = 0;
+	ms->newres[0] = ms->newres[1] = ms->newres[2] = 40;
+
+	ms->varr_alloc_size = ms->varr_size = 0;
+	ms->varr = ms->narr = 0;
 
 	ms->dx = ms->dy = ms->dz = 0.001;
-	ms->flip = 0;
+	ms->flags = 0;
 
 	return 0;
 }
 
-void msurf_set_user_data(struct metasurface *ms, void *udata)
+void msurf_enable(struct metasurface *ms, unsigned int opt)
 {
-	ms->udata = udata;
+	ms->flags |= opt;
 }
 
-void *msurf_get_user_data(struct metasurface *ms)
+void msurf_disable(struct metasurface *ms, unsigned int opt)
 {
-	return ms->udata;
+	ms->flags &= ~opt;
+}
+
+int msurf_is_enabled(struct metasurface *ms, unsigned int opt)
+{
+	return ms->flags & opt;
 }
 
 void msurf_set_inside(struct metasurface *ms, int inside)
 {
 	switch(inside) {
 	case MSURF_GREATER:
-		ms->flip = 0;
+		msurf_disable(ms, MSURF_FLIP);
 		break;
 
 	case MSURF_LESS:
-		ms->flip = 1;
+		msurf_enable(ms, MSURF_FLIP);
 		break;
 
 	default:
@@ -121,22 +99,7 @@ void msurf_set_inside(struct metasurface *ms, int inside)
 
 int msurf_get_inside(struct metasurface *ms)
 {
-	return ms->flip ? MSURF_LESS : MSURF_GREATER;
-}
-
-void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)
-{
-	ms->eval = func;
-}
-
-void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)
-{
-	ms->vertex = func;
-}
-
-void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)
-{
-	ms->normal = func;
+	return msurf_is_enabled(ms, MSURF_FLIP) ? MSURF_LESS : MSURF_GREATER;
 }
 
 void msurf_set_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)
@@ -161,9 +124,9 @@ void msurf_get_bounds(struct metasurface *ms, float *xmin, float *ymin, float *z
 
 void msurf_set_resolution(struct metasurface *ms, int xres, int yres, int zres)
 {
-	ms->res[0] = xres;
-	ms->res[1] = yres;
-	ms->res[2] = zres;
+	ms->newres[0] = xres;
+	ms->newres[1] = yres;
+	ms->newres[2] = zres;
 }
 
 void msurf_get_resolution(struct metasurface *ms, int *xres, int *yres, int *zres)
@@ -184,112 +147,155 @@ float msurf_get_threshold(struct metasurface *ms)
 }
 
 
+float *msurf_voxels(struct metasurface *ms)
+{
+	if(ms->res[0] != ms->newres[0] || ms->res[1] != ms->newres[1] || ms->res[2] != ms->newres[2]) {
+		int count;
+		ms->res[0] = ms->newres[0];
+		ms->res[1] = ms->newres[1];
+		ms->res[2] = ms->newres[2];
+		count = ms->res[0] * ms->res[1] * ms->res[2];
+		free(ms->voxels);
+		if(!(ms->voxels = malloc(count * sizeof *ms->voxels))) {
+			return 0;
+		}
+	}
+	return ms->voxels;
+}
+
+float *msurf_slice(struct metasurface *ms, int idx)
+{
+	float *vox = msurf_voxels(ms);
+	if(!vox) return 0;
+
+	return vox + ms->res[0] * ms->res[1] * idx;
+}
+
 int msurf_polygonize(struct metasurface *ms)
 {
 	int i, j, k;
 	vec3 pos, delta;
 
-	if(!ms->eval || !ms->vertex) {
-		fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");
-		return -1;
-	}
+	if(!ms->voxels) return -1;
+
+	ms->varr_size = 0;
 
 	for(i=0; i<3; i++) {
 		delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];
 	}
 
-	pos[0] = ms->min[0];
 	for(i=0; i<ms->res[0] - 1; i++) {
-		pos[1] = ms->min[1];
 		for(j=0; j<ms->res[1] - 1; j++) {
-
-			pos[2] = ms->min[2];
 			for(k=0; k<ms->res[2] - 1; k++) {
 
-				process_cell(ms, pos, delta);
+				pos[0] = ms->min[0] + i * delta[0];
+				pos[1] = ms->min[1] + j * delta[1];
+				pos[2] = ms->min[2] + k * delta[2];
 
-				pos[2] += delta[2];
+				process_cell(ms, i, j, k, pos, delta);
 			}
-			pos[1] += delta[1];
 		}
-		pos[0] += delta[0];
 	}
 	return 0;
 }
 
-
-static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)
+int msurf_vertex_count(struct metasurface *ms)
 {
-	int i;
-	vec3 p[8];
-	float val[8];
-
-#ifdef USE_MTETRA
-	static int tetra[][4] = {
-		{0, 2, 3, 7},
-		{0, 2, 6, 7},
-		{0, 4, 6, 7},
-		{0, 6, 1, 2},
-		{0, 6, 1, 4},
-		{5, 6, 1, 4}
-	};
-#endif
-
-	static const float offs[][3] = {
-		{0.0f, 0.0f, 0.0f},
-		{1.0f, 0.0f, 0.0f},
-		{1.0f, 1.0f, 0.0f},
-		{0.0f, 1.0f, 0.0f},
-		{0.0f, 0.0f, 1.0f},
-		{1.0f, 0.0f, 1.0f},
-		{1.0f, 1.0f, 1.0f},
-		{0.0f, 1.0f, 1.0f}
-	};
-
-	for(i=0; i<8; i++) {
-		p[i][0] = pos[0] + sz[0] * offs[i][2];
-		p[i][1] = pos[1] + sz[1] * offs[i][1];
-		p[i][2] = pos[2] + sz[2] * offs[i][0];
-
-		val[i] = ms->eval(ms, p[i][0], p[i][1], p[i][2]);
-	}
-
-#ifdef USE_MTETRA
-	for(i=0; i<6; i++) {
-		process_tetra(ms, tetra[i], p, val);
-	}
-#endif
-#ifdef USE_MCUBES
-	process_cube(ms, p, val);
-#endif
+	return ms->varr_size / 3;
 }
 
+float *msurf_vertices(struct metasurface *ms)
+{
+	return ms->varr;
+}
 
-/* ---- marching cubes implementation ---- */
-#ifdef USE_MCUBES
+float *msurf_normals(struct metasurface *ms)
+{
+	return ms->narr;
+}
 
 static unsigned int mc_bitcode(float *val, float thres);
 
-static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
+static void process_cell(struct metasurface *ms, int xcell, int ycell, int zcell, vec3 cellpos, vec3 cellsz)
 {
+	int i, j, k, slice_size;
+	vec3 pos[8];
+	float dfdx[8], dfdy[8], dfdz[8];
+	vec3 vert[12], norm[12];
+	float val[8];
+	float *cellptr;
+	unsigned int code;
+
+	static const int offs[][3] = {
+		{0, 0, 0},
+		{1, 0, 0},
+		{1, 1, 0},
+		{0, 1, 0},
+		{0, 0, 1},
+		{1, 0, 1},
+		{1, 1, 1},
+		{0, 1, 1}
+	};
+
 	static const int pidx[12][2] = {
 		{0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},
 		{6, 7},	{7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
 	};
-	int i, j;
-	vec3 vert[12];
-	unsigned int code = mc_bitcode(val, ms->thres);
 
-	if(ms->flip) {
+	slice_size = ms->res[0] * ms->res[1];
+	cellptr = ms->voxels + slice_size * zcell + ms->res[0] * ycell + xcell;
+
+	val[0] = *cellptr;
+	val[1] = *(cellptr + 1);
+	val[2] = *(cellptr + 1 + ms->res[0]);
+	val[3] = *(cellptr + ms->res[0]);
+	cellptr += slice_size;
+	val[4] = *cellptr;
+	val[5] = *(cellptr + 1);
+	val[6] = *(cellptr + 1 + ms->res[0]);
+	val[7] = *(cellptr + ms->res[0]);
+	cellptr -= slice_size;
+
+	code = mc_bitcode(val, ms->thres);
+	if(ms->flags & MSURF_FLIP) {
 		code = ~code & 0xff;
 	}
-
 	if(mc_edge_table[code] == 0) {
 		return;
 	}
 
+	/* calculate normals at the 8 corners */
+	for(i=0; i<8; i++) {
+		float *ptr = cellptr + offs[i][2] * slice_size + offs[i][1] * ms->res[0] + offs[i][0];
+
+		if(xcell < ms->res[0] - 1) {
+			dfdx[i] = *(ptr + 1) - *ptr;
+		} else {
+			dfdx[i] = *ptr - *(ptr - 1);
+		}
+		if(ycell < ms->res[1] - 1) {
+			dfdy[i] = *(ptr + ms->res[0]) - *ptr;
+		} else {
+			dfdy[i] = *ptr - *(ptr - ms->res[0]);
+		}
+		if(zcell < ms->res[2] - 1) {
+			dfdz[i] = *(ptr + slice_size) - *ptr;
+		} else {
+			dfdz[i] = *ptr - *(ptr - slice_size);
+		}
+	}
+
+	/* calculate the world-space position of each corner */
+	for(i=0; i<8; i++) {
+		pos[i][0] = cellpos[0] + cellsz[0] * offs[i][0];
+		pos[i][1] = cellpos[1] + cellsz[1] * offs[i][1];
+		pos[i][2] = cellpos[2] + cellsz[2] * offs[i][2];
+	}
+
+	/* generate up to a max of 12 vertices per cube. interpolate positions and normals for each one */
 	for(i=0; i<12; i++) {
 		if(mc_edge_table[code] & (1 << i)) {
+			float nx, ny, nz;
 			int p0 = pidx[i][0];
 			int p1 = pidx[i][1];
 
@@ -297,29 +303,61 @@ static void process_cube(struct metasurface *ms, vec3 *pos, float *val)
 			vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;
 			vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;
 			vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;
+
+			nx = dfdx[p0] + (dfdx[p1] - dfdx[p0]) * t;
+			ny = dfdy[p0] + (dfdy[p1] - dfdy[p0]) * t;
+			nz = dfdz[p0] + (dfdz[p1] - dfdz[p0]) * t;
+
+			if(ms->flags & MSURF_FLIP) {
+				nx = -nx;
+				ny = -ny;
+				nz = -nz;
+			}
+
+			if(ms->flags & MSURF_NORMALIZE) {
+				float len = sqrt(nx * nx + ny * ny + nz * nz);
+				if(len != 0.0) {
+					float s = 1.0f / len;
+					nx *= s;
+					ny *= s;
+					nz *= s;
+				}
+			}
+
+			norm[i][0] = nx;
+			norm[i][1] = ny;
+			norm[i][2] = nz;
 		}
 	}
 
+	/* for each triangle of the cube add the appropriate vertices to the vertex buffer */
 	for(i=0; mc_tri_table[code][i] != -1; i+=3) {
 		for(j=0; j<3; j++) {
-			float *v = vert[mc_tri_table[code][i + j]];
-
-			if(ms->normal) {
-				float dfdx, dfdy, dfdz;
-				dfdx = ms->eval(ms, v[0] - ms->dx, v[1], v[2]) - ms->eval(ms, v[0] + ms->dx, v[1], v[2]);
-				dfdy = ms->eval(ms, v[0], v[1] - ms->dy, v[2]) - ms->eval(ms, v[0], v[1] + ms->dy, v[2]);
-				dfdz = ms->eval(ms, v[0], v[1], v[2] - ms->dz) - ms->eval(ms, v[0], v[1], v[2] + ms->dz);
-
-				if(ms->flip) {
-					dfdx = -dfdx;
-					dfdy = -dfdy;
-					dfdz = -dfdz;
+			int idx = mc_tri_table[code][i + j];
+			float *v = vert[idx];
+			float *n = norm[idx];
+
+			/* TODO multithreadied polygon emit */
+			if(ms->varr_size + 3 > ms->varr_alloc_size) {
+				int newsz = ms->varr_alloc_size ? ms->varr_alloc_size * 2 : 1024;
+				float *new_varr, *new_narr;
+
+				if(!(new_varr = realloc(ms->varr, newsz * sizeof *new_varr)) ||
+						!(new_narr = realloc(ms->narr, newsz * sizeof *new_narr))) {
+					free(new_varr);
+					error_log("msurf_polygonize: failed to grow vertex buffers to %d elements\n", newsz);
+					return;
 				}
-				ms->normal(ms, dfdx, dfdy, dfdz);
+				ms->varr = new_varr;
+				ms->narr = new_narr;
+				ms->varr_alloc_size = newsz;
 			}
 
-			/* TODO multithreadied polygon emmit */
-			ms->vertex(ms, v[0], v[1], v[2]);
+			for(k=0; k<3; k++) {
+				ms->varr[ms->varr_size] = v[k];
+				ms->narr[ms->varr_size] = n[k];
+				++ms->varr_size;
+			}
 		}
 	}
 }
@@ -335,118 +373,3 @@ static unsigned int mc_bitcode(float *val, float thres)
 	}
 	return res;
 }
-#endif	/* USE_MCUBES */
-
-
-/* ---- marching tetrahedra implementation (incomplete) ---- */
-#ifdef USE_MTETRA
-
-static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);
-static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
-
-
-#define REVBIT(x)	((x) & 8)
-#define INV(x)		(~(x) & 0xf)
-#define EDGE(a, b)	emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))
-static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)
-{
-	unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);
-
-	switch(code) {
-	case 1:
-	case INV(1):
-		EDGE(0, 1);
-		EDGE(0, 2);
-		EDGE(0, 3);
-		break;
-
-	case 2:
-	case INV(2):
-		EDGE(1, 0);
-		EDGE(1, 3);
-		EDGE(1, 2);
-		break;
-
-	case 3:
-	case INV(3):
-		EDGE(0, 3);
-		EDGE(0, 2);
-		EDGE(1, 3);
-
-		EDGE(1, 3);
-		EDGE(1, 2);
-		EDGE(0, 2);
-		break;
-
-	case 4:
-	case INV(4):
-		EDGE(2, 0);
-		EDGE(2, 1);
-		EDGE(2, 3);
-		break;
-
-	case 5:
-	case INV(5):
-		EDGE(0, 1);
-		EDGE(2, 3);
-		EDGE(0, 3);
-
-		EDGE(0, 1);
-		EDGE(1, 2);
-		EDGE(2, 3);
-		break;
-
-	case 6:
-	case INV(6):
-		EDGE(0, 1);
-		EDGE(1, 3);
-		EDGE(2, 3);
-
-		EDGE(0, 1);
-		EDGE(0, 2);
-		EDGE(2, 3);
-		break;
-
-	case 7:
-	case INV(7):
-		EDGE(3, 0);
-		EDGE(3, 2);
-		EDGE(3, 1);
-		break;
-
-	default:
-		break;	/* cases 0 and 15 */
-	}
-}
-
-#define BIT(i)	((v##i > thres) ? (1 << i) : 0)
-static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)
-{
-	return BIT(0) | BIT(1) | BIT(2) | BIT(3);
-}
-
-static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)
-{
-	int i;
-	float t = (ms->thres - v0) / (v1 - v0);
-
-	vec3 p;
-	for(i=0; i<3; i++) {
-		p[i] = p0[i] + (p1[i] - p0[i]) * t;
-	}
-	ms->vertex(ms, p[0], p[1], p[2]);
-
-	/*for(i=0; i<3; i++) {
-		ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;
-	}
-
-	if(++ms->nverts >= 3) {
-		ms->nverts = 0;
-
-		for(i=0; i<3; i++) {
-			int idx = rev ? (2 - i) : i;
-			ms->vertex(ms, ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);
-		}
-	}*/
-}
-#endif	/* USE_MTETRA */