clray: e7f79c6ad246 src/scene.cc

clray

view src/scene.cc @ 44:e7f79c6ad246

added maximum node face limit by default

author	John Tsiombikas <nuclear@member.fsf.org>
date	Sat, 28 Aug 2010 21:50:17 +0100
parents	4bcf78e572d6
children	54a96b738afe

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 void flatten_kdtree(const KDNode *node, KDNodeGPU *kdbuf, int *count, std::map<const KDNode*, int> *flatmap);

18 static void fix_child_links(const KDNode *node, KDNodeGPU *kdbuf, const std::map<const KDNode*, int> &flatmap);

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

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

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

22 static void free_kdtree(KDNode *node);

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

26 static int accel_param[NUM_ACCEL_PARAMS] = {

27 64, // max tree depth

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

29 5, // estimated traversal cost

30 15 // estimated interseciton cost

31 };

34 void set_accel_param(int p, int v)

35 {

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

37 accel_param[p] = v;

38 }

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

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

43 {

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

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

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

47 return false;

48 }

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

50 return false;

51 }

52 }

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

54 return false;

55 }

56 }

57 return true;

58 }

60 float AABBox::calc_surface_area() const

61 {

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

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

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

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

67 }

69 KDNode::KDNode()

70 {

71 left = right = 0;

72 cost = 0.0;

73 }

76 Scene::Scene()

77 {

78 facebuf = 0;

79 num_faces = -1;

80 kdtree = 0;

81 kdbuf = 0;

82 }

84 Scene::~Scene()

85 {

86 delete [] facebuf;

87 delete [] kdbuf;

88 free_kdtree(kdtree);

89 }

91 bool Scene::add_mesh(Mesh *m)

92 {

93 // make sure triangles have material ids

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

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

96 }

98 try {

99 meshes.push_back(m);

100 }

101 catch(...) {

102 return false;

103 }

104

105 // invalidate facebuffer and count

106 delete [] facebuf;

107 facebuf = 0;

108 num_faces = -1;

109

110 return true;

111 }

112

113 int Scene::get_num_meshes() const

114 {

115 return (int)meshes.size();

116 }

117

118 int Scene::get_num_faces() const

119 {

120 if(num_faces >= 0) {

121 return num_faces;

122 }

123

124 num_faces = 0;

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

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

127 }

128 return num_faces;

129 }

130

131 int Scene::get_num_materials() const

132 {

133 return (int)matlib.size();

134 }

135

136 int Scene::get_num_kdnodes() const

137 {

138 return kdtree_nodes(kdtree);

139 }

140

141 Material *Scene::get_materials()

142 {

143 if(matlib.empty()) {

144 return 0;

145 }

146 return &matlib[0];

147 }

148

149 const Material *Scene::get_materials() const

150 {

151 if(matlib.empty()) {

152 return 0;

153 }

154 return &matlib[0];

155 }

156

157 const Face *Scene::get_face_buffer() const

158 {

159 if(facebuf) {

160 return facebuf;

161 }

162

163 int num_meshes = get_num_meshes();

164

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

166 facebuf = new Face[num_faces];

167 Face *fptr = facebuf;

168

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

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

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

172 }

173 }

174 return facebuf;

175 }

176

177 static int find_node_index(const KDNode *node, const std::map<const KDNode*, int> &flatmap);

178

179 static bool validate_flat_tree(const KDNode *node, const KDNodeGPU *kdbuf, int count, const std::map<const KDNode*, int> &flatmap)

180 {

181 if(!node) return true;

182

183 int idx = find_node_index(node, flatmap);

184 int left_idx = find_node_index(node->left, flatmap);

185 int right_idx = find_node_index(node->right, flatmap);

186

187 if(kdbuf[idx].left != left_idx) {

188 fprintf(stderr, "%d's left should be %d. it's: %d\n", idx, left_idx, kdbuf[idx].left);

189 return false;

190 }

191 if(kdbuf[idx].right != right_idx) {

192 fprintf(stderr, "%d's right should be %d. it's: %d\n", idx, right_idx, kdbuf[idx].right);

193 return false;

194 }

195 return validate_flat_tree(node->left, kdbuf, count, flatmap) && validate_flat_tree(node->right, kdbuf, count, flatmap);

196 }

197

198 const KDNodeGPU *Scene::get_kdtree_buffer() const

199 {

200 if(kdbuf) {

201 return kdbuf;

202 }

203

204 if(!kdtree) {

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

206 }

207

208 /* we'll associate kdnodes with flattened array indices through this map as

209 * we add them so that the second child-index pass, will be able to find

210 * the children indices of each node.

211 */

212 std::map<const KDNode*, int> flatmap;

213 flatmap[0] = -1;

214

215 int num_nodes = get_num_kdnodes();

216 kdbuf = new KDNodeGPU[num_nodes];

217

218 int count = 0;

219

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

221 flatten_kdtree(kdtree, kdbuf, &count, &flatmap);

222

223 // then fix all the left/right links to point to the correct nodes

224 fix_child_links(kdtree, kdbuf, flatmap);

225

226 assert(validate_flat_tree(kdtree, kdbuf, count, flatmap));

227

228 return kdbuf;

229 }

230

231 static void flatten_kdtree(const KDNode *node, KDNodeGPU *kdbuf, int *count, std::map<const KDNode*, int> *flatmap)

232 {

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

234 int idx = (*count)++;

235

236 // copy the node

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

238 kdbuf[idx].num_faces = 0;

239

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

241 if(i >= max_node_items) {

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

243 break;

244 }

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

246 kdbuf[idx].num_faces++;

247 }

248

249 // update the node* -> array-position mapping

250 (*flatmap)[node] = idx;

251

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

253 if(node->left) {

254 assert(node->right);

255

256 flatten_kdtree(node->left, kdbuf, count, flatmap);

257 flatten_kdtree(node->right, kdbuf, count, flatmap);

258 }

259 }

260

261 static int find_node_index(const KDNode *node, const std::map<const KDNode*, int> &flatmap)

262 {

263 std::map<const KDNode*, int>::const_iterator it = flatmap.find(node);

264 assert(it != flatmap.end());

265 return it->second;

266 }

267

268 static void fix_child_links(const KDNode *node, KDNodeGPU *kdbuf, const std::map<const KDNode*, int> &flatmap)

269 {

270 if(!node) return;

271

272 int idx = find_node_index(node, flatmap);

273

274 kdbuf[idx].left = find_node_index(node->left, flatmap);

275 if(!node->left && kdbuf[idx].left != -1) abort();

276 kdbuf[idx].right = find_node_index(node->right, flatmap);

277 if(!node->right && kdbuf[idx].right != -1) abort();

278

279 fix_child_links(node->left, kdbuf, flatmap);

280 fix_child_links(node->right, kdbuf, flatmap);

281 }

282

283 void Scene::draw_kdtree() const

284 {

285 glPushAttrib(GL_ENABLE_BIT);

286 glDisable(GL_LIGHTING);

287 glDepthMask(0);

288

289 glBegin(GL_LINES);

290 ::draw_kdtree(kdtree, 0);

291 glEnd();

292

293 glDepthMask(1);

294 glPopAttrib();

295 }

296

297 static float palette[][3] = {

298 {0, 1, 0},

299 {1, 0, 0},

300 {0, 0, 1},

301 {1, 1, 0},

302 {0, 0, 1},

303 {1, 0, 1}

304 };

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

306

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

308 {

309 if(!node) return;

310

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

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

313

314 glColor3fv(palette[level % pal_size]);

315

316 glVertex3fv(node->aabb.min);

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

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

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

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

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

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

323 glVertex3fv(node->aabb.min);

324

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

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

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

328 glVertex3fv(node->aabb.max);

329 glVertex3fv(node->aabb.max);

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

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

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

333

334 glVertex3fv(node->aabb.min);

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

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

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

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

339 glVertex3fv(node->aabb.max);

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

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

342 /*if(!node->left) return;

343

344 AABBox *bleft = &node->left->aabb;

345

346 int axis = level % 3;

347 switch(axis) {

348 case 0:

349 glVertex3f(bleft->max[0], bleft->min[1], bleft->min[2]);

350 glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]);

351 glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]);

352 glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]);

353 glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]);

354 glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]);

355 glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]);

356 glVertex3f(bleft->max[0], bleft->min[1], bleft->min[2]);

357 break;

358

359 case 1:

360 glVertex3f(bleft->min[0], bleft->min[1], bleft->max[2]);

361 glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]);

362 glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]);

363 glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]);

364 glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]);

365 glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]);

366 glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]);

367 glVertex3f(bleft->min[0], bleft->min[1], bleft->max[2]);

368 break;

369

370 case 2:

371 glVertex3f(bleft->min[0], bleft->max[1], bleft->min[2]);

372 glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]);

373 glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]);

374 glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]);

375 glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]);

376 glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]);

377 glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]);

378 glVertex3f(bleft->min[0], bleft->max[1], bleft->min[2]);

379 break;

380

381 default:

382 break;

383 }*/

384 }

385

386 bool Scene::build_kdtree()

387 {

388 assert(kdtree == 0);

389

390 const Face *faces = get_face_buffer();

391 int num_faces = get_num_faces();

392

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

394

395 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

396 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

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

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

399

400 free_kdtree(kdtree);

401 kdtree = new KDNode;

402

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

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

405 */

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

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

408

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

410 const Face *face = faces + i;

411

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

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

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

415

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

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

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

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

420 }

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

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

423 }

424 }

425 }

426

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

428 }

429

430 CHECK_AABB(kdtree->aabb);

431

432 // calculate the heuristic for the root

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

434

435 // now proceed splitting the root recursively

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

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

438 return false;

439 }

440

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

442 print_item_counts(kdtree, 0);

443 return true;

444 }

445

446 struct Split {

447 int axis;

448 float pos;

449 float sum_cost;

450 float cost_left, cost_right;

451 };

452

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

454 {

455 Split best_split;

456 best_split.sum_cost = FLT_MAX;

457

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

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

460

461 float splitpt[2];

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

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

464

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

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

467 continue;

468 }

469

470 AABBox aabb_left, aabb_right;

471 aabb_left = aabb_right = node->aabb;

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

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

474

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

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

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

478

479 if(sum_cost < best_split.sum_cost) {

480 best_split.cost_left = left_cost;

481 best_split.cost_right = right_cost;

482 best_split.sum_cost = sum_cost;

483 best_split.pos = splitpt[j];

484 }

485 }

486 }

487

488 assert(split);

489 *split = best_split;

490 split->axis = axis;

491 }

492

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

494 {

495 int opt_max_depth = accel_param[ACCEL_PARAM_MAX_TREE_DEPTH];

496 int opt_max_items = accel_param[ACCEL_PARAM_MAX_NODE_ITEMS];

497

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

499 return true;

500 }

501

502 Split best_split;

503 best_split.axis = -1;

504 best_split.sum_cost = FLT_MAX;

505

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

507 Split split;

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

509

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

511 best_split = split;

512 }

513 }

514

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

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

517 }

518

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

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

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

522 }

523

524 kd->axis = best_split.axis;

525

526 // create the two children

527 KDNode *kdleft, *kdright;

528 kdleft = new KDNode;

529 kdright = new KDNode;

530

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

532

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

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

535

536 kdleft->cost = best_split.cost_left;

537 kdright->cost = best_split.cost_right;

538

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

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

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

542 const Face *face = faces + fidx;

543

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

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

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

547 kdleft->face_idx.push_back(fidx);

548 }

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

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

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

552 kdright->face_idx.push_back(fidx);

553 }

554 }

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

556

557 kd->left = kdleft;

558 kd->right = kdright;

559

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

561 }

562

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

564 {

565 int num_inside = 0;

566 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

567 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

568

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

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

571

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

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

574 num_inside++;

575 break;

576 }

577 }

578 }

579

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

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

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

583

584 if(dx < 0.0) {

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

586 abort();

587 }

588 if(dy < 0.0) {

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

590 abort();

591 }

592 if(dz < 0.0) {

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

594 abort();

595 }

596

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

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

599 }

600

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

602 return tcost + sarea * num_inside * icost;

603 }

604

605 static void free_kdtree(KDNode *node)

606 {

607 if(node) {

608 free_kdtree(node->left);

609 free_kdtree(node->right);

610 delete node;

611 }

612 }

613

614 int kdtree_depth(const KDNode *node)

615 {

616 if(!node) return 0;

617

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

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

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

621 }

622

623 int kdtree_nodes(const KDNode *node)

624 {

625 if(!node) return 0;

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

627 }

628

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

630 {

631 if(!node) return;

632

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

634 fputs(" ", stdout);

635 }

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

637

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

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

640 }