qvolray: f40e4edfee7e src/volray.c

qvolray

view src/volray.c @ 7:f40e4edfee7e

shading

author	John Tsiombikas <nuclear@member.fsf.org>
date	Thu, 05 Apr 2012 00:52:11 +0300
parents	0c3874aa717a
children	a6765984e057

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.7, xfer_sdev = 0.1;

58 int xfertex_needs_recalc = 1;

60 static int uimode;

61 static float cur_z = 0.0;

62 static float ray_step = 0.01;

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

65 {

66 glutInit(&argc, argv);

67 glutInitWindowSize(1280, 720);

68 glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);

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

71 return 1;

72 }

74 glutCreateWindow("Volume Raytracer");

76 glutDisplayFunc(disp);

77 glutReshapeFunc(reshape);

78 glutKeyboardFunc(keyb);

79 glutKeyboardUpFunc(keyb_up);

80 glutMouseFunc(mouse);

81 glutMotionFunc(motion);

83 glewInit();

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

86 return 1;

87 }

89 glutMainLoop();

90 return 0;

91 }

93 int init(void)

94 {

95 int i, vol_xsz, vol_ysz;

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

98 return -1;

99 }

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

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

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

103 set_uniform_float(vol_sdr, "ray_step", ray_step);

104 set_uniform_float(vol_sdr, "zclip", cur_z);

105

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

107 return -1;

108 }

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

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

111

112 glGenTextures(1, &vol_tex);

113 glBindTexture(GL_TEXTURE_3D, vol_tex);

114 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

115 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

116 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

117 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

118 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);

119

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

121 int xsz, ysz;

122 void *pix;

123 struct slice_file *sfile = flist;

124 flist = flist->next;

125

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

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

128 return -1;

129 }

130

131 if(i == 0) {

132 /* allocate storage for the texture */

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

134

135 vol_xsz = xsz;

136 vol_ysz = ysz;

137

138 } else {

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

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

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

142 img_free_pixels(pix);

143 return -1;

144 }

145 }

146 free(sfile);

147

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

149 img_free_pixels(pix);

150 }

151 return 0;

152 }

153

154 void disp(void)

155 {

156 /* recalculate primary ray texture if needed */

157 if(raytex_needs_recalc) {

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

159 }

160 /* recalculate transfer function texture if needed */

161 if(xfertex_needs_recalc) {

162 create_transfer_map(xfer_mean, xfer_sdev);

163 }

164

165 render_volume();

166 draw_slice();

167 draw_xfer_func();

168

169 glutSwapBuffers();

170 assert(glGetError() == GL_NO_ERROR);

171 }

172

173 void render_volume(void)

174 {

175 /* set the camera transformation */

176 glMatrixMode(GL_MODELVIEW);

177 glPushMatrix();

178 glLoadIdentity();

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

180 glTranslatef(cam_x, cam_y, -cam_z);

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

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

183 glTranslatef(0, 0, -cam_dist);

184

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

186 glMatrixMode(GL_TEXTURE);

187 glPushMatrix();

188 glLoadIdentity();

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

190

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

192 glActiveTexture(GL_TEXTURE0);

193 glBindTexture(GL_TEXTURE_3D, vol_tex);

194 glEnable(GL_TEXTURE_3D);

195

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

197 glActiveTexture(GL_TEXTURE1);

198 glBindTexture(GL_TEXTURE_2D, ray_tex);

199 glEnable(GL_TEXTURE_2D);

200

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

202 glActiveTexture(GL_TEXTURE2);

203 glBindTexture(GL_TEXTURE_1D, xfer_tex);

204 glEnable(GL_TEXTURE_1D);

205

206 bind_program(vol_sdr);

207 glBegin(GL_QUADS);

208 glColor3f(1, 1, 1);

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

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

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

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

213 glEnd();

214 bind_program(0);

215

216 glActiveTexture(GL_TEXTURE2);

217 glDisable(GL_TEXTURE_1D);

218 glActiveTexture(GL_TEXTURE1);

219 glDisable(GL_TEXTURE_2D);

220 glActiveTexture(GL_TEXTURE0);

221 glDisable(GL_TEXTURE_3D);

222

223 glMatrixMode(GL_TEXTURE);

224 glPopMatrix();

225 glMatrixMode(GL_MODELVIEW);

226 glPopMatrix();

227 }

228

229 void draw_slice(void)

230 {

231 glMatrixMode(GL_MODELVIEW);

232 glPushMatrix();

233 glTranslatef(0.9, 0.9, 0);

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

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

236

237 glActiveTexture(GL_TEXTURE0);

238 glBindTexture(GL_TEXTURE_3D, vol_tex);

239 glEnable(GL_TEXTURE_3D);

240

241 glActiveTexture(GL_TEXTURE1);

242 glBindTexture(GL_TEXTURE_1D, xfer_tex);

243 glEnable(GL_TEXTURE_1D);

244

245 bind_program(slice_sdr);

246

247 glBegin(GL_QUADS);

248 glColor3f(1, 1, 1);

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

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

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

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

253 glEnd();

254

255 bind_program(0);

256

257 glActiveTexture(GL_TEXTURE1);

258 glDisable(GL_TEXTURE_1D);

259 glActiveTexture(GL_TEXTURE0);

260 glDisable(GL_TEXTURE_3D);

261 glPopMatrix();

262 }

263

264 void draw_xfer_func(void)

265 {

266 glMatrixMode(GL_MODELVIEW);

267 glPushMatrix();

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

269 glScalef(0.5, 0.1, 1);

270

271 glBindTexture(GL_TEXTURE_1D, xfer_tex);

272 glEnable(GL_TEXTURE_1D);

273

274 glBegin(GL_QUADS);

275 glColor3f(1, 1, 1);

276 glTexCoord1f(1);

277 glVertex2f(1, 0);

278 glVertex2f(1, 1);

279 glTexCoord1f(0);

280 glVertex2f(0, 1);

281 glVertex2f(0, 0);

282 glEnd();

283

284 glDisable(GL_TEXTURE_1D);

285

286 glLineWidth(2.0);

287 glBegin(GL_LINE_LOOP);

288 if(uimode == UIMODE_XFER) {

289 glColor3f(1, 0, 0);

290 } else {

291 glColor3f(0, 0, 1);

292 }

293 glVertex2f(0, 0);

294 glVertex2f(1, 0);

295 glVertex2f(1, 1);

296 glVertex2f(0, 1);

297 glEnd();

298

299 glPopMatrix();

300 }

301

302 void reshape(int x, int y)

303 {

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

305

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

307 raytex_needs_recalc = 1;

308 win_xsz = x;

309 win_ysz = y;

310 }

311 }

312

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

314 {

315 switch(key) {

316 case 27:

317 exit(0);

318

319 case 'x':

320 uimode = UIMODE_XFER;

321 glutPostRedisplay();

322 break;

323

324 case 'c':

325 uimode = UIMODE_CURSOR;

326 glutPostRedisplay();

327 break;

328

329 default:

330 break;

331 }

332 }

333

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

335 {

336 switch(key) {

337 case 'x':

338 if(uimode == UIMODE_XFER) {

339 uimode = UIMODE_DEFAULT;

340 glutPostRedisplay();

341 }

342 break;

343

344 case 'c':

345 if(uimode == UIMODE_CURSOR) {

346 uimode = UIMODE_DEFAULT;

347 glutPostRedisplay();

348 }

349 break;

350

351 default:

352 break;

353 }

354 }

355

356 static int bnstate[32];

357 static int prev_x, prev_y;

358

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

360 {

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

362 prev_x = x;

363 prev_y = y;

364 }

365

366 void motion(int x, int y)

367 {

368 int dx = x - prev_x;

369 int dy = y - prev_y;

370 prev_x = x;

371 prev_y = y;

372

373 switch(uimode) {

374 case UIMODE_XFER:

375 if(dx || dy) {

376 xfer_mean += dx / (float)win_xsz;

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

378

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

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

381

382 xfertex_needs_recalc = 1;

383 glutPostRedisplay();

384 }

385 break;

386

387 case UIMODE_CURSOR:

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

389

390 if(cur_z < 0.0)

391 cur_z = 0.0;

392 if(cur_z > 1.0)

393 cur_z = 1.0;

394

395 set_uniform_float(vol_sdr, "zclip", cur_z);

396 glutPostRedisplay();

397 break;

398

399 default:

400 /* view control */

401 if(bnstate[0]) {

402 cam_theta += dx * 0.5;

403 cam_phi += dy * 0.5;

404

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

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

407

408 glutPostRedisplay();

409 }

410

411 if(bnstate[1]) {

412 cam_x += dx * 0.025;

413 cam_y += dy * 0.025;

414 glutPostRedisplay();

415 }

416

417 if(bnstate[2]) {

418 cam_dist += dy * 0.025;

419 if(cam_dist < 0.0) cam_dist = 0.0;

420 glutPostRedisplay();

421 }

422 }

423 }

424

425

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

427 {

428 int i;

429 struct slice_file *tail;

430 char *endp;

431

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

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

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

435 case 'm':

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

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

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

439 return -1;

440 }

441 break;

442

443 case 'd':

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

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

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

447 return -1;

448 }

449 break;

450

451 default:

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

453 return -1;

454 }

455 } else {

456 struct slice_file *sfile;

457

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

459 perror("failed to allocate memory");

460 return -1;

461 }

462 sfile->name = argv[i];

463 sfile->next = 0;

464

465 if(!flist) {

466 flist = tail = sfile;

467 } else {

468 tail->next = sfile;

469 tail = sfile;

470 }

471 nslices++;

472 }

473 }

474

475 if(!nslices) {

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

477 return -1;

478 }

479 return 0;

480 }

481

482

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

484 {

485 int i, j;

486 int cur_tex_xsz, cur_tex_ysz;

487 int tex_xsz = round_pow2(xsz);

488 int tex_ysz = round_pow2(ysz);

489 float *teximg, *dir;

490

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

492 return;

493 }

494 dir = teximg;

495

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

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

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

499 *dir++ = rdir.x;

500 *dir++ = rdir.y;

501 *dir++ = rdir.z;

502 }

503 }

504

505 if(!ray_tex) {

506 glGenTextures(1, &ray_tex);

507 }

508

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

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

511

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

513 glBindTexture(GL_TEXTURE_2D, ray_tex);

514 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);

515 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);

516 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);

517 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);

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

519 }

520

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

522 free(teximg);

523

524 if(tex_scale) {

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

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

527 }

528 raytex_needs_recalc = 0;

529 }

530

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

532 {

533 float vfov = M_PI * vfov_deg / 180.0;

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

535

536 float ysz = 2.0;

537 float xsz = aspect * ysz;

538

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

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

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

542

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

544

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

546 }

547

548 static int round_pow2(int x)

549 {

550 x--;

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

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

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

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

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

556 return x + 1;

557 }

558

559 static void create_transfer_map(float mean, float sdev)

560 {

561 static float map[XFER_MAP_SZ];

562 int i;

563

564 if(!xfer_tex) {

565 glGenTextures(1, &xfer_tex);

566 glBindTexture(GL_TEXTURE_1D, xfer_tex);

567 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

568 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

569 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

570 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

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

572 }

573

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

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

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

577 }

578 putchar('\n');

579

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

581 xfertex_needs_recalc = 0;

582 }