qvolray: 6f275934717b src/volray.c

qvolray

view src/volray.c @ 3:6f275934717b

foon

author	John Tsiombikas <nuclear@member.fsf.org>
date	Mon, 02 Apr 2012 14:42:03 +0300
parents	0b73aa7317e1
children	3e53a16d4667

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 reshape(int x, int y);

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

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

28 void motion(int x, int y);

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

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

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

33 static int round_pow2(int x);

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

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

37 float cam_x, cam_y, cam_z;

39 vec2_t tex_scale;

40 struct slice_file *flist;

41 int nslices;

42 unsigned int sdr, vol_tex, ray_tex;

43 int win_xsz, win_ysz;

44 int raytex_needs_recalc = 1;

46 unsigned int xfer_tex;

47 float xfer_mean = 0.5, xfer_sdev = 1.0;

48 int xfertex_needs_recalc = 1;

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

51 {

52 glutInit(&argc, argv);

53 glutInitWindowSize(1280, 720);

54 glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);

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

57 return 1;

58 }

60 glutCreateWindow("Volume Raytracer");

62 glutDisplayFunc(disp);

63 glutReshapeFunc(reshape);

64 glutKeyboardFunc(keyb);

65 glutMouseFunc(mouse);

66 glutMotionFunc(motion);

68 glewInit();

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

71 return 1;

72 }

74 glutMainLoop();

75 return 0;

76 }

78 int init(void)

79 {

80 int i, vol_xsz, vol_ysz;

82 if(!(sdr = create_program_load("volray.v.glsl", "volray.p.glsl"))) {

83 return -1;

84 }

85 set_uniform_int(sdr, "volume", 0);

86 set_uniform_int(sdr, "ray_tex", 1);

87 set_uniform_int(sdr, "xfer_tex", 2);

89 glGenTextures(1, &vol_tex);

90 glBindTexture(GL_TEXTURE_3D, vol_tex);

91 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

92 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

93 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);

94 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);

95 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);

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

98 int xsz, ysz;

99 void *pix;

100 struct slice_file *sfile = flist;

101 flist = flist->next;

102

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

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

105 return -1;

106 }

107

108 if(i == 0) {

109 /* allocate storage for the texture */

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

111

112 vol_xsz = xsz;

113 vol_ysz = ysz;

114

115 } else {

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

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

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

119 img_free_pixels(pix);

120 return -1;

121 }

122 }

123 free(sfile);

124

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

126 img_free_pixels(pix);

127 }

128 return 0;

129 }

130

131 void disp(void)

132 {

133 if(raytex_needs_recalc) {

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

135 }

136 if(xfertex_needs_recalc) {

137 create_transfer_map(xfer_mean, xfer_sdev);

138 }

139

140 glMatrixMode(GL_MODELVIEW);

141 glLoadIdentity();

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

143 glTranslatef(cam_x, cam_y, -cam_z);

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

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

146 glTranslatef(0, 0, -cam_dist);

147

148 glMatrixMode(GL_TEXTURE);

149 glLoadIdentity();

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

151

152 glActiveTexture(GL_TEXTURE0);

153 glBindTexture(GL_TEXTURE_3D, vol_tex);

154 glEnable(GL_TEXTURE_3D);

155

156 glActiveTexture(GL_TEXTURE1);

157 glBindTexture(GL_TEXTURE_2D, ray_tex);

158 glEnable(GL_TEXTURE_2D);

159

160 glActiveTexture(GL_TEXTURE2);

161 glBindTexture(GL_TEXTURE_1D, xfer_tex);

162 glEnable(GL_TEXTURE_1D);

163

164 bind_program(sdr);

165 glBegin(GL_QUADS);

166 glColor3f(1, 1, 1);

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

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

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

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

171 glEnd();

172 bind_program(0);

173

174 glActiveTexture(GL_TEXTURE2);

175 glDisable(GL_TEXTURE_1D);

176 glActiveTexture(GL_TEXTURE1);

177 glDisable(GL_TEXTURE_2D);

178 glActiveTexture(GL_TEXTURE0);

179 glDisable(GL_TEXTURE_3D);

180

181 glMatrixMode(GL_TEXTURE);

182 glLoadIdentity();

183

184 glutSwapBuffers();

185 assert(glGetError() == GL_NO_ERROR);

186 }

187

188 void reshape(int x, int y)

189 {

190 printf("reshape: %dx%d\n", x, y);

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

192

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

194 raytex_needs_recalc = 1;

195 win_xsz = x;

196 win_ysz = y;

197 }

198 }

199

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

201 {

202 switch(key) {

203 case 27:

204 exit(0);

205 }

206 }

207

208 static int bnstate[32];

209 static int prev_x, prev_y;

210

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

212 {

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

214 prev_x = x;

215 prev_y = y;

216 }

217

218 void motion(int x, int y)

219 {

220 int dx = x - prev_x;

221 int dy = y - prev_y;

222

223 prev_x = x;

224 prev_y = y;

225

226 if(bnstate[0]) {

227 cam_theta += dx * 0.5;

228 cam_phi += dy * 0.5;

229

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

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

232

233 glutPostRedisplay();

234 }

235

236 if(bnstate[1]) {

237 cam_x += dx * 0.025;

238 cam_y += dy * 0.025;

239 glutPostRedisplay();

240 }

241

242 if(bnstate[2]) {

243 cam_dist += dy * 0.025;

244 if(cam_dist < 0.0) cam_dist = 0.0;

245 glutPostRedisplay();

246 }

247 }

248

249

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

251 {

252 int i;

253 struct slice_file *tail;

254 char *endp;

255

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

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

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

259 case 'm':

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

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

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

263 return -1;

264 }

265 break;

266

267 case 'v':

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

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

270 fprintf(stderr, "-v must be followed by the transfer function sdeviance\n");

271 return -1;

272 }

273 break;

274

275 default:

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

277 return -1;

278 }

279 } else {

280 struct slice_file *sfile;

281

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

283 perror("failed to allocate memory");

284 return -1;

285 }

286 sfile->name = argv[i];

287 sfile->next = 0;

288

289 if(!flist) {

290 flist = tail = sfile;

291 } else {

292 tail->next = sfile;

293 tail = sfile;

294 }

295 nslices++;

296 }

297 }

298

299 if(!nslices) {

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

301 return -1;

302 }

303 return 0;

304 }

305

306

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

308 {

309 int i, j;

310 int cur_tex_xsz, cur_tex_ysz;

311 int tex_xsz = round_pow2(xsz);

312 int tex_ysz = round_pow2(ysz);

313 float *teximg, *dir;

314

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

316 return;

317 }

318 dir = teximg;

319

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

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

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

323 *dir++ = rdir.x;

324 *dir++ = rdir.y;

325 *dir++ = rdir.z;

326 }

327 }

328

329 if(!ray_tex) {

330 glGenTextures(1, &ray_tex);

331 }

332

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

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

335

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

337 glBindTexture(GL_TEXTURE_2D, ray_tex);

338 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);

339 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);

340 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);

341 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);

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

343 }

344

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

346 free(teximg);

347

348 if(tex_scale) {

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

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

351 }

352 raytex_needs_recalc = 0;

353 }

354

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

356 {

357 float vfov = M_PI * vfov_deg / 180.0;

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

359

360 float ysz = 2.0;

361 float xsz = aspect * ysz;

362

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

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

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

366

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

368

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

370 }

371

372 static int round_pow2(int x)

373 {

374 x--;

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

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

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

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

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

380 return x + 1;

381 }

382

383 static void create_transfer_map(float mean, float sdev)

384 {

385 static float map[XFER_MAP_SZ];

386 int i;

387

388 if(!xfer_tex) {

389 glGenTextures(1, &xfer_tex);

390 glBindTexture(GL_TEXTURE_1D, xfer_tex);

391 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

392 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

393 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

394 glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

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

396 }

397

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

399 float x = (float)i / (float)XFER_MAP_SZ;

400 map[i] = gaussian(x, mean, sdev);

401 }

402

403 glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, XFER_MAP_SZ, 1, GL_LUMINANCE, GL_FLOAT, map);

404 xfertex_needs_recalc = 0;

405 }