src/kdtree.cc
branchpyrit
changeset 15 a0a3e334744f
parent 14 fc18ac4833f2
child 16 20bceb605f48
equal deleted inserted replaced
14:fc18ac4833f2 15:a0a3e334744f
    21 };
    21 };
    22 
    22 
    23 class SortableShapeList: public vector<SortableShape>
    23 class SortableShapeList: public vector<SortableShape>
    24 {
    24 {
    25 public:
    25 public:
    26 	SortableShapeList(ShapeList &shapes, short axis)
    26 	SortableShapeList(ShapeList *shapes, short axis)
    27 	{
    27 	{
    28 		ShapeList::iterator shape;
    28 		ShapeList::iterator shape;
    29 		for (shape = shapes.begin(); shape != shapes.end(); shape++)
    29 		for (shape = shapes->begin(); shape != shapes->end(); shape++)
    30 			push_back(SortableShape(*shape, axis));
    30 			push_back(SortableShape(*shape, axis));
    31 	};
    31 	};
    32 };
    32 };
    33 
    33 
    34 class SplitPos
    34 class SplitPos
    87 	return nearest_shape;
    87 	return nearest_shape;
    88 }
    88 }
    89 
    89 
    90 void KdNode::subdivide(BBox bbox, int depth)
    90 void KdNode::subdivide(BBox bbox, int depth)
    91 {
    91 {
    92 	if (depth >= 20 || shapes.size() <= 4)
    92 	if (depth >= 20 || shapes->size() <= 4)
    93 	{
    93 	{
    94 		leaf = true;
    94 		setLeaf();
    95 		return;
    95 		return;
    96 	}
    96 	}
    97 
    97 
    98 	// choose split axis
    98 	// choose split axis
    99 	axis = 0;
    99 	axis = 0;
   171 	}
   171 	}
   172 
   172 
   173 	// choose best split pos
   173 	// choose best split pos
   174 	const float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection
   174 	const float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection
   175 	float SAV = 2*(bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node
   175 	float SAV = 2*(bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node
   176 	float cost = SAV * (K + shapes.size()); // initial cost = non-split cost
   176 	float cost = SAV * (K + shapes->size()); // initial cost = non-split cost
   177 	SplitPos *splitpos = NULL;
   177 	SplitPos *splitpos = NULL;
   178 	leaf = true;
   178 	bool leaf = true;
   179 	BBox lbb = bbox;
   179 	BBox lbb = bbox;
   180 	BBox rbb = bbox;
   180 	BBox rbb = bbox;
   181 	for (spl = splitlist.begin(); spl != splitlist.end(); spl++)
   181 	for (spl = splitlist.begin(); spl != splitlist.end(); spl++)
   182 	{
   182 	{
   183 		// calculate SAH cost of this split
   183 		// calculate SAH cost of this split
   194 			splitpos = &*spl;
   194 			splitpos = &*spl;
   195 		}
   195 		}
   196 	}
   196 	}
   197 
   197 
   198 	if (leaf)
   198 	if (leaf)
       
   199 	{
       
   200 		setLeaf();
   199 		return;
   201 		return;
   200 
   202 	}
   201 #if 1
   203 
       
   204 #if 0
   202 // export kd-tree as .obj for visualization
   205 // export kd-tree as .obj for visualization
   203 // this would be hard to reconstruct later
   206 // this would be hard to reconstruct later
   204 	static ofstream F("kdtree.obj");
   207 	static ofstream F("kdtree.obj");
   205 	Vector3 v;
   208 	Vector3 v;
   206 	static int f=0;
   209 	static int f=0;
   287 }
   290 }
   288 
   291 
   289 void KdTree::build()
   292 void KdTree::build()
   290 {
   293 {
   291 	root = new KdNode();
   294 	root = new KdNode();
   292 	root->shapes = shapes;
   295 	root->shapes = &shapes;
   293 	root->subdivide(bbox, 0);
   296 	root->subdivide(bbox, 0);
   294 	built = true;
   297 	built = true;
   295 }
   298 }
   296 
   299 
   297 /* algorithm by Vlastimil Havran, Heuristic Ray Shooting Algorithms, appendix C */
   300 /* algorithm by Vlastimil Havran, Heuristic Ray Shooting Algorithms, appendix C */
   318 	enPt->t = a; /* set the signed distance */
   321 	enPt->t = a; /* set the signed distance */
   319 	enPt->node = NULL;
   322 	enPt->node = NULL;
   320 
   323 
   321 	/* distinguish between internal and external origin */
   324 	/* distinguish between internal and external origin */
   322 	if (a >= 0.0) /* a ray with external origin */
   325 	if (a >= 0.0) /* a ray with external origin */
   323 		enPt->pb = ray.a + ray.dir * a;
   326 		enPt->pb = ray.o + ray.dir * a;
   324 	else /* a ray with internal origin */
   327 	else /* a ray with internal origin */
   325 		enPt->pb = ray.a;
   328 		enPt->pb = ray.o;
   326 
   329 
   327 	/* setup initial exit point in the stack */
   330 	/* setup initial exit point in the stack */
   328 	StackElem *exPt = new StackElem();
   331 	StackElem *exPt = new StackElem();
   329 	exPt->t = b;
   332 	exPt->t = b;
   330 	exPt->pb = ray.a + ray.dir * b;
   333 	exPt->pb = ray.o + ray.dir * b;
   331 	exPt->node = NULL; /* set termination flag */
   334 	exPt->node = NULL; /* set termination flag */
   332 
   335 
   333 	st.push(enPt);
   336 	st.push(enPt);
   334 
   337 
   335 	/* loop, traverse through the whole kd-tree, until an object is intersected or ray leaves the scene */
   338 	/* loop, traverse through the whole kd-tree, until an object is intersected or ray leaves the scene */
   372 			}
   375 			}
   373 
   376 
   374 			/* case P4 or N4 . . . traverse both children */
   377 			/* case P4 or N4 . . . traverse both children */
   375 
   378 
   376 			/* signed distance to the splitting plane */
   379 			/* signed distance to the splitting plane */
   377 			t = (splitVal - ray.a.cell[axis]) / ray.dir.cell[axis];
   380 			t = (splitVal - ray.o.cell[axis]) / ray.dir.cell[axis];
   378 
   381 
   379 			/* setup the new exit point */
   382 			/* setup the new exit point */
   380 			st.push(exPt);
   383 			st.push(exPt);
   381 			exPt = new StackElem();
   384 			exPt = new StackElem();
   382 
   385 
   383 			/* push values onto the stack */
   386 			/* push values onto the stack */
   384 			exPt->t = t;
   387 			exPt->t = t;
   385 			exPt->node = farchild;
   388 			exPt->node = farchild;
   386 			exPt->pb[axis] = splitVal;
   389 			exPt->pb[axis] = splitVal;
   387 			exPt->pb[(axis+1)%3] = ray.a.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
   390 			exPt->pb[(axis+1)%3] = ray.o.cell[(axis+1)%3] + t * ray.dir.cell[(axis+1)%3];
   388 			exPt->pb[(axis+2)%3] = ray.a.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
   391 			exPt->pb[(axis+2)%3] = ray.o.cell[(axis+2)%3] + t * ray.dir.cell[(axis+2)%3];
   389 		} /* while */
   392 		} /* while */
   390 
   393 
   391 		/* current node is the leaf . . . empty or full */
   394 		/* current node is the leaf . . . empty or full */
   392 		/* "intersect ray with each object in the object list, discarding "
   395 		/* "intersect ray with each object in the object list, discarding "
   393 		"those lying before stack[enPt].t or farther than stack[exPt].t" */
   396 		"those lying before stack[enPt].t or farther than stack[exPt].t" */
   394 		Shape *nearest_shape = NULL;
   397 		Shape *nearest_shape = NULL;
   395 		float dist = exPt->t;
   398 		float dist = exPt->t;
   396 		ShapeList::iterator shape;
   399 		ShapeList::iterator shape;
   397 		for (shape = node->shapes.begin(); shape != node->shapes.end(); shape++)
   400 		for (shape = node->shapes->begin(); shape != node->shapes->end(); shape++)
   398 			if (*shape != origin_shape && (*shape)->intersect(ray, dist)
   401 			if (*shape != origin_shape && (*shape)->intersect(ray, dist)
   399 			&& dist >= enPt->t)
   402 			&& dist >= enPt->t)
   400 				nearest_shape = *shape;
   403 				nearest_shape = *shape;
   401 
   404 
   402 		if (nearest_shape)
   405 		if (nearest_shape)
   429 	if (!built)
   432 	if (!built)
   430 		return;
   433 		return;
   431 	if (node == NULL)
   434 	if (node == NULL)
   432 		node = root;
   435 		node = root;
   433 	if (node->isLeaf())
   436 	if (node->isLeaf())
   434 		str << "(leaf: " << node->shapes.size() << " shapes)";
   437 		str << "(leaf: " << node->shapes->size() << " shapes)";
   435 	else
   438 	else
   436 	{
   439 	{
   437 		str << "(split " << (char)('x'+node->getAxis()) << " at " << node->getSplit() << "; L=";
   440 		str << "(split " << (char)('x'+node->getAxis()) << " at " << node->getSplit() << "; L=";
   438 		save(str, node->getLeftChild());
   441 		save(str, node->getLeftChild());
   439 		str << "; R=";
   442 		str << "; R=";