# HG changeset patch
# User John Tsiombikas <nuclear@member.fsf.org>
# Date 1282073700 -3600
# Node ID 8b2f2ad14ae771b7f6235b572ba0ef1415705032
# Parent  c740ae431d51d05baf298843f6089654f5bedbeb
semi-fixed the kdtree construction

diff -r c740ae431d51 -r 8b2f2ad14ae7 src/clray.cc
--- a/src/clray.cc	Tue Aug 17 01:19:43 2010 +0100
+++ b/src/clray.cc	Tue Aug 17 20:35:00 2010 +0100
@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <ctype.h>
 #include <errno.h>
 #ifndef __APPLE__
 #include <GL/glut.h>
@@ -25,7 +26,9 @@
 static float cam_theta, cam_phi = 25.0;
 static float cam_dist = 10.0;
 
-static bool dbg_glrender = false;
+static bool dbg_glrender = true;
+static bool dbg_show_kdtree = false;
+static bool dbg_show_obj = true;
 
 static Scene scn;
 
@@ -36,11 +39,46 @@
 
 	int loaded = 0;
 	for(int i=1; i<argc; i++) {
-		if(!scn.load(argv[i])) {
-			fprintf(stderr, "failed to load scene: %s\n", argv[i]);
-			return false;
+		if(argv[i][0] == '-' && argv[i][2] == 0) {
+			switch(argv[i][1]) {
+			case 'i':
+				if(!argv[++i] || !isdigit(argv[i][0])) {
+					fprintf(stderr, "-i must be followed by the intersection cost\n");
+					return 1;
+				}
+
+				set_accel_param(ACCEL_PARAM_COST_INTERSECT, atoi(argv[i]));
+				break;
+
+			case 't':
+				if(!argv[++i] || !isdigit(argv[i][0])) {
+					fprintf(stderr, "-t must be followed by the traversal cost\n");
+					return 1;
+				}
+
+				set_accel_param(ACCEL_PARAM_COST_TRAVERSE, atoi(argv[i]));
+				break;
+
+			case 'c':
+				if(!argv[++i] || !isdigit(argv[i][0])) {
+					fprintf(stderr, "-c must be followed by the max number of items per leaf node\n");
+					return 1;
+				}
+
+				set_accel_param(ACCEL_PARAM_MAX_NODE_ITEMS, atoi(argv[i]));
+				break;
+
+			default:
+				fprintf(stderr, "unrecognized option: %s\n", argv[i]);
+				return 1;
+			}
+		} else {
+			if(!scn.load(argv[i])) {
+				fprintf(stderr, "failed to load scene: %s\n", argv[i]);
+				return false;
+			}
+			loaded++;
 		}
-		loaded++;
 	}
 
 	if(!loaded) {
@@ -69,6 +107,11 @@
 	}
 	atexit(cleanup);
 
+	if(!scn.build_kdtree()) {
+		return 1;
+	}
+
+
 	/*glGenTextures(1, &tex);
 	glBindTexture(GL_TEXTURE_2D, tex);*/
 	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
@@ -125,7 +168,7 @@
 	if(dbg_glrender) {
 		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 		glLoadMatrixf(inv_mat.m);
-		dbg_render_gl(&scn);
+		dbg_render_gl(&scn, dbg_show_kdtree, dbg_show_obj);
 	} else {
 		glEnable(GL_TEXTURE_2D);
 		glDisable(GL_LIGHTING);
@@ -175,6 +218,20 @@
 		glutPostRedisplay();
 		break;
 
+	case 'k':
+		dbg_show_kdtree = !dbg_show_kdtree;
+		if(dbg_glrender) {
+			glutPostRedisplay();
+		}
+		break;
+
+	case 'o':
+		dbg_show_obj = !dbg_show_obj;
+		if(dbg_glrender) {
+			glutPostRedisplay();
+		}
+		break;
+
 	default:
 		break;
 	}
diff -r c740ae431d51 -r 8b2f2ad14ae7 src/rt.cc
--- a/src/rt.cc	Tue Aug 17 01:19:43 2010 +0100
+++ b/src/rt.cc	Tue Aug 17 20:35:00 2010 +0100
@@ -135,6 +135,7 @@
 	return true;
 }
 
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
 static void dbg_set_gl_material(Material *mat)
 {
 	static Material def_mat = {{0.7, 0.7, 0.7, 1}, {0, 0, 0, 0}, 0, 0, 0};
@@ -143,10 +144,10 @@
 
 	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat->kd);
 	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat->ks);
-	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, mat->spow);
+	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, MIN(mat->spow, 128.0f));
 }
 
-void dbg_render_gl(Scene *scn)
+void dbg_render_gl(Scene *scn, bool show_tree, bool show_obj)
 {
 	glPushAttrib(GL_ENABLE_BIT | GL_TRANSFORM_BIT | GL_LIGHTING_BIT);
 
@@ -170,27 +171,33 @@
 	glLoadIdentity();
 	gluPerspective(45.0, (float)rinf.xsz / (float)rinf.ysz, 0.5, 1000.0);
 
-	Material *materials = scn->get_materials();
+	if(show_obj) {
+		Material *materials = scn->get_materials();
 
-	int num_faces = scn->get_num_faces();
-	int cur_mat = -1;
+		int num_faces = scn->get_num_faces();
+		int cur_mat = -1;
 
-	for(int i=0; i<num_faces; i++) {
-		if(faces[i].matid != cur_mat) {
-			if(cur_mat != -1) {
-				glEnd();
+		for(int i=0; i<num_faces; i++) {
+			if(faces[i].matid != cur_mat) {
+				if(cur_mat != -1) {
+					glEnd();
+				}
+				dbg_set_gl_material(materials ? materials + faces[i].matid : 0);
+				cur_mat = faces[i].matid;
+				glBegin(GL_TRIANGLES);
 			}
-			dbg_set_gl_material(materials ? materials + faces[i].matid : 0);
-			cur_mat = faces[i].matid;
-			glBegin(GL_TRIANGLES);
+
+			for(int j=0; j<3; j++) {
+				glNormal3fv(faces[i].v[j].normal);
+				glVertex3fv(faces[i].v[j].pos);
+			}
 		}
+		glEnd();
+	}
 
-		for(int j=0; j<3; j++) {
-			glNormal3fv(faces[i].v[j].normal);
-			glVertex3fv(faces[i].v[j].pos);
-		}
+	if(show_tree) {
+		scn->draw_kdtree();
 	}
-	glEnd();
 
 	glPopMatrix();
 	glPopAttrib();
diff -r c740ae431d51 -r 8b2f2ad14ae7 src/rt.h
--- a/src/rt.h	Tue Aug 17 01:19:43 2010 +0100
+++ b/src/rt.h	Tue Aug 17 20:35:00 2010 +0100
@@ -8,6 +8,6 @@
 bool render();
 void set_xform(float *matrix, float *invtrans);
 
-void dbg_render_gl(Scene *scn);
+void dbg_render_gl(Scene *scn, bool show_tree = false, bool show_obj = true);
 
 #endif	/* RT_H_ */
diff -r c740ae431d51 -r 8b2f2ad14ae7 src/scene.cc
--- a/src/scene.cc	Tue Aug 17 01:19:43 2010 +0100
+++ b/src/scene.cc	Tue Aug 17 20:35:00 2010 +0100
@@ -2,11 +2,15 @@
 #include <float.h>
 #include <assert.h>
 #include "scene.h"
+#include "ogl.h"
 
 
-static int build_kdtree(KDNode *kd);
-static float eval_cost(const std::list<const Face*> &faces, const AABBox &aabb, int axis);
+static void draw_kdtree(const KDNode *node, int level = 0);
+static bool build_kdtree(KDNode *kd, int level = 0);
+static float eval_cost(const std::list<const Face*> &faces, const AABBox &aabb, int axis, float par_sarea = 1.0);
 static void free_kdtree(KDNode *node);
+static int kdtree_depth(const KDNode *node);
+static void print_item_counts(const KDNode *node, int level);
 
 
 static int accel_param[NUM_ACCEL_PARAMS] = {
@@ -155,13 +159,91 @@
 }
 
 
-void Scene::build_kdtree()
+void Scene::draw_kdtree() const
+{
+	glPushAttrib(GL_ENABLE_BIT);
+	glDisable(GL_LIGHTING);
+	glDepthMask(0);
+
+	glBegin(GL_LINES);
+	::draw_kdtree(kdtree, 0);
+	glEnd();
+
+	glDepthMask(1);
+	glPopAttrib();
+}
+
+static float palette[][3] = {
+	{0, 1, 0},
+	{0, 1, 0},
+	{0, 1, 0},
+	{1, 0, 0},
+	{1, 0, 0},
+	{1, 0, 0},
+	{0, 0, 1},
+	{0, 0, 1},
+	{0, 0, 1},
+	{1, 1, 0},
+	{1, 1, 0},
+	{1, 1, 0},
+	{0, 0, 1},
+	{0, 0, 1},
+	{0, 0, 1},
+	{1, 0, 1},
+	{1, 0, 1},
+	{1, 0, 1}
+};
+static int pal_size = sizeof palette / sizeof *palette;
+
+static void draw_kdtree(const KDNode *node, int level)
+{
+	if(!node) return;
+
+	draw_kdtree(node->left, level + 1);
+	draw_kdtree(node->right, level + 1);
+
+	glColor3fv(palette[level % pal_size]);
+
+	glVertex3fv(node->aabb.min);
+	glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]);
+	glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]);
+	glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]);
+	glVertex3fv(node->aabb.min);
+
+	glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]);
+	glVertex3fv(node->aabb.max);
+	glVertex3fv(node->aabb.max);
+	glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]);
+	glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]);
+	glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]);
+
+	glVertex3fv(node->aabb.min);
+	glVertex3f(node->aabb.min[0], node->aabb.min[1], node->aabb.max[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.min[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.min[1], node->aabb.max[2]);
+	glVertex3f(node->aabb.max[0], node->aabb.max[1], node->aabb.min[2]);
+	glVertex3fv(node->aabb.max);
+	glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.min[2]);
+	glVertex3f(node->aabb.min[0], node->aabb.max[1], node->aabb.max[2]);
+}
+
+bool Scene::build_kdtree()
 {
 	const Face *faces = get_face_buffer();
 	int num_faces = get_num_faces();
 
 	printf("Constructing kd-tree out of %d faces ...\n", num_faces);
 
+	int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];
+	int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];
+	printf("  max items per leaf: %d\n", accel_param[ACCEL_PARAM_MAX_NODE_ITEMS]);
+	printf("  SAH parameters - tcost: %d - icost: %d\n", tcost, icost);
+
 	free_kdtree(kdtree);
 	kdtree = new KDNode;
 
@@ -197,17 +279,27 @@
 	kdtree->cost = eval_cost(kdtree->faces, kdtree->aabb, kdtree->axis);
 
 	// now proceed splitting the root recursively
-	::build_kdtree(kdtree);
+	if(!::build_kdtree(kdtree)) {
+		fprintf(stderr, "failed to build kdtree\n");
+		return false;
+	}
+
+	printf("  tree depth: %d\n", kdtree_depth(kdtree));
+	print_item_counts(kdtree, 0);
+	return true;
 }
 
-static int build_kdtree(KDNode *kd)
+static bool build_kdtree(KDNode *kd, int level)
 {
 	int opt_max_items = accel_param[ACCEL_PARAM_MAX_NODE_ITEMS];
-	if(kd->num_faces == 0 || (opt_max_items > 1 && kd->num_faces < opt_max_items)) {
-		return 0;
+	int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];
+
+	if(kd->num_faces == 0) {
+		return true;
 	}
 
-	int axis = (kd->axis + 1) % 3;
+	int axis = level % 3;
+	//float parent_sa = kd->aabb.calc_surface_area();
 
 	float best_cost[2], best_sum_cost = FLT_MAX;
 	float best_split;
@@ -226,7 +318,7 @@
 
 			float left_cost = eval_cost(kd->faces, aabb_left, axis);
 			float right_cost = eval_cost(kd->faces, aabb_right, axis);
-			float sum_cost = left_cost + right_cost;
+			float sum_cost = left_cost + right_cost - tcost;	// tcost is added twice
 
 			if(sum_cost < best_sum_cost) {
 				best_cost[0] = left_cost;
@@ -237,8 +329,9 @@
 		}
 	}
 
-	if(best_sum_cost >= kd->cost) {
-		return 0;	// stop splitting if it doesn't reduce the cost
+	printf("current cost: %f,   best_cost: %f\n", kd->cost, best_sum_cost);
+	if(best_sum_cost > kd->cost && (opt_max_items == 0 || kd->num_faces <= opt_max_items)) {
+		return true;	// stop splitting if it doesn't reduce the cost
 	}
 	kd->pt = best_split;
 
@@ -247,7 +340,6 @@
 	kdleft = new KDNode;
 	kdright = new KDNode;
 
-	kdleft->axis = kdright->axis = axis;
 	kdleft->aabb = kdright->aabb = kd->aabb;
 
 	kdleft->aabb.max[axis] = best_split;
@@ -256,6 +348,8 @@
 	kdleft->cost = best_cost[0];
 	kdright->cost = best_cost[1];
 
+	//kdleft->axis = kdright->axis = (axis + 1) % 3;
+
 	it = kd->faces.begin();
 	while(it != kd->faces.end()) {
 		const Face *face = *it++;
@@ -273,13 +367,15 @@
 			kdright->num_faces++;
 		}
 	}
+	kd->faces.clear();	// only leaves have faces
 
 	kd->left = kdleft;
 	kd->right = kdright;
-	return 0;
+
+	return build_kdtree(kd->left, level + 1) && build_kdtree(kd->right, level + 1);
 }
 
-static float eval_cost(const std::list<const Face*> &faces, const AABBox &aabb, int axis)
+static float eval_cost(const std::list<const Face*> &faces, const AABBox &aabb, int axis, float par_sarea)
 {
 	int num_inside = 0;
 	int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];
@@ -297,7 +393,12 @@
 		}
 	}
 
-	return tcost + aabb.calc_surface_area() * num_inside * icost;
+	float sarea = aabb.calc_surface_area();
+	if(sarea < 1e-8) {
+		return FLT_MAX;	// heavily penalize 0-area voxels
+	}
+
+	return tcost + (sarea / par_sarea) * num_inside * icost;
 }
 
 static void free_kdtree(KDNode *node)
@@ -308,3 +409,25 @@
 		delete node;
 	}
 }
+
+static int kdtree_depth(const KDNode *node)
+{
+	if(!node) return 0;
+
+	int left = kdtree_depth(node->left);
+	int right = kdtree_depth(node->right);
+	return (left > right ? left : right) + 1;
+}
+
+static void print_item_counts(const KDNode *node, int level)
+{
+	if(!node) return;
+
+	for(int i=0; i<level; i++) {
+		fputs("   ", stdout);
+	}
+	printf("- %d (cost: %f)\n", node->num_faces, node->cost);
+
+	print_item_counts(node->left, level + 1);
+	print_item_counts(node->right, level + 1);
+}
diff -r c740ae431d51 -r 8b2f2ad14ae7 src/scene.h
--- a/src/scene.h	Tue Aug 17 01:19:43 2010 +0100
+++ b/src/scene.h	Tue Aug 17 20:35:00 2010 +0100
@@ -98,7 +98,9 @@
 	bool load(FILE *fp);
 
 	const Face *get_face_buffer() const;
-	void build_kdtree();
+
+	void draw_kdtree() const;
+	bool build_kdtree();
 };
 
 enum {