clray: 54a96b738afe src/scene.cc

clray

view src/scene.cc @ 53:54a96b738afe

removed the unnecessary second pass for the kdtree flattening

author	John Tsiombikas <nuclear@member.fsf.org>
date	Sun, 05 Sep 2010 16:43:55 +0100
parents	e7f79c6ad246
children	6a30f27fa1e6

line source

1 #include <stdlib.h>

2 #include <math.h>

3 #include <float.h>

4 #include <assert.h>

5 #include <map>

6 #include "scene.h"

7 #include "ogl.h"

9 #define CHECK_AABB(aabb) \

10 assert(aabb.max[0] >= aabb.min[0] && aabb.max[1] >= aabb.min[1] && aabb.max[2] >= aabb.min[2])

13 #define MIN(a, b) ((a) < (b) ? (a) : (b))

14 #define MAX(a, b) ((a) > (b) ? (a) : (b))

17 static int flatten_kdtree(const KDNode *node, KDNodeGPU *kdbuf, int *count);

18 static void draw_kdtree(const KDNode *node, int level = 0);

19 static bool build_kdtree(KDNode *kd, const Face *faces, int level = 0);

20 static float eval_cost(const Face *faces, const int *face_idx, int num_faces, const AABBox &aabb, int axis);

21 static void free_kdtree(KDNode *node);

22 static void print_item_counts(const KDNode *node, int level);

25 static int accel_param[NUM_ACCEL_PARAMS] = {

26 64, // max tree depth

27 MAX_NODE_FACES, // max items per node (0 means ignore limit)

28 5, // estimated traversal cost

29 15 // estimated interseciton cost

30 };

33 void set_accel_param(int p, int v)

34 {

35 assert(p >= 0 && p < NUM_ACCEL_PARAMS);

36 accel_param[p] = v;

37 }

40 #define FEQ(a, b) (fabs((a) - (b)) < 1e-8)

41 bool Face::operator ==(const Face &f) const

42 {

43 for(int i=0; i<3; i++) {

44 for(int j=0; j<3; j++) {

45 if(!FEQ(v[i].pos[j], f.v[i].pos[j])) {

46 return false;

47 }

48 if(!FEQ(v[i].normal[j], f.v[i].normal[j])) {

49 return false;

50 }

51 }

52 if(!FEQ(normal[i], f.normal[i])) {

53 return false;

54 }

55 }

56 return true;

57 }

59 float AABBox::calc_surface_area() const

60 {

61 float area1 = (max[0] - min[0]) * (max[1] - min[1]);

62 float area2 = (max[3] - min[3]) * (max[1] - min[1]);

63 float area3 = (max[0] - min[0]) * (max[3] - min[3]);

65 return 2.0f * (area1 + area2 + area3);

66 }

68 KDNode::KDNode()

69 {

70 left = right = 0;

71 cost = 0.0;

72 }

75 Scene::Scene()

76 {

77 facebuf = 0;

78 num_faces = -1;

79 kdtree = 0;

80 kdbuf = 0;

81 }

83 Scene::~Scene()

84 {

85 delete [] facebuf;

86 delete [] kdbuf;

87 free_kdtree(kdtree);

88 }

90 bool Scene::add_mesh(Mesh *m)

91 {

92 // make sure triangles have material ids

93 for(size_t i=0; i<m->faces.size(); i++) {

94 m->faces[i].matid = m->matid;

95 }

97 try {

98 meshes.push_back(m);

99 }

100 catch(...) {

101 return false;

102 }

103

104 // invalidate facebuffer and count

105 delete [] facebuf;

106 facebuf = 0;

107 num_faces = -1;

108

109 return true;

110 }

111

112 int Scene::get_num_meshes() const

113 {

114 return (int)meshes.size();

115 }

116

117 int Scene::get_num_faces() const

118 {

119 if(num_faces >= 0) {

120 return num_faces;

121 }

122

123 num_faces = 0;

124 for(size_t i=0; i<meshes.size(); i++) {

125 num_faces += meshes[i]->faces.size();

126 }

127 return num_faces;

128 }

129

130 int Scene::get_num_materials() const

131 {

132 return (int)matlib.size();

133 }

134

135 int Scene::get_num_kdnodes() const

136 {

137 return kdtree_nodes(kdtree);

138 }

139

140 Material *Scene::get_materials()

141 {

142 if(matlib.empty()) {

143 return 0;

144 }

145 return &matlib[0];

146 }

147

148 const Material *Scene::get_materials() const

149 {

150 if(matlib.empty()) {

151 return 0;

152 }

153 return &matlib[0];

154 }

155

156 const Face *Scene::get_face_buffer() const

157 {

158 if(facebuf) {

159 return facebuf;

160 }

161

162 int num_meshes = get_num_meshes();

163

164 printf("constructing face buffer with %d faces (out of %d meshes)\n", get_num_faces(), num_meshes);

165 facebuf = new Face[num_faces];

166 Face *fptr = facebuf;

167

168 for(int i=0; i<num_meshes; i++) {

169 for(size_t j=0; j<meshes[i]->faces.size(); j++) {

170 *fptr++ = meshes[i]->faces[j];

171 }

172 }

173 return facebuf;

174 }

175

176 const KDNodeGPU *Scene::get_kdtree_buffer() const

177 {

178 if(kdbuf) {

179 return kdbuf;

180 }

181

182 if(!kdtree) {

183 ((Scene*)this)->build_kdtree();

184 }

185

186 int num_nodes = get_num_kdnodes();

187 kdbuf = new KDNodeGPU[num_nodes];

188

189 int count = 0;

190

191 // first arrange the kdnodes into an array (flatten)

192 flatten_kdtree(kdtree, kdbuf, &count);

193

194 return kdbuf;

195 }

196

197 static int flatten_kdtree(const KDNode *node, KDNodeGPU *kdbuf, int *count)

198 {

199 const size_t max_node_items = sizeof kdbuf[0].face_idx / sizeof kdbuf[0].face_idx[0];

200 int idx = (*count)++;

201

202 // copy the node

203 kdbuf[idx].aabb = node->aabb;

204 kdbuf[idx].num_faces = 0;

205

206 for(size_t i=0; i<node->face_idx.size(); i++) {

207 if(i >= max_node_items) {

208 fprintf(stderr, "WARNING too many faces per leaf node!\n");

209 break;

210 }

211 kdbuf[idx].face_idx[i] = node->face_idx[i];

212 kdbuf[idx].num_faces++;

213 }

214

215 // recurse to the left/right (if we're not in a leaf node)

216 if(node->left) {

217 assert(node->right);

218

219 kdbuf[idx].left = flatten_kdtree(node->left, kdbuf, count);

220 kdbuf[idx].right = flatten_kdtree(node->right, kdbuf, count);

221 } else {

222 kdbuf[idx].left = kdbuf[idx].right = -1;

223 }

224

225 return idx;

226 }

227

228 void Scene::draw_kdtree() const

229 {

230 glPushAttrib(GL_ENABLE_BIT);

231 glDisable(GL_LIGHTING);

232 glDepthMask(0);

233

234 glBegin(GL_LINES);

235 ::draw_kdtree(kdtree, 0);

236 glEnd();

237

238 glDepthMask(1);

239 glPopAttrib();

240 }

241

242 static float palette[][3] = {

243 {0, 1, 0},

244 {1, 0, 0},

245 {0, 0, 1},

246 {1, 1, 0},

247 {0, 0, 1},

248 {1, 0, 1}

249 };

250 static int pal_size = sizeof palette / sizeof *palette;

251

252 static void draw_kdtree(const KDNode *node, int level)

253 {

254 if(!node) return;

255

256 draw_kdtree(node->left, level + 1);

257 draw_kdtree(node->right, level + 1);

258

259 glColor3fv(palette[level % pal_size]);

260

261 glVertex3fv(node->aabb.min);

262 glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]);

263 glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]);

264 glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]);

265 glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]);

266 glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]);

267 glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]);

268 glVertex3fv(node->aabb.min);

269

270 glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]);

271 glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]);

272 glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]);

273 glVertex3fv(node->aabb.max);

274 glVertex3fv(node->aabb.max);

275 glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]);

276 glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]);

277 glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]);

278

279 glVertex3fv(node->aabb.min);

280 glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]);

281 glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]);

282 glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]);

283 glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]);

284 glVertex3fv(node->aabb.max);

285 glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]);

286 glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]);

287 }

288

289 bool Scene::build_kdtree()

290 {

291 assert(kdtree == 0);

292

293 const Face *faces = get_face_buffer();

294 int num_faces = get_num_faces();

295

296 printf("Constructing kd-tree out of %d faces ...\n", num_faces);

297

298 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

299 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

300 printf(" max items per leaf: %d\n", accel_param[ACCEL_PARAM_MAX_NODE_ITEMS]);

301 printf(" SAH parameters - tcost: %d - icost: %d\n", tcost, icost);

302

303 free_kdtree(kdtree);

304 kdtree = new KDNode;

305

306 /* Start the construction of the kdtree by adding all faces of the scene

307 * to the new root node. At the same time calculate the root's AABB.

308 */

309 kdtree->aabb.min[0] = kdtree->aabb.min[1] = kdtree->aabb.min[2] = FLT_MAX;

310 kdtree->aabb.max[0] = kdtree->aabb.max[1] = kdtree->aabb.max[2] = -FLT_MAX;

311

312 for(int i=0; i<num_faces; i++) {

313 const Face *face = faces + i;

314

315 // for each vertex of the face ...

316 for(int j=0; j<3; j++) {

317 const float *pos = face->v[j].pos;

318

319 // for each element (xyz) of the position vector ...

320 for(int k=0; k<3; k++) {

321 if(pos[k] < kdtree->aabb.min[k]) {

322 kdtree->aabb.min[k] = pos[k];

323 }

324 if(pos[k] > kdtree->aabb.max[k]) {

325 kdtree->aabb.max[k] = pos[k];

326 }

327 }

328 }

329

330 kdtree->face_idx.push_back(i); // add the face

331 }

332

333 CHECK_AABB(kdtree->aabb);

334

335 // calculate the heuristic for the root

336 kdtree->cost = eval_cost(faces, &kdtree->face_idx[0], kdtree->face_idx.size(), kdtree->aabb, 0);

337

338 // now proceed splitting the root recursively

339 if(!::build_kdtree(kdtree, faces)) {

340 fprintf(stderr, "failed to build kdtree\n");

341 return false;

342 }

343

344 printf(" tree depth: %d\n", kdtree_depth(kdtree));

345 print_item_counts(kdtree, 0);

346 return true;

347 }

348

349 struct Split {

350 int axis;

351 float pos;

352 float sum_cost;

353 float cost_left, cost_right;

354 };

355

356 static void find_best_split(const KDNode *node, int axis, const Face *faces, Split *split)

357 {

358 Split best_split;

359 best_split.sum_cost = FLT_MAX;

360

361 for(size_t i=0; i<node->face_idx.size(); i++) {

362 const Face *face = faces + node->face_idx[i];

363

364 float splitpt[2];

365 splitpt[0] = MIN(face->v[0].pos[axis], MIN(face->v[1].pos[axis], face->v[2].pos[axis]));

366 splitpt[1] = MAX(face->v[0].pos[axis], MAX(face->v[1].pos[axis], face->v[2].pos[axis]));

367

368 for(int j=0; j<2; j++) {

369 if(splitpt[j] <= node->aabb.min[axis] || splitpt[j] >= node->aabb.max[axis]) {

370 continue;

371 }

372

373 AABBox aabb_left, aabb_right;

374 aabb_left = aabb_right = node->aabb;

375 aabb_left.max[axis] = splitpt[j];

376 aabb_right.min[axis] = splitpt[j];

377

378 float left_cost = eval_cost(faces, &node->face_idx[0], node->face_idx.size(), aabb_left, axis);

379 float right_cost = eval_cost(faces, &node->face_idx[0], node->face_idx.size(), aabb_right, axis);

380 float sum_cost = left_cost + right_cost - accel_param[ACCEL_PARAM_COST_TRAVERSE]; // tcost is added twice

381

382 if(sum_cost < best_split.sum_cost) {

383 best_split.cost_left = left_cost;

384 best_split.cost_right = right_cost;

385 best_split.sum_cost = sum_cost;

386 best_split.pos = splitpt[j];

387 }

388 }

389 }

390

391 assert(split);

392 *split = best_split;

393 split->axis = axis;

394 }

395

396 static bool build_kdtree(KDNode *kd, const Face *faces, int level)

397 {

398 int opt_max_depth = accel_param[ACCEL_PARAM_MAX_TREE_DEPTH];

399 int opt_max_items = accel_param[ACCEL_PARAM_MAX_NODE_ITEMS];

400

401 if(kd->face_idx.empty() || level >= opt_max_depth) {

402 return true;

403 }

404

405 Split best_split;

406 best_split.axis = -1;

407 best_split.sum_cost = FLT_MAX;

408

409 for(int i=0; i<3; i++) {

410 Split split;

411 find_best_split(kd, i, faces, &split);

412

413 if(split.sum_cost < best_split.sum_cost) {

414 best_split = split;

415 }

416 }

417

418 if(best_split.axis == -1) {

419 return true; // can't split any more, only 0-area splits available

420 }

421

422 //printf("current cost: %f, best_cost: %f\n", kd->cost, best_sum_cost);

423 if(best_split.sum_cost > kd->cost && (opt_max_items == 0 || (int)kd->face_idx.size() <= opt_max_items)) {

424 return true; // stop splitting if it doesn't reduce the cost

425 }

426

427 kd->axis = best_split.axis;

428

429 // create the two children

430 KDNode *kdleft, *kdright;

431 kdleft = new KDNode;

432 kdright = new KDNode;

433

434 kdleft->aabb = kdright->aabb = kd->aabb;

435

436 kdleft->aabb.max[kd->axis] = best_split.pos;

437 kdright->aabb.min[kd->axis] = best_split.pos;

438

439 kdleft->cost = best_split.cost_left;

440 kdright->cost = best_split.cost_right;

441

442 // TODO would it be much better if we actually split faces that straddle the splitting plane?

443 for(size_t i=0; i<kd->face_idx.size(); i++) {

444 int fidx = kd->face_idx[i];

445 const Face *face = faces + fidx;

446

447 if(face->v[0].pos[kd->axis] < best_split.pos ||

448 face->v[1].pos[kd->axis] < best_split.pos ||

449 face->v[2].pos[kd->axis] < best_split.pos) {

450 kdleft->face_idx.push_back(fidx);

451 }

452 if(face->v[0].pos[kd->axis] >= best_split.pos ||

453 face->v[1].pos[kd->axis] >= best_split.pos ||

454 face->v[2].pos[kd->axis] >= best_split.pos) {

455 kdright->face_idx.push_back(fidx);

456 }

457 }

458 kd->face_idx.clear(); // only leaves have faces

459

460 kd->left = kdleft;

461 kd->right = kdright;

462

463 return build_kdtree(kd->left, faces, level + 1) && build_kdtree(kd->right, faces, level + 1);

464 }

465

466 static float eval_cost(const Face *faces, const int *face_idx, int num_faces, const AABBox &aabb, int axis)

467 {

468 int num_inside = 0;

469 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

470 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

471

472 for(int i=0; i<num_faces; i++) {

473 const Face *face = faces + face_idx[i];

474

475 for(int j=0; j<3; j++) {

476 if(face->v[j].pos[axis] >= aabb.min[axis] && face->v[j].pos[axis] < aabb.max[axis]) {

477 num_inside++;

478 break;

479 }

480 }

481 }

482

483 float dx = aabb.max[0] - aabb.min[0];

484 float dy = aabb.max[1] - aabb.min[1];

485 float dz = aabb.max[2] - aabb.min[2];

486

487 if(dx < 0.0) {

488 fprintf(stderr, "FOO DX = %f\n", dx);

489 abort();

490 }

491 if(dy < 0.0) {

492 fprintf(stderr, "FOO DX = %f\n", dy);

493 abort();

494 }

495 if(dz < 0.0) {

496 fprintf(stderr, "FOO DX = %f\n", dz);

497 abort();

498 }

499

500 if(dx < 1e-6 || dy < 1e-6 || dz < 1e-6) {

501 return FLT_MAX; // heavily penalize 0-area voxels

502 }

503

504 float sarea = 2.0 * (dx + dy + dz);//aabb.calc_surface_area();

505 return tcost + sarea * num_inside * icost;

506 }

507

508 static void free_kdtree(KDNode *node)

509 {

510 if(node) {

511 free_kdtree(node->left);

512 free_kdtree(node->right);

513 delete node;

514 }

515 }

516

517 int kdtree_depth(const KDNode *node)

518 {

519 if(!node) return 0;

520

521 int left = kdtree_depth(node->left);

522 int right = kdtree_depth(node->right);

523 return (left > right ? left : right) + 1;

524 }

525

526 int kdtree_nodes(const KDNode *node)

527 {

528 if(!node) return 0;

529 return kdtree_nodes(node->left) + kdtree_nodes(node->right) + 1;

530 }

531

532 static void print_item_counts(const KDNode *node, int level)

533 {

534 if(!node) return;

535

536 for(int i=0; i<level; i++) {

537 fputs(" ", stdout);

538 }

539 printf("- %d (cost: %f)\n", (int)node->face_idx.size(), node->cost);

540

541 print_item_counts(node->left, level + 1);

542 print_item_counts(node->right, level + 1);

543 }