1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/volray.cc	Mon Apr 09 23:42:57 2012 +0300
     1.3 @@ -0,0 +1,490 @@
     1.4 +#include <stdio.h>
     1.5 +#include <stdlib.h>
     1.6 +#include <assert.h>
     1.7 +
     1.8 +#include <GL/glew.h>
     1.9 +#ifndef __APPLE__
    1.10 +#include <GL/glut.h>
    1.11 +#else
    1.12 +#include <GLUT/glut.h>
    1.13 +#endif
    1.14 +
    1.15 +#include <vmath/vmath.h>
    1.16 +#include <imago2.h>
    1.17 +#include "sdr.h"
    1.18 +#include "volume.h"
    1.19 +
    1.20 +#define XFER_MAP_SZ	512
    1.21 +
    1.22 +static void render_volume();
    1.23 +static void draw_slice();
    1.24 +static void draw_xfer_func();
    1.25 +
    1.26 +/*
    1.27 +void keyb(unsigned char key, int x, int y);
    1.28 +void keyb_up(unsigned char key, int x, int y);
    1.29 +void mouse(int bn, int state, int x, int y);
    1.30 +void motion(int x, int y);
    1.31 +int parse_args(int argc, char **argv);
    1.32 +*/
    1.33 +
    1.34 +static void create_ray_texture(int xsz, int ysz, float vfov, Vector2 *tex_scale);
    1.35 +static Vector3 get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg);
    1.36 +static int round_pow2(int x);
    1.37 +static void create_transfer_map(float mean, float sdev);
    1.38 +
    1.39 +static float cam_theta = 0, cam_phi = 0, cam_dist = 4.0;
    1.40 +static float cam_x, cam_y, cam_z;
    1.41 +
    1.42 +static vec2_t tex_scale;
    1.43 +static unsigned int vol_sdr, slice_sdr, ray_tex;
    1.44 +static int win_xsz, win_ysz;
    1.45 +static bool raytex_needs_recalc = true;
    1.46 +
    1.47 +static unsigned int xfer_tex;
    1.48 +static float xfer_mean = 0.7, xfer_sdev = 0.1;
    1.49 +static bool xfertex_needs_recalc = true;
    1.50 +
    1.51 +static float cur_z = 0.0;
    1.52 +static float ray_step = 0.01;
    1.53 +
    1.54 +static Volume volume;
    1.55 +
    1.56 +
    1.57 +bool volray_init()
    1.58 +{
    1.59 +	glewInit();
    1.60 +
    1.61 +	if(!(vol_sdr = create_program_load("sdr/volray.v.glsl", "sdr/volray.p.glsl"))) {
    1.62 +		return false;
    1.63 +	}
    1.64 +	set_uniform_int(vol_sdr, "volume", 0);
    1.65 +	set_uniform_int(vol_sdr, "ray_tex", 1);
    1.66 +	set_uniform_int(vol_sdr, "xfer_tex", 2);
    1.67 +	set_uniform_float(vol_sdr, "ray_step", ray_step);
    1.68 +	set_uniform_float(vol_sdr, "zclip", cur_z);
    1.69 +
    1.70 +	if(!(slice_sdr = create_program_load(0, "sdr/slice.p.glsl"))) {
    1.71 +		return false;
    1.72 +	}
    1.73 +	set_uniform_int(slice_sdr, "volume", 0);
    1.74 +	set_uniform_int(slice_sdr, "xfer_tex", 1);
    1.75 +
    1.76 +	if(!volume.load(fname)) {
    1.77 +		return false;
    1.78 +	}
    1.79 +	return true;
    1.80 +}
    1.81 +
    1.82 +void volray_draw(void)
    1.83 +{
    1.84 +	/* recalculate primary ray texture if needed */
    1.85 +	if(raytex_needs_recalc) {
    1.86 +		create_ray_texture(win_xsz, win_ysz, 50.0, &tex_scale);
    1.87 +	}
    1.88 +	/* recalculate transfer function texture if needed */
    1.89 +	if(xfertex_needs_recalc) {
    1.90 +		create_transfer_map(xfer_mean, xfer_sdev);
    1.91 +	}
    1.92 +
    1.93 +	render_volume();
    1.94 +	draw_slice();
    1.95 +	draw_xfer_func();
    1.96 +
    1.97 +	assert(glGetError() == GL_NO_ERROR);
    1.98 +}
    1.99 +
   1.100 +static void render_volume(void)
   1.101 +{
   1.102 +	/* set the camera transformation */
   1.103 +	glMatrixMode(GL_MODELVIEW);
   1.104 +	glPushMatrix();
   1.105 +	glLoadIdentity();
   1.106 +	glRotatef(-90, 1, 0, 0);
   1.107 +	glTranslatef(cam_x, cam_y, -cam_z);
   1.108 +	glRotatef(cam_theta, 0, 1, 0);
   1.109 +	glRotatef(cam_phi, 1, 0, 0);
   1.110 +	glTranslatef(0, 0, -cam_dist);
   1.111 +
   1.112 +	/* setup the texture matrix to map the useful part of the ray texture to [0,1] */
   1.113 +	glMatrixMode(GL_TEXTURE);
   1.114 +	glPushMatrix();
   1.115 +	glLoadIdentity();
   1.116 +	glScalef(tex_scale.x, tex_scale.y, 1.0);
   1.117 +
   1.118 +	/* tex unit0: volume data 3D texture */
   1.119 +	glActiveTexture(GL_TEXTURE0);
   1.120 +	glBindTexture(GL_TEXTURE_3D, volume.get_texture());
   1.121 +	glEnable(GL_TEXTURE_3D);
   1.122 +
   1.123 +	/* tex unit1: primary rays in view space */
   1.124 +	glActiveTexture(GL_TEXTURE1);
   1.125 +	glBindTexture(GL_TEXTURE_2D, ray_tex);
   1.126 +	glEnable(GL_TEXTURE_2D);
   1.127 +
   1.128 +	/* tex unit2: transfer function (1d) */
   1.129 +	glActiveTexture(GL_TEXTURE2);
   1.130 +	glBindTexture(GL_TEXTURE_1D, xfer_tex);
   1.131 +	glEnable(GL_TEXTURE_1D);
   1.132 +
   1.133 +	bind_program(vol_sdr);
   1.134 +	glBegin(GL_QUADS);
   1.135 +	glColor3f(1, 1, 1);
   1.136 +	glTexCoord2f(0, 1); glVertex2f(-1, -1);
   1.137 +	glTexCoord2f(1, 1); glVertex2f(1, -1);
   1.138 +	glTexCoord2f(1, 0); glVertex2f(1, 1);
   1.139 +	glTexCoord2f(0, 0); glVertex2f(-1, 1);
   1.140 +	glEnd();
   1.141 +	bind_program(0);
   1.142 +
   1.143 +	glActiveTexture(GL_TEXTURE2);
   1.144 +	glDisable(GL_TEXTURE_1D);
   1.145 +	glActiveTexture(GL_TEXTURE1);
   1.146 +	glDisable(GL_TEXTURE_2D);
   1.147 +	glActiveTexture(GL_TEXTURE0);
   1.148 +	glDisable(GL_TEXTURE_3D);
   1.149 +
   1.150 +	glMatrixMode(GL_TEXTURE);
   1.151 +	glPopMatrix();
   1.152 +	glMatrixMode(GL_MODELVIEW);
   1.153 +	glPopMatrix();
   1.154 +}
   1.155 +
   1.156 +static void draw_slice(void)
   1.157 +{
   1.158 +	glMatrixMode(GL_MODELVIEW);
   1.159 +	glPushMatrix();
   1.160 +	glTranslatef(0.9, 0.9, 0);
   1.161 +	glScalef(0.3, 0.3 * ((float)win_xsz / win_ysz), 1);
   1.162 +	glTranslatef(-1, -1, 0);
   1.163 +
   1.164 +	glActiveTexture(GL_TEXTURE0);
   1.165 +	glBindTexture(GL_TEXTURE_3D, volume.get_texture());
   1.166 +	glEnable(GL_TEXTURE_3D);
   1.167 +
   1.168 +	glActiveTexture(GL_TEXTURE1);
   1.169 +	glBindTexture(GL_TEXTURE_1D, xfer_tex);
   1.170 +	glEnable(GL_TEXTURE_1D);
   1.171 +
   1.172 +	bind_program(slice_sdr);
   1.173 +
   1.174 +	glBegin(GL_QUADS);
   1.175 +	glColor3f(1, 1, 1);
   1.176 +	glTexCoord3f(0, 1, cur_z); glVertex2f(-1, -1);
   1.177 +	glTexCoord3f(1, 1, cur_z); glVertex2f(1, -1);
   1.178 +	glTexCoord3f(1, 0, cur_z); glVertex2f(1, 1);
   1.179 +	glTexCoord3f(0, 0, cur_z); glVertex2f(-1, 1);
   1.180 +	glEnd();
   1.181 +
   1.182 +	bind_program(0);
   1.183 +
   1.184 +	glActiveTexture(GL_TEXTURE1);
   1.185 +	glDisable(GL_TEXTURE_1D);
   1.186 +	glActiveTexture(GL_TEXTURE0);
   1.187 +	glDisable(GL_TEXTURE_3D);
   1.188 +	glPopMatrix();
   1.189 +}
   1.190 +
   1.191 +static void draw_xfer_func(void)
   1.192 +{
   1.193 +	glMatrixMode(GL_MODELVIEW);
   1.194 +	glPushMatrix();
   1.195 +	glTranslatef(-0.9, -0.9, 0);
   1.196 +	glScalef(0.5, 0.1, 1);
   1.197 +
   1.198 +	glBindTexture(GL_TEXTURE_1D, xfer_tex);
   1.199 +	glEnable(GL_TEXTURE_1D);
   1.200 +
   1.201 +	glBegin(GL_QUADS);
   1.202 +	glColor3f(1, 1, 1);
   1.203 +	glTexCoord1f(1);
   1.204 +	glVertex2f(1, 0);
   1.205 +	glVertex2f(1, 1);
   1.206 +	glTexCoord1f(0);
   1.207 +	glVertex2f(0, 1);
   1.208 +	glVertex2f(0, 0);
   1.209 +	glEnd();
   1.210 +
   1.211 +	glDisable(GL_TEXTURE_1D);
   1.212 +
   1.213 +	glLineWidth(2.0);
   1.214 +	glBegin(GL_LINE_LOOP);
   1.215 +	if(uimode == UIMODE_XFER) {
   1.216 +		glColor3f(1, 0, 0);
   1.217 +	} else {
   1.218 +		glColor3f(0, 0, 1);
   1.219 +	}
   1.220 +	glVertex2f(0, 0);
   1.221 +	glVertex2f(1, 0);
   1.222 +	glVertex2f(1, 1);
   1.223 +	glVertex2f(0, 1);
   1.224 +	glEnd();
   1.225 +
   1.226 +	glPopMatrix();
   1.227 +}
   1.228 +
   1.229 +void volray_resize(int x, int y)
   1.230 +{
   1.231 +	glViewport(0, 0, x, y);
   1.232 +
   1.233 +	if(x != win_xsz || y != win_ysz) {
   1.234 +		raytex_needs_recalc = true;
   1.235 +		win_xsz = x;
   1.236 +		win_ysz = y;
   1.237 +	}
   1.238 +}
   1.239 +
   1.240 +#if 0
   1.241 +void keyb(unsigned char key, int x, int y)
   1.242 +{
   1.243 +	switch(key) {
   1.244 +	case 27:
   1.245 +		exit(0);
   1.246 +
   1.247 +	case 'x':
   1.248 +		uimode = UIMODE_XFER;
   1.249 +		glutPostRedisplay();
   1.250 +		break;
   1.251 +
   1.252 +	case 'c':
   1.253 +		uimode = UIMODE_CURSOR;
   1.254 +		glutPostRedisplay();
   1.255 +		break;
   1.256 +
   1.257 +	default:
   1.258 +		break;
   1.259 +	}
   1.260 +}
   1.261 +
   1.262 +void keyb_up(unsigned char key, int x, int y)
   1.263 +{
   1.264 +	switch(key) {
   1.265 +	case 'x':
   1.266 +		if(uimode == UIMODE_XFER) {
   1.267 +			uimode = UIMODE_DEFAULT;
   1.268 +			glutPostRedisplay();
   1.269 +		}
   1.270 +		break;
   1.271 +
   1.272 +	case 'c':
   1.273 +		if(uimode == UIMODE_CURSOR) {
   1.274 +			uimode = UIMODE_DEFAULT;
   1.275 +			glutPostRedisplay();
   1.276 +		}
   1.277 +		break;
   1.278 +
   1.279 +	default:
   1.280 +		break;
   1.281 +	}
   1.282 +}
   1.283 +
   1.284 +static int bnstate[32];
   1.285 +static int prev_x, prev_y;
   1.286 +
   1.287 +void mouse(int bn, int state, int x, int y)
   1.288 +{
   1.289 +	bnstate[bn - GLUT_LEFT_BUTTON] = state == GLUT_DOWN;
   1.290 +	prev_x = x;
   1.291 +	prev_y = y;
   1.292 +}
   1.293 +
   1.294 +void motion(int x, int y)
   1.295 +{
   1.296 +	int dx = x - prev_x;
   1.297 +	int dy = y - prev_y;
   1.298 +	prev_x = x;
   1.299 +	prev_y = y;
   1.300 +
   1.301 +	switch(uimode) {
   1.302 +	case UIMODE_XFER:
   1.303 +		if(dx || dy) {
   1.304 +			xfer_mean += dx / (float)win_xsz;
   1.305 +			xfer_sdev += 0.5 * dy / (float)win_ysz;
   1.306 +
   1.307 +			xfer_mean = xfer_mean < 0.0 ? 0.0 : (xfer_mean > 1.0 ? 1.0 : xfer_mean);
   1.308 +			xfer_sdev = xfer_sdev < 0.0 ? 0.0 : (xfer_sdev > 1.0 ? 1.0 : xfer_sdev);
   1.309 +
   1.310 +			xfertex_needs_recalc = true;
   1.311 +			glutPostRedisplay();
   1.312 +		}
   1.313 +		break;
   1.314 +
   1.315 +	case UIMODE_CURSOR:
   1.316 +		cur_z += 0.5 * dy / (float)win_ysz;
   1.317 +
   1.318 +		if(cur_z < 0.0)
   1.319 +			cur_z = 0.0;
   1.320 +		if(cur_z > 1.0)
   1.321 +			cur_z = 1.0;
   1.322 +
   1.323 +		set_uniform_float(vol_sdr, "zclip", cur_z);
   1.324 +		glutPostRedisplay();
   1.325 +		break;
   1.326 +
   1.327 +	default:
   1.328 +		/* view control */
   1.329 +		if(bnstate[0]) {
   1.330 +			cam_theta += dx * 0.5;
   1.331 +			cam_phi += dy * 0.5;
   1.332 +
   1.333 +			if(cam_phi <= -90) cam_phi = -89;
   1.334 +			if(cam_phi >= 90) cam_phi = 89;
   1.335 +
   1.336 +			glutPostRedisplay();
   1.337 +		}
   1.338 +
   1.339 +		if(bnstate[1]) {
   1.340 +			cam_x += dx * 0.025;
   1.341 +			cam_y += dy * 0.025;
   1.342 +			glutPostRedisplay();
   1.343 +		}
   1.344 +
   1.345 +		if(bnstate[2]) {
   1.346 +			cam_dist += dy * 0.025;
   1.347 +			if(cam_dist < 0.0) cam_dist = 0.0;
   1.348 +			glutPostRedisplay();
   1.349 +		}
   1.350 +	}
   1.351 +}
   1.352 +
   1.353 +
   1.354 +int parse_args(int argc, char **argv)
   1.355 +{
   1.356 +	int i;
   1.357 +	char *endp;
   1.358 +
   1.359 +	for(i=1; i<argc; i++) {
   1.360 +		if(argv[i][0] == '-' && argv[i][2] == 0) {
   1.361 +			switch(argv[i][1]) {
   1.362 +			case 'm':
   1.363 +				xfer_mean = strtod(argv[++i], &endp);
   1.364 +				if(endp == argv[i]) {
   1.365 +					fprintf(stderr, "-m must be followed by the transfer function mean\n");
   1.366 +					return -1;
   1.367 +				}
   1.368 +				break;
   1.369 +
   1.370 +			case 'd':
   1.371 +				xfer_sdev = strtod(argv[++i], &endp);
   1.372 +				if(endp == argv[i]) {
   1.373 +					fprintf(stderr, "-d must be followed by the transfer function std.deviation\n");
   1.374 +					return -1;
   1.375 +				}
   1.376 +				break;
   1.377 +
   1.378 +			default:
   1.379 +				fprintf(stderr, "unrecognized option: %s\n", argv[i]);
   1.380 +				return -1;
   1.381 +			}
   1.382 +		} else {
   1.383 +			if(fname) {
   1.384 +				fprintf(stderr, "unexpected argument: %s\n", argv[i]);
   1.385 +				return -1;
   1.386 +			}
   1.387 +			fname = argv[i];
   1.388 +		}
   1.389 +	}
   1.390 +
   1.391 +	if(!fname) {
   1.392 +		fprintf(stderr, "pass the volume descriptor filename\n");
   1.393 +		return -1;
   1.394 +	}
   1.395 +	return 0;
   1.396 +}
   1.397 +#endif
   1.398 +
   1.399 +
   1.400 +static void create_ray_texture(int xsz, int ysz, float vfov, Vector2 *tex_scale)
   1.401 +{
   1.402 +	int cur_tex_xsz, cur_tex_ysz;
   1.403 +	int tex_xsz = round_pow2(xsz);
   1.404 +	int tex_ysz = round_pow2(ysz);
   1.405 +	float *teximg, *dir;
   1.406 +
   1.407 +	teximg = new float[3 * xsz * ysz];
   1.408 +	dir = teximg;
   1.409 +
   1.410 +	for(int i=0; i<ysz; i++) {
   1.411 +		for(int j=0; j<xsz; j++) {
   1.412 +			Vector3 rdir = get_primary_ray_dir(j, i, xsz, ysz, vfov);
   1.413 +			*dir++ = rdir.x;
   1.414 +			*dir++ = rdir.y;
   1.415 +			*dir++ = rdir.z;
   1.416 +		}
   1.417 +	}
   1.418 +
   1.419 +	if(!ray_tex) {
   1.420 +		glGenTextures(1, &ray_tex);
   1.421 +	}
   1.422 +
   1.423 +	glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_WIDTH, &cur_tex_xsz);
   1.424 +	glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_HEIGHT, &cur_tex_ysz);
   1.425 +
   1.426 +	if(tex_xsz > cur_tex_xsz || tex_ysz > cur_tex_ysz) {
   1.427 +		glBindTexture(GL_TEXTURE_2D, ray_tex);
   1.428 +		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
   1.429 +		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
   1.430 +		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
   1.431 +		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
   1.432 +		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F_ARB, tex_xsz, tex_ysz, 0, GL_RGB, GL_FLOAT, 0);
   1.433 +	}
   1.434 +
   1.435 +	glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, xsz, ysz, GL_RGB, GL_FLOAT, teximg);
   1.436 +	delete [] teximg;
   1.437 +
   1.438 +	if(tex_scale) {
   1.439 +		tex_scale->x = (float)xsz / (float)tex_xsz;
   1.440 +		tex_scale->y = (float)ysz / (float)tex_ysz;
   1.441 +	}
   1.442 +	raytex_needs_recalc = false;
   1.443 +}
   1.444 +
   1.445 +static Vector3 get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg)
   1.446 +{
   1.447 +	float vfov = M_PI * vfov_deg / 180.0;
   1.448 +	float aspect = (float)w / (float)h;
   1.449 +
   1.450 +	float ysz = 2.0;
   1.451 +	float xsz = aspect * ysz;
   1.452 +
   1.453 +	float px = ((float)x / (float)w) * xsz - xsz / 2.0;
   1.454 +	float py = 1.0 - ((float)y / (float)h) * ysz;
   1.455 +	float pz = 1.0 / tan(0.5 * vfov);
   1.456 +
   1.457 +	float mag = sqrt(px * px + py * py + pz * pz);
   1.458 +	return Vector3(px / mag, py / mag, pz / mag);
   1.459 +}
   1.460 +
   1.461 +static int round_pow2(int x)
   1.462 +{
   1.463 +	x--;
   1.464 +	x = (x >> 1) | x;
   1.465 +	x = (x >> 2) | x;
   1.466 +	x = (x >> 4) | x;
   1.467 +	x = (x >> 8) | x;
   1.468 +	x = (x >> 16) | x;
   1.469 +	return x + 1;
   1.470 +}
   1.471 +
   1.472 +static void create_transfer_map(float mean, float sdev)
   1.473 +{
   1.474 +	static float map[XFER_MAP_SZ];
   1.475 +
   1.476 +	if(!xfer_tex) {
   1.477 +		glGenTextures(1, &xfer_tex);
   1.478 +		glBindTexture(GL_TEXTURE_1D, xfer_tex);
   1.479 +		glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
   1.480 +		glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
   1.481 +		glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
   1.482 +		glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
   1.483 +		glTexImage1D(GL_TEXTURE_1D, 0, GL_LUMINANCE32F_ARB, XFER_MAP_SZ, 0, GL_LUMINANCE, GL_FLOAT, 0);
   1.484 +	}
   1.485 +
   1.486 +	for(int i=0; i<XFER_MAP_SZ; i++) {
   1.487 +		float x = (float)i / (float)(XFER_MAP_SZ - 1);
   1.488 +		map[i] = gaussian(x, mean, sdev) - 1.0;
   1.489 +	}
   1.490 +
   1.491 +	glTexSubImage1D(GL_TEXTURE_1D, 0, 0, XFER_MAP_SZ, GL_LUMINANCE, GL_FLOAT, map);
   1.492 +	xfertex_needs_recalc = false;
   1.493 +}