nuclear@35: #include <stdlib.h>
John@15: #include <math.h>
nuclear@25: #include <float.h>
nuclear@26: #include <assert.h>
nuclear@35: #include <map>
nuclear@22: #include "scene.h"
nuclear@27: #include "ogl.h"
nuclear@6: 
nuclear@38: #define CHECK_AABB(aabb)	\
nuclear@38: 	assert(aabb.max[0] >= aabb.min[0] && aabb.max[1] >= aabb.min[1] && aabb.max[2] >= aabb.min[2])
nuclear@38: 
nuclear@26: 
nuclear@37: #define MIN(a, b)	((a) < (b) ? (a) : (b))
nuclear@37: #define MAX(a, b)	((a) > (b) ? (a) : (b))
nuclear@37: 
nuclear@37: 
nuclear@53: static int flatten_kdtree(const KDNode *node, KDNodeGPU *kdbuf, int *count);
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@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@44: 	64,	// max tree depth
nuclear@44: 	MAX_NODE_FACES,	// 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; i<m->faces.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: 
nuclear@54: bool Scene::add_light(const Light &lt)
nuclear@54: {
nuclear@54: 	try {
nuclear@54: 		lights.push_back(lt);
nuclear@54: 	}
nuclear@54: 	catch(...) {
nuclear@54: 		return false;
nuclear@54: 	}
nuclear@54: 	return true;
nuclear@54: }
nuclear@54: 
John@14: int Scene::get_num_meshes() const
John@14: {
John@14: 	return (int)meshes.size();
John@14: }
John@14: 
nuclear@54: int Scene::get_num_lights() const
nuclear@54: {
nuclear@54: 	return (int)lights.size();
nuclear@54: }
nuclear@54: 
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; i<meshes.size(); i++) {
nuclear@13: 		num_faces += meshes[i]->faces.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: 
nuclear@35: int Scene::get_num_kdnodes() const
nuclear@35: {
nuclear@35: 	return kdtree_nodes(kdtree);
nuclear@35: }
nuclear@35: 
nuclear@54: Mesh **Scene::get_meshes()
nuclear@54: {
nuclear@54: 	if(meshes.empty()) {
nuclear@54: 		return 0;
nuclear@54: 	}
nuclear@54: 	return &meshes[0];
nuclear@54: }
nuclear@54: 
nuclear@54: const Mesh * const *Scene::get_meshes() const
nuclear@54: {
nuclear@54: 	if(meshes.empty()) {
nuclear@54: 		return 0;
nuclear@54: 	}
nuclear@54: 	return &meshes[0];
nuclear@54: }
nuclear@54: 
nuclear@54: Light *Scene::get_lights()
nuclear@54: {
nuclear@54: 	if(lights.empty()) {
nuclear@54: 		return 0;
nuclear@54: 	}
nuclear@54: 	return &lights[0];
nuclear@54: }
nuclear@54: 
nuclear@54: const Light *Scene::get_lights() const
nuclear@54: {
nuclear@54: 	if(lights.empty()) {
nuclear@54: 		return 0;
nuclear@54: 	}
nuclear@54: 	return &lights[0];
nuclear@54: }
nuclear@54: 
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; i<num_meshes; i++) {
nuclear@24: 		for(size_t j=0; j<meshes[i]->faces.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@35: 	int num_nodes = get_num_kdnodes();
nuclear@35: 	kdbuf = new KDNodeGPU[num_nodes];
nuclear@35: 
nuclear@35: 	int count = 0;
nuclear@35: 
nuclear@35: 	// first arrange the kdnodes into an array (flatten)
nuclear@53: 	flatten_kdtree(kdtree, kdbuf, &count);
nuclear@35: 
nuclear@28: 	return kdbuf;
nuclear@28: }
nuclear@28: 
nuclear@53: static int flatten_kdtree(const KDNode *node, KDNodeGPU *kdbuf, int *count)
nuclear@28: {
nuclear@38: 	const size_t max_node_items = sizeof kdbuf[0].face_idx / sizeof kdbuf[0].face_idx[0];
nuclear@35: 	int idx = (*count)++;
nuclear@29: 
nuclear@35: 	// copy the node
nuclear@35: 	kdbuf[idx].aabb = node->aabb;
nuclear@38: 	kdbuf[idx].num_faces = 0;
nuclear@38: 
nuclear@35: 	for(size_t i=0; i<node->face_idx.size(); i++) {
nuclear@38: 		if(i >= max_node_items) {
nuclear@38: 			fprintf(stderr, "WARNING too many faces per leaf node!\n");
nuclear@38: 			break;
nuclear@38: 		}
nuclear@35: 		kdbuf[idx].face_idx[i] = node->face_idx[i];
nuclear@38: 		kdbuf[idx].num_faces++;
nuclear@28: 	}
nuclear@35: 
nuclear@35: 	// recurse to the left/right (if we're not in a leaf node)
nuclear@35: 	if(node->left) {
nuclear@35: 		assert(node->right);
nuclear@35: 
nuclear@53: 		kdbuf[idx].left = flatten_kdtree(node->left, kdbuf, count);
nuclear@53: 		kdbuf[idx].right = flatten_kdtree(node->right, kdbuf, count);
nuclear@53: 	} else {
nuclear@53: 		kdbuf[idx].left = kdbuf[idx].right = -1;
nuclear@35: 	}
nuclear@28: 
nuclear@53: 	return idx;
nuclear@29: }
nuclear@24: 
nuclear@27: void Scene::draw_kdtree() const
nuclear@27: {
nuclear@27: 	glPushAttrib(GL_ENABLE_BIT);
nuclear@27: 	glDisable(GL_LIGHTING);
nuclear@27: 	glDepthMask(0);
nuclear@27: 
nuclear@27: 	glBegin(GL_LINES);
nuclear@27: 	::draw_kdtree(kdtree, 0);
nuclear@27: 	glEnd();
nuclear@27: 
nuclear@27: 	glDepthMask(1);
nuclear@27: 	glPopAttrib();
nuclear@27: }
nuclear@27: 
nuclear@27: static float palette[][3] = {
nuclear@27: 	{0, 1, 0},
nuclear@27: 	{1, 0, 0},
nuclear@27: 	{0, 0, 1},
nuclear@27: 	{1, 1, 0},
nuclear@27: 	{0, 0, 1},
nuclear@27: 	{1, 0, 1}
nuclear@27: };
nuclear@27: static int pal_size = sizeof palette / sizeof *palette;
nuclear@27: 
nuclear@27: static void draw_kdtree(const KDNode *node, int level)
nuclear@27: {
nuclear@27: 	if(!node) return;
nuclear@27: 
nuclear@27: 	draw_kdtree(node->left, 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@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; i<num_faces; i++) {
nuclear@25: 		const Face *face = faces + i;
nuclear@25: 
nuclear@25: 		// for each vertex of the face ...
nuclear@25: 		for(int j=0; j<3; j++) {
nuclear@25: 			const float *pos = face->v[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@38: 	CHECK_AABB(kdtree->aabb);
nuclear@38: 
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@37: struct Split {
nuclear@37: 	int axis;
nuclear@37: 	float pos;
nuclear@37: 	float sum_cost;
nuclear@37: 	float cost_left, cost_right;
nuclear@37: };
nuclear@37: 
nuclear@37: static void find_best_split(const KDNode *node, int axis, const Face *faces, Split *split)
nuclear@37: {
nuclear@37: 	Split best_split;
nuclear@37: 	best_split.sum_cost = FLT_MAX;
nuclear@37: 
nuclear@37: 	for(size_t i=0; i<node->face_idx.size(); i++) {
nuclear@37: 		const Face *face = faces + node->face_idx[i];
nuclear@37: 
nuclear@37: 		float splitpt[2];
nuclear@37: 		splitpt[0] = MIN(face->v[0].pos[axis], MIN(face->v[1].pos[axis], face->v[2].pos[axis]));
nuclear@37: 		splitpt[1] = MAX(face->v[0].pos[axis], MAX(face->v[1].pos[axis], face->v[2].pos[axis]));
nuclear@37: 
nuclear@37: 		for(int j=0; j<2; j++) {
nuclear@38: 			if(splitpt[j] <= node->aabb.min[axis] || splitpt[j] >= node->aabb.max[axis]) {
nuclear@38: 				continue;
nuclear@38: 			}
nuclear@38: 
nuclear@37: 			AABBox aabb_left, aabb_right;
nuclear@37: 			aabb_left = aabb_right = node->aabb;
nuclear@37: 			aabb_left.max[axis] = splitpt[j];
nuclear@37: 			aabb_right.min[axis] = splitpt[j];
nuclear@37: 
nuclear@37: 			float left_cost = eval_cost(faces, &node->face_idx[0], node->face_idx.size(), aabb_left, axis);
nuclear@37: 			float right_cost = eval_cost(faces, &node->face_idx[0], node->face_idx.size(), aabb_right, axis);
nuclear@37: 			float sum_cost = left_cost + right_cost - accel_param[ACCEL_PARAM_COST_TRAVERSE]; // tcost is added twice
nuclear@37: 
nuclear@37: 			if(sum_cost < best_split.sum_cost) {
nuclear@37: 				best_split.cost_left = left_cost;
nuclear@37: 				best_split.cost_right = right_cost;
nuclear@37: 				best_split.sum_cost = sum_cost;
nuclear@37: 				best_split.pos = splitpt[j];
nuclear@37: 			}
nuclear@37: 		}
nuclear@37: 	}
nuclear@37: 
nuclear@37: 	assert(split);
nuclear@37: 	*split = best_split;
nuclear@37: 	split->axis = axis;
nuclear@37: }
nuclear@37: 
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: 
nuclear@32: 	if(kd->face_idx.empty() || level >= opt_max_depth) {
nuclear@27: 		return true;
nuclear@25: 	}
nuclear@25: 
nuclear@37: 	Split best_split;
nuclear@38: 	best_split.axis = -1;
nuclear@37: 	best_split.sum_cost = FLT_MAX;
nuclear@26: 
nuclear@38: 	for(int i=0; i<3; i++) {
nuclear@37: 		Split split;
nuclear@37: 		find_best_split(kd, i, faces, &split);
nuclear@26: 
nuclear@37: 		if(split.sum_cost < best_split.sum_cost) {
nuclear@37: 			best_split = split;
nuclear@26: 		}
nuclear@26: 	}
nuclear@26: 
nuclear@38: 	if(best_split.axis == -1) {
nuclear@38: 		return true;	// can't split any more, only 0-area splits available
nuclear@38: 	}
nuclear@37: 
nuclear@29: 	//printf("current cost: %f,   best_cost: %f\n", kd->cost, best_sum_cost);
nuclear@37: 	if(best_split.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@38: 	kd->axis = best_split.axis;
nuclear@38: 
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@37: 	kdleft->aabb.max[kd->axis] = best_split.pos;
nuclear@37: 	kdright->aabb.min[kd->axis] = best_split.pos;
nuclear@26: 
nuclear@37: 	kdleft->cost = best_split.cost_left;
nuclear@37: 	kdright->cost = best_split.cost_right;
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; i<kd->face_idx.size(); i++) {
nuclear@32: 		int fidx = kd->face_idx[i];
nuclear@32: 		const Face *face = faces + fidx;
nuclear@26: 
nuclear@37: 		if(face->v[0].pos[kd->axis] < best_split.pos ||
nuclear@37: 				face->v[1].pos[kd->axis] < best_split.pos ||
nuclear@37: 				face->v[2].pos[kd->axis] < best_split.pos) {
nuclear@32: 			kdleft->face_idx.push_back(fidx);
nuclear@26: 		}
nuclear@37: 		if(face->v[0].pos[kd->axis] >= best_split.pos ||
nuclear@37: 				face->v[1].pos[kd->axis] >= best_split.pos ||
nuclear@37: 				face->v[2].pos[kd->axis] >= best_split.pos) {
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; i<num_faces; i++) {
nuclear@32: 		const Face *face = faces + face_idx[i];
nuclear@26: 
nuclear@32: 		for(int j=0; j<3; j++) {
nuclear@32: 			if(face->v[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@38: 	float dx = aabb.max[0] - aabb.min[0];
nuclear@38: 	float dy = aabb.max[1] - aabb.min[1];
nuclear@38: 	float dz = aabb.max[2] - aabb.min[2];
nuclear@38: 
nuclear@38: 	if(dx < 0.0) {
nuclear@38: 		fprintf(stderr, "FOO DX = %f\n", dx);
nuclear@38: 		abort();
nuclear@38: 	}
nuclear@38: 	if(dy < 0.0) {
nuclear@38: 		fprintf(stderr, "FOO DX = %f\n", dy);
nuclear@38: 		abort();
nuclear@38: 	}
nuclear@38: 	if(dz < 0.0) {
nuclear@38: 		fprintf(stderr, "FOO DX = %f\n", dz);
nuclear@38: 		abort();
nuclear@38: 	}
nuclear@38: 
nuclear@38: 	if(dx < 1e-6 || dy < 1e-6 || dz < 1e-6) {
nuclear@27: 		return FLT_MAX;	// heavily penalize 0-area voxels
nuclear@27: 	}
nuclear@27: 
nuclear@38: 	float sarea = 2.0 * (dx + dy + dz);//aabb.calc_surface_area();
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@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; i<level; i++) {
nuclear@27: 		fputs("   ", stdout);
nuclear@27: 	}
nuclear@32: 	printf("- %d (cost: %f)\n", (int)node->face_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: }