1.1 --- a/src/scene.cc	Sat Aug 14 07:00:04 2010 +0100
     1.2 +++ b/src/scene.cc	Tue Aug 17 01:19:43 2010 +0100
     1.3 @@ -1,7 +1,29 @@
     1.4  #include <math.h>
     1.5  #include <float.h>
     1.6 +#include <assert.h>
     1.7  #include "scene.h"
     1.8  
     1.9 +
    1.10 +static int build_kdtree(KDNode *kd);
    1.11 +static float eval_cost(const std::list<const Face*> &faces, const AABBox &aabb, int axis);
    1.12 +static void free_kdtree(KDNode *node);
    1.13 +
    1.14 +
    1.15 +static int accel_param[NUM_ACCEL_PARAMS] = {
    1.16 +	75,	// max tree depth
    1.17 +	0,	// max items per node (0 means ignore limit)
    1.18 +	5,	// estimated traversal cost
    1.19 +	15	// estimated interseciton cost
    1.20 +};
    1.21 +
    1.22 +
    1.23 +void set_accel_param(int p, int v)
    1.24 +{
    1.25 +	assert(p >= 0 && p < NUM_ACCEL_PARAMS);
    1.26 +	accel_param[p] = v;
    1.27 +}
    1.28 +
    1.29 +
    1.30  #define FEQ(a, b)	(fabs((a) - (b)) < 1e-8)
    1.31  bool Face::operator ==(const Face &f) const
    1.32  {
    1.33 @@ -30,6 +52,14 @@
    1.34  	return 2.0f * (area1 + area2 + area3);
    1.35  }
    1.36  
    1.37 +KDNode::KDNode()
    1.38 +{
    1.39 +	axis = 0;
    1.40 +	pt = 0.0;
    1.41 +	left = right = 0;
    1.42 +	num_faces = 0;
    1.43 +}
    1.44 +
    1.45  
    1.46  Scene::Scene()
    1.47  {
    1.48 @@ -124,8 +154,6 @@
    1.49  	return facebuf;
    1.50  }
    1.51  
    1.52 -static int build_kdtree(KDNode *kd);
    1.53 -static void free_kdtree(KDNode *node);
    1.54  
    1.55  void Scene::build_kdtree()
    1.56  {
    1.57 @@ -136,7 +164,6 @@
    1.58  
    1.59  	free_kdtree(kdtree);
    1.60  	kdtree = new KDNode;
    1.61 -	kdtree->left = kdtree->right = 0;
    1.62  
    1.63  	/* Start the construction of the kdtree by adding all faces of the scene
    1.64  	 * to the new root node. At the same time calculate the root's AABB.
    1.65 @@ -163,20 +190,114 @@
    1.66  		}
    1.67  
    1.68  		kdtree->faces.push_back(face);	// add the face
    1.69 +		kdtree->num_faces++;
    1.70  	}
    1.71  
    1.72 +	// calculate the heuristic for the root
    1.73 +	kdtree->cost = eval_cost(kdtree->faces, kdtree->aabb, kdtree->axis);
    1.74 +
    1.75  	// now proceed splitting the root recursively
    1.76  	::build_kdtree(kdtree);
    1.77  }
    1.78  
    1.79  static int build_kdtree(KDNode *kd)
    1.80  {
    1.81 -	if(kd->faces.empty()) {
    1.82 +	int opt_max_items = accel_param[ACCEL_PARAM_MAX_NODE_ITEMS];
    1.83 +	if(kd->num_faces == 0 || (opt_max_items > 1 && kd->num_faces < opt_max_items)) {
    1.84  		return 0;
    1.85  	}
    1.86  
    1.87 -	// XXX cont.
    1.88 -	return -1;
    1.89 +	int axis = (kd->axis + 1) % 3;
    1.90 +
    1.91 +	float best_cost[2], best_sum_cost = FLT_MAX;
    1.92 +	float best_split;
    1.93 +
    1.94 +	std::list<const Face*>::iterator it = kd->faces.begin();
    1.95 +	while(it != kd->faces.end()) {
    1.96 +		const Face *face = *it++;
    1.97 +
    1.98 +		for(int i=0; i<3; i++) {
    1.99 +			AABBox aabb_left, aabb_right;
   1.100 +			const float *split = face->v[i].pos;
   1.101 +
   1.102 +			aabb_left = aabb_right = kd->aabb;
   1.103 +			aabb_left.max[axis] = split[axis];
   1.104 +			aabb_right.min[axis] = split[axis];
   1.105 +
   1.106 +			float left_cost = eval_cost(kd->faces, aabb_left, axis);
   1.107 +			float right_cost = eval_cost(kd->faces, aabb_right, axis);
   1.108 +			float sum_cost = left_cost + right_cost;
   1.109 +
   1.110 +			if(sum_cost < best_sum_cost) {
   1.111 +				best_cost[0] = left_cost;
   1.112 +				best_cost[1] = right_cost;
   1.113 +				best_sum_cost = sum_cost;
   1.114 +				best_split = split[axis];
   1.115 +			}
   1.116 +		}
   1.117 +	}
   1.118 +
   1.119 +	if(best_sum_cost >= kd->cost) {
   1.120 +		return 0;	// stop splitting if it doesn't reduce the cost
   1.121 +	}
   1.122 +	kd->pt = best_split;
   1.123 +
   1.124 +	// create the two children
   1.125 +	KDNode *kdleft, *kdright;
   1.126 +	kdleft = new KDNode;
   1.127 +	kdright = new KDNode;
   1.128 +
   1.129 +	kdleft->axis = kdright->axis = axis;
   1.130 +	kdleft->aabb = kdright->aabb = kd->aabb;
   1.131 +
   1.132 +	kdleft->aabb.max[axis] = best_split;
   1.133 +	kdright->aabb.min[axis] = best_split;
   1.134 +
   1.135 +	kdleft->cost = best_cost[0];
   1.136 +	kdright->cost = best_cost[1];
   1.137 +
   1.138 +	it = kd->faces.begin();
   1.139 +	while(it != kd->faces.end()) {
   1.140 +		const Face *face = *it++;
   1.141 +
   1.142 +		if(face->v[0].pos[axis] < best_split ||
   1.143 +				face->v[1].pos[axis] < best_split ||
   1.144 +				face->v[2].pos[axis] < best_split) {
   1.145 +			kdleft->faces.push_back(face);
   1.146 +			kdleft->num_faces++;
   1.147 +		}
   1.148 +		if(face->v[0].pos[axis] >= best_split ||
   1.149 +				face->v[1].pos[axis] >= best_split ||
   1.150 +				face->v[2].pos[axis] >= best_split) {
   1.151 +			kdright->faces.push_back(face);
   1.152 +			kdright->num_faces++;
   1.153 +		}
   1.154 +	}
   1.155 +
   1.156 +	kd->left = kdleft;
   1.157 +	kd->right = kdright;
   1.158 +	return 0;
   1.159 +}
   1.160 +
   1.161 +static float eval_cost(const std::list<const Face*> &faces, const AABBox &aabb, int axis)
   1.162 +{
   1.163 +	int num_inside = 0;
   1.164 +	int tcost = accel_param[ACCEL_PARAM_COST_TRAVERSE];
   1.165 +	int icost = accel_param[ACCEL_PARAM_COST_INTERSECT];
   1.166 +
   1.167 +	std::list<const Face*>::const_iterator it = faces.begin();
   1.168 +	while(it != faces.end()) {
   1.169 +		const Face *face = *it++;
   1.170 +
   1.171 +		for(int i=0; i<3; i++) {
   1.172 +			if(face->v[i].pos[axis] >= aabb.min[axis] && face->v[i].pos[axis] < aabb.max[axis]) {
   1.173 +				num_inside++;
   1.174 +				break;
   1.175 +			}
   1.176 +		}
   1.177 +	}
   1.178 +
   1.179 +	return tcost + aabb.calc_surface_area() * num_inside * icost;
   1.180  }
   1.181  
   1.182  static void free_kdtree(KDNode *node)