qvolray

view src/volray.c @ 4:3e53a16d4667

foo
author John Tsiombikas <nuclear@member.fsf.org>
date Mon, 02 Apr 2012 17:59:46 +0300
parents 6f275934717b
children 0c3874aa717a
line source
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <assert.h>
5 #include <GL/glew.h>
6 #ifndef __APPLE__
7 #include <GL/glut.h>
8 #else
9 #include <GLUT/glut.h>
10 #endif
12 #include <vmath/vmath.h>
13 #include <imago2.h>
14 #include "sdr.h"
16 #define XFER_MAP_SZ 512
18 struct slice_file {
19 char *name;
20 struct slice_file *next;
21 };
23 int init(void);
24 void disp(void);
25 void render_volume(void);
26 void draw_xfer_func(void);
27 void reshape(int x, int y);
28 void keyb(unsigned char key, int x, int y);
29 void keyb_up(unsigned char key, int x, int y);
30 void mouse(int bn, int state, int x, int y);
31 void motion(int x, int y);
32 int parse_args(int argc, char **argv);
34 static void create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale);
35 static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg);
36 static int round_pow2(int x);
37 static void create_transfer_map(float mean, float sdev);
39 float cam_theta = 0, cam_phi = 0, cam_dist = 4.0;
40 float cam_x, cam_y, cam_z;
42 vec2_t tex_scale;
43 struct slice_file *flist;
44 int nslices;
45 unsigned int sdr, vol_tex, ray_tex;
46 int win_xsz, win_ysz;
47 int raytex_needs_recalc = 1;
49 unsigned int xfer_tex;
50 float xfer_mean = 0.5, xfer_sdev = 1.0;
51 int xfertex_needs_recalc = 1;
53 static int uimode_xfer;
55 int main(int argc, char **argv)
56 {
57 glutInit(&argc, argv);
58 glutInitWindowSize(1280, 720);
59 glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
61 if(parse_args(argc, argv) == -1) {
62 return 1;
63 }
65 glutCreateWindow("Volume Raytracer");
67 glutDisplayFunc(disp);
68 glutReshapeFunc(reshape);
69 glutKeyboardFunc(keyb);
70 glutKeyboardUpFunc(keyb_up);
71 glutMouseFunc(mouse);
72 glutMotionFunc(motion);
74 glewInit();
76 if(init() == -1) {
77 return 1;
78 }
80 glutMainLoop();
81 return 0;
82 }
84 int init(void)
85 {
86 int i, vol_xsz, vol_ysz;
88 if(!(sdr = create_program_load("volray.v.glsl", "volray.p.glsl"))) {
89 return -1;
90 }
91 set_uniform_int(sdr, "volume", 0);
92 set_uniform_int(sdr, "ray_tex", 1);
93 set_uniform_int(sdr, "xfer_tex", 2);
95 glGenTextures(1, &vol_tex);
96 glBindTexture(GL_TEXTURE_3D, vol_tex);
97 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
98 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
99 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
100 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
101 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
103 for(i=0; i<nslices; i++) {
104 int xsz, ysz;
105 void *pix;
106 struct slice_file *sfile = flist;
107 flist = flist->next;
109 if(!(pix = img_load_pixels(sfile->name, &xsz, &ysz, IMG_FMT_RGBA32))) {
110 fprintf(stderr, "failed to load image: %s\n", sfile->name);
111 return -1;
112 }
114 if(i == 0) {
115 /* allocate storage for the texture */
116 glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA, xsz, ysz, nslices, 0, GL_RGBA, GL_UNSIGNED_BYTE, 0);
118 vol_xsz = xsz;
119 vol_ysz = ysz;
121 } else {
122 if(xsz != vol_xsz || ysz != vol_ysz) {
123 fprintf(stderr, "%s: inconsistent slice size: %dx%d. expected: %dx%d\n",
124 sfile->name, xsz, ysz, vol_xsz, vol_ysz);
125 img_free_pixels(pix);
126 return -1;
127 }
128 }
129 free(sfile);
131 glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, i, xsz, ysz, 1, GL_RGBA, GL_UNSIGNED_BYTE, pix);
132 img_free_pixels(pix);
133 }
134 return 0;
135 }
137 void disp(void)
138 {
139 /* recalculate primary ray texture if needed */
140 if(raytex_needs_recalc) {
141 create_ray_texture(win_xsz, win_ysz, 50.0, &tex_scale);
142 }
143 /* recalculate transfer function texture if needed */
144 if(xfertex_needs_recalc) {
145 create_transfer_map(xfer_mean, xfer_sdev);
146 }
148 render_volume();
149 draw_xfer_func();
151 glutSwapBuffers();
152 assert(glGetError() == GL_NO_ERROR);
153 }
155 void render_volume(void)
156 {
157 /* set the camera transformation */
158 glMatrixMode(GL_MODELVIEW);
159 glPushMatrix();
160 glLoadIdentity();
161 glRotatef(-90, 1, 0, 0);
162 glTranslatef(cam_x, cam_y, -cam_z);
163 glRotatef(cam_theta, 0, 1, 0);
164 glRotatef(cam_phi, 1, 0, 0);
165 glTranslatef(0, 0, -cam_dist);
167 /* setup the texture matrix to map the useful part of the ray texture to [0,1] */
168 glMatrixMode(GL_TEXTURE);
169 glPushMatrix();
170 glLoadIdentity();
171 glScalef(tex_scale.x, tex_scale.y, 1.0);
173 /* tex unit0: volume data 3D texture */
174 glActiveTexture(GL_TEXTURE0);
175 glBindTexture(GL_TEXTURE_3D, vol_tex);
176 glEnable(GL_TEXTURE_3D);
178 /* tex unit1: primary rays in view space */
179 glActiveTexture(GL_TEXTURE1);
180 glBindTexture(GL_TEXTURE_2D, ray_tex);
181 glEnable(GL_TEXTURE_2D);
183 /* tex unit2: transfer function (1d) */
184 glActiveTexture(GL_TEXTURE2);
185 glBindTexture(GL_TEXTURE_1D, xfer_tex);
186 glEnable(GL_TEXTURE_1D);
188 bind_program(sdr);
189 glBegin(GL_QUADS);
190 glColor3f(1, 1, 1);
191 glTexCoord2f(0, 1); glVertex2f(-1, -1);
192 glTexCoord2f(1, 1); glVertex2f(1, -1);
193 glTexCoord2f(1, 0); glVertex2f(1, 1);
194 glTexCoord2f(0, 0); glVertex2f(-1, 1);
195 glEnd();
196 bind_program(0);
198 glActiveTexture(GL_TEXTURE2);
199 glDisable(GL_TEXTURE_1D);
200 glActiveTexture(GL_TEXTURE1);
201 glDisable(GL_TEXTURE_2D);
202 glActiveTexture(GL_TEXTURE0);
203 glDisable(GL_TEXTURE_3D);
205 glMatrixMode(GL_TEXTURE);
206 glPopMatrix();
207 glMatrixMode(GL_MODELVIEW);
208 glPopMatrix();
209 }
211 void draw_xfer_func(void)
212 {
213 glMatrixMode(GL_MODELVIEW);
214 glPushMatrix();
215 glTranslatef(-0.9, -0.9, 0);
216 glScalef(0.5, 0.1, 1);
218 glBindTexture(GL_TEXTURE_1D, xfer_tex);
219 glEnable(GL_TEXTURE_1D);
221 glBegin(GL_QUADS);
222 glColor3f(1, 1, 1);
223 glTexCoord1f(1);
224 glVertex2f(1, 0);
225 glVertex2f(1, 1);
226 glTexCoord1f(0);
227 glVertex2f(0, 1);
228 glVertex2f(0, 0);
229 glEnd();
231 glDisable(GL_TEXTURE_1D);
233 glLineWidth(2.0);
234 glBegin(GL_LINE_LOOP);
235 glColor3f(uimode_xfer ? 1 : 0, 0, uimode_xfer ? 0 : 1);
236 glVertex2f(0, 0);
237 glVertex2f(1, 0);
238 glVertex2f(1, 1);
239 glVertex2f(0, 1);
240 glEnd();
242 glPopMatrix();
243 }
245 void reshape(int x, int y)
246 {
247 glViewport(0, 0, x, y);
249 if(x != win_xsz || y != win_ysz) {
250 raytex_needs_recalc = 1;
251 win_xsz = x;
252 win_ysz = y;
253 }
254 }
256 void keyb(unsigned char key, int x, int y)
257 {
258 switch(key) {
259 case 27:
260 exit(0);
262 case 'x':
263 uimode_xfer = 1;
264 glutPostRedisplay();
265 break;
266 }
267 }
269 void keyb_up(unsigned char key, int x, int y)
270 {
271 switch(key) {
272 case 'x':
273 uimode_xfer = 0;
274 glutPostRedisplay();
275 break;
276 }
277 }
279 static int bnstate[32];
280 static int prev_x, prev_y;
282 void mouse(int bn, int state, int x, int y)
283 {
284 bnstate[bn - GLUT_LEFT_BUTTON] = state == GLUT_DOWN;
285 prev_x = x;
286 prev_y = y;
287 }
289 void motion(int x, int y)
290 {
291 int dx = x - prev_x;
292 int dy = y - prev_y;
293 prev_x = x;
294 prev_y = y;
296 if(uimode_xfer) {
297 if(dx || dy) {
298 xfer_mean += dx / (float)win_xsz;
299 xfer_sdev += 0.5 * dy / (float)win_ysz;
301 xfer_mean = xfer_mean < 0.0 ? 0.0 : (xfer_mean > 1.0 ? 1.0 : xfer_mean);
302 xfer_sdev = xfer_sdev < 0.0 ? 0.0 : (xfer_sdev > 1.0 ? 1.0 : xfer_sdev);
304 xfertex_needs_recalc = 1;
305 glutPostRedisplay();
306 }
307 } else {
309 if(bnstate[0]) {
310 cam_theta += dx * 0.5;
311 cam_phi += dy * 0.5;
313 if(cam_phi <= -90) cam_phi = -89;
314 if(cam_phi >= 90) cam_phi = 89;
316 glutPostRedisplay();
317 }
319 if(bnstate[1]) {
320 cam_x += dx * 0.025;
321 cam_y += dy * 0.025;
322 glutPostRedisplay();
323 }
325 if(bnstate[2]) {
326 cam_dist += dy * 0.025;
327 if(cam_dist < 0.0) cam_dist = 0.0;
328 glutPostRedisplay();
329 }
330 }
331 }
334 int parse_args(int argc, char **argv)
335 {
336 int i;
337 struct slice_file *tail;
338 char *endp;
340 for(i=1; i<argc; i++) {
341 if(argv[i][0] == '-' && argv[i][2] == 0) {
342 switch(argv[i][1]) {
343 case 'm':
344 xfer_mean = strtod(argv[++i], &endp);
345 if(endp == argv[i]) {
346 fprintf(stderr, "-m must be followed by the transfer function mean\n");
347 return -1;
348 }
349 break;
351 case 'd':
352 xfer_sdev = strtod(argv[++i], &endp);
353 if(endp == argv[i]) {
354 fprintf(stderr, "-d must be followed by the transfer function std.deviation\n");
355 return -1;
356 }
357 break;
359 default:
360 fprintf(stderr, "unrecognized option: %s\n", argv[i]);
361 return -1;
362 }
363 } else {
364 struct slice_file *sfile;
366 if(!(sfile = malloc(sizeof *sfile))) {
367 perror("failed to allocate memory");
368 return -1;
369 }
370 sfile->name = argv[i];
371 sfile->next = 0;
373 if(!flist) {
374 flist = tail = sfile;
375 } else {
376 tail->next = sfile;
377 tail = sfile;
378 }
379 nslices++;
380 }
381 }
383 if(!nslices) {
384 fprintf(stderr, "pass the slice filenames\n");
385 return -1;
386 }
387 return 0;
388 }
391 static void create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale)
392 {
393 int i, j;
394 int cur_tex_xsz, cur_tex_ysz;
395 int tex_xsz = round_pow2(xsz);
396 int tex_ysz = round_pow2(ysz);
397 float *teximg, *dir;
399 if(!(teximg = malloc(3 * xsz * ysz * sizeof *teximg))) {
400 return;
401 }
402 dir = teximg;
404 for(i=0; i<ysz; i++) {
405 for(j=0; j<xsz; j++) {
406 vec3_t rdir = get_primary_ray_dir(j, i, xsz, ysz, vfov);
407 *dir++ = rdir.x;
408 *dir++ = rdir.y;
409 *dir++ = rdir.z;
410 }
411 }
413 if(!ray_tex) {
414 glGenTextures(1, &ray_tex);
415 }
417 glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_WIDTH, &cur_tex_xsz);
418 glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_HEIGHT, &cur_tex_ysz);
420 if(tex_xsz > cur_tex_xsz || tex_ysz > cur_tex_ysz) {
421 glBindTexture(GL_TEXTURE_2D, ray_tex);
422 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
423 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
424 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
425 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
426 glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F_ARB, tex_xsz, tex_ysz, 0, GL_RGB, GL_FLOAT, 0);
427 }
429 glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, xsz, ysz, GL_RGB, GL_FLOAT, teximg);
430 free(teximg);
432 if(tex_scale) {
433 tex_scale->x = (float)xsz / (float)tex_xsz;
434 tex_scale->y = (float)ysz / (float)tex_ysz;
435 }
436 raytex_needs_recalc = 0;
437 }
439 static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg)
440 {
441 float vfov = M_PI * vfov_deg / 180.0;
442 float aspect = (float)w / (float)h;
444 float ysz = 2.0;
445 float xsz = aspect * ysz;
447 float px = ((float)x / (float)w) * xsz - xsz / 2.0;
448 float py = 1.0 - ((float)y / (float)h) * ysz;
449 float pz = 1.0 / tan(0.5 * vfov);
451 float mag = sqrt(px * px + py * py + pz * pz);
453 return v3_cons(px / mag, py / mag, pz / mag);
454 }
456 static int round_pow2(int x)
457 {
458 x--;
459 x = (x >> 1) | x;
460 x = (x >> 2) | x;
461 x = (x >> 4) | x;
462 x = (x >> 8) | x;
463 x = (x >> 16) | x;
464 return x + 1;
465 }
467 static void create_transfer_map(float mean, float sdev)
468 {
469 static float map[XFER_MAP_SZ];
470 int i;
472 if(!xfer_tex) {
473 glGenTextures(1, &xfer_tex);
474 glBindTexture(GL_TEXTURE_1D, xfer_tex);
475 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
476 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
477 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
478 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
479 glTexImage1D(GL_TEXTURE_1D, 0, GL_LUMINANCE32F_ARB, XFER_MAP_SZ, 0, GL_LUMINANCE, GL_FLOAT, 0);
480 }
482 for(i=0; i<XFER_MAP_SZ; i++) {
483 float x = (float)i / (float)(XFER_MAP_SZ - 1);
484 map[i] = gaussian(x, mean, sdev) - 1.0;
485 }
486 putchar('\n');
488 glTexSubImage1D(GL_TEXTURE_1D, 0, 0, XFER_MAP_SZ, GL_LUMINANCE, GL_FLOAT, map);
489 xfertex_needs_recalc = 0;
490 }