1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/volray.c	Sat Mar 31 02:08:42 2012 +0300
     1.3 @@ -0,0 +1,332 @@
     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 +
    1.19 +struct slice_file {
    1.20 +	char *name;
    1.21 +	struct slice_file *next;
    1.22 +};
    1.23 +
    1.24 +int init(void);
    1.25 +void disp(void);
    1.26 +void reshape(int x, int y);
    1.27 +void keyb(unsigned char key, int x, int y);
    1.28 +void mouse(int bn, int state, int x, int y);
    1.29 +void motion(int x, int y);
    1.30 +int parse_args(int argc, char **argv);
    1.31 +
    1.32 +unsigned int create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale);
    1.33 +static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg);
    1.34 +static int round_pow2(int x);
    1.35 +
    1.36 +float cam_theta = 0, cam_phi = 0, cam_dist = 4.0;
    1.37 +float cam_x, cam_y, cam_z;
    1.38 +
    1.39 +vec2_t tex_scale;
    1.40 +struct slice_file *flist;
    1.41 +int nslices;
    1.42 +unsigned int sdr, vol_tex, ray_tex;
    1.43 +int win_xsz, win_ysz;
    1.44 +
    1.45 +int main(int argc, char **argv)
    1.46 +{
    1.47 +	glutInit(&argc, argv);
    1.48 +	glutInitWindowSize(1280, 720);
    1.49 +	glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    1.50 +
    1.51 +	if(parse_args(argc, argv) == -1) {
    1.52 +		return 1;
    1.53 +	}
    1.54 +
    1.55 +	glutCreateWindow("Volume Raytracer");
    1.56 +
    1.57 +	glutDisplayFunc(disp);
    1.58 +	glutReshapeFunc(reshape);
    1.59 +	glutKeyboardFunc(keyb);
    1.60 +	glutMouseFunc(mouse);
    1.61 +	glutMotionFunc(motion);
    1.62 +
    1.63 +	glewInit();
    1.64 +
    1.65 +	if(init() == -1) {
    1.66 +		return 1;
    1.67 +	}
    1.68 +
    1.69 +	glutMainLoop();
    1.70 +	return 0;
    1.71 +}
    1.72 +
    1.73 +int init(void)
    1.74 +{
    1.75 +	int i, vol_xsz, vol_ysz;
    1.76 +
    1.77 +	if(!(sdr = create_program_load("volray.v.glsl", "volray.p.glsl"))) {
    1.78 +		return 1;
    1.79 +	}
    1.80 +
    1.81 +	glGenTextures(1, &vol_tex);
    1.82 +	glBindTexture(GL_TEXTURE_3D, vol_tex);
    1.83 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    1.84 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    1.85 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    1.86 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
    1.87 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
    1.88 +
    1.89 +	for(i=0; i<nslices; i++) {
    1.90 +		int xsz, ysz;
    1.91 +		void *pix;
    1.92 +		struct slice_file *sfile = flist;
    1.93 +		flist = flist->next;
    1.94 +
    1.95 +		if(!(pix = img_load_pixels(sfile->name, &xsz, &ysz, IMG_FMT_RGBA32))) {
    1.96 +			fprintf(stderr, "failed to load image: %s\n", sfile->name);
    1.97 +			return -1;
    1.98 +		}
    1.99 +
   1.100 +		if(i == 0) {
   1.101 +			/* allocate storage for the texture */
   1.102 +			glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA, xsz, ysz, nslices, 0, GL_RGBA, GL_UNSIGNED_BYTE, 0);
   1.103 +
   1.104 +			vol_xsz = xsz;
   1.105 +			vol_ysz = ysz;
   1.106 +
   1.107 +		} else {
   1.108 +			if(xsz != vol_xsz || ysz != vol_ysz) {
   1.109 +				fprintf(stderr, "%s: inconsistent slice size: %dx%d. expected: %dx%d\n",
   1.110 +						sfile->name, xsz, ysz, vol_xsz, vol_ysz);
   1.111 +				img_free_pixels(pix);
   1.112 +				return -1;
   1.113 +			}
   1.114 +		}
   1.115 +		free(sfile);
   1.116 +
   1.117 +		glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, i, xsz, ysz, 1, GL_RGBA, GL_UNSIGNED_BYTE, pix);
   1.118 +		img_free_pixels(pix);
   1.119 +	}
   1.120 +
   1.121 +	return 0;
   1.122 +}
   1.123 +
   1.124 +void disp(void)
   1.125 +{
   1.126 +	glMatrixMode(GL_MODELVIEW);
   1.127 +	glLoadIdentity();
   1.128 +	glTranslatef(cam_x, cam_y, -cam_z);
   1.129 +	glRotatef(cam_theta, 0, 1, 0);
   1.130 +	glRotatef(cam_phi, 1, 0, 0);
   1.131 +	glTranslatef(0, 0, -cam_dist);
   1.132 +
   1.133 +	glMatrixMode(GL_TEXTURE);
   1.134 +	glPushMatrix();
   1.135 +	glScalef(tex_scale.x, tex_scale.y, 1.0);
   1.136 +
   1.137 +	glActiveTexture(GL_TEXTURE0);
   1.138 +	glEnable(GL_TEXTURE_3D);
   1.139 +	glBindTexture(GL_TEXTURE_3D, vol_tex);
   1.140 +
   1.141 +	glActiveTexture(GL_TEXTURE1);
   1.142 +	glEnable(GL_TEXTURE_2D);
   1.143 +	glBindTexture(GL_TEXTURE_2D, ray_tex);
   1.144 +
   1.145 +	bind_program(sdr);
   1.146 +	glBegin(GL_QUADS);
   1.147 +	glColor3f(1, 1, 1);
   1.148 +	glTexCoord2f(0, 1);
   1.149 +	glVertex2f(-1, -1);
   1.150 +	glTexCoord2f(1, 1);
   1.151 +	glVertex2f(1, -1);
   1.152 +	glTexCoord2f(1, 0);
   1.153 +	glVertex2f(1, 1);
   1.154 +	glTexCoord2f(0, 0);
   1.155 +	glVertex2f(-1, 1);
   1.156 +	glEnd();
   1.157 +	bind_program(0);
   1.158 +
   1.159 +	glDisable(GL_TEXTURE_2D);
   1.160 +	glActiveTexture(GL_TEXTURE0);
   1.161 +	glDisable(GL_TEXTURE_3D);
   1.162 +
   1.163 +	glMatrixMode(GL_TEXTURE);
   1.164 +	glPopMatrix();
   1.165 +
   1.166 +	glutSwapBuffers();
   1.167 +	assert(glGetError() == GL_NO_ERROR);
   1.168 +}
   1.169 +
   1.170 +void reshape(int x, int y)
   1.171 +{
   1.172 +	glViewport(0, 0, x, y);
   1.173 +
   1.174 +	if(x != win_xsz || y != win_ysz) {
   1.175 +		ray_tex = create_ray_texture(x, y, 50.0, &tex_scale);
   1.176 +		win_xsz = x;
   1.177 +		win_ysz = y;
   1.178 +	}
   1.179 +}
   1.180 +
   1.181 +void keyb(unsigned char key, int x, int y)
   1.182 +{
   1.183 +	switch(key) {
   1.184 +	case 27:
   1.185 +		exit(0);
   1.186 +	}
   1.187 +}
   1.188 +
   1.189 +static int bnstate[32];
   1.190 +static int prev_x, prev_y;
   1.191 +
   1.192 +void mouse(int bn, int state, int x, int y)
   1.193 +{
   1.194 +	bnstate[bn - GLUT_LEFT_BUTTON] = state == GLUT_DOWN;
   1.195 +	prev_x = x;
   1.196 +	prev_y = y;
   1.197 +}
   1.198 +
   1.199 +void motion(int x, int y)
   1.200 +{
   1.201 +	int dx = x - prev_x;
   1.202 +	int dy = y - prev_y;
   1.203 +
   1.204 +	prev_x = x;
   1.205 +	prev_y = y;
   1.206 +
   1.207 +	if(bnstate[0]) {
   1.208 +		cam_theta += dx * 0.5;
   1.209 +		cam_phi += dy * 0.5;
   1.210 +
   1.211 +		if(cam_phi < -90) cam_phi = -90;
   1.212 +		if(cam_phi > 90) cam_phi = 90;
   1.213 +
   1.214 +		glutPostRedisplay();
   1.215 +	}
   1.216 +
   1.217 +	if(bnstate[1]) {
   1.218 +		cam_y += dy * 0.1;
   1.219 +		glutPostRedisplay();
   1.220 +	}
   1.221 +
   1.222 +	if(bnstate[2]) {
   1.223 +		cam_dist += dy * 0.1;
   1.224 +		if(cam_dist < 0.0) cam_dist = 0.0;
   1.225 +		glutPostRedisplay();
   1.226 +	}
   1.227 +}
   1.228 +
   1.229 +
   1.230 +int parse_args(int argc, char **argv)
   1.231 +{
   1.232 +	int i;
   1.233 +	struct slice_file *tail;
   1.234 +
   1.235 +	for(i=1; i<argc; i++) {
   1.236 +		struct slice_file *sfile;
   1.237 +
   1.238 +		if(!(sfile = malloc(sizeof *sfile))) {
   1.239 +			perror("failed to allocate memory");
   1.240 +			return -1;
   1.241 +		}
   1.242 +		sfile->name = argv[i];
   1.243 +		sfile->next = 0;
   1.244 +
   1.245 +		if(!flist) {
   1.246 +			flist = tail = sfile;
   1.247 +		} else {
   1.248 +			tail->next = sfile;
   1.249 +			tail = sfile;
   1.250 +		}
   1.251 +		nslices++;
   1.252 +	}
   1.253 +
   1.254 +	if(!nslices) {
   1.255 +		fprintf(stderr, "pass the slice filenames\n");
   1.256 +		return -1;
   1.257 +	}
   1.258 +
   1.259 +	return 0;
   1.260 +}
   1.261 +
   1.262 +
   1.263 +unsigned int create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale)
   1.264 +{
   1.265 +	int i, j;
   1.266 +	unsigned int tex;
   1.267 +	int tex_xsz = round_pow2(xsz);
   1.268 +	int tex_ysz = round_pow2(ysz);
   1.269 +	float *teximg, *dir;
   1.270 +
   1.271 +	if(!(teximg = malloc(3 * tex_xsz * tex_ysz * sizeof *teximg))) {
   1.272 +		return 0;
   1.273 +	}
   1.274 +	dir = teximg;
   1.275 +
   1.276 +	for(i=0; i<tex_ysz; i++) {
   1.277 +		for(j=0; j<tex_xsz; j++) {
   1.278 +			if(j < xsz && i < ysz) {
   1.279 +				vec3_t rdir = get_primary_ray_dir(j, i, xsz, ysz, vfov);
   1.280 +				dir[0] = rdir.x;
   1.281 +				dir[1] = rdir.y;
   1.282 +				dir[2] = rdir.z;
   1.283 +			} else {
   1.284 +				dir[0] = dir[1] = 0.0f;
   1.285 +				dir[2] = 1.0f;
   1.286 +			}
   1.287 +
   1.288 +			dir += 3;
   1.289 +		}
   1.290 +	}
   1.291 +
   1.292 +	glGenTextures(1, &tex);
   1.293 +	glBindTexture(GL_TEXTURE_2D, tex);
   1.294 +	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
   1.295 +	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
   1.296 +	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
   1.297 +	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
   1.298 +	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F_ARB, tex_xsz, tex_ysz, 0, GL_RGB, GL_FLOAT, teximg);
   1.299 +	free(teximg);
   1.300 +
   1.301 +	if(tex_scale) {
   1.302 +		tex_scale->x = (float)xsz / (float)tex_xsz;
   1.303 +		tex_scale->y = (float)ysz / (float)tex_ysz;
   1.304 +	}
   1.305 +	return tex;
   1.306 +}
   1.307 +
   1.308 +static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg)
   1.309 +{
   1.310 +	float vfov = M_PI * vfov_deg / 180.0;
   1.311 +	float aspect = (float)w / (float)h;
   1.312 +
   1.313 +	float ysz = 2.0;
   1.314 +	float xsz = aspect * ysz;
   1.315 +
   1.316 +	float px = ((float)x / (float)w) * xsz - xsz / 2.0;
   1.317 +	float py = 1.0 - ((float)y / (float)h) * ysz;
   1.318 +	float pz = 1.0 / tan(0.5 * vfov);
   1.319 +
   1.320 +	float mag = sqrt(px * px + py * py + pz * pz);
   1.321 +
   1.322 +	return v3_cons(px / mag, py / mag, pz / mag);
   1.323 +}
   1.324 +
   1.325 +static int round_pow2(int x)
   1.326 +{
   1.327 +	x--;
   1.328 +	x = (x >> 1) | x;
   1.329 +	x = (x >> 2) | x;
   1.330 +	x = (x >> 4) | x;
   1.331 +	x = (x >> 8) | x;
   1.332 +	x = (x >> 16) | x;
   1.333 +	return x + 1;
   1.334 +}
   1.335 +