John@15: #include nuclear@25: #include nuclear@26: #include nuclear@22: #include "scene.h" nuclear@27: #include "ogl.h" nuclear@6: nuclear@26: nuclear@27: static void draw_kdtree(const KDNode *node, int level = 0); nuclear@32: static bool build_kdtree(KDNode *kd, const Face *faces, int level = 0); nuclear@32: static float eval_cost(const Face *faces, const int *face_idx, int num_faces, const AABBox &aabb, int axis); nuclear@26: static void free_kdtree(KDNode *node); nuclear@32: static void kdtree_gpu_flatten(KDNodeGPU *kdbuf, int idx, const KDNode *node); nuclear@27: static void print_item_counts(const KDNode *node, int level); nuclear@26: nuclear@26: nuclear@26: static int accel_param[NUM_ACCEL_PARAMS] = { nuclear@28: 40, // max tree depth nuclear@26: 0, // max items per node (0 means ignore limit) nuclear@26: 5, // estimated traversal cost nuclear@26: 15 // estimated interseciton cost nuclear@26: }; nuclear@26: nuclear@26: nuclear@26: void set_accel_param(int p, int v) nuclear@26: { nuclear@26: assert(p >= 0 && p < NUM_ACCEL_PARAMS); nuclear@26: accel_param[p] = v; nuclear@26: } nuclear@26: nuclear@26: John@15: #define FEQ(a, b) (fabs((a) - (b)) < 1e-8) John@15: bool Face::operator ==(const Face &f) const John@15: { John@15: for(int i=0; i<3; i++) { John@15: for(int j=0; j<3; j++) { John@15: if(!FEQ(v[i].pos[j], f.v[i].pos[j])) { John@15: return false; John@15: } John@15: if(!FEQ(v[i].normal[j], f.v[i].normal[j])) { John@15: return false; John@15: } John@15: } John@15: if(!FEQ(normal[i], f.normal[i])) { John@15: return false; John@15: } John@15: } John@15: return true; John@15: } John@15: nuclear@25: float AABBox::calc_surface_area() const nuclear@25: { nuclear@25: float area1 = (max[0] - min[0]) * (max[1] - min[1]); nuclear@25: float area2 = (max[3] - min[3]) * (max[1] - min[1]); nuclear@25: float area3 = (max[0] - min[0]) * (max[3] - min[3]); nuclear@25: nuclear@25: return 2.0f * (area1 + area2 + area3); nuclear@25: } nuclear@25: nuclear@26: KDNode::KDNode() nuclear@26: { nuclear@26: left = right = 0; nuclear@32: cost = 0.0; nuclear@26: } nuclear@26: nuclear@25: nuclear@24: Scene::Scene() nuclear@24: { nuclear@24: facebuf = 0; nuclear@24: num_faces = -1; nuclear@24: kdtree = 0; nuclear@28: kdbuf = 0; nuclear@24: } nuclear@24: nuclear@24: Scene::~Scene() nuclear@24: { nuclear@24: delete [] facebuf; nuclear@28: delete [] kdbuf; nuclear@28: free_kdtree(kdtree); nuclear@24: } nuclear@24: nuclear@13: bool Scene::add_mesh(Mesh *m) nuclear@13: { nuclear@13: // make sure triangles have material ids nuclear@13: for(size_t i=0; ifaces.size(); i++) { nuclear@13: m->faces[i].matid = m->matid; nuclear@13: } nuclear@24: nuclear@24: try { nuclear@24: meshes.push_back(m); nuclear@24: } nuclear@24: catch(...) { nuclear@24: return false; nuclear@24: } nuclear@24: nuclear@24: // invalidate facebuffer and count nuclear@24: delete [] facebuf; nuclear@24: facebuf = 0; nuclear@24: num_faces = -1; nuclear@24: nuclear@13: return true; nuclear@13: } nuclear@13: John@14: int Scene::get_num_meshes() const John@14: { John@14: return (int)meshes.size(); John@14: } John@14: nuclear@13: int Scene::get_num_faces() const nuclear@13: { nuclear@24: if(num_faces >= 0) { nuclear@24: return num_faces; nuclear@24: } nuclear@24: nuclear@24: num_faces = 0; nuclear@13: for(size_t i=0; ifaces.size(); nuclear@13: } nuclear@13: return num_faces; nuclear@13: } nuclear@13: John@14: int Scene::get_num_materials() const John@14: { John@14: return (int)matlib.size(); John@14: } John@14: John@14: Material *Scene::get_materials() John@14: { John@14: if(matlib.empty()) { John@14: return 0; John@14: } John@14: return &matlib[0]; John@14: } John@14: John@14: const Material *Scene::get_materials() const John@14: { John@14: if(matlib.empty()) { John@14: return 0; John@14: } John@14: return &matlib[0]; John@14: } nuclear@24: nuclear@24: const Face *Scene::get_face_buffer() const nuclear@24: { nuclear@24: if(facebuf) { nuclear@24: return facebuf; nuclear@24: } nuclear@24: nuclear@24: int num_meshes = get_num_meshes(); nuclear@24: nuclear@24: printf("constructing face buffer with %d faces (out of %d meshes)\n", get_num_faces(), num_meshes); nuclear@24: facebuf = new Face[num_faces]; nuclear@24: Face *fptr = facebuf; nuclear@24: nuclear@24: for(int i=0; ifaces.size(); j++) { nuclear@24: *fptr++ = meshes[i]->faces[j]; nuclear@24: } nuclear@24: } nuclear@24: return facebuf; nuclear@24: } nuclear@24: nuclear@28: const KDNodeGPU *Scene::get_kdtree_buffer() const nuclear@28: { nuclear@28: if(kdbuf) { nuclear@28: return kdbuf; nuclear@28: } nuclear@28: nuclear@28: if(!kdtree) { nuclear@28: ((Scene*)this)->build_kdtree(); nuclear@28: } nuclear@28: nuclear@29: int max_nodes = (int)pow(2, kdtree_depth(kdtree)) - 1; nuclear@29: printf("allocating storage for the complete tree (%d)\n", max_nodes); nuclear@28: nuclear@29: kdbuf = new KDNodeGPU[max_nodes + 1]; nuclear@32: kdtree_gpu_flatten(kdbuf, 1, kdtree); nuclear@28: return kdbuf; nuclear@28: } nuclear@28: nuclear@29: static int ipow(int x, int n) nuclear@28: { nuclear@29: assert(n >= 0); nuclear@29: nuclear@29: int res = 1; nuclear@29: for(int i=0; ileft, level + 1); nuclear@27: draw_kdtree(node->right, level + 1); nuclear@27: nuclear@27: glColor3fv(palette[level % pal_size]); nuclear@27: nuclear@27: glVertex3fv(node->aabb.min); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]); nuclear@27: glVertex3fv(node->aabb.min); nuclear@27: nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]); nuclear@27: glVertex3fv(node->aabb.max); nuclear@27: glVertex3fv(node->aabb.max); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]); nuclear@27: nuclear@27: glVertex3fv(node->aabb.min); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]); nuclear@27: glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]); nuclear@27: glVertex3fv(node->aabb.max); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]); nuclear@27: glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]); nuclear@34: /*if(!node->left) return; nuclear@34: nuclear@34: AABBox *bleft = &node->left->aabb; nuclear@34: nuclear@34: int axis = level % 3; nuclear@34: switch(axis) { nuclear@34: case 0: nuclear@34: glVertex3f(bleft->max[0], bleft->min[1], bleft->min[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->min[1], bleft->min[2]); nuclear@34: break; nuclear@34: nuclear@34: case 1: nuclear@34: glVertex3f(bleft->min[0], bleft->min[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->min[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->min[0], bleft->min[1], bleft->max[2]); nuclear@34: break; nuclear@34: nuclear@34: case 2: nuclear@34: glVertex3f(bleft->min[0], bleft->max[1], bleft->min[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->min[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->max[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->min[0], bleft->max[1], bleft->max[2]); nuclear@34: glVertex3f(bleft->min[0], bleft->max[1], bleft->min[2]); nuclear@34: break; nuclear@34: nuclear@34: default: nuclear@34: break; nuclear@34: }*/ nuclear@27: } nuclear@27: nuclear@27: bool Scene::build_kdtree() nuclear@24: { nuclear@29: assert(kdtree == 0); nuclear@29: nuclear@24: const Face *faces = get_face_buffer(); nuclear@24: int num_faces = get_num_faces(); nuclear@24: nuclear@25: printf("Constructing kd-tree out of %d faces ...\n", num_faces); nuclear@25: nuclear@27: int icost = accel_param[ACCEL_PARAM_COST_INTERSECT]; nuclear@27: int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE]; nuclear@27: printf(" max items per leaf: %d\n", accel_param[ACCEL_PARAM_MAX_NODE_ITEMS]); nuclear@27: printf(" SAH parameters - tcost: %d - icost: %d\n", tcost, icost); nuclear@27: nuclear@25: free_kdtree(kdtree); nuclear@25: kdtree = new KDNode; nuclear@25: nuclear@25: /* Start the construction of the kdtree by adding all faces of the scene nuclear@25: * to the new root node. At the same time calculate the root's AABB. nuclear@25: */ nuclear@25: kdtree->aabb.min[0] = kdtree->aabb.min[1] = kdtree->aabb.min[2] = FLT_MAX; nuclear@25: kdtree->aabb.max[0] = kdtree->aabb.max[1] = kdtree->aabb.max[2] = -FLT_MAX; nuclear@25: nuclear@24: for(int i=0; iv[j].pos; nuclear@25: nuclear@25: // for each element (xyz) of the position vector ... nuclear@25: for(int k=0; k<3; k++) { nuclear@25: if(pos[k] < kdtree->aabb.min[k]) { nuclear@25: kdtree->aabb.min[k] = pos[k]; nuclear@25: } nuclear@25: if(pos[k] > kdtree->aabb.max[k]) { nuclear@25: kdtree->aabb.max[k] = pos[k]; nuclear@25: } nuclear@25: } nuclear@25: } nuclear@25: nuclear@32: kdtree->face_idx.push_back(i); // add the face nuclear@24: } nuclear@24: nuclear@26: // calculate the heuristic for the root nuclear@32: kdtree->cost = eval_cost(faces, &kdtree->face_idx[0], kdtree->face_idx.size(), kdtree->aabb, 0); nuclear@26: nuclear@25: // now proceed splitting the root recursively nuclear@32: if(!::build_kdtree(kdtree, faces)) { nuclear@27: fprintf(stderr, "failed to build kdtree\n"); nuclear@27: return false; nuclear@27: } nuclear@27: nuclear@27: printf(" tree depth: %d\n", kdtree_depth(kdtree)); nuclear@27: print_item_counts(kdtree, 0); nuclear@27: return true; nuclear@24: } nuclear@24: nuclear@32: static bool build_kdtree(KDNode *kd, const Face *faces, int level) nuclear@24: { nuclear@28: int opt_max_depth = accel_param[ACCEL_PARAM_MAX_TREE_DEPTH]; nuclear@26: int opt_max_items = accel_param[ACCEL_PARAM_MAX_NODE_ITEMS]; nuclear@27: int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE]; nuclear@27: nuclear@32: if(kd->face_idx.empty() || level >= opt_max_depth) { nuclear@27: return true; nuclear@25: } nuclear@25: nuclear@27: int axis = level % 3; nuclear@26: nuclear@26: float best_cost[2], best_sum_cost = FLT_MAX; nuclear@26: float best_split; nuclear@26: nuclear@32: for(size_t i=0; iface_idx.size(); i++) { nuclear@32: const Face *face = faces + kd->face_idx[i]; nuclear@26: nuclear@32: for(int j=0; j<3; j++) { nuclear@26: AABBox aabb_left, aabb_right; nuclear@32: const float *split = face->v[j].pos; nuclear@26: nuclear@26: aabb_left = aabb_right = kd->aabb; nuclear@26: aabb_left.max[axis] = split[axis]; nuclear@26: aabb_right.min[axis] = split[axis]; nuclear@26: nuclear@32: float left_cost = eval_cost(faces, &kd->face_idx[0], kd->face_idx.size(), aabb_left, axis); nuclear@32: float right_cost = eval_cost(faces, &kd->face_idx[0], kd->face_idx.size(), aabb_right, axis); nuclear@27: float sum_cost = left_cost + right_cost - tcost; // tcost is added twice nuclear@26: nuclear@26: if(sum_cost < best_sum_cost) { nuclear@26: best_cost[0] = left_cost; nuclear@26: best_cost[1] = right_cost; nuclear@26: best_sum_cost = sum_cost; nuclear@26: best_split = split[axis]; nuclear@26: } nuclear@26: } nuclear@26: } nuclear@26: nuclear@29: //printf("current cost: %f, best_cost: %f\n", kd->cost, best_sum_cost); nuclear@32: if(best_sum_cost > kd->cost && (opt_max_items == 0 || (int)kd->face_idx.size() <= opt_max_items)) { nuclear@27: return true; // stop splitting if it doesn't reduce the cost nuclear@26: } nuclear@26: nuclear@26: // create the two children nuclear@26: KDNode *kdleft, *kdright; nuclear@26: kdleft = new KDNode; nuclear@26: kdright = new KDNode; nuclear@26: nuclear@26: kdleft->aabb = kdright->aabb = kd->aabb; nuclear@26: nuclear@26: kdleft->aabb.max[axis] = best_split; nuclear@26: kdright->aabb.min[axis] = best_split; nuclear@26: nuclear@26: kdleft->cost = best_cost[0]; nuclear@26: kdright->cost = best_cost[1]; nuclear@26: nuclear@34: // TODO would it be much better if we actually split faces that straddle the splitting plane? nuclear@32: for(size_t i=0; iface_idx.size(); i++) { nuclear@32: int fidx = kd->face_idx[i]; nuclear@32: const Face *face = faces + fidx; nuclear@26: nuclear@26: if(face->v[0].pos[axis] < best_split || nuclear@26: face->v[1].pos[axis] < best_split || nuclear@26: face->v[2].pos[axis] < best_split) { nuclear@32: kdleft->face_idx.push_back(fidx); nuclear@26: } nuclear@26: if(face->v[0].pos[axis] >= best_split || nuclear@26: face->v[1].pos[axis] >= best_split || nuclear@26: face->v[2].pos[axis] >= best_split) { nuclear@32: kdright->face_idx.push_back(fidx); nuclear@26: } nuclear@26: } nuclear@32: kd->face_idx.clear(); // only leaves have faces nuclear@26: nuclear@26: kd->left = kdleft; nuclear@26: kd->right = kdright; nuclear@27: nuclear@32: return build_kdtree(kd->left, faces, level + 1) && build_kdtree(kd->right, faces, level + 1); nuclear@26: } nuclear@26: nuclear@32: static float eval_cost(const Face *faces, const int *face_idx, int num_faces, const AABBox &aabb, int axis) nuclear@26: { nuclear@26: int num_inside = 0; nuclear@26: int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE]; nuclear@26: int icost = accel_param[ACCEL_PARAM_COST_INTERSECT]; nuclear@26: nuclear@32: for(int i=0; iv[j].pos[axis] >= aabb.min[axis] && face->v[j].pos[axis] < aabb.max[axis]) { nuclear@26: num_inside++; nuclear@26: break; nuclear@26: } nuclear@26: } nuclear@26: } nuclear@26: nuclear@27: float sarea = aabb.calc_surface_area(); nuclear@34: if(sarea < 1e-6) { nuclear@27: return FLT_MAX; // heavily penalize 0-area voxels nuclear@27: } nuclear@27: nuclear@32: return tcost + sarea * num_inside * icost; nuclear@24: } nuclear@25: nuclear@25: static void free_kdtree(KDNode *node) nuclear@25: { nuclear@25: if(node) { nuclear@25: free_kdtree(node->left); nuclear@25: free_kdtree(node->right); nuclear@25: delete node; nuclear@25: } nuclear@25: } nuclear@27: nuclear@28: int kdtree_depth(const KDNode *node) nuclear@27: { nuclear@27: if(!node) return 0; nuclear@27: nuclear@27: int left = kdtree_depth(node->left); nuclear@27: int right = kdtree_depth(node->right); nuclear@27: return (left > right ? left : right) + 1; nuclear@27: } nuclear@27: nuclear@28: int kdtree_nodes(const KDNode *node) nuclear@28: { nuclear@28: if(!node) return 0; nuclear@28: return kdtree_nodes(node->left) + kdtree_nodes(node->right) + 1; nuclear@28: } nuclear@28: nuclear@28: #define MAX_FACES (sizeof dest->face_idx / sizeof *dest->face_idx) nuclear@32: static void kdtree_gpu_flatten(KDNodeGPU *kdbuf, int idx, const KDNode *node) nuclear@28: { nuclear@28: KDNodeGPU *dest = kdbuf + idx; nuclear@28: nuclear@28: dest->aabb = node->aabb; nuclear@28: dest->num_faces = 0; nuclear@28: nuclear@32: for(size_t i=0; iface_idx.size(); i++) { nuclear@28: if(dest->num_faces >= (int)MAX_FACES) { nuclear@28: fprintf(stderr, "kdtree_gpu_flatten WARNING: more than %d faces in node, skipping!\n", (int)MAX_FACES); nuclear@28: break; nuclear@28: } nuclear@32: dest->face_idx[dest->num_faces++] = node->face_idx[i]; nuclear@28: } nuclear@28: nuclear@28: if(node->left) { nuclear@28: assert(node->right); nuclear@31: assert(!dest->num_faces); nuclear@31: nuclear@31: dest->num_faces = -1; nuclear@31: nuclear@32: kdtree_gpu_flatten(kdbuf, idx * 2, node->left); nuclear@32: kdtree_gpu_flatten(kdbuf, idx * 2 + 1, node->right); nuclear@28: } nuclear@28: } nuclear@28: nuclear@27: static void print_item_counts(const KDNode *node, int level) nuclear@27: { nuclear@27: if(!node) return; nuclear@27: nuclear@30: for(int i=0; iface_idx.size(), node->cost); nuclear@27: nuclear@27: print_item_counts(node->left, level + 1); nuclear@27: print_item_counts(node->right, level + 1); nuclear@27: }