qvolray: 0c3874aa717a src/volray.c

qvolray

view src/volray.c @ 5:0c3874aa717a

slice

author	John Tsiombikas <nuclear@member.fsf.org>
date	Wed, 04 Apr 2012 01:37:33 +0300
parents	3e53a16d4667
children	f40e4edfee7e

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 enum {

24 UIMODE_DEFAULT,

25 UIMODE_XFER,

26 UIMODE_CURSOR

27 };

29 int init(void);

30 void disp(void);

31 void render_volume(void);

32 void draw_slice(void);

33 void draw_xfer_func(void);

34 void reshape(int x, int y);

35 void keyb(unsigned char key, int x, int y);

36 void keyb_up(unsigned char key, int x, int y);

37 void mouse(int bn, int state, int x, int y);

38 void motion(int x, int y);

39 int parse_args(int argc, char **argv);

41 static void create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale);

42 static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg);

43 static int round_pow2(int x);

44 static void create_transfer_map(float mean, float sdev);

46 float cam_theta = 0, cam_phi = 0, cam_dist = 4.0;

47 float cam_x, cam_y, cam_z;

49 vec2_t tex_scale;

50 struct slice_file *flist;

51 int nslices;

52 unsigned int vol_sdr, slice_sdr, vol_tex, ray_tex;

53 int win_xsz, win_ysz;

54 int raytex_needs_recalc = 1;

56 unsigned int xfer_tex;

57 float xfer_mean = 0.5, xfer_sdev = 1.0;

58 int xfertex_needs_recalc = 1;

60 static int uimode;

61 static float cur_z = 0.5;

63 int main(int argc, char **argv)

64 {

65 glutInit(&argc, argv);

66 glutInitWindowSize(1280, 720);

67 glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);

69 if(parse_args(argc, argv) == -1) {

70 return 1;

71 }

73 glutCreateWindow("Volume Raytracer");

75 glutDisplayFunc(disp);

76 glutReshapeFunc(reshape);

77 glutKeyboardFunc(keyb);

78 glutKeyboardUpFunc(keyb_up);

79 glutMouseFunc(mouse);

80 glutMotionFunc(motion);

82 glewInit();

84 if(init() == -1) {

85 return 1;

86 }

88 glutMainLoop();

89 return 0;

90 }

92 int init(void)

93 {

94 int i, vol_xsz, vol_ysz;

96 if(!(vol_sdr = create_program_load("volray.v.glsl", "volray.p.glsl"))) {

97 return -1;

98 }

99 set_uniform_int(vol_sdr, "volume", 0);

100 set_uniform_int(vol_sdr, "ray_tex", 1);

101 set_uniform_int(vol_sdr, "xfer_tex", 2);

102

103 if(!(slice_sdr = create_program_load(0, "slice.p.glsl"))) {

104 return -1;

105 }

106 set_uniform_int(slice_sdr, "volume", 0);

107 set_uniform_int(slice_sdr, "xfer_tex", 1);

108

109 glGenTextures(1, &vol_tex);

110 glBindTexture(GL_TEXTURE_3D, vol_tex);

111 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

112 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

113 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

114 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

115 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);

116

117 for(i=0; i<nslices; i++) {

118 int xsz, ysz;

119 void *pix;

120 struct slice_file *sfile = flist;

121 flist = flist->next;

122

123 if(!(pix = img_load_pixels(sfile->name, &xsz, &ysz, IMG_FMT_RGBA32))) {

124 fprintf(stderr, "failed to load image: %s\n", sfile->name);

125 return -1;

126 }

127

128 if(i == 0) {

129 /* allocate storage for the texture */

130 glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA, xsz, ysz, nslices, 0, GL_RGBA, GL_UNSIGNED_BYTE, 0);

131

132 vol_xsz = xsz;

133 vol_ysz = ysz;

134

135 } else {

136 if(xsz != vol_xsz || ysz != vol_ysz) {

137 fprintf(stderr, "%s: inconsistent slice size: %dx%d. expected: %dx%d\n",

138 sfile->name, xsz, ysz, vol_xsz, vol_ysz);

139 img_free_pixels(pix);

140 return -1;

141 }

142 }

143 free(sfile);

144

145 glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, i, xsz, ysz, 1, GL_RGBA, GL_UNSIGNED_BYTE, pix);

146 img_free_pixels(pix);

147 }

148 return 0;

149 }

150

151 void disp(void)

152 {

153 /* recalculate primary ray texture if needed */

154 if(raytex_needs_recalc) {

155 create_ray_texture(win_xsz, win_ysz, 50.0, &tex_scale);

156 }

157 /* recalculate transfer function texture if needed */

158 if(xfertex_needs_recalc) {

159 create_transfer_map(xfer_mean, xfer_sdev);

160 }

161

162 render_volume();

163 draw_slice();

164 draw_xfer_func();

165

166 glutSwapBuffers();

167 assert(glGetError() == GL_NO_ERROR);

168 }

169

170 void render_volume(void)

171 {

172 /* set the camera transformation */

173 glMatrixMode(GL_MODELVIEW);

174 glPushMatrix();

175 glLoadIdentity();

176 glRotatef(-90, 1, 0, 0);

177 glTranslatef(cam_x, cam_y, -cam_z);

178 glRotatef(cam_theta, 0, 1, 0);

179 glRotatef(cam_phi, 1, 0, 0);

180 glTranslatef(0, 0, -cam_dist);

181

182 /* setup the texture matrix to map the useful part of the ray texture to [0,1] */

183 glMatrixMode(GL_TEXTURE);

184 glPushMatrix();

185 glLoadIdentity();

186 glScalef(tex_scale.x, tex_scale.y, 1.0);

187

188 /* tex unit0: volume data 3D texture */

189 glActiveTexture(GL_TEXTURE0);

190 glBindTexture(GL_TEXTURE_3D, vol_tex);

191 glEnable(GL_TEXTURE_3D);

192

193 /* tex unit1: primary rays in view space */

194 glActiveTexture(GL_TEXTURE1);

195 glBindTexture(GL_TEXTURE_2D, ray_tex);

196 glEnable(GL_TEXTURE_2D);

197

198 /* tex unit2: transfer function (1d) */

199 glActiveTexture(GL_TEXTURE2);

200 glBindTexture(GL_TEXTURE_1D, xfer_tex);

201 glEnable(GL_TEXTURE_1D);

202

203 bind_program(vol_sdr);

204 glBegin(GL_QUADS);

205 glColor3f(1, 1, 1);

206 glTexCoord2f(0, 1); glVertex2f(-1, -1);

207 glTexCoord2f(1, 1); glVertex2f(1, -1);

208 glTexCoord2f(1, 0); glVertex2f(1, 1);

209 glTexCoord2f(0, 0); glVertex2f(-1, 1);

210 glEnd();

211 bind_program(0);

212

213 glActiveTexture(GL_TEXTURE2);

214 glDisable(GL_TEXTURE_1D);

215 glActiveTexture(GL_TEXTURE1);

216 glDisable(GL_TEXTURE_2D);

217 glActiveTexture(GL_TEXTURE0);

218 glDisable(GL_TEXTURE_3D);

219

220 glMatrixMode(GL_TEXTURE);

221 glPopMatrix();

222 glMatrixMode(GL_MODELVIEW);

223 glPopMatrix();

224 }

225

226 void draw_slice(void)

227 {

228 glMatrixMode(GL_MODELVIEW);

229 glPushMatrix();

230 glTranslatef(0.9, 0.9, 0);

231 glScalef(0.3, 0.3 * ((float)win_xsz / win_ysz), 1);

232 glTranslatef(-1, -1, 0);

233

234 glActiveTexture(GL_TEXTURE0);

235 glBindTexture(GL_TEXTURE_3D, vol_tex);

236 glEnable(GL_TEXTURE_3D);

237

238 glActiveTexture(GL_TEXTURE1);

239 glBindTexture(GL_TEXTURE_1D, xfer_tex);

240 glEnable(GL_TEXTURE_1D);

241

242 bind_program(slice_sdr);

243

244 glBegin(GL_QUADS);

245 glColor3f(1, 1, 1);

246 glTexCoord3f(0, 1, cur_z); glVertex2f(-1, -1);

247 glTexCoord3f(1, 1, cur_z); glVertex2f(1, -1);

248 glTexCoord3f(1, 0, cur_z); glVertex2f(1, 1);

249 glTexCoord3f(0, 0, cur_z); glVertex2f(-1, 1);

250 glEnd();

251

252 bind_program(0);

253

254 glActiveTexture(GL_TEXTURE1);

255 glDisable(GL_TEXTURE_1D);

256 glActiveTexture(GL_TEXTURE0);

257 glDisable(GL_TEXTURE_3D);

258 glPopMatrix();

259 }

260

261 void draw_xfer_func(void)

262 {

263 glMatrixMode(GL_MODELVIEW);

264 glPushMatrix();

265 glTranslatef(-0.9, -0.9, 0);

266 glScalef(0.5, 0.1, 1);

267

268 glBindTexture(GL_TEXTURE_1D, xfer_tex);

269 glEnable(GL_TEXTURE_1D);

270

271 glBegin(GL_QUADS);

272 glColor3f(1, 1, 1);

273 glTexCoord1f(1);

274 glVertex2f(1, 0);

275 glVertex2f(1, 1);

276 glTexCoord1f(0);

277 glVertex2f(0, 1);

278 glVertex2f(0, 0);

279 glEnd();

280

281 glDisable(GL_TEXTURE_1D);

282

283 glLineWidth(2.0);

284 glBegin(GL_LINE_LOOP);

285 if(uimode == UIMODE_XFER) {

286 glColor3f(1, 0, 0);

287 } else {

288 glColor3f(0, 0, 1);

289 }

290 glVertex2f(0, 0);

291 glVertex2f(1, 0);

292 glVertex2f(1, 1);

293 glVertex2f(0, 1);

294 glEnd();

295

296 glPopMatrix();

297 }

298

299 void reshape(int x, int y)

300 {

301 glViewport(0, 0, x, y);

302

303 if(x != win_xsz || y != win_ysz) {

304 raytex_needs_recalc = 1;

305 win_xsz = x;

306 win_ysz = y;

307 }

308 }

309

310 void keyb(unsigned char key, int x, int y)

311 {

312 switch(key) {

313 case 27:

314 exit(0);

315

316 case 'x':

317 uimode = UIMODE_XFER;

318 glutPostRedisplay();

319 break;

320

321 case 'c':

322 uimode = UIMODE_CURSOR;

323 glutPostRedisplay();

324 break;

325

326 default:

327 break;

328 }

329 }

330

331 void keyb_up(unsigned char key, int x, int y)

332 {

333 switch(key) {

334 case 'x':

335 if(uimode == UIMODE_XFER) {

336 uimode = UIMODE_DEFAULT;

337 glutPostRedisplay();

338 }

339 break;

340

341 case 'c':

342 if(uimode == UIMODE_CURSOR) {

343 uimode = UIMODE_DEFAULT;

344 glutPostRedisplay();

345 }

346 break;

347

348 default:

349 break;

350 }

351 }

352

353 static int bnstate[32];

354 static int prev_x, prev_y;

355

356 void mouse(int bn, int state, int x, int y)

357 {

358 bnstate[bn - GLUT_LEFT_BUTTON] = state == GLUT_DOWN;

359 prev_x = x;

360 prev_y = y;

361 }

362

363 void motion(int x, int y)

364 {

365 int dx = x - prev_x;

366 int dy = y - prev_y;

367 prev_x = x;

368 prev_y = y;

369

370 switch(uimode) {

371 case UIMODE_XFER:

372 if(dx || dy) {

373 xfer_mean += dx / (float)win_xsz;

374 xfer_sdev += 0.5 * dy / (float)win_ysz;

375

376 xfer_mean = xfer_mean < 0.0 ? 0.0 : (xfer_mean > 1.0 ? 1.0 : xfer_mean);

377 xfer_sdev = xfer_sdev < 0.0 ? 0.0 : (xfer_sdev > 1.0 ? 1.0 : xfer_sdev);

378

379 xfertex_needs_recalc = 1;

380 glutPostRedisplay();

381 }

382 break;

383

384 case UIMODE_CURSOR:

385 cur_z += 0.5 * dy / (float)win_ysz;

386

387 if(cur_z < 0.0)

388 cur_z = 0.0;

389 if(cur_z > 1.0)

390 cur_z = 1.0;

391 glutPostRedisplay();

392 break;

393

394 default:

395 /* view control */

396 if(bnstate[0]) {

397 cam_theta += dx * 0.5;

398 cam_phi += dy * 0.5;

399

400 if(cam_phi <= -90) cam_phi = -89;

401 if(cam_phi >= 90) cam_phi = 89;

402

403 glutPostRedisplay();

404 }

405

406 if(bnstate[1]) {

407 cam_x += dx * 0.025;

408 cam_y += dy * 0.025;

409 glutPostRedisplay();

410 }

411

412 if(bnstate[2]) {

413 cam_dist += dy * 0.025;

414 if(cam_dist < 0.0) cam_dist = 0.0;

415 glutPostRedisplay();

416 }

417 }

418 }

419

420

421 int parse_args(int argc, char **argv)

422 {

423 int i;

424 struct slice_file *tail;

425 char *endp;

426

427 for(i=1; i<argc; i++) {

428 if(argv[i][0] == '-' && argv[i][2] == 0) {

429 switch(argv[i][1]) {

430 case 'm':

431 xfer_mean = strtod(argv[++i], &endp);

432 if(endp == argv[i]) {

433 fprintf(stderr, "-m must be followed by the transfer function mean\n");

434 return -1;

435 }

436 break;

437

438 case 'd':

439 xfer_sdev = strtod(argv[++i], &endp);

440 if(endp == argv[i]) {

441 fprintf(stderr, "-d must be followed by the transfer function std.deviation\n");

442 return -1;

443 }

444 break;

445

446 default:

447 fprintf(stderr, "unrecognized option: %s\n", argv[i]);

448 return -1;

449 }

450 } else {

451 struct slice_file *sfile;

452

453 if(!(sfile = malloc(sizeof *sfile))) {

454 perror("failed to allocate memory");

455 return -1;

456 }

457 sfile->name = argv[i];

458 sfile->next = 0;

459

460 if(!flist) {

461 flist = tail = sfile;

462 } else {

463 tail->next = sfile;

464 tail = sfile;

465 }

466 nslices++;

467 }

468 }

469

470 if(!nslices) {

471 fprintf(stderr, "pass the slice filenames\n");

472 return -1;

473 }

474 return 0;

475 }

476

477

478 static void create_ray_texture(int xsz, int ysz, float vfov, vec2_t *tex_scale)

479 {

480 int i, j;

481 int cur_tex_xsz, cur_tex_ysz;

482 int tex_xsz = round_pow2(xsz);

483 int tex_ysz = round_pow2(ysz);

484 float *teximg, *dir;

485

486 if(!(teximg = malloc(3 * xsz * ysz * sizeof *teximg))) {

487 return;

488 }

489 dir = teximg;

490

491 for(i=0; i<ysz; i++) {

492 for(j=0; j<xsz; j++) {

493 vec3_t rdir = get_primary_ray_dir(j, i, xsz, ysz, vfov);

494 *dir++ = rdir.x;

495 *dir++ = rdir.y;

496 *dir++ = rdir.z;

497 }

498 }

499

500 if(!ray_tex) {

501 glGenTextures(1, &ray_tex);

502 }

503

504 glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_WIDTH, &cur_tex_xsz);

505 glGetTexLevelParameteriv(GL_TEXTURE_2D, 0, GL_TEXTURE_HEIGHT, &cur_tex_ysz);

506

507 if(tex_xsz > cur_tex_xsz || tex_ysz > cur_tex_ysz) {

508 glBindTexture(GL_TEXTURE_2D, ray_tex);

509 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);

510 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);

511 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);

512 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);

513 glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F_ARB, tex_xsz, tex_ysz, 0, GL_RGB, GL_FLOAT, 0);

514 }

515

516 glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, xsz, ysz, GL_RGB, GL_FLOAT, teximg);

517 free(teximg);

518

519 if(tex_scale) {

520 tex_scale->x = (float)xsz / (float)tex_xsz;

521 tex_scale->y = (float)ysz / (float)tex_ysz;

522 }

523 raytex_needs_recalc = 0;

524 }

525

526 static vec3_t get_primary_ray_dir(int x, int y, int w, int h, float vfov_deg)

527 {

528 float vfov = M_PI * vfov_deg / 180.0;

529 float aspect = (float)w / (float)h;

530

531 float ysz = 2.0;

532 float xsz = aspect * ysz;

533

534 float px = ((float)x / (float)w) * xsz - xsz / 2.0;

535 float py = 1.0 - ((float)y / (float)h) * ysz;

536 float pz = 1.0 / tan(0.5 * vfov);

537

538 float mag = sqrt(px * px + py * py + pz * pz);

539

540 return v3_cons(px / mag, py / mag, pz / mag);

541 }

542

543 static int round_pow2(int x)

544 {

545 x--;

546 x = (x >> 1) | x;

547 x = (x >> 2) | x;

548 x = (x >> 4) | x;

549 x = (x >> 8) | x;

550 x = (x >> 16) | x;

551 return x + 1;

552 }

553

554 static void create_transfer_map(float mean, float sdev)

555 {

556 static float map[XFER_MAP_SZ];

557 int i;

558

559 if(!xfer_tex) {

560 glGenTextures(1, &xfer_tex);

561 glBindTexture(GL_TEXTURE_1D, xfer_tex);

562 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

563 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

564 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

565 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

566 glTexImage1D(GL_TEXTURE_1D, 0, GL_LUMINANCE32F_ARB, XFER_MAP_SZ, 0, GL_LUMINANCE, GL_FLOAT, 0);

567 }

568

569 for(i=0; i<XFER_MAP_SZ; i++) {

570 float x = (float)i / (float)(XFER_MAP_SZ - 1);

571 map[i] = gaussian(x, mean, sdev) - 1.0;

572 }

573 putchar('\n');

574

575 glTexSubImage1D(GL_TEXTURE_1D, 0, 0, XFER_MAP_SZ, GL_LUMINANCE, GL_FLOAT, map);

576 xfertex_needs_recalc = 0;

577 }