# HG changeset patch
# User John Tsiombikas <nuclear@member.fsf.org>
# Date 1319538603 -10800
# Node ID 52664d3451ad8dc382bdff0241948f5c594ba429
# Parent  9ab057fba0c51d50180471fb7a821e75c247302d
added volume rendering example

diff -r 9ab057fba0c5 -r 52664d3451ad examples/metaballs/src/metaballs.c
--- a/examples/metaballs/src/metaballs.c	Tue Oct 25 08:59:07 2011 +0300
+++ b/examples/metaballs/src/metaballs.c	Tue Oct 25 13:30:03 2011 +0300
@@ -31,7 +31,6 @@
 int num_mballs = sizeof mball / sizeof *mball;
 
 float eval(float x, float y, float z);
-float eval_cached(float x, float y, float z);
 void vertex(float x, float y, float z);
 void render(void);
 void disp(void);
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/Makefile
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/Makefile	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,34 @@
+src = $(wildcard src/*.c)
+obj = $(src:.c=.o)
+dep = $(obj:.o=.d)
+bin = volume
+
+# uncomment the following line to disable shaders (for old systems like SGI or whatever)
+#nosdr = -DNO_SHADERS
+
+CC = gcc
+CFLAGS = -pedantic -Wall -g -I../../src $(nosdr)
+LDFLAGS = -L../.. -lmetasurf $(libgl) -limago
+
+ifeq ($(shell uname -s), Darwin)
+	libgl = -framework OpenGL -framework GLUT -lGLEW
+else
+	libgl = -lGL -lGLU -lglut -lGLEW
+endif
+
+$(bin): $(obj)
+	$(CC) -o $@ $(obj) $(LDFLAGS)
+
+-include $(dep)
+
+%.d: %.c
+	@$(CPP) $(CFLAGS) $< -MM -MT $(@:.d=.o) >$@
+
+.PHONY: clean
+clean:
+	rm -f $(obj) $(bin)
+
+
+.PHONY: cleandep
+cleandep:
+	rm -f $(dep)
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/sdr/frag.glsl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/sdr/frag.glsl	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,19 @@
+varying vec3 vpos, ldir, norm;
+
+void main()
+{
+	vec3 n = normalize(norm);
+	vec3 l = normalize(ldir);
+	vec3 v = -normalize(vpos);
+	vec3 h = normalize(v + l);
+
+	const vec3 kd = vec3(0.87, 0.82, 0.74);
+
+	float diff = abs(dot(n, l));
+	float spec = pow(abs(dot(n, h)), 60.0);
+
+	vec3 dcol = kd * diff;
+	vec3 scol = vec3(0.8, 0.8, 0.8) * spec;
+
+	gl_FragColor = vec4(dcol + scol, 1.0);
+}
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/sdr/vert.glsl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/sdr/vert.glsl	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,10 @@
+varying vec3 vpos, ldir, norm;
+
+void main()
+{
+	gl_Position = ftransform();
+
+	vpos = (gl_ModelViewMatrix * gl_Vertex).xyz;
+	norm = gl_NormalMatrix * gl_Normal;
+	ldir = gl_LightSource[0].position.xyz;
+}
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/src/cam.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/src/cam.c	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,216 @@
+#include <math.h>
+#include "cam.h"
+
+#if defined(__APPLE__) && defined(__MACH__)
+#include <OpenGL/gl.h>
+#else
+#include <GL/gl.h>
+#endif
+
+#define DR				(M_PI / 180.0)
+#define DEG_TO_RAD(x)	((x) * DR)
+
+typedef struct vec3 {
+	float x, y, z;
+} vec3_t;
+
+/* viewing parameters */
+#define DEF_THETA	0
+#define DEF_PHI		0
+#define DEF_DIST	0
+#define DEF_X		0
+#define DEF_Y		0
+#define DEF_Z		0
+
+static float cam_theta = DEF_THETA, cam_phi = DEF_PHI;
+static float cam_dist = DEF_DIST;
+static vec3_t cam_pos = {DEF_X, DEF_Y, DEF_Z};
+
+/* projection parameters */
+#define DEF_VFOV	45.0
+#define DEF_ASPECT	1.3333333
+#define DEF_NEAR	1.0
+#define DEF_FAR		1000.0
+
+static float vfov = DEF_VFOV;
+static float aspect = DEF_ASPECT;
+static float nearclip = DEF_NEAR, farclip = DEF_FAR;
+
+/* stereo parameters */
+#define DEF_EYE_SEP		0.1
+#define DEF_FOCUS_DIST	1.0
+
+static float eye_sep = DEF_EYE_SEP;
+static float focus_dist = DEF_FOCUS_DIST;
+
+static float pan_speed = 0.001;
+static float rot_speed = 0.5;
+static float zoom_speed = 0.1;
+
+static float vmin_deg = -90.0, vmax_deg = 90.0;
+
+void cam_reset(void)
+{
+	cam_reset_view();
+	cam_reset_proj();
+	cam_reset_stereo();
+}
+
+void cam_reset_view(void)
+{
+	cam_theta = DEF_THETA;
+	cam_phi = DEF_PHI;
+	cam_dist = DEF_DIST;
+	cam_pos.x = DEF_X;
+	cam_pos.y = DEF_Y;
+	cam_pos.z = DEF_Z;
+}
+
+void cam_reset_proj(void)
+{
+	vfov = DEF_VFOV;
+	aspect = DEF_ASPECT;
+	nearclip = DEF_NEAR;
+	farclip = DEF_FAR;
+}
+
+void cam_reset_stereo(void)
+{
+	eye_sep = DEF_EYE_SEP;
+	focus_dist = DEF_FOCUS_DIST;
+}
+
+void cam_set_vrange(float min_deg, float max_deg)
+{
+	vmin_deg = min_deg;
+	vmax_deg = max_deg;
+}
+
+void cam_move(float x, float y, float z)
+{
+	cam_pos.x += x;
+	cam_pos.y += y;
+	cam_pos.z += z;
+}
+
+void cam_rotate(float theta, float phi)
+{
+	cam_phi += phi;
+	cam_theta += theta;
+}
+
+void cam_dolly(float dist)
+{
+	cam_dist += dist;
+}
+
+void cam_inp_pan_speed(float speed)
+{
+	pan_speed = speed;
+}
+
+void cam_inp_rotate_speed(float speed)
+{
+	rot_speed = speed;
+}
+
+void cam_inp_zoom_speed(float speed)
+{
+	zoom_speed = speed;
+}
+
+
+void cam_inp_pan(int dx, int dy)
+{
+	float dxf = dx * pan_speed;
+	float dyf = dy * pan_speed;
+	float angle = -DEG_TO_RAD(cam_theta);
+
+	cam_pos.x += cos(angle) * dxf + sin(angle) * dyf;
+	cam_pos.z += -sin(angle) * dxf + cos(angle) * dyf;
+}
+
+void cam_inp_height(int dh)
+{
+	cam_pos.y += dh * pan_speed;
+}
+
+void cam_inp_rotate(int dx, int dy)
+{
+	cam_theta += dx * rot_speed;
+	cam_phi += dy * rot_speed;
+
+	if(cam_phi < vmin_deg) cam_phi = vmin_deg;
+	if(cam_phi > vmax_deg) cam_phi = vmax_deg;
+}
+
+void cam_inp_zoom(int dz)
+{
+	cam_dist += dz * zoom_speed;
+	if(cam_dist < 0.001) {
+		cam_dist = 0.001;
+	}
+}
+
+void cam_clip(float n, float f)
+{
+	nearclip = n;
+	farclip = f;
+}
+
+void cam_fov(float f)
+{
+	vfov = f;
+}
+
+void cam_aspect(float a)
+{
+	aspect = a;
+}
+
+void cam_separation(float s)
+{
+	eye_sep = s;
+}
+
+void cam_focus_dist(float d)
+{
+	focus_dist = d;
+
+	cam_separation(d / 30.0);
+}
+
+void cam_view_matrix(void)
+{
+	cam_stereo_view_matrix(CAM_CENTER);
+}
+
+void cam_stereo_view_matrix(int eye)
+{
+	static const float offs_sign[] = {0.0f, 0.5f, -0.5f};	/* center, left, right */
+	float offs = eye_sep * offs_sign[eye];
+
+	glTranslatef(offs, 0, 0);
+
+	glTranslatef(0, 0, -cam_dist);
+	glRotatef(cam_phi, 1, 0, 0);
+	glRotatef(cam_theta, 0, 1, 0);
+	glTranslatef(-cam_pos.x, -cam_pos.y, -cam_pos.z);
+}
+
+void cam_proj_matrix(void)
+{
+	cam_stereo_proj_matrix(CAM_CENTER);
+}
+
+void cam_stereo_proj_matrix(int eye)
+{
+	float vfov_rad = M_PI * vfov / 180.0;
+	float top = nearclip * tan(vfov_rad * 0.5);
+	float right = top * aspect;
+
+	static const float offs_sign[] = {0.0f, 1.0, -1.0};	/* center, left, right */
+	float frust_shift = offs_sign[eye] * (eye_sep * 0.5 * nearclip / focus_dist);
+
+	glFrustum(-right + frust_shift, right + frust_shift, -top, top, nearclip, farclip);
+}
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/src/cam.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/src/cam.h	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,54 @@
+#ifndef CAM_H_
+#define CAM_H_
+
+enum {
+	CAM_CENTER,
+	CAM_LEFT,
+	CAM_RIGHT
+};
+
+/* reset to the initial state */
+void cam_reset(void);			/* all */
+void cam_reset_view(void);		/* view parameters */
+void cam_reset_proj(void);		/* projection parameters */
+void cam_reset_stereo(void);	/* stereo parameters */
+
+void cam_set_vrange(float min_deg, float max_deg);
+
+void cam_move(float x, float y, float z);
+void cam_rotate(float theta, float phi);
+void cam_dolly(float dist);
+
+/* camera input handling */
+void cam_inp_pan_speed(float speed);
+void cam_inp_rotate_speed(float speed);
+void cam_inp_zoom_speed(float speed);
+
+void cam_inp_pan(int dx, int dy);		/* pan across X/Z plane */
+void cam_inp_height(int dh);			/* move verticaly */
+void cam_inp_rotate(int dx, int dy);	/* rotate around local Y and X axis */
+void cam_inp_zoom(int dz);				/* dolly the camera fwd/back */
+
+/* camera projection parameters */
+void cam_clip(float n, float f);	/* set clipping planes */
+void cam_fov(float f);				/* vertical field of view in degrees */
+void cam_aspect(float a);			/* aspect ratio (width / height) */
+
+/* stereo parameters */
+void cam_separation(float s);
+void cam_focus_dist(float d);
+
+
+/* multiply the camera view matrix on top of the current matrix stack
+ * (which should be GL_MODELVIEW)
+ */
+void cam_view_matrix(void);
+void cam_stereo_view_matrix(int eye);
+
+/* multiply the camera projection matrix on top of the current matrix stack
+ * (which should be GL_PROJECTION)
+ */
+void cam_proj_matrix(void);
+void cam_stereo_proj_matrix(int eye);
+
+#endif	/* CAM_H_ */
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/src/sdr.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/src/sdr.c	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,333 @@
+#ifndef NO_SHADERS
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <errno.h>
+#include <assert.h>
+#include <GL/glew.h>
+
+#if defined(unix) || defined(__unix__)
+#include <unistd.h>
+#include <sys/stat.h>
+#endif	/* unix */
+
+#include "sdr.h"
+
+unsigned int create_vertex_shader(const char *src)
+{
+	return create_shader(src, GL_VERTEX_SHADER);
+}
+
+unsigned int create_pixel_shader(const char *src)
+{
+	return create_shader(src, GL_FRAGMENT_SHADER);
+}
+
+unsigned int create_shader(const char *src, unsigned int sdr_type)
+{
+	unsigned int sdr;
+	int success, info_len;
+	char *info_str = 0;
+	GLenum err;
+
+	sdr = glCreateShader(sdr_type);
+	assert(glGetError() == GL_NO_ERROR);
+	glShaderSource(sdr, 1, &src, 0);
+	err = glGetError();
+	assert(err == GL_NO_ERROR);
+	glCompileShader(sdr);
+	assert(glGetError() == GL_NO_ERROR);
+
+	glGetShaderiv(sdr, GL_COMPILE_STATUS, &success);
+	assert(glGetError() == GL_NO_ERROR);
+	glGetShaderiv(sdr, GL_INFO_LOG_LENGTH, &info_len);
+	assert(glGetError() == GL_NO_ERROR);
+
+	if(info_len) {
+		if((info_str = malloc(info_len + 1))) {
+			glGetShaderInfoLog(sdr, info_len, 0, info_str);
+			assert(glGetError() == GL_NO_ERROR);
+		}
+	}
+
+	if(success) {
+		fprintf(stderr, info_str ? "done: %s\n" : "done\n", info_str);
+	} else {
+		fprintf(stderr, info_str ? "failed: %s\n" : "failed\n", info_str);
+		glDeleteShader(sdr);
+		sdr = 0;
+	}
+
+	free(info_str);
+	return sdr;
+}
+
+void free_shader(unsigned int sdr)
+{
+	glDeleteShader(sdr);
+}
+
+unsigned int load_vertex_shader(const char *fname)
+{
+	return load_shader(fname, GL_VERTEX_SHADER);
+}
+
+unsigned int load_pixel_shader(const char *fname)
+{
+	return load_shader(fname, GL_FRAGMENT_SHADER);
+}
+
+unsigned int load_shader(const char *fname, unsigned int sdr_type)
+{
+#if defined(unix) || defined(__unix__)
+	struct stat st;
+#endif
+	unsigned int sdr;
+	size_t filesize;
+	FILE *fp;
+	char *src;
+
+	if(!(fp = fopen(fname, "r"))) {
+		fprintf(stderr, "failed to open shader %s: %s\n", fname, strerror(errno));
+		return 0;
+	}
+
+#if defined(unix) || defined(__unix__)
+	fstat(fileno(fp), &st);
+	filesize = st.st_size;
+#else
+	fseek(fp, 0, SEEK_END);
+	filesize = ftell(fp);
+	fseek(fp, 0, SEEK_SET);
+#endif	/* unix */
+
+	if(!(src = malloc(filesize + 1))) {
+		fclose(fp);
+		return 0;
+	}
+	fread(src, 1, filesize, fp);
+	src[filesize] = 0;
+	fclose(fp);
+
+	fprintf(stderr, "compiling %s shader: %s... ", (sdr_type == GL_VERTEX_SHADER ? "vertex" : "pixel"), fname);
+	sdr = create_shader(src, sdr_type);
+
+	free(src);
+	return sdr;
+}
+
+
+unsigned int get_vertex_shader(const char *fname)
+{
+	return get_shader(fname, GL_VERTEX_SHADER);
+}
+
+unsigned int get_pixel_shader(const char *fname)
+{
+	return get_shader(fname, GL_FRAGMENT_SHADER);
+}
+
+unsigned int get_shader(const char *fname, unsigned int sdr_type)
+{
+	unsigned int sdr;
+
+	if(!(sdr = load_shader(fname, sdr_type))) {
+		return 0;
+	}
+	return sdr;
+}
+
+
+/* ---- gpu programs ---- */
+
+unsigned int create_program(void)
+{
+	unsigned int prog = glCreateProgram();
+	assert(glGetError() == GL_NO_ERROR);
+	return prog;
+}
+
+unsigned int create_program_link(unsigned int vs, unsigned int ps)
+{
+	unsigned int prog;
+
+	if(!(prog = create_program())) {
+		return 0;
+	}
+
+	if(vs) {
+		attach_shader(prog, vs);
+		assert(glGetError() == GL_NO_ERROR);
+	}
+	if(ps) {
+		attach_shader(prog, ps);
+		assert(glGetError() == GL_NO_ERROR);
+	}
+
+	if(link_program(prog) == -1) {
+		free_program(prog);
+		return 0;
+	}
+	return prog;
+}
+
+unsigned int create_program_load(const char *vfile, const char *pfile)
+{
+	unsigned int vs = 0, ps = 0;
+
+	if(vfile && !(vs = get_vertex_shader(vfile))) {
+		return 0;
+	}
+	if(pfile && !(ps = get_pixel_shader(pfile))) {
+		return 0;
+	}
+	return create_program_link(vs, ps);
+}
+
+void free_program(unsigned int sdr)
+{
+	glDeleteProgram(sdr);
+}
+
+void attach_shader(unsigned int prog, unsigned int sdr)
+{
+	glAttachShader(prog, sdr);
+	assert(glGetError() == GL_NO_ERROR);
+}
+
+int link_program(unsigned int prog)
+{
+	int linked, info_len, retval = 0;
+	char *info_str = 0;
+
+	glLinkProgram(prog);
+	assert(glGetError() == GL_NO_ERROR);
+	glGetProgramiv(prog, GL_LINK_STATUS, &linked);
+	assert(glGetError() == GL_NO_ERROR);
+	glGetProgramiv(prog, GL_INFO_LOG_LENGTH, &info_len);
+	assert(glGetError() == GL_NO_ERROR);
+
+	if(info_len) {
+		if((info_str = malloc(info_len + 1))) {
+			glGetProgramInfoLog(prog, info_len, 0, info_str);
+			assert(glGetError() == GL_NO_ERROR);
+		}
+	}
+
+	if(linked) {
+		fprintf(stderr, info_str ? "linking done: %s\n" : "linking done\n", info_str);
+	} else {
+		fprintf(stderr, info_str ? "linking failed: %s\n" : "linking failed\n", info_str);
+		retval = -1;
+	}
+
+	free(info_str);
+	return retval;
+}
+
+int bind_program(unsigned int prog)
+{
+	GLenum err;
+
+	glUseProgram(prog);
+	if(prog && (err = glGetError()) != GL_NO_ERROR) {
+		/* maybe the program is not linked, try linking first */
+		if(err == GL_INVALID_OPERATION) {
+			if(link_program(prog) == -1) {
+				return -1;
+			}
+			glUseProgram(prog);
+			return glGetError() == GL_NO_ERROR ? 0 : -1;
+		}
+		return -1;
+	}
+	return 0;
+}
+
+/* ugly but I'm not going to write the same bloody code over and over */
+#define BEGIN_UNIFORM_CODE \
+	int loc, curr_prog; \
+	glGetIntegerv(GL_CURRENT_PROGRAM, &curr_prog); \
+	if(curr_prog != prog && bind_program(prog) == -1) { \
+		return -1; \
+	} \
+	if((loc = glGetUniformLocation(prog, name)) != -1)
+
+#define END_UNIFORM_CODE \
+	if(curr_prog != prog) { \
+		bind_program(curr_prog); \
+	} \
+	return loc == -1 ? -1 : 0
+
+int set_uniform_int(unsigned int prog, const char *name, int val)
+{
+	BEGIN_UNIFORM_CODE {
+		glUniform1i(loc, val);
+	}
+	END_UNIFORM_CODE;
+}
+
+int set_uniform_float(unsigned int prog, const char *name, float val)
+{
+	BEGIN_UNIFORM_CODE {
+		glUniform1f(loc, val);
+	}
+	END_UNIFORM_CODE;
+}
+
+int set_uniform_float3(unsigned int prog, const char *name, float x, float y, float z)
+{
+	BEGIN_UNIFORM_CODE {
+		glUniform3f(loc, x, y, z);
+	}
+	END_UNIFORM_CODE;
+}
+
+int set_uniform_float4(unsigned int prog, const char *name, float x, float y, float z, float w)
+{
+	BEGIN_UNIFORM_CODE {
+		glUniform4f(loc, x, y, z, w);
+	}
+	END_UNIFORM_CODE;
+}
+
+int set_uniform_matrix4(unsigned int prog, const char *name, float *mat)
+{
+	BEGIN_UNIFORM_CODE {
+		glUniformMatrix4fv(loc, 1, GL_FALSE, mat);
+	}
+	END_UNIFORM_CODE;
+}
+
+int set_uniform_matrix4_transposed(unsigned int prog, const char *name, float *mat)
+{
+	BEGIN_UNIFORM_CODE {
+		glUniformMatrix4fv(loc, 1, GL_TRUE, mat);
+	}
+	END_UNIFORM_CODE;
+}
+
+int get_attrib_loc(unsigned int prog, const char *name)
+{
+	int loc, curr_prog;
+
+	glGetIntegerv(GL_CURRENT_PROGRAM, &curr_prog);
+	if(curr_prog != prog && bind_program(prog) == -1) {
+		return -1;
+	}
+
+	loc = glGetAttribLocation(prog, (char*)name);
+
+	if(curr_prog != prog) {
+		bind_program(curr_prog);
+	}
+	return loc;
+}
+
+void set_attrib_float3(int attr_loc, float x, float y, float z)
+{
+	glVertexAttrib3f(attr_loc, x, y, z);
+}
+
+#endif	/* !NO_SHADERS */
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/src/sdr.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/src/sdr.h	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,49 @@
+#ifndef SDR_H_
+#define SDR_H_
+
+#ifdef __cplusplus
+extern "C" {
+#endif	/* __cplusplus */
+
+/* ---- shaders ---- */
+unsigned int create_vertex_shader(const char *src);
+unsigned int create_pixel_shader(const char *src);
+unsigned int create_shader(const char *src, unsigned int sdr_type);
+void free_shader(unsigned int sdr);
+
+unsigned int load_vertex_shader(const char *fname);
+unsigned int load_pixel_shader(const char *fname);
+unsigned int load_shader(const char *src, unsigned int sdr_type);
+
+unsigned int get_vertex_shader(const char *fname);
+unsigned int get_pixel_shader(const char *fname);
+unsigned int get_shader(const char *fname, unsigned int sdr_type);
+
+int add_shader(const char *fname, unsigned int sdr);
+int remove_shader(const char *fname);
+
+/* ---- gpu programs ---- */
+unsigned int create_program(void);
+unsigned int create_program_link(unsigned int vs, unsigned int ps);
+unsigned int create_program_load(const char *vfile, const char *pfile);
+void free_program(unsigned int sdr);
+
+void attach_shader(unsigned int prog, unsigned int sdr);
+int link_program(unsigned int prog);
+int bind_program(unsigned int prog);
+
+int set_uniform_int(unsigned int prog, const char *name, int val);
+int set_uniform_float(unsigned int prog, const char *name, float val);
+int set_uniform_float3(unsigned int prog, const char *name, float x, float y, float z);
+int set_uniform_float4(unsigned int prog, const char *name, float x, float y, float z, float w);
+int set_uniform_matrix4(unsigned int prog, const char *name, float *mat);
+int set_uniform_matrix4_transposed(unsigned int prog, const char *name, float *mat);
+
+int get_attrib_loc(unsigned int prog, const char *name);
+void set_attrib_float3(int attr_loc, float x, float y, float z);
+
+#ifdef __cplusplus
+}
+#endif	/* __cplusplus */
+
+#endif	/* SDR_H_ */
diff -r 9ab057fba0c5 -r 52664d3451ad examples/volume/src/volume.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/volume/src/volume.c	Tue Oct 25 13:30:03 2011 +0300
@@ -0,0 +1,391 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>
+
+#ifndef NO_SHADERS
+#include <GL/glew.h>
+#include "sdr.h"
+#endif
+
+#ifndef __APPLE__
+#include <GL/glut.h>
+#else
+#include <GLUT/glut.h>
+#endif
+
+#include <imago2.h>
+
+#include "cam.h"
+#include "metasurf.h"
+
+float eval(float x, float y, float z);
+void vertex(float x, float y, float z);
+void render(void);
+void disp(void);
+void reshape(int x, int y);
+void keyb(unsigned char key, int x, int y);
+void mouse(int bn, int state, int x, int y);
+void motion(int x, int y);
+int parse_args(int argc, char **argv);
+
+int stereo, fullscreen;
+int orig_xsz, orig_ysz;
+
+struct metasurface *msurf;
+float threshold = 0.5;
+#ifndef NO_SHADERS
+unsigned int sdr;
+#endif
+
+struct img_pixmap *volume;
+int xres, yres, num_slices;
+
+int dlist, need_update = 1;
+
+int main(int argc, char **argv)
+{
+	float amb[] = {0, 0, 0, 0};
+	float lpos[] = {-0.2, 0.2, 1, 0};
+
+	glutInitWindowSize(1280, 720);
+	glutInit(&argc, argv);
+
+	if(parse_args(argc, argv) == -1) {
+		return 1;
+	}
+	glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE | (stereo ? GLUT_STEREO : 0));
+	glutCreateWindow("metasurf - volume rendering");
+
+	orig_xsz = glutGet(GLUT_WINDOW_WIDTH);
+	orig_ysz = glutGet(GLUT_WINDOW_HEIGHT);
+
+	if(fullscreen) {
+		glutFullScreen();
+	}
+
+	glutDisplayFunc(disp);
+	glutReshapeFunc(reshape);
+	glutKeyboardFunc(keyb);
+	glutMouseFunc(mouse);
+	glutMotionFunc(motion);
+
+#ifndef NO_SHADERS
+	glewInit();
+	if(!(sdr = create_program_load("sdr/vert.glsl", "sdr/frag.glsl"))) {
+		return 1;
+	}
+#endif
+
+	glEnable(GL_CULL_FACE);
+	glEnable(GL_DEPTH_TEST);
+
+	glLightModelfv(GL_LIGHT_MODEL_AMBIENT, amb);
+
+	glEnable(GL_LIGHTING);
+	glEnable(GL_LIGHT0);
+	glLightfv(GL_LIGHT0, GL_POSITION, lpos);
+
+	glEnable(GL_NORMALIZE);
+
+	cam_focus_dist(3.0);
+	cam_clip(0.1, 200.0);
+	cam_rotate(0, 0);
+	cam_dolly(2);
+
+	msurf = msurf_create();
+	msurf_eval_func(msurf, eval);
+	msurf_vertex_func(msurf, vertex);
+	msurf_threshold(msurf, threshold);
+	msurf_resolution(msurf, xres, yres, num_slices);
+	msurf_bounds(msurf, -1, -1, -1, 1, 1, 1);
+
+	glClearColor(0.6, 0.6, 0.6, 1.0);
+
+	dlist = glGenLists(1);
+
+	glutMainLoop();
+	return 0;
+}
+
+float eval(float x, float y, float z)
+{
+	int px, py, slice;
+	struct img_pixmap *img;
+
+	px = round((x * 0.5 + 0.5) * xres);
+	py = round((y * 0.5 + 0.5) * yres);
+	slice = round((z * 0.5 + 0.5) * num_slices);
+
+	if(px < 0) px = 0;
+	if(px >= xres) px = xres - 1;
+
+	if(py < 0) py = 0;
+	if(py >= yres) py = yres - 1;
+
+	if(slice < 0) slice = 0;
+	if(slice >= num_slices) slice = num_slices - 1;
+
+	img = volume + slice;
+	return *((unsigned char*)img->pixels + py * img->width + px) / 255.0;
+}
+
+void vertex(float x, float y, float z)
+{
+	float dx = 1.0 / xres;
+	float dy = 1.0 / yres;
+	float dz = 1.0 / num_slices;
+	float dfdx = eval(x - dx, y, z) - eval(x + dx, y, z);
+	float dfdy = eval(x, y - dy, z) - eval(x, y + dy, z);
+	float dfdz = eval(x, y, z - dz) - eval(x, y, z + dz);
+
+	glNormal3f(dfdx, dfdy, dfdz);
+	glVertex3f(x, y, z);
+}
+
+void render(void)
+{
+	float kd[] = {0.87, 0.82, 0.74, 1.0};
+	float ks[] = {0.9, 0.9, 0.9, 1.0};
+
+	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, kd);
+	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, ks);
+	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 60.0);
+
+#ifndef NO_SHADERS
+	bind_program(sdr);
+#endif
+
+	glMatrixMode(GL_MODELVIEW);
+	glPushMatrix();
+	glRotatef(90, 1, 0, 0);
+
+	if(need_update) {
+		glNewList(dlist, GL_COMPILE);
+		glBegin(GL_TRIANGLES);
+		printf("generating mesh... ");
+		fflush(stdout);
+		msurf_polygonize(msurf);
+		glEnd();
+		glEndList();
+		need_update = 0;
+		printf("done\n");
+	}
+	glCallList(dlist);
+
+	glPopMatrix();
+
+	assert(glGetError() == GL_NO_ERROR);
+}
+
+void disp(void)
+{
+	if(stereo) {
+		glDrawBuffer(GL_BACK_LEFT);
+	}
+
+	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+	glMatrixMode(GL_PROJECTION);
+	glLoadIdentity();
+	cam_stereo_proj_matrix(stereo ? CAM_LEFT : CAM_CENTER);
+
+	glMatrixMode(GL_MODELVIEW);
+	glLoadIdentity();
+	cam_stereo_view_matrix(stereo ? CAM_LEFT : CAM_CENTER);
+
+	render();
+
+	if(stereo) {
+		glDrawBuffer(GL_BACK_RIGHT);
+		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+		glMatrixMode(GL_PROJECTION);
+		glLoadIdentity();
+		cam_stereo_proj_matrix(CAM_RIGHT);
+
+		glMatrixMode(GL_MODELVIEW);
+		glLoadIdentity();
+		cam_stereo_view_matrix(CAM_RIGHT);
+
+		render();
+	}
+	glutSwapBuffers();
+}
+
+void reshape(int x, int y)
+{
+	glViewport(0, 0, x, y);
+	cam_aspect((float)x / (float)y);
+}
+
+void keyb(unsigned char key, int x, int y)
+{
+	static int wire;
+
+	switch(key) {
+	case 27:
+		exit(0);
+
+	case 'f':
+		fullscreen = !fullscreen;
+		if(fullscreen) {
+			glutFullScreen();
+		} else {
+			glutReshapeWindow(orig_xsz, orig_ysz);
+		}
+		break;
+
+	case 's':
+		stereo = !stereo;
+		glutPostRedisplay();
+		break;
+
+	case 'w':
+		wire = !wire;
+		glPolygonMode(GL_FRONT_AND_BACK, wire ? GL_LINE : GL_FILL);
+		glutPostRedisplay();
+		break;
+
+	case '=':
+		threshold += 0.05;
+		msurf_threshold(msurf, threshold);
+		printf("threshold: %f\n", threshold);
+		glutPostRedisplay();
+		need_update = 1;
+		break;
+
+	case '-':
+		threshold -= 0.05;
+		msurf_threshold(msurf, threshold);
+		printf("threshold: %f\n", threshold);
+		glutPostRedisplay();
+		need_update = 1;
+		break;
+
+	default:
+		break;
+	}
+}
+
+int bnstate[32];
+int prev_x, prev_y;
+
+void mouse(int bn, int state, int x, int y)
+{
+	bnstate[bn] = state == GLUT_DOWN;
+	prev_x = x;
+	prev_y = y;
+}
+
+void motion(int x, int y)
+{
+	int dx, dy;
+
+	dx = x - prev_x;
+	dy = y - prev_y;
+	prev_x = x;
+	prev_y = y;
+
+	if(bnstate[GLUT_LEFT_BUTTON]) {
+		cam_inp_rotate(dx, dy);
+	}
+	if(bnstate[GLUT_RIGHT_BUTTON]) {
+		cam_inp_zoom(dy);
+	}
+	glutPostRedisplay();
+}
+
+struct list_node {
+	struct img_pixmap img;
+	struct list_node *next;
+};
+
+int parse_args(int argc, char **argv)
+{
+	int i;
+	struct list_node *head = 0, *tail = 0;
+
+	for(i=1; i<argc; i++) {
+		if(argv[i][0] == '-' && argv[i][2] == 0) {
+			switch(argv[i][1]) {
+			case 'f':
+				fullscreen = !fullscreen;
+				break;
+
+			case 's':
+				stereo = !stereo;
+				break;
+
+			case 'h':
+				printf("usage: %s [opt]\n", argv[0]);
+				printf("options:\n");
+				printf("  -f    start in fullscreen\n");
+				printf("  -s    enable stereoscopic rendering\n");
+				printf("  -h    print usage and exit\n");
+				exit(0);
+
+			default:
+				fprintf(stderr, "unrecognized option: %s\n", argv[i]);
+				return -1;
+			}
+		} else {
+			struct list_node *slice;
+
+			if(!(slice = malloc(sizeof *slice))) {
+				fprintf(stderr, "failed to allocate volume slice: %d\n", num_slices);
+				return -1;
+			}
+			slice->next = 0;
+
+			img_init(&slice->img);
+			if(img_load(&slice->img, argv[i]) == -1) {
+				fprintf(stderr, "failed to load volume slice %d: %s\n", num_slices, argv[i]);
+				free(slice);
+				return -1;
+			}
+			img_convert(&slice->img, IMG_FMT_GREY8);
+
+			if(num_slices > 0 && (xres != slice->img.width || yres != slice->img.height)) {
+				fprintf(stderr, "error: slice %d (%s) is %dx%d, up to now we had %dx%d images\n", num_slices, argv[i],
+						slice->img.width, slice->img.height, xres, yres);
+				img_destroy(&slice->img);
+				free(slice);
+				return -1;
+			}
+			xres = slice->img.width;
+			yres = slice->img.height;
+
+			if(head) {
+				tail->next = slice;
+				tail = slice;
+			} else {
+				head = tail = slice;
+			}
+			printf("loaded volume slice %d: %s\n", num_slices++, argv[i]);
+		}
+	}
+
+	if(!head) {
+		fprintf(stderr, "you must specify a list of images for the volume data slices\n");
+		return -1;
+	}
+
+	if(!(volume = malloc(num_slices * sizeof *volume))) {
+		fprintf(stderr, "failed to allocate volume data (%d slices)\n", num_slices);
+		return -1;
+	}
+
+	for(i=0; i<num_slices; i++) {
+		void *tmp;
+
+		assert(head);
+		volume[i] = head->img;
+
+		tmp = head;
+		head = head->next;
+		free(tmp);
+	}
+
+	return 0;
+}