clray: 3d13924b22e6 src/scene.cc

clray

view src/scene.cc @ 58:3d13924b22e6

implementing polygon split

author	John Tsiombikas <nuclear@member.fsf.org>
date	Sun, 12 Sep 2010 00:19:04 +0100
parents	6a30f27fa1e6
children	eb97f9c92e1d

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);

23 static int split_face(const Face *inface, int axis, Face *face1, Face *face2);

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 bool Scene::add_light(const Light &lt)

114 {

115 try {

116 lights.push_back(lt);

117 }

118 catch(...) {

119 return false;

120 }

121 return true;

122 }

123

124 int Scene::get_num_meshes() const

125 {

126 return (int)meshes.size();

127 }

128

129 int Scene::get_num_lights() const

130 {

131 return (int)lights.size();

132 }

133

134 int Scene::get_num_faces() const

135 {

136 if(num_faces >= 0) {

137 return num_faces;

138 }

139

140 num_faces = 0;

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

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

143 }

144 return num_faces;

145 }

146

147 int Scene::get_num_materials() const

148 {

149 return (int)matlib.size();

150 }

151

152 int Scene::get_num_kdnodes() const

153 {

154 return kdtree_nodes(kdtree);

155 }

156

157 Mesh **Scene::get_meshes()

158 {

159 if(meshes.empty()) {

160 return 0;

161 }

162 return &meshes[0];

163 }

164

165 const Mesh * const *Scene::get_meshes() const

166 {

167 if(meshes.empty()) {

168 return 0;

169 }

170 return &meshes[0];

171 }

172

173 Light *Scene::get_lights()

174 {

175 if(lights.empty()) {

176 return 0;

177 }

178 return &lights[0];

179 }

180

181 const Light *Scene::get_lights() const

182 {

183 if(lights.empty()) {

184 return 0;

185 }

186 return &lights[0];

187 }

188

189 Material *Scene::get_materials()

190 {

191 if(matlib.empty()) {

192 return 0;

193 }

194 return &matlib[0];

195 }

196

197 const Material *Scene::get_materials() const

198 {

199 if(matlib.empty()) {

200 return 0;

201 }

202 return &matlib[0];

203 }

204

205 const Face *Scene::get_face_buffer() const

206 {

207 if(facebuf) {

208 return facebuf;

209 }

210

211 int num_meshes = get_num_meshes();

212

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

214 facebuf = new Face[num_faces];

215 Face *fptr = facebuf;

216

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

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

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

220 }

221 }

222 return facebuf;

223 }

224

225 const KDNodeGPU *Scene::get_kdtree_buffer() const

226 {

227 if(kdbuf) {

228 return kdbuf;

229 }

230

231 if(!kdtree) {

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

233 }

234

235 int num_nodes = get_num_kdnodes();

236 kdbuf = new KDNodeGPU[num_nodes];

237

238 int count = 0;

239

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

241 flatten_kdtree(kdtree, kdbuf, &count);

242

243 return kdbuf;

244 }

245

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

247 {

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

249 int idx = (*count)++;

250

251 // copy the node

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

253 kdbuf[idx].num_faces = 0;

254

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

256 if(i >= max_node_items) {

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

258 break;

259 }

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

261 kdbuf[idx].num_faces++;

262 }

263

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

265 if(node->left) {

266 assert(node->right);

267

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

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

270 } else {

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

272 }

273

274 return idx;

275 }

276

277 void Scene::draw_kdtree() const

278 {

279 glPushAttrib(GL_ENABLE_BIT);

280 glDisable(GL_LIGHTING);

281 glDepthMask(0);

282

283 glBegin(GL_LINES);

284 ::draw_kdtree(kdtree, 0);

285 glEnd();

286

287 glDepthMask(1);

288 glPopAttrib();

289 }

290

291 static float palette[][3] = {

292 {0, 1, 0},

293 {1, 0, 0},

294 {0, 0, 1},

295 {1, 1, 0},

296 {0, 0, 1},

297 {1, 0, 1}

298 };

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

300

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

302 {

303 if(!node) return;

304

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

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

307

308 glColor3fv(palette[level % pal_size]);

309

310 glVertex3fv(node->aabb.min);

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

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

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

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

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

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

317 glVertex3fv(node->aabb.min);

318

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

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

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

322 glVertex3fv(node->aabb.max);

323 glVertex3fv(node->aabb.max);

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

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

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

327

328 glVertex3fv(node->aabb.min);

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

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

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

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

333 glVertex3fv(node->aabb.max);

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

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

336 }

337

338 bool Scene::build_kdtree()

339 {

340 assert(kdtree == 0);

341

342 const Face *faces = get_face_buffer();

343 int num_faces = get_num_faces();

344

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

346

347 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

348 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

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

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

351

352 free_kdtree(kdtree);

353 kdtree = new KDNode;

354

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

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

357 */

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

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

360

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

362 const Face *face = faces + i;

363

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

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

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

367

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

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

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

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

372 }

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

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

375 }

376 }

377 }

378

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

380 }

381

382 CHECK_AABB(kdtree->aabb);

383

384 // calculate the heuristic for the root

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

386

387 // now proceed splitting the root recursively

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

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

390 return false;

391 }

392

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

394 print_item_counts(kdtree, 0);

395 return true;

396 }

397

398 struct Split {

399 int axis;

400 float pos;

401 float sum_cost;

402 float cost_left, cost_right;

403 };

404

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

406 {

407 Split best_split;

408 best_split.sum_cost = FLT_MAX;

409

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

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

412

413 float splitpt[2];

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

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

416

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

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

419 continue;

420 }

421

422 AABBox aabb_left, aabb_right;

423 aabb_left = aabb_right = node->aabb;

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

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

426

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

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

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

430

431 if(sum_cost < best_split.sum_cost) {

432 best_split.cost_left = left_cost;

433 best_split.cost_right = right_cost;

434 best_split.sum_cost = sum_cost;

435 best_split.pos = splitpt[j];

436 }

437 }

438 }

439

440 assert(split);

441 *split = best_split;

442 split->axis = axis;

443 }

444

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

446 {

447 int opt_max_depth = accel_param[ACCEL_PARAM_MAX_TREE_DEPTH];

448 int opt_max_items = accel_param[ACCEL_PARAM_MAX_NODE_ITEMS];

449

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

451 return true;

452 }

453

454 Split best_split;

455 best_split.axis = -1;

456 best_split.sum_cost = FLT_MAX;

457

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

459 Split split;

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

461

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

463 best_split = split;

464 }

465 }

466

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

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

469 }

470

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

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

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

474 }

475

476 kd->axis = best_split.axis;

477

478 // create the two children

479 KDNode *kdleft, *kdright;

480 kdleft = new KDNode;

481 kdright = new KDNode;

482

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

484

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

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

487

488 kdleft->cost = best_split.cost_left;

489 kdright->cost = best_split.cost_right;

490

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

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

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

494 const Face *face = faces + fidx;

495

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

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

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

499 kdleft->face_idx.push_back(fidx);

500 }

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

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

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

504 kdright->face_idx.push_back(fidx);

505 }

506 }

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

508

509 kd->left = kdleft;

510 kd->right = kdright;

511

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

513 }

514

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

516 {

517 int num_inside = 0;

518 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

519 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

520

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

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

523

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

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

526 num_inside++;

527 break;

528 }

529 }

530 }

531

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

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

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

535

536 if(dx < 0.0) {

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

538 abort();

539 }

540 if(dy < 0.0) {

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

542 abort();

543 }

544 if(dz < 0.0) {

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

546 abort();

547 }

548

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

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

551 }

552

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

554 return tcost + sarea * num_inside * icost;

555 }

556

557 static void free_kdtree(KDNode *node)

558 {

559 if(node) {

560 free_kdtree(node->left);

561 free_kdtree(node->right);

562 delete node;

563 }

564 }

565

566 int kdtree_depth(const KDNode *node)

567 {

568 if(!node) return 0;

569

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

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

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

573 }

574

575 int kdtree_nodes(const KDNode *node)

576 {

577 if(!node) return 0;

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

579 }

580

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

582 {

583 if(!node) return;

584

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

586 fputs(" ", stdout);

587 }

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

589

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

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

592 }

593 /*

594 #define SWAP(a, b, type) \

595 do { \

596 type tmp = a; \

597 a = b; \

598 b = tmp; \

599 } while(0)

600

601 static bool clip_face(const Face *inface, float splitpos, int axis, Face *face1, Face *face2)

602 {

603 assert(inface && face1 && face2);

604 assert(inface != face1 && inface != face2);

605 assert(axis >= 0 && axis < 3);

606

607 std::vector<Vertex> verts;

608

609 // find the edges that must be split

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

611 verts.push_back(inface->v[i]);

612

613 float start = inface->v[i].pos[axis];

614 float end = inface->v[(i + 1) % 3].pos[axis];

615

616 if((splitpos >= start && splitpos < end) || (splitpos >= end && splitpos < start)) {

617

618 }

619 }

620 }*/