qvolray

annotate src/volray.c @ 3:6f275934717b

find changesets by author, revision, files, or words in the commit message
foon
author John Tsiombikas <nuclear@member.fsf.org>
date Mon, 02 Apr 2012 14:42:03 +0300
parents 0b73aa7317e1
children 3e53a16d4667
rev   line source
nuclear@0 1 #include <stdio.h>
nuclear@0 2 #include <stdlib.h>
nuclear@0 3 #include <assert.h>
nuclear@0 4
nuclear@0 5 #include <GL/glew.h>
nuclear@0 6 #ifndef __APPLE__
nuclear@0 7 #include <GL/glut.h>
nuclear@0 8 #else
nuclear@0 9 #include <GLUT/glut.h>
nuclear@0 10 #endif
nuclear@0 11
nuclear@0 12 #include <vmath/vmath.h>
nuclear@0 13 #include <imago2.h>
nuclear@0 14 #include "sdr.h"
nuclear@0 15
nuclear@3 16 #define XFER_MAP_SZ 512
nuclear@3 17
nuclear@0 18 struct slice_file {
nuclear@0 19 char *name;
nuclear@0 20 struct slice_file *next;
nuclear@0 21 };
nuclear@0 22
nuclear@0 23 int init(void);
nuclear@0 24 void disp(void);
nuclear@0 25 void reshape(int x, int y);
nuclear@0 26 void keyb(unsigned char key, int x, int y);
nuclear@0 27 void mouse(int bn, int state, int x, int y);
nuclear@0 28 void motion(int x, int y);
nuclear@0 29 int parse_args(int argc, char **argv);
nuclear@0 30
nuclear@1 31 static void create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale);
nuclear@0 32 static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg);
nuclear@0 33 static int round_pow2(int x);
nuclear@3 34 static void create_transfer_map(float mean, float sdev);
nuclear@0 35
nuclear@0 36 float cam_theta = 0, cam_phi = 0, cam_dist = 4.0;
nuclear@0 37 float cam_x, cam_y, cam_z;
nuclear@0 38
nuclear@0 39 vec2_t tex_scale;
nuclear@0 40 struct slice_file *flist;
nuclear@0 41 int nslices;
nuclear@0 42 unsigned int sdr, vol_tex, ray_tex;
nuclear@0 43 int win_xsz, win_ysz;
nuclear@1 44 int raytex_needs_recalc = 1;
nuclear@0 45
nuclear@3 46 unsigned int xfer_tex;
nuclear@3 47 float xfer_mean = 0.5, xfer_sdev = 1.0;
nuclear@3 48 int xfertex_needs_recalc = 1;
nuclear@3 49
nuclear@0 50 int main(int argc, char **argv)
nuclear@0 51 {
nuclear@0 52 glutInit(&argc, argv);
nuclear@0 53 glutInitWindowSize(1280, 720);
nuclear@0 54 glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
nuclear@0 55
nuclear@0 56 if(parse_args(argc, argv) == -1) {
nuclear@0 57 return 1;
nuclear@0 58 }
nuclear@0 59
nuclear@0 60 glutCreateWindow("Volume Raytracer");
nuclear@0 61
nuclear@0 62 glutDisplayFunc(disp);
nuclear@0 63 glutReshapeFunc(reshape);
nuclear@0 64 glutKeyboardFunc(keyb);
nuclear@0 65 glutMouseFunc(mouse);
nuclear@0 66 glutMotionFunc(motion);
nuclear@0 67
nuclear@0 68 glewInit();
nuclear@0 69
nuclear@0 70 if(init() == -1) {
nuclear@0 71 return 1;
nuclear@0 72 }
nuclear@0 73
nuclear@0 74 glutMainLoop();
nuclear@0 75 return 0;
nuclear@0 76 }
nuclear@0 77
nuclear@0 78 int init(void)
nuclear@0 79 {
nuclear@0 80 int i, vol_xsz, vol_ysz;
nuclear@0 81
nuclear@0 82 if(!(sdr = create_program_load("volray.v.glsl", "volray.p.glsl"))) {
nuclear@1 83 return -1;
nuclear@0 84 }
nuclear@1 85 set_uniform_int(sdr, "volume", 0);
nuclear@1 86 set_uniform_int(sdr, "ray_tex", 1);
nuclear@3 87 set_uniform_int(sdr, "xfer_tex", 2);
nuclear@0 88
nuclear@0 89 glGenTextures(1, &vol_tex);
nuclear@0 90 glBindTexture(GL_TEXTURE_3D, vol_tex);
nuclear@0 91 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
nuclear@0 92 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
nuclear@1 93 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
nuclear@1 94 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
nuclear@1 95 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
nuclear@0 96
nuclear@0 97 for(i=0; i<nslices; i++) {
nuclear@0 98 int xsz, ysz;
nuclear@0 99 void *pix;
nuclear@0 100 struct slice_file *sfile = flist;
nuclear@0 101 flist = flist->next;
nuclear@0 102
nuclear@0 103 if(!(pix = img_load_pixels(sfile->name, &xsz, &ysz, IMG_FMT_RGBA32))) {
nuclear@0 104 fprintf(stderr, "failed to load image: %s\n", sfile->name);
nuclear@0 105 return -1;
nuclear@0 106 }
nuclear@0 107
nuclear@0 108 if(i == 0) {
nuclear@0 109 /* allocate storage for the texture */
nuclear@0 110 glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA, xsz, ysz, nslices, 0, GL_RGBA, GL_UNSIGNED_BYTE, 0);
nuclear@0 111
nuclear@0 112 vol_xsz = xsz;
nuclear@0 113 vol_ysz = ysz;
nuclear@0 114
nuclear@0 115 } else {
nuclear@0 116 if(xsz != vol_xsz || ysz != vol_ysz) {
nuclear@0 117 fprintf(stderr, "%s: inconsistent slice size: %dx%d. expected: %dx%d\n",
nuclear@0 118 sfile->name, xsz, ysz, vol_xsz, vol_ysz);
nuclear@0 119 img_free_pixels(pix);
nuclear@0 120 return -1;
nuclear@0 121 }
nuclear@0 122 }
nuclear@0 123 free(sfile);
nuclear@0 124
nuclear@0 125 glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, i, xsz, ysz, 1, GL_RGBA, GL_UNSIGNED_BYTE, pix);
nuclear@0 126 img_free_pixels(pix);
nuclear@0 127 }
nuclear@0 128 return 0;
nuclear@0 129 }
nuclear@0 130
nuclear@0 131 void disp(void)
nuclear@0 132 {
nuclear@1 133 if(raytex_needs_recalc) {
nuclear@1 134 create_ray_texture(win_xsz, win_ysz, 50.0, &tex_scale);
nuclear@1 135 }
nuclear@3 136 if(xfertex_needs_recalc) {
nuclear@3 137 create_transfer_map(xfer_mean, xfer_sdev);
nuclear@3 138 }
nuclear@1 139
nuclear@0 140 glMatrixMode(GL_MODELVIEW);
nuclear@0 141 glLoadIdentity();
nuclear@1 142 glRotatef(-90, 1, 0, 0);
nuclear@0 143 glTranslatef(cam_x, cam_y, -cam_z);
nuclear@0 144 glRotatef(cam_theta, 0, 1, 0);
nuclear@0 145 glRotatef(cam_phi, 1, 0, 0);
nuclear@0 146 glTranslatef(0, 0, -cam_dist);
nuclear@0 147
nuclear@0 148 glMatrixMode(GL_TEXTURE);
nuclear@1 149 glLoadIdentity();
nuclear@0 150 glScalef(tex_scale.x, tex_scale.y, 1.0);
nuclear@0 151
nuclear@0 152 glActiveTexture(GL_TEXTURE0);
nuclear@1 153 glBindTexture(GL_TEXTURE_3D, vol_tex);
nuclear@0 154 glEnable(GL_TEXTURE_3D);
nuclear@0 155
nuclear@0 156 glActiveTexture(GL_TEXTURE1);
nuclear@1 157 glBindTexture(GL_TEXTURE_2D, ray_tex);
nuclear@0 158 glEnable(GL_TEXTURE_2D);
nuclear@0 159
nuclear@3 160 glActiveTexture(GL_TEXTURE2);
nuclear@3 161 glBindTexture(GL_TEXTURE_1D, xfer_tex);
nuclear@3 162 glEnable(GL_TEXTURE_1D);
nuclear@3 163
nuclear@0 164 bind_program(sdr);
nuclear@0 165 glBegin(GL_QUADS);
nuclear@0 166 glColor3f(1, 1, 1);
nuclear@1 167 glTexCoord2f(0, 1); glVertex2f(-1, -1);
nuclear@1 168 glTexCoord2f(1, 1); glVertex2f(1, -1);
nuclear@1 169 glTexCoord2f(1, 0); glVertex2f(1, 1);
nuclear@1 170 glTexCoord2f(0, 0); glVertex2f(-1, 1);
nuclear@0 171 glEnd();
nuclear@0 172 bind_program(0);
nuclear@0 173
nuclear@3 174 glActiveTexture(GL_TEXTURE2);
nuclear@3 175 glDisable(GL_TEXTURE_1D);
nuclear@1 176 glActiveTexture(GL_TEXTURE1);
nuclear@0 177 glDisable(GL_TEXTURE_2D);
nuclear@0 178 glActiveTexture(GL_TEXTURE0);
nuclear@0 179 glDisable(GL_TEXTURE_3D);
nuclear@0 180
nuclear@0 181 glMatrixMode(GL_TEXTURE);
nuclear@1 182 glLoadIdentity();
nuclear@0 183
nuclear@0 184 glutSwapBuffers();
nuclear@0 185 assert(glGetError() == GL_NO_ERROR);
nuclear@0 186 }
nuclear@0 187
nuclear@0 188 void reshape(int x, int y)
nuclear@0 189 {
nuclear@1 190 printf("reshape: %dx%d\n", x, y);
nuclear@0 191 glViewport(0, 0, x, y);
nuclear@0 192
nuclear@0 193 if(x != win_xsz || y != win_ysz) {
nuclear@1 194 raytex_needs_recalc = 1;
nuclear@0 195 win_xsz = x;
nuclear@0 196 win_ysz = y;
nuclear@0 197 }
nuclear@0 198 }
nuclear@0 199
nuclear@0 200 void keyb(unsigned char key, int x, int y)
nuclear@0 201 {
nuclear@0 202 switch(key) {
nuclear@0 203 case 27:
nuclear@0 204 exit(0);
nuclear@0 205 }
nuclear@0 206 }
nuclear@0 207
nuclear@0 208 static int bnstate[32];
nuclear@0 209 static int prev_x, prev_y;
nuclear@0 210
nuclear@0 211 void mouse(int bn, int state, int x, int y)
nuclear@0 212 {
nuclear@0 213 bnstate[bn - GLUT_LEFT_BUTTON] = state == GLUT_DOWN;
nuclear@0 214 prev_x = x;
nuclear@0 215 prev_y = y;
nuclear@0 216 }
nuclear@0 217
nuclear@0 218 void motion(int x, int y)
nuclear@0 219 {
nuclear@0 220 int dx = x - prev_x;
nuclear@0 221 int dy = y - prev_y;
nuclear@0 222
nuclear@0 223 prev_x = x;
nuclear@0 224 prev_y = y;
nuclear@0 225
nuclear@0 226 if(bnstate[0]) {
nuclear@0 227 cam_theta += dx * 0.5;
nuclear@0 228 cam_phi += dy * 0.5;
nuclear@0 229
nuclear@1 230 if(cam_phi <= -90) cam_phi = -89;
nuclear@1 231 if(cam_phi >= 90) cam_phi = 89;
nuclear@0 232
nuclear@0 233 glutPostRedisplay();
nuclear@0 234 }
nuclear@0 235
nuclear@0 236 if(bnstate[1]) {
nuclear@1 237 cam_x += dx * 0.025;
nuclear@1 238 cam_y += dy * 0.025;
nuclear@0 239 glutPostRedisplay();
nuclear@0 240 }
nuclear@0 241
nuclear@0 242 if(bnstate[2]) {
nuclear@1 243 cam_dist += dy * 0.025;
nuclear@0 244 if(cam_dist < 0.0) cam_dist = 0.0;
nuclear@0 245 glutPostRedisplay();
nuclear@0 246 }
nuclear@0 247 }
nuclear@0 248
nuclear@0 249
nuclear@0 250 int parse_args(int argc, char **argv)
nuclear@0 251 {
nuclear@0 252 int i;
nuclear@0 253 struct slice_file *tail;
nuclear@3 254 char *endp;
nuclear@0 255
nuclear@0 256 for(i=1; i<argc; i++) {
nuclear@2 257 if(argv[i][0] == '-' && argv[i][2] == 0) {
nuclear@2 258 switch(argv[i][1]) {
nuclear@3 259 case 'm':
nuclear@3 260 xfer_mean = strtod(argv[++i], &endp);
nuclear@3 261 if(endp == argv[i]) {
nuclear@3 262 fprintf(stderr, "-m must be followed by the transfer function mean\n");
nuclear@3 263 return -1;
nuclear@3 264 }
nuclear@2 265 break;
nuclear@3 266
nuclear@3 267 case 'v':
nuclear@3 268 xfer_sdev = strtod(argv[++i], &endp);
nuclear@3 269 if(endp == argv[i]) {
nuclear@3 270 fprintf(stderr, "-v must be followed by the transfer function sdeviance\n");
nuclear@3 271 return -1;
nuclear@3 272 }
nuclear@3 273 break;
nuclear@3 274
nuclear@3 275 default:
nuclear@3 276 fprintf(stderr, "unrecognized option: %s\n", argv[i]);
nuclear@3 277 return -1;
nuclear@2 278 }
nuclear@2 279 } else {
nuclear@2 280 struct slice_file *sfile;
nuclear@0 281
nuclear@2 282 if(!(sfile = malloc(sizeof *sfile))) {
nuclear@2 283 perror("failed to allocate memory");
nuclear@2 284 return -1;
nuclear@2 285 }
nuclear@2 286 sfile->name = argv[i];
nuclear@2 287 sfile->next = 0;
nuclear@2 288
nuclear@2 289 if(!flist) {
nuclear@2 290 flist = tail = sfile;
nuclear@2 291 } else {
nuclear@2 292 tail->next = sfile;
nuclear@2 293 tail = sfile;
nuclear@2 294 }
nuclear@2 295 nslices++;
nuclear@0 296 }
nuclear@0 297 }
nuclear@0 298
nuclear@0 299 if(!nslices) {
nuclear@0 300 fprintf(stderr, "pass the slice filenames\n");
nuclear@0 301 return -1;
nuclear@0 302 }
nuclear@0 303 return 0;
nuclear@0 304 }
nuclear@0 305
nuclear@0 306
nuclear@1 307 static void create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale)
nuclear@0 308 {
nuclear@0 309 int i, j;
nuclear@1 310 int cur_tex_xsz, cur_tex_ysz;
nuclear@0 311 int tex_xsz = round_pow2(xsz);
nuclear@0 312 int tex_ysz = round_pow2(ysz);
nuclear@0 313 float *teximg, *dir;
nuclear@0 314
nuclear@1 315 if(!(teximg = malloc(3 * xsz * ysz * sizeof *teximg))) {
nuclear@1 316 return;
nuclear@0 317 }
nuclear@0 318 dir = teximg;
nuclear@0 319
nuclear@1 320 for(i=0; i<ysz; i++) {
nuclear@1 321 for(j=0; j<xsz; j++) {
nuclear@1 322 vec3_t rdir = get_primary_ray_dir(j, i, xsz, ysz, vfov);
nuclear@1 323 *dir++ = rdir.x;
nuclear@1 324 *dir++ = rdir.y;
nuclear@1 325 *dir++ = rdir.z;
nuclear@0 326 }
nuclear@0 327 }
nuclear@0 328
nuclear@1 329 if(!ray_tex) {
nuclear@1 330 glGenTextures(1, &ray_tex);
nuclear@1 331 }
nuclear@1 332
nuclear@1 333 glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_WIDTH, &cur_tex_xsz);
nuclear@1 334 glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_HEIGHT, &cur_tex_ysz);
nuclear@1 335
nuclear@1 336 if(tex_xsz > cur_tex_xsz || tex_ysz > cur_tex_ysz) {
nuclear@1 337 glBindTexture(GL_TEXTURE_2D, ray_tex);
nuclear@1 338 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
nuclear@1 339 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
nuclear@1 340 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
nuclear@1 341 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
nuclear@1 342 glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F_ARB, tex_xsz, tex_ysz, 0, GL_RGB, GL_FLOAT, 0);
nuclear@1 343 }
nuclear@1 344
nuclear@1 345 glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, xsz, ysz, GL_RGB, GL_FLOAT, teximg);
nuclear@0 346 free(teximg);
nuclear@0 347
nuclear@0 348 if(tex_scale) {
nuclear@0 349 tex_scale->x = (float)xsz / (float)tex_xsz;
nuclear@0 350 tex_scale->y = (float)ysz / (float)tex_ysz;
nuclear@0 351 }
nuclear@1 352 raytex_needs_recalc = 0;
nuclear@0 353 }
nuclear@0 354
nuclear@0 355 static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg)
nuclear@0 356 {
nuclear@0 357 float vfov = M_PI * vfov_deg / 180.0;
nuclear@0 358 float aspect = (float)w / (float)h;
nuclear@0 359
nuclear@0 360 float ysz = 2.0;
nuclear@0 361 float xsz = aspect * ysz;
nuclear@0 362
nuclear@0 363 float px = ((float)x / (float)w) * xsz - xsz / 2.0;
nuclear@0 364 float py = 1.0 - ((float)y / (float)h) * ysz;
nuclear@0 365 float pz = 1.0 / tan(0.5 * vfov);
nuclear@0 366
nuclear@0 367 float mag = sqrt(px * px + py * py + pz * pz);
nuclear@0 368
nuclear@0 369 return v3_cons(px / mag, py / mag, pz / mag);
nuclear@0 370 }
nuclear@0 371
nuclear@0 372 static int round_pow2(int x)
nuclear@0 373 {
nuclear@0 374 x--;
nuclear@0 375 x = (x >> 1) | x;
nuclear@0 376 x = (x >> 2) | x;
nuclear@0 377 x = (x >> 4) | x;
nuclear@0 378 x = (x >> 8) | x;
nuclear@0 379 x = (x >> 16) | x;
nuclear@0 380 return x + 1;
nuclear@0 381 }
nuclear@0 382
nuclear@3 383 static void create_transfer_map(float mean, float sdev)
nuclear@3 384 {
nuclear@3 385 static float map[XFER_MAP_SZ];
nuclear@3 386 int i;
nuclear@3 387
nuclear@3 388 if(!xfer_tex) {
nuclear@3 389 glGenTextures(1, &xfer_tex);
nuclear@3 390 glBindTexture(GL_TEXTURE_1D, xfer_tex);
nuclear@3 391 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
nuclear@3 392 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
nuclear@3 393 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
nuclear@3 394 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
nuclear@3 395 glTexImage1D(GL_TEXTURE_1D, 0, GL_LUMINANCE32F_ARB, XFER_MAP_SZ, 0, GL_LUMINANCE, GL_FLOAT, 0);
nuclear@3 396 }
nuclear@3 397
nuclear@3 398 for(i=0; i<XFER_MAP_SZ; i++) {
nuclear@3 399 float x = (float)i / (float)XFER_MAP_SZ;
nuclear@3 400 map[i] = gaussian(x, mean, sdev);
nuclear@3 401 }
nuclear@3 402
nuclear@3 403 glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, XFER_MAP_SZ, 1, GL_LUMINANCE, GL_FLOAT, map);
nuclear@3 404 xfertex_needs_recalc = 0;
nuclear@3 405 }