clray: d9520da6b801 src/scene.cc

clray

view src/scene.cc @ 62:d9520da6b801

minor readme fix

author	John Tsiombikas <nuclear@member.fsf.org>
date	Mon, 28 Dec 2015 10:31:58 +0200
parents	3d13924b22e6
children

line source

1 #include <stdlib.h>

2 #include <string.h>

3 #include <math.h>

4 #include <float.h>

5 #include <assert.h>

6 #include <map>

7 #include "scene.h"

8 #include "ogl.h"

9 #include "vector.h"

11 #define CHECK_AABB(aabb) \

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

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

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

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

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

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

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

23 static void free_kdtree(KDNode *node);

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

25 static int clip_face(const Face &inface, float splitpos, int axis, int sign, Face *faces);

26 static float calc_sq_area(const Vector3 &a, const Vector3 &b, const Vector3 &c);

29 static int accel_param[NUM_ACCEL_PARAMS] = {

30 64, // max tree depth

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

32 5, // estimated traversal cost

33 15 // estimated interseciton cost

34 };

37 void set_accel_param(int p, int v)

38 {

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

40 accel_param[p] = v;

41 }

44 float AABBox::calc_surface_area() const

45 {

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

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

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

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

51 }

53 KDNode::KDNode()

54 {

55 left = right = 0;

56 cost = 0.0;

57 }

60 Scene::Scene()

61 {

62 facebuf = 0;

63 num_faces = -1;

64 kdtree = 0;

65 kdbuf = 0;

66 }

68 Scene::~Scene()

69 {

70 delete [] facebuf;

71 delete [] kdbuf;

72 free_kdtree(kdtree);

73 }

75 bool Scene::add_mesh(Mesh *m)

76 {

77 // make sure triangles have material ids

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

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

80 }

82 try {

83 meshes.push_back(m);

84 }

85 catch(...) {

86 return false;

87 }

89 // invalidate facebuffer and count

90 delete [] facebuf;

91 facebuf = 0;

92 num_faces = -1;

94 return true;

95 }

97 bool Scene::add_light(const Light &lt)

98 {

99 try {

100 lights.push_back(lt);

101 }

102 catch(...) {

103 return false;

104 }

105 return true;

106 }

107

108 int Scene::get_num_meshes() const

109 {

110 return (int)meshes.size();

111 }

112

113 int Scene::get_num_lights() const

114 {

115 return (int)lights.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 Mesh **Scene::get_meshes()

142 {

143 if(meshes.empty()) {

144 return 0;

145 }

146 return &meshes[0];

147 }

148

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

150 {

151 if(meshes.empty()) {

152 return 0;

153 }

154 return &meshes[0];

155 }

156

157 Light *Scene::get_lights()

158 {

159 if(lights.empty()) {

160 return 0;

161 }

162 return &lights[0];

163 }

164

165 const Light *Scene::get_lights() const

166 {

167 if(lights.empty()) {

168 return 0;

169 }

170 return &lights[0];

171 }

172

173 Material *Scene::get_materials()

174 {

175 if(matlib.empty()) {

176 return 0;

177 }

178 return &matlib[0];

179 }

180

181 const Material *Scene::get_materials() const

182 {

183 if(matlib.empty()) {

184 return 0;

185 }

186 return &matlib[0];

187 }

188

189 const Face *Scene::get_face_buffer() const

190 {

191 if(facebuf) {

192 return facebuf;

193 }

194

195 int num_meshes = get_num_meshes();

196

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

198 facebuf = new Face[num_faces];

199 Face *fptr = facebuf;

200

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

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

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

204 }

205 }

206 return facebuf;

207 }

208

209 const KDNodeGPU *Scene::get_kdtree_buffer() const

210 {

211 if(kdbuf) {

212 return kdbuf;

213 }

214

215 if(!kdtree) {

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

217 }

218

219 int num_nodes = get_num_kdnodes();

220 kdbuf = new KDNodeGPU[num_nodes];

221

222 int count = 0;

223

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

225 flatten_kdtree(kdtree, kdbuf, &count);

226

227 return kdbuf;

228 }

229

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

231 {

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

233 int idx = (*count)++;

234

235 // copy the node

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

237 kdbuf[idx].num_faces = 0;

238

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

240 if(i >= max_node_items) {

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

242 break;

243 }

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

245 kdbuf[idx].num_faces++;

246 }

247

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

249 if(node->left) {

250 assert(node->right);

251

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

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

254 } else {

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

256 }

257

258 return idx;

259 }

260

261 void Scene::draw_kdtree() const

262 {

263 glPushAttrib(GL_ENABLE_BIT);

264 glDisable(GL_LIGHTING);

265 glDepthMask(0);

266

267 glBegin(GL_LINES);

268 ::draw_kdtree(kdtree, 0);

269 glEnd();

270

271 glDepthMask(1);

272 glPopAttrib();

273 }

274

275 static float palette[][3] = {

276 {0, 1, 0},

277 {1, 0, 0},

278 {0, 0, 1},

279 {1, 1, 0},

280 {0, 0, 1},

281 {1, 0, 1}

282 };

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

284

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

286 {

287 if(!node) return;

288

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

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

291

292 glColor3fv(palette[level % pal_size]);

293

294 glVertex3fv(node->aabb.min);

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

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

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

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

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

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

301 glVertex3fv(node->aabb.min);

302

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

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

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

306 glVertex3fv(node->aabb.max);

307 glVertex3fv(node->aabb.max);

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

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

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

311

312 glVertex3fv(node->aabb.min);

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

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

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

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

317 glVertex3fv(node->aabb.max);

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

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

320 }

321

322 bool Scene::build_kdtree()

323 {

324 assert(kdtree == 0);

325

326 const Face *faces = get_face_buffer();

327 int num_faces = get_num_faces();

328

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

330

331 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

332 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

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

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

335

336 free_kdtree(kdtree);

337 kdtree = new KDNode;

338

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

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

341 */

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

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

344

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

346 const Face *face = faces + i;

347

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

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

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

351

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

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

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

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

356 }

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

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

359 }

360 }

361 }

362

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

364 }

365

366 CHECK_AABB(kdtree->aabb);

367

368 // calculate the heuristic for the root

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

370

371 // now proceed splitting the root recursively

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

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

374 return false;

375 }

376

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

378 print_item_counts(kdtree, 0);

379 return true;

380 }

381

382 struct Split {

383 int axis;

384 float pos;

385 float sum_cost;

386 float cost_left, cost_right;

387 };

388

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

390 {

391 Split best_split;

392 best_split.sum_cost = FLT_MAX;

393

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

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

396

397 float splitpt[2];

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

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

400

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

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

403 continue;

404 }

405

406 AABBox aabb_left, aabb_right;

407 aabb_left = aabb_right = node->aabb;

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

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

410

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

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

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

414

415 if(sum_cost < best_split.sum_cost) {

416 best_split.cost_left = left_cost;

417 best_split.cost_right = right_cost;

418 best_split.sum_cost = sum_cost;

419 best_split.pos = splitpt[j];

420 }

421 }

422 }

423

424 assert(split);

425 *split = best_split;

426 split->axis = axis;

427 }

428

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

430 {

431 int opt_max_depth = accel_param[ACCEL_PARAM_MAX_TREE_DEPTH];

432 int opt_max_items = accel_param[ACCEL_PARAM_MAX_NODE_ITEMS];

433

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

435 return true;

436 }

437

438 Split best_split;

439 best_split.axis = -1;

440 best_split.sum_cost = FLT_MAX;

441

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

443 Split split;

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

445

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

447 best_split = split;

448 }

449 }

450

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

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

453 }

454

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

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

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

458 }

459

460 kd->axis = best_split.axis;

461

462 // create the two children

463 KDNode *kdleft, *kdright;

464 kdleft = new KDNode;

465 kdright = new KDNode;

466

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

468

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

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

471

472 kdleft->cost = best_split.cost_left;

473 kdright->cost = best_split.cost_right;

474

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

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

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

478 const Face *face = faces + fidx;

479

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

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

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

483 kdleft->face_idx.push_back(fidx);

484 }

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

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

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

488 kdright->face_idx.push_back(fidx);

489 }

490 }

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

492

493 kd->left = kdleft;

494 kd->right = kdright;

495

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

497 }

498

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

500 {

501 int num_inside = 0;

502 int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];

503 int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];

504

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

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

507

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

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

510 num_inside++;

511 break;

512 }

513 }

514 }

515

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

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

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

519

520 if(dx < 0.0) {

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

522 abort();

523 }

524 if(dy < 0.0) {

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

526 abort();

527 }

528 if(dz < 0.0) {

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

530 abort();

531 }

532

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

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

535 }

536

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

538 return tcost + sarea * num_inside * icost;

539 }

540

541 static void free_kdtree(KDNode *node)

542 {

543 if(node) {

544 free_kdtree(node->left);

545 free_kdtree(node->right);

546 delete node;

547 }

548 }

549

550 int kdtree_depth(const KDNode *node)

551 {

552 if(!node) return 0;

553

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

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

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

557 }

558

559 int kdtree_nodes(const KDNode *node)

560 {

561 if(!node) return 0;

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

563 }

564

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

566 {

567 if(!node) return;

568

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

570 fputs(" ", stdout);

571 }

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

573

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

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

576 }

577

578 #define SGN(x) ((x) >= 0 ? 1 : -1)

579 #define INSIDE(x) (SGN((x) - (splitpos)) == sign)

580 #define OUTSIDE(x) (!INSIDE(x))

581

582 #define LERPV3(res, a, b, t) \

583 do { \

584 (res)[0] = (a)[0] + ((b)[0] - (a)[0]) * (t); \

585 (res)[1] = (a)[1] + ((b)[1] - (a)[1]) * (t); \

586 (res)[2] = (a)[2] + ((b)[2] - (a)[2]) * (t); \

587 } while(0)

588

589 #define NORMALIZE(v) \

590 do { \

591 float mag = (float)sqrt((v)[0] * (v)[0] + (v)[1] * (v)[1] + (v)[2] * (v)[2]); \

592 (v)[0] /= mag; \

593 (v)[1] /= mag; \

594 (v)[2] /= mag; \

595 } while(0)

596

597 static int clip_face(const Face &inface, float splitpos, int axis, int sign, Face *faces)

598 {

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

600

601 std::vector<Vertex> verts;

602 bool clipped = false;

603

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

605 const Vertex *vstart = inface.v + i;

606 const Vertex *vend = inface.v + ((i + 1) % 3);

607

608 float start = vstart->pos[axis];

609 float end = vend->pos[axis];

610

611 if(OUTSIDE(start) && INSIDE(end)) {

612 float t = (splitpos - start) / (end - start);

613

614 Vertex newv;

615 memset(&newv, 0, sizeof newv);

616 LERPV3(newv.pos, vstart->pos, vend->pos, t);

617 LERPV3(newv.normal, vstart->normal, vend->normal, t);

618 LERPV3(newv.tex, vstart->tex, vend->tex, t);

619 NORMALIZE(newv.normal);

620

621 verts.push_back(newv);

622 clipped = true;

623

624 } else if(INSIDE(start) && INSIDE(end)) {

625 verts.push_back(inface.v[i]);

626 } else if(INSIDE(start) && OUTSIDE(end)) {

627 verts.push_back(inface.v[i]);

628

629 float t = (splitpos - start) / (end - start);

630

631 Vertex newv;

632 memset(&newv, 0, sizeof newv);

633 LERPV3(newv.pos, vstart->pos, vend->pos, t);

634 LERPV3(newv.normal, vstart->normal, vend->normal, t);

635 LERPV3(newv.tex, vstart->tex, vend->tex, t);

636 NORMALIZE(newv.normal);

637

638 verts.push_back(newv);

639 clipped = true;

640 }

641 }

642

643 if(!clipped) {

644 return 0;

645 }

646

647 assert(verts.size() < 5);

648 bool quad = verts.size() > 3;

649

650 if(!quad) {

651 faces[0] = inface;

652 faces[0].v[0] = verts[0];

653 faces[0].v[1] = verts[1];

654 faces[0].v[2] = verts[2];

655 return 1;

656 }

657

658 /* calculate triangle areas for both possible splits and pick the one

659 * with the smallest absolute difference to avoid slivers.

660 */

661 float area1, area2;

662

663 area1 = calc_sq_area(verts[0].pos, verts[1].pos, verts[2].pos);

664 area2 = calc_sq_area(verts[0].pos, verts[2].pos, verts[3].pos);

665 float s1diff = fabs(area1 - area2);

666

667 area1 = calc_sq_area(verts[0].pos, verts[1].pos, verts[3].pos);

668 area2 = calc_sq_area(verts[1].pos, verts[2].pos, verts[3].pos);

669 float s2diff = fabs(area1 - area2);

670

671 faces[0] = faces[1] = inface;

672 if(s1diff < s2diff) {

673 faces[0].v[0] = verts[0];

674 faces[0].v[1] = verts[1];

675 faces[0].v[2] = verts[2];

676 faces[1].v[0] = verts[0];

677 faces[1].v[1] = verts[2];

678 faces[1].v[2] = verts[3];

679 } else {

680 faces[0].v[0] = verts[0];

681 faces[0].v[1] = verts[1];

682 faces[0].v[2] = verts[3];

683 faces[1].v[0] = verts[1];

684 faces[1].v[1] = verts[2];

685 faces[1].v[2] = verts[3];

686 }

687 return 2;

688 }

689

690 static float calc_sq_area(const Vector3 &a, const Vector3 &b, const Vector3 &c)

691 {

692 Vector3 v1 = b - a;

693 Vector3 v2 = c - a;

694 return cross(v1, v2).lengthsq();

695 }