Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/VssTree.cpp @ 1359

Revision 1359, 47.2 KB checked in by mattausch, 18 years ago (diff)
worked on global object sorting

Rev	Line
[372]	1
	2	// ================================================================
	3	// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
	4	// ****************************************************************
	5	/**
	6	The KD tree based LSDS
	7	*/
	8	// Initial coding by
	9	/**
	10	@author Jiri Bittner
	11	*/
	12
	13	// Standard headers
	14	#include <stack>
	15	#include <queue>
	16	#include <algorithm>
	17	#include <fstream>
	18	#include <string>
	19
	20	#include "VssTree.h"
	21
	22	#include "Environment.h"
	23	#include "VssRay.h"
[382]	24	#include "Intersectable.h"
[401]	25	#include "Ray.h"
[372]	26
[863]	27	namespace GtpVisibilityPreprocessor {
[860]	28
[438]	29	#define DEBUG_SPLIT_COST 0
[372]	30
	31	// Static variables
	32	int
	33	VssTreeLeaf::mailID = 0;
	34
[382]	35	inline void
	36	AddObject2Pvs(Intersectable *object,
[434]	37	const int side,
	38	int &pvsBack,
	39	int &pvsFront)
[382]	40	{
[434]	41
	42	if (!object)
	43	return;
	44
	45	if (side <= 0) {
	46	if (!object->Mailed() && !object->Mailed(2)) {
	47	pvsBack++;
	48	if (object->Mailed(1))
	49	object->Mail(2);
	50	else
	51	object->Mail();
[382]	52	}
[434]	53	}
	54
	55	if (side >= 0) {
	56	if (!object->Mailed(1) && !object->Mailed(2)) {
	57	pvsFront++;
	58	if (object->Mailed())
	59	object->Mail(2);
	60	else
	61	object->Mail(1);
[382]	62	}
[434]	63	}
[382]	64	}
	65
[372]	66	// Constructor
	67	VssTree::VssTree()
	68	{
[1004]	69	Environment::GetSingleton()->GetIntValue("VssTree.maxDepth", termMaxDepth);
	70	Environment::GetSingleton()->GetIntValue("VssTree.minPvs", termMinPvs);
	71	Environment::GetSingleton()->GetIntValue("VssTree.minRays", termMinRays);
	72	Environment::GetSingleton()->GetFloatValue("VssTree.maxRayContribution", termMaxRayContribution);
	73	Environment::GetSingleton()->GetFloatValue("VssTree.maxCostRatio", termMaxCostRatio);
[372]	74
[1004]	75	Environment::GetSingleton()->GetFloatValue("VssTree.minSize", termMinSize);
[372]	76	termMinSize = sqr(termMinSize);
[382]	77
[1004]	78	Environment::GetSingleton()->GetFloatValue("VssTree.refDirBoxMaxSize", refDirBoxMaxSize);
[372]	79	refDirBoxMaxSize = sqr(refDirBoxMaxSize);
	80
[1004]	81	Environment::GetSingleton()->GetFloatValue("VssTree.epsilon", epsilon);
	82	Environment::GetSingleton()->GetFloatValue("VssTree.ct_div_ci", ct_div_ci);
[382]	83
[1004]	84	Environment::GetSingleton()->GetFloatValue("VssTree.maxTotalMemory", maxTotalMemory);
	85	Environment::GetSingleton()->GetFloatValue("VssTree.maxStaticMemory", maxStaticMemory);
[372]	86
	87
	88
	89
	90	float refDirAngle;
[1004]	91	Environment::GetSingleton()->GetFloatValue("VssTree.refDirAngle", refDirAngle);
[434]	92
[1004]	93	Environment::GetSingleton()->GetIntValue("VssTree.accessTimeThreshold", accessTimeThreshold);
[372]	94	//= 1000;
[1004]	95	Environment::GetSingleton()->GetIntValue("VssTree.minCollapseDepth", minCollapseDepth);
[372]	96	// int minCollapseDepth = 4;
	97
	98	// pRefDirThresh = cos(0.5M_PI - M_PIrefDirAngle/180.0);
[434]	99	// cosRefDir = cos(M_PI*refDirAngle/180.0);
	100	// sinRefDir = sin(M_PI*refDirAngle/180.0);
[372]	101
	102
	103	// split type
[434]	104	char sname[128];
[1004]	105	Environment::GetSingleton()->GetStringValue("VssTree.splitType", sname);
[372]	106	string name(sname);
	107
	108	if (name.compare("regular") == 0)
	109	splitType = ESplitRegular;
	110	else
[386]	111	if (name.compare("heuristic") == 0)
	112	splitType = ESplitHeuristic;
[434]	113	else
	114	if (name.compare("hybrid") == 0)
	115	splitType = ESplitHybrid;
	116	else {
	117	cerr<<"Invalid VssTree split type "<<name<<endl;
	118	exit(1);
	119	}
[372]	120
[1004]	121	Environment::GetSingleton()->GetBoolValue("VssTree.randomize", randomize);
	122	Environment::GetSingleton()->GetBoolValue("VssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
	123	Environment::GetSingleton()->GetBoolValue("VssTree.useRss", mUseRss);
[438]	124
[1004]	125	Environment::GetSingleton()->GetBoolValue("VssTree.interleaveDirSplits", mInterleaveDirSplits);
	126	Environment::GetSingleton()->GetIntValue("VssTree.dirSplitDepth", mDirSplitDepth);
[438]	127
[372]	128	root = NULL;
	129
[382]	130	splitCandidates = new vector<SortableEntry>;
[372]	131	}
	132
	133
	134	VssTree::~VssTree()
	135	{
	136	if (root)
	137	delete root;
	138	}
	139
	140
	141
	142
	143	void
	144	VssStatistics::Print(ostream &app) const
	145	{
	146	app << "===== VssTree statistics ===============\n";
	147
[382]	148	app << "#N_RAYS ( Number of rays )\n"
[372]	149	<< rays <<endl;
[382]	150
[434]	151	app << "#N_INITPVS ( Initial PVS size )\n"
[382]	152	<< initialPvsSize <<endl;
[372]	153
	154	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
	155
	156	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
	157
[667]	158	app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
[372]	159	for (int i=0; i<7; i++)
	160	app << splits[i] <<" ";
	161	app <<endl;
	162
	163	app << "#N_RAYREFS ( Number of rayRefs )\n" <<
	164	rayRefs << "\n";
	165
	166	app << "#N_RAYRAYREFS ( Number of rayRefs / ray )\n" <<
	167	rayRefs/(double)rays << "\n";
	168
	169	app << "#N_LEAFRAYREFS ( Number of rayRefs / leaf )\n" <<
	170	rayRefs/(double)Leaves() << "\n";
	171
	172	app << "#N_MAXRAYREFS ( Max number of rayRefs / leaf )\n" <<
	173	maxRayRefs << "\n";
	174
	175
	176	// app << setprecision(4);
	177
	178	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
	179	maxDepthNodes*100/(double)Leaves()<<endl;
	180
[551]	181	app << "#N_PMINPVSLEAVES ( Percentage of leaves with minPvs )\n"<<
[403]	182	minPvsNodes*100/(double)Leaves()<<endl;
	183
[551]	184	app << "#N_PMINRAYSLEAVES ( Percentage of leaves with minRays )\n"<<
[403]	185	minRaysNodes*100/(double)Leaves()<<endl;
[382]	186
[434]	187	app << "#N_PMINSIZELEAVES ( Percentage of leaves with minSize )\n"<<
[382]	188	minSizeNodes*100/(double)Leaves()<<endl;
[372]	189
[434]	190	app << "#N_PMAXRAYCONTRIBLEAVES ( Percentage of leaves with maximal ray contribution )\n"<<
[382]	191	maxRayContribNodes*100/(double)Leaves()<<endl;
	192
[434]	193	app << "#N_PMAXCOSTRATIOLEAVES ( Percentage of leaves with max cost ratio )\n"<<
[427]	194	maxCostRatioNodes*100/(double)Leaves()<<endl;
	195
[372]	196	app << "#N_ADDED_RAYREFS (Number of dynamically added ray references )\n"<<
	197	addedRayRefs<<endl;
	198
	199	app << "#N_REMOVED_RAYREFS (Number of dynamically removed ray references )\n"<<
	200	removedRayRefs<<endl;
	201
	202	// app << setprecision(4);
	203
	204	app << "#N_CTIME ( Construction time [s] )\n"
	205	<< Time() << " \n";
	206
	207	app << "===== END OF VssTree statistics ==========\n";
	208
	209	}
	210
	211
[395]	212	void
[401]	213	VssTreeLeaf::UpdatePvsSize()
[395]	214	{
[434]	215	if (!mValidPvs) {
	216	Intersectable::NewMail();
	217	int pvsSize = 0;
	218	for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
	219	ri != rays.end();
	220	ri++)
	221	if ((*ri).mRay->IsActive()) {
	222	Intersectable *object;
[395]	223	#if BIDIRECTIONAL_RAY
[434]	224	object = (*ri).mRay->mOriginObject;
	225	if (object && !object->Mailed()) {
	226	pvsSize++;
	227	object->Mail();
	228	}
[395]	229	#endif
[434]	230	object = (*ri).mRay->mTerminationObject;
	231	if (object && !object->Mailed()) {
	232	pvsSize++;
	233	object->Mail();
	234	}
	235	}
	236	mPvsSize = pvsSize;
	237	mValidPvs = true;
[438]	238
[434]	239	}
[395]	240	}
[372]	241
[427]	242	bool
	243	VssTree::ClipRay(
[434]	244	VssTreeNode::RayInfo &rayInfo,
	245	const AxisAlignedBox3 &box
	246	)
[427]	247	{
[434]	248	float tmin, tmax;
	249	static Ray ray;
	250	ray.Init(rayInfo.mRay->GetOrigin(), rayInfo.mRay->GetDir(), Ray::LINE_SEGMENT);
	251	box.ComputeMinMaxT(ray, &tmin, &tmax);
	252	if (tmin >= tmax)
	253	return false;
[427]	254
[434]	255	if (tmin > 0.0f)
	256	rayInfo.SetMinT(tmin);
[372]	257
[434]	258	if (tmax < 1.0f)
	259	rayInfo.SetMaxT(tmax);
[427]	260
[434]	261	// vssRay->SetupEndPoints(
	262	// origin,
	263	// termination
	264	// );
	265	return true;
[427]	266	}
	267
	268
[372]	269	void
	270	VssTree::Construct(
[434]	271	VssRayContainer &rays,
	272	AxisAlignedBox3 *forcedBoundingBox
	273	)
[372]	274	{
	275	stat.Start();
	276
[434]	277	maxMemory = maxStaticMemory;
[372]	278
	279	if (root)
	280	delete root;
[434]	281
[1359]	282	root = new VssTreeLeaf(NULL, (int)rays.size());
[372]	283	// first construct a leaf that will get subdivide
	284	VssTreeLeaf leaf = (VssTreeLeaf ) root;
	285
	286	stat.nodes = 1;
	287
	288	bbox.Initialize();
	289	dirBBox.Initialize();
[434]	290
	291	if (mUseRss)
	292	forcedBoundingBox = NULL;
	293
[372]	294	for(VssRayContainer::const_iterator ri = rays.begin();
	295	ri != rays.end();
	296	ri++) {
[427]	297
[434]	298	VssTreeNode::RayInfo info(*ri);
	299	if (forcedBoundingBox)
	300	if (!ClipRay(info, *forcedBoundingBox))
	301	continue;
	302	leaf->AddRay(info);
	303
[372]	304	bbox.Include((*ri)->GetOrigin());
	305	bbox.Include((*ri)->GetTermination());
	306
[382]	307
[434]	308	dirBBox.Include(Vector3(
	309	(*ri)->GetDirParametrization(0),
	310	(*ri)->GetDirParametrization(1),
	311	0
	312	)
	313	);
	314	}
[382]	315
	316
[434]	317	if ( forcedBoundingBox )
	318	bbox = *forcedBoundingBox;
[427]	319
[434]	320	cout<<"Bbox = "<<bbox<<endl;
	321	cout<<"Dirr Bbox = "<<dirBBox<<endl;
[382]	322
[1359]	323	stat.rays = (int)leaf->rays.size();
[434]	324	leaf->UpdatePvsSize();
[438]	325	leaf->ComputeEntropyImportance();
	326
[395]	327	stat.initialPvsSize = leaf->GetPvsSize();
[372]	328	// Subdivide();
	329	root = Subdivide(TraversalData(leaf, bbox, 0));
	330
	331	if (splitCandidates) {
	332	// force realease of this vector
	333	delete splitCandidates;
[382]	334	splitCandidates = new vector<SortableEntry>;
[372]	335	}
	336
	337	stat.Stop();
	338
[434]	339	stat.Print(cout);
	340	cout<<"#Total memory="<<GetMemUsage()<<endl;
[382]	341
[372]	342	}
	343
[427]	344	int
	345	VssTree::UpdateSubdivision()
	346	{
[434]	347	priority_queue<TraversalData> tStack;
[427]	348	// stack<TraversalData> tStack;
	349
[434]	350	tStack.push(TraversalData(root, bbox, 0));
[427]	351
	352	AxisAlignedBox3 backBox;
	353	AxisAlignedBox3 frontBox;
[372]	354
[434]	355	maxMemory = maxTotalMemory;
	356	int subdivided = 0;
	357	int lastMem = 0;
[427]	358	while (!tStack.empty()) {
	359
[434]	360	float mem = GetMemUsage();
[427]	361
[434]	362	if ( lastMem/10 != ((int)mem)/10) {
	363	cout<<mem<<" MB"<<endl;
	364	}
	365	lastMem = (int)mem;
[427]	366
[434]	367	if ( mem > maxMemory ) {
[427]	368	// count statistics on unprocessed leafs
	369	while (!tStack.empty()) {
[434]	370	// EvaluateLeafStats(tStack.top());
	371	tStack.pop();
[427]	372	}
	373	break;
	374	}
	375
	376	TraversalData data = tStack.top();
	377	tStack.pop();
[382]	378
[434]	379	if (data.node->IsLeaf()) {
	380	VssTreeNode node = SubdivideNode((VssTreeLeaf ) data.node,
	381	data.bbox,
	382	backBox,
	383	frontBox
	384	);
	385	if (!node->IsLeaf()) {
	386	subdivided++;
	387	VssTreeInterior interior = (VssTreeInterior ) node;
	388	// push the children on the stack
	389	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
	390	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
	391	} else {
	392	// EvaluateLeafStats(data);
	393	}
	394	} else {
	395	VssTreeInterior interior = (VssTreeInterior ) data.node;
	396	tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
	397	tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
	398	}
[427]	399	}
[434]	400	return subdivided;
[427]	401	}
	402
	403
[372]	404	VssTreeNode *
	405	VssTree::Subdivide(const TraversalData &tdata)
	406	{
	407	VssTreeNode *result = NULL;
	408
	409	priority_queue<TraversalData> tStack;
	410	// stack<TraversalData> tStack;
	411
	412	tStack.push(tdata);
	413
	414	AxisAlignedBox3 backBox;
	415	AxisAlignedBox3 frontBox;
	416
	417
[434]	418	int lastMem = 0;
[372]	419	while (!tStack.empty()) {
	420
[434]	421	float mem = GetMemUsage();
[386]	422
[434]	423	if ( lastMem/10 != ((int)mem)/10) {
	424	cout<<mem<<" MB"<<endl;
	425	}
	426	lastMem = (int)mem;
[386]	427
[434]	428	if ( mem > maxMemory ) {
[372]	429	// count statistics on unprocessed leafs
	430	while (!tStack.empty()) {
[434]	431	EvaluateLeafStats(tStack.top());
	432	tStack.pop();
[372]	433	}
	434	break;
	435	}
	436
	437	TraversalData data = tStack.top();
	438	tStack.pop();
	439
	440	VssTreeNode node = SubdivideNode((VssTreeLeaf ) data.node,
[434]	441	data.bbox,
	442	backBox,
	443	frontBox
	444	);
[372]	445	if (result == NULL)
	446	result = node;
	447
	448	if (!node->IsLeaf()) {
	449
	450	VssTreeInterior interior = (VssTreeInterior ) node;
	451	// push the children on the stack
	452	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
	453	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
	454
	455	} else {
	456	EvaluateLeafStats(data);
	457	}
	458	}
	459
	460	return result;
	461	}
	462
	463
	464	// returns selected plane for subdivision
	465	int
	466	VssTree::SelectPlane(
[434]	467	VssTreeLeaf *leaf,
	468	const AxisAlignedBox3 &box,
	469	float &position,
	470	int &raysBack,
	471	int &raysFront,
	472	int &pvsBack,
	473	int &pvsFront
	474	)
[372]	475	{
	476
	477	int minDirDepth = 6;
[382]	478	int axis;
[434]	479	float costRatio;
[382]	480
[434]	481	costRatio = BestCostRatio(leaf,
	482	axis,
	483	position,
	484	raysBack,
	485	raysFront,
	486	pvsBack,
	487	pvsFront
	488	);
[438]	489	#if DEBUG_SPLIT_COST
[435]	490	cout<<axis<<" r="<<costRatio<<endl;
[438]	491	#endif
[434]	492	if (costRatio > termMaxCostRatio) {
	493	// cout<<"Too big cost ratio "<<costRatio<<endl;
	494	stat.maxCostRatioNodes++;
	495	return -1;
[382]	496	}
[434]	497
[382]	498	#if 0
[434]	499	cout<<
	500	"pvs="<<leaf->mPvsSize<<
	501	" rays="<<leaf->rays.size()<<
	502	" rc="<<leaf->GetAvgRayContribution()<<
	503	" axis="<<axis<<endl;
[382]	504	#endif
	505
[434]	506	return axis;
[382]	507	}
	508
	509
[434]	510	float
	511	VssTree::GetCostRatio(
	512	VssTreeLeaf *leaf,
	513	const int axis,
	514	const float position,
	515	const int raysBack,
	516	const int raysFront,
	517	const int pvsBack,
	518	const int pvsFront
	519	)
	520	{
	521	bool costImportance = true;
	522
	523	float ratio;
	524	AxisAlignedBox3 box;
	525	float minBox, maxBox;
	526
	527	if (axis < 3) {
	528	box = GetBBox(leaf);
	529	minBox = box.Min(axis);
	530	maxBox = box.Max(axis);
	531	} else {
	532	box = GetDirBBox(leaf);
	533	minBox = box.Min(axis-3);
	534	maxBox = box.Max(axis-3);
	535	}
	536
	537	float sizeBox = maxBox - minBox;
	538
	539	int pvsSize = leaf->GetPvsSize();
	540
	541	if (!costImportance) {
	542	// float sum = raysBack(position - minBox) + raysFront(maxBox - position);
	543	float sum = pvsBack(position - minBox) + pvsFront(maxBox - position);
	544	float newCost = ct_div_ci + sum/sizeBox;
[1297]	545	float oldCost = (float)pvsSize;
[434]	546	ratio = newCost/oldCost;
[1297]	547
[434]	548	} else {
	549	// importance based cost
	550	#if 0
	551	float newContrib =
	552	((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
	553	(maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
	554
	555	// float newContrib =
	556	// sqr(pvsBack/(raysBack + Limits::Small)) +
	557	// sqr(pvsFront/(raysFront + Limits::Small));
	558	float oldContrib = sqr(leaf->GetAvgRayContribution());
	559	ratio = oldContrib/newContrib;
	560	#else
	561	#if 1
[1297]	562	float newCost = float(raysBackpvsBack + raysFrontpvsFront);
[485]	563	float oldCost = (float)leaf->rays.size()*pvsSize;
[434]	564	ratio = newCost/oldCost;
	565	#else
	566	float newCost = (pvsBack + pvsFront)*0.5f;
	567	float oldCost = pvsSize;
	568	ratio = newCost/oldCost;
	569	#endif
	570	#endif
	571	}
	572
	573	return ratio;
	574	}
[382]	575
	576
	577	float
	578	VssTree::EvalCostRatio(
[434]	579	VssTreeLeaf *leaf,
	580	const int axis,
	581	const float position,
	582	int &raysBack,
	583	int &raysFront,
	584	int &pvsBack,
	585	int &pvsFront
	586	)
[382]	587	{
[434]	588	raysBack = 0;
	589	raysFront = 0;
	590	pvsFront = 0;
	591	pvsBack = 0;
[382]	592
[434]	593
	594	Intersectable::NewMail(3);
	595
	596	if (axis <= VssTreeNode::SPLIT_Z) {
[382]	597	// this is the main ray classification loop!
	598	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	599	ri != leaf->rays.end();
	600	ri++)
[382]	601	if ((*ri).mRay->IsActive()) {
[463]	602	float t;
[434]	603	// determine the side of this ray with respect to the plane
[463]	604	int side = (*ri).ComputeRayIntersection(axis, position, t);
[434]	605	// (*ri).mRay->mSide = side;
[382]	606
[434]	607	if (side <= 0)
	608	raysBack++;
[382]	609
[434]	610	if (side >= 0)
	611	raysFront++;
[382]	612
[434]	613	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
[382]	614	}
	615
	616	} else {
	617
[434]	618	// directional split
[382]	619	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	620	ri != leaf->rays.end();
	621	ri++)
[382]	622	if ((*ri).mRay->IsActive()) {
	623
[434]	624	// determine the side of this ray with respect to the plane
	625	int side;
	626	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
	627	side = 1;
	628	else
	629	side = -1;
[382]	630
[434]	631	if (side <= 0)
	632	raysBack++;
[382]	633
[434]	634	if (side >= 0)
	635	raysFront++;
[382]	636
[434]	637	// (*ri).mRay->mSide = side;
	638	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
[427]	639
[382]	640	}
[434]	641	}
[382]	642
[434]	643	float ratio = GetCostRatio(
	644	leaf,
	645	axis,
	646	position,
	647	raysBack,
	648	raysFront,
	649	pvsBack,
	650	pvsFront);
[427]	651
[434]	652	// cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
	653	// float oldCost = leaf->rays.size();
	654
	655	// cout<<"ratio="<<ratio<<endl;
	656
	657	return ratio;
[382]	658	}
	659
	660	float
[434]	661	VssTree::BestCostRatio(
	662	VssTreeLeaf *leaf,
	663	int &axis,
	664	float &position,
	665	int &raysBack,
	666	int &raysFront,
	667	int &pvsBack,
	668	int &pvsFront
	669	)
	670	{
	671	int nRaysBack[6], nRaysFront[6];
	672	int nPvsBack[6], nPvsFront[6];
	673	float nPosition[6];
	674	float nCostRatio[6];
	675	int bestAxis = -1;
[386]	676
[434]	677	AxisAlignedBox3 sBox = GetBBox(leaf);
	678	AxisAlignedBox3 dBox = GetDirBBox(leaf);
	679	// int sAxis = box.Size().DrivingAxis();
	680	int sAxis = sBox.Size().DrivingAxis();
	681	int dAxis = dBox.Size().DrivingAxis() + 3;
[387]	682
	683
[434]	684	float dirSplitBoxSize = 0.01f;
	685	bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
[427]	686
[438]	687
	688	for (axis = 0; axis < 5; axis++)
	689	if (mInterleaveDirSplits \|\|
	690	(axis < 3 && leaf->depth < mDirSplitDepth) \|\|
	691	(axis >= 3 && leaf->depth >= mDirSplitDepth)
	692	) {
	693	if (!mSplitUseOnlyDrivingAxis \|\| axis == sAxis \|\| axis == dAxis) {
	694
	695
	696	if (splitType == ESplitRegular) {
	697	if (axis < 3)
	698	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
	699	else
	700	nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
	701
	702	nCostRatio[axis] = EvalCostRatio(leaf,
	703	axis,
	704	nPosition[axis],
	705	nRaysBack[axis],
	706	nRaysFront[axis],
	707	nPvsBack[axis],
	708	nPvsFront[axis]
	709	);
	710	} else
	711	if (splitType == ESplitHeuristic) {
	712	nCostRatio[axis] = EvalCostRatioHeuristic(
	713	leaf,
	714	axis,
	715	nPosition[axis],
	716	nRaysBack[axis],
	717	nRaysFront[axis],
	718	nPvsBack[axis],
	719	nPvsFront[axis]);
	720	} else
	721	if (splitType == ESplitHybrid) {
	722	if (leaf->depth > 7)
	723	nCostRatio[axis] = EvalCostRatioHeuristic(
	724	leaf,
	725	axis,
	726	nPosition[axis],
	727	nRaysBack[axis],
	728	nRaysFront[axis],
	729	nPvsBack[axis],
	730	nPvsFront[axis]);
	731	else {
	732	if (axis < 3)
	733	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
	734	else
	735	nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
[434]	736
[438]	737	nCostRatio[axis] = EvalCostRatio(leaf,
	738	axis,
	739	nPosition[axis],
	740	nRaysBack[axis],
	741	nRaysFront[axis],
	742	nPvsBack[axis],
	743	nPvsFront[axis]
	744	);
	745	}
	746	} else {
	747	cerr<<"VssTree: Unknown split heuristics\n";
	748	exit(1);
[434]	749	}
[438]	750
	751
	752	if ( bestAxis == -1)
[434]	753	bestAxis = axis;
[438]	754	else
	755	if ( nCostRatio[axis] < nCostRatio[bestAxis] )
	756	bestAxis = axis;
	757	}
[387]	758	}
[438]	759
[434]	760	axis = bestAxis;
	761	position = nPosition[bestAxis];
[387]	762
[434]	763	raysBack = nRaysBack[bestAxis];
	764	raysFront = nRaysFront[bestAxis];
[387]	765
[434]	766	pvsBack = nPvsBack[bestAxis];
	767	pvsFront = nPvsFront[bestAxis];
[387]	768
[434]	769	return nCostRatio[bestAxis];
[386]	770	}
	771
[434]	772
[386]	773	float
[434]	774	VssTree::EvalCostRatioHeuristic(
	775	VssTreeLeaf *leaf,
	776	const int axis,
	777	float &bestPosition,
	778	int &raysBack,
	779	int &raysFront,
	780	int &pvsBack,
	781	int &pvsFront
	782	)
[386]	783	{
[434]	784	AxisAlignedBox3 box;
	785	float minBox, maxBox;
[387]	786
[434]	787	if (axis < 3) {
	788	box = GetBBox(leaf);
	789	minBox = box.Min(axis);
	790	maxBox = box.Max(axis);
[438]	791	} else {
[434]	792	box = GetDirBBox(leaf);
	793	minBox = box.Min(axis-3);
	794	maxBox = box.Max(axis-3);
	795	}
	796
[1233]	797	SortSubdivisionCandidates(leaf, axis);
[372]	798
[434]	799	// go through the lists, count the number of objects left and right
	800	// and evaluate the following cost funcion:
	801	// C = ct_div_ci + (qlrl + qrrr)/queries
[382]	802
[1302]	803	int rl=0, rr = (int)leaf->rays.size();
	804	int pl=0, pr = (int)leaf->GetPvsSize();
[434]	805	float sizeBox = maxBox - minBox;
[387]	806
[1302]	807	float minBand = minBox + 0.1f*(maxBox - minBox);
	808	float maxBand = minBox + 0.9f*(maxBox - minBox);
[387]	809
[1297]	810	float minRatio = 1e20f;
[387]	811
[434]	812	Intersectable::NewMail();
	813	// set all object as belonging to the fron pvs
	814	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	815	ri != leaf->rays.end();
	816	ri++)
	817	if ((*ri).mRay->IsActive()) {
	818	Intersectable object = (ri).mRay->mTerminationObject;
	819	if (object)
	820	if (!object->Mailed()) {
	821	object->Mail();
	822	object->mCounter = 1;
	823	} else
	824	object->mCounter++;
	825	}
[387]	826
[434]	827	Intersectable::NewMail();
[387]	828
[434]	829	for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
	830	ci < splitCandidates->end();
	831	ci++) {
	832	VssRay *ray;
	833	switch ((*ci).type) {
	834	case SortableEntry::ERayMin: {
	835	rl++;
	836	ray = (VssRay ) (ci).data;
	837	Intersectable *object = ray->mTerminationObject;
	838	if (object && !object->Mailed()) {
	839	object->Mail();
	840	pl++;
	841	}
	842	break;
	843	}
	844	case SortableEntry::ERayMax: {
	845	rr--;
	846	ray = (VssRay ) (ci).data;
	847	Intersectable *object = ray->mTerminationObject;
	848	if (object) {
	849	if (--object->mCounter == 0)
	850	pr--;
	851	}
	852	break;
	853	}
	854	}
	855
	856	float position = (*ci).value;
	857
	858	if (position > minBand && position < maxBand) {
[387]	859
[434]	860	float ratio = GetCostRatio(
	861	leaf,
	862	axis,
	863	position,
	864	rl,
	865	rr,
	866	pl,
	867	pr);
[387]	868
	869
[434]	870	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
	871	// cout<<"cost= "<<sum<<endl;
	872
	873	if (ratio < minRatio) {
	874	minRatio = ratio;
	875	bestPosition = position;
[382]	876
[434]	877	raysBack = rl;
	878	raysFront = rr;
[382]	879
[434]	880	pvsBack = pl;
	881	pvsFront = pr;
[387]	882
[382]	883	}
	884	}
[372]	885	}
[434]	886
[372]	887
[382]	888	// cout<<"===================="<<endl;
	889	// cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
	890	// <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
[434]	891	return minRatio;
[382]	892	}
	893
	894	void
[1233]	895	VssTree::SortSubdivisionCandidates(
[434]	896	VssTreeLeaf *node,
	897	const int axis
	898	)
[382]	899	{
[372]	900
[382]	901	splitCandidates->clear();
	902
[1302]	903	int requestedSize = 2*((int)node->rays.size());
[382]	904	// creates a sorted split candidates array
	905	if (splitCandidates->capacity() > 500000 &&
	906	requestedSize < (int)(splitCandidates->capacity()/10) ) {
	907
	908	delete splitCandidates;
	909	splitCandidates = new vector<SortableEntry>;
[372]	910	}
	911
[382]	912	splitCandidates->reserve(requestedSize);
	913
	914	// insert all queries
	915	for(VssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
	916	ri < node->rays.end();
	917	ri++) {
[434]	918	if ((*ri).mRay->IsActive()) {
	919	if (axis < 3) {
	920	bool positive = (*ri).mRay->HasPosDir(axis);
	921	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
	922	SortableEntry::ERayMax,
	923	(*ri).ExtrapOrigin(axis),
	924	(void )(ri).mRay)
	925	);
	926
	927	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
	928	SortableEntry::ERayMin,
	929	(*ri).ExtrapTermination(axis),
	930	(void )(ri).mRay)
	931	);
	932	} else {
	933	float pos = (*ri).mRay->GetDirParametrization(axis-3);
	934	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
	935	pos - Limits::Small,
	936	(void )(ri).mRay)
	937	);
	938
	939	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMax,
	940	pos + Limits::Small,
	941	(void )(ri).mRay)
	942	);
	943	}
	944	}
[382]	945	}
[434]	946
[382]	947	stable_sort(splitCandidates->begin(), splitCandidates->end());
[372]	948	}
	949
	950
	951	void
	952	VssTree::EvaluateLeafStats(const TraversalData &data)
	953	{
	954
	955	// the node became a leaf -> evaluate stats for leafs
	956	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
	957
[382]	958	if (data.depth >= termMaxDepth)
[372]	959	stat.maxDepthNodes++;
	960
[434]	961	// if ( (int)(leaf->rays.size()) < termMinCost)
	962	// stat.minCostNodes++;
	963	if ( leaf->GetPvsSize() < termMinPvs)
	964	stat.minPvsNodes++;
[403]	965
[434]	966	if ( leaf->GetPvsSize() < termMinRays)
	967	stat.minRaysNodes++;
[403]	968
[434]	969	if (0 && leaf->GetAvgRayContribution() > termMaxRayContribution )
	970	stat.maxRayContribNodes++;
[382]	971
[434]	972	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
	973	stat.minSizeNodes++;
	974	}
[372]	975
[434]	976	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
[469]	977	stat.maxRayRefs = (int)leaf->rays.size();
[382]	978
[372]	979	}
	980
[427]	981	bool
	982	VssTree::TerminationCriteriaSatisfied(VssTreeLeaf *leaf)
	983	{
[434]	984	return ( (leaf->GetPvsSize() < termMinPvs) \|\|
	985	(leaf->rays.size() < termMinRays) \|\|
	986	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) \|\|
	987	(leaf->depth >= termMaxDepth) \|\|
[446]	988	(SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
	989	// \|\|
	990	// (mUseRss && leaf->mPassingRays == leaf->rays.size())
[434]	991	);
[427]	992	}
[372]	993
	994
	995	VssTreeNode *
	996	VssTree::SubdivideNode(
[434]	997	VssTreeLeaf *leaf,
	998	const AxisAlignedBox3 &box,
	999	AxisAlignedBox3 &backBBox,
	1000	AxisAlignedBox3 &frontBBox
	1001	)
[372]	1002	{
	1003
[434]	1004	if (TerminationCriteriaSatisfied(leaf)) {
[382]	1005	#if 0
[434]	1006	if (leaf->depth >= termMaxDepth) {
	1007	cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
	1008	cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
	1009	}
[382]	1010	#endif
	1011
[434]	1012	return leaf;
	1013	}
[382]	1014
[434]	1015	float position;
[382]	1016
[434]	1017	// first count ray sides
[382]	1018	int raysBack;
	1019	int raysFront;
[434]	1020	int pvsBack;
	1021	int pvsFront;
[382]	1022
[372]	1023	// select subdivision axis
[382]	1024	int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
[434]	1025	// cout<<axis<<" ";
[427]	1026
[434]	1027	// cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
	1028	// pvsFront<<endl;
[382]	1029
[372]	1030	if (axis == -1) {
	1031	return leaf;
	1032	}
	1033
	1034	stat.nodes+=2;
	1035	stat.splits[axis]++;
	1036
	1037	// add the new nodes to the tree
	1038	VssTreeInterior *node = new VssTreeInterior(leaf->parent);
	1039
	1040	node->axis = axis;
	1041	node->position = position;
	1042	node->bbox = box;
	1043	node->dirBBox = GetDirBBox(leaf);
	1044
	1045	backBBox = box;
	1046	frontBBox = box;
	1047
	1048	VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
	1049	VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);
	1050
	1051	// replace a link from node's parent
	1052	if ( leaf->parent )
	1053	leaf->parent->ReplaceChildLink(leaf, node);
	1054	// and setup child links
	1055	node->SetupChildLinks(back, front);
[382]	1056
[372]	1057	if (axis <= VssTreeNode::SPLIT_Z) {
[434]	1058	backBBox.SetMax(axis, position);
[382]	1059	frontBBox.SetMin(axis, position);
	1060
[434]	1061	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1062	ri != leaf->rays.end();
	1063	ri++) {
[372]	1064	if ((*ri).mRay->IsActive()) {
	1065
[434]	1066	// first unref ray from the former leaf
	1067	(*ri).mRay->Unref();
[382]	1068
[434]	1069	// Debug << "computed t: " << (*ri).mRay->mT << endl;
	1070	// determine the side of this ray with respect to the plane
[463]	1071	float t;
	1072	int side = node->ComputeRayIntersection(*ri, t);
[434]	1073
	1074	if (side == 0) {
	1075	if ((*ri).mRay->HasPosDir(axis)) {
	1076	back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
	1077	(*ri).mMinT,
[463]	1078	t)
[434]	1079	);
	1080	front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
[463]	1081	t,
[434]	1082	(*ri).mMaxT));
	1083	} else {
	1084	back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
[463]	1085	t,
[434]	1086	(*ri).mMaxT));
	1087	front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
	1088	(*ri).mMinT,
[463]	1089	t));
[434]	1090	}
	1091	} else
	1092	if (side == 1)
	1093	front->AddRay(*ri);
	1094	else
	1095	back->AddRay(*ri);
[372]	1096	} else
[434]	1097	(*ri).mRay->Unref();
[372]	1098	}
	1099	} else {
	1100	// rays front/back
	1101
	1102
	1103	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	1104	ri != leaf->rays.end();
	1105	ri++) {
[372]	1106	if ((*ri).mRay->IsActive()) {
[434]	1107	// first unref ray from the former leaf
	1108	(*ri).mRay->Unref();
[382]	1109
[434]	1110	int side;
	1111	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
	1112	side = 1;
	1113	else
	1114	side = -1;
[372]	1115
[434]	1116	if (side == 1)
	1117	front->AddRay(*ri);
	1118	else
	1119	back->AddRay(*ri);
[372]	1120
	1121	} else
[434]	1122	(*ri).mRay->Unref();
[372]	1123	}
	1124	}
[434]	1125
	1126	front->SetPvsSize(pvsFront);
	1127	back->SetPvsSize(pvsBack);
[438]	1128	// compute entropy as well
	1129	front->ComputeEntropyImportance();
	1130	back->ComputeEntropyImportance();
	1131
[372]	1132	// update stats
[469]	1133	stat.rayRefs -= (int)leaf->rays.size();
[372]	1134	stat.rayRefs += raysBack + raysFront;
	1135
	1136
	1137	delete leaf;
	1138	return node;
	1139	}
	1140
	1141
	1142
	1143
	1144
	1145
	1146	int
	1147	VssTree::ReleaseMemory(const int time)
	1148	{
	1149	stack<VssTreeNode *> tstack;
	1150
	1151	// find a node in the tree which subtree will be collapsed
	1152	int maxAccessTime = time - accessTimeThreshold;
	1153	int released;
	1154
	1155	tstack.push(root);
	1156
	1157	while (!tstack.empty()) {
	1158	VssTreeNode *node = tstack.top();
	1159	tstack.pop();
	1160
	1161
	1162	if (!node->IsLeaf()) {
	1163	VssTreeInterior in = (VssTreeInterior )node;
	1164	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
	1165	if (in->depth >= minCollapseDepth &&
[434]	1166	in->lastAccessTime <= maxAccessTime) {
	1167	released = CollapseSubtree(node, time);
	1168	break;
[372]	1169	}
	1170
	1171	if (in->back->GetAccessTime() <
[434]	1172	in->front->GetAccessTime()) {
	1173	tstack.push(in->front);
	1174	tstack.push(in->back);
[372]	1175	} else {
[434]	1176	tstack.push(in->back);
	1177	tstack.push(in->front);
[372]	1178	}
	1179	}
	1180	}
	1181
	1182	while (tstack.empty()) {
	1183	// could find node to collaps...
	1184	// cout<<"Could not find a node to release "<<endl;
	1185	break;
	1186	}
	1187
	1188	return released;
	1189	}
	1190
	1191
	1192
	1193
	1194	VssTreeNode *
	1195	VssTree::SubdivideLeaf(
[434]	1196	VssTreeLeaf *leaf
	1197	)
[372]	1198	{
	1199	VssTreeNode *node = leaf;
[427]	1200
[372]	1201	AxisAlignedBox3 leafBBox = GetBBox(leaf);
	1202
[434]	1203	static int pass = 0;
	1204	pass ++;
[372]	1205
	1206	// check if we should perform a dynamic subdivision of the leaf
[434]	1207	if (!TerminationCriteriaSatisfied(leaf)) {
[372]	1208
[434]	1209	// memory check and realese...
[372]	1210	if (GetMemUsage() > maxTotalMemory) {
	1211	ReleaseMemory( pass );
	1212	}
	1213
	1214	AxisAlignedBox3 backBBox, frontBBox;
	1215
	1216	// subdivide the node
	1217	node =
	1218	SubdivideNode(leaf,
[434]	1219	leafBBox,
	1220	backBBox,
	1221	frontBBox
	1222	);
[372]	1223	}
[427]	1224
[372]	1225	return node;
	1226	}
	1227
	1228
	1229
	1230	void
	1231	VssTree::UpdateRays(VssRayContainer &remove,
[434]	1232	VssRayContainer &add
	1233	)
[372]	1234	{
	1235	VssTreeLeaf::NewMail();
	1236
	1237	// schedule rays for removal
	1238	for(VssRayContainer::const_iterator ri = remove.begin();
	1239	ri != remove.end();
	1240	ri++) {
	1241	(*ri)->ScheduleForRemoval();
	1242	}
	1243
	1244	int inactive=0;
	1245
	1246	for(VssRayContainer::const_iterator ri = remove.begin();
	1247	ri != remove.end();
	1248	ri++) {
	1249	if ((*ri)->ScheduledForRemoval())
[434]	1250	// RemoveRay(*ri, NULL, false);
	1251	// !!! BUG - with true it does not work correctly - aggreated delete
[372]	1252	RemoveRay(*ri, NULL, true);
	1253	else
	1254	inactive++;
	1255	}
	1256
	1257
	1258	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
	1259
	1260	for(VssRayContainer::const_iterator ri = add.begin();
	1261	ri != add.end();
	1262	ri++) {
[434]	1263	VssTreeNode::RayInfo info(*ri);
	1264	if (ClipRay(info, bbox))
	1265	AddRay(info);
[372]	1266	}
	1267	}
	1268
	1269
	1270	void
	1271	VssTree::RemoveRay(VssRay *ray,
[434]	1272	vector<VssTreeLeaf > affectedLeaves,
	1273	const bool removeAllScheduledRays
	1274	)
[372]	1275	{
	1276
	1277	stack<RayTraversalData> tstack;
	1278
	1279	tstack.push(RayTraversalData(root, VssTreeNode::RayInfo(ray)));
	1280
	1281	RayTraversalData data;
	1282
	1283	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
	1284
	1285	while (!tstack.empty()) {
	1286	data = tstack.top();
	1287	tstack.pop();
	1288
	1289	if (!data.node->IsLeaf()) {
	1290	// split the set of rays in two groups intersecting the
	1291	// two subtrees
	1292
	1293	TraverseInternalNode(data, tstack);
	1294
	1295	} else {
	1296	// remove the ray from the leaf
	1297	// find the ray in the leaf and swap it with the last ray...
	1298	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
	1299
	1300	if (!leaf->Mailed()) {
[434]	1301	leaf->Mail();
	1302	if (affectedLeaves)
	1303	affectedLeaves->push_back(leaf);
[372]	1304
[434]	1305	if (removeAllScheduledRays) {
[1302]	1306	int tail = (int)leaf->rays.size()-1;
[372]	1307
[434]	1308	for (int i=0; i < (int)(leaf->rays.size()); i++) {
	1309	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
	1310	// find a ray to replace it with
	1311	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
	1312	stat.removedRayRefs++;
	1313	leaf->rays[tail].mRay->Unref();
	1314	leaf->rays.pop_back();
	1315	tail--;
	1316	}
[372]	1317
[434]	1318	if (tail < i)
	1319	break;
[372]	1320
[434]	1321	stat.removedRayRefs++;
	1322	leaf->rays[i].mRay->Unref();
	1323	leaf->rays[i] = leaf->rays[tail];
	1324	leaf->rays.pop_back();
	1325	tail--;
	1326	}
	1327	}
	1328	}
[372]	1329	}
	1330
	1331	if (!removeAllScheduledRays)
[434]	1332	for (int i=0; i < (int)leaf->rays.size(); i++) {
	1333	if (leaf->rays[i].mRay == ray) {
	1334	stat.removedRayRefs++;
	1335	ray->Unref();
	1336	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
	1337	leaf->rays.pop_back();
	1338	// check this ray again
	1339	break;
	1340	}
	1341	}
[372]	1342
	1343	}
	1344	}
	1345
	1346	if (ray->RefCount() != 0) {
	1347	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
	1348	exit(1);
	1349	}
	1350
	1351	}
	1352
	1353
	1354	void
[427]	1355	VssTree::AddRay(VssTreeNode::RayInfo &info)
[372]	1356	{
	1357
	1358	stack<RayTraversalData> tstack;
	1359
[427]	1360	tstack.push(RayTraversalData(root, info));
[372]	1361
	1362	RayTraversalData data;
	1363
	1364	while (!tstack.empty()) {
	1365	data = tstack.top();
	1366	tstack.pop();
	1367
	1368	if (!data.node->IsLeaf()) {
	1369	TraverseInternalNode(data, tstack);
	1370	} else {
	1371	// remove the ray from the leaf
	1372	// find the ray in the leaf and swap it with the last ray...
	1373	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
	1374	leaf->AddRay(data.rayData);
	1375	stat.addedRayRefs++;
	1376	}
	1377	}
	1378	}
	1379
	1380	void
	1381	VssTree::TraverseInternalNode(
[434]	1382	RayTraversalData &data,
	1383	stack<RayTraversalData> &tstack)
[372]	1384	{
	1385	VssTreeInterior in = (VssTreeInterior ) data.node;
	1386
	1387	if (in->axis <= VssTreeNode::SPLIT_Z) {
[463]	1388	float t;
[372]	1389	// determine the side of this ray with respect to the plane
	1390	int side = in->ComputeRayIntersection(data.rayData,
[463]	1391	t);
[372]	1392
	1393
	1394	if (side == 0) {
	1395	if (data.rayData.mRay->HasPosDir(in->axis)) {
[434]	1396	tstack.push(RayTraversalData(in->back,
	1397	VssTreeNode::RayInfo(data.rayData.mRay,
	1398	data.rayData.mMinT,
[463]	1399	t))
[434]	1400	);
[372]	1401
[434]	1402	tstack.push(RayTraversalData(in->front,
	1403	VssTreeNode::RayInfo(data.rayData.mRay,
[463]	1404	t,
[434]	1405	data.rayData.mMaxT
	1406	))
	1407	);
[372]	1408
	1409	} else {
[434]	1410	tstack.push(RayTraversalData(in->back,
	1411	VssTreeNode::RayInfo(data.rayData.mRay,
[463]	1412	t,
[434]	1413	data.rayData.mMaxT
	1414	))
	1415	);
[372]	1416
[434]	1417	tstack.push(RayTraversalData(in->front,
	1418	VssTreeNode::RayInfo(data.rayData.mRay,
	1419	data.rayData.mMinT,
[463]	1420	t))
[434]	1421	);
[372]	1422
	1423
	1424	}
	1425	} else
	1426	if (side == 1)
[434]	1427	tstack.push(RayTraversalData(in->front, data.rayData));
[372]	1428	else
[434]	1429	tstack.push(RayTraversalData(in->back, data.rayData));
[372]	1430	}
	1431	else {
	1432	// directional split
[434]	1433	if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position)
	1434	tstack.push(RayTraversalData(in->front, data.rayData));
	1435	else
	1436	tstack.push(RayTraversalData(in->back, data.rayData));
[372]	1437	}
	1438	}
	1439
	1440
	1441	int
	1442	VssTree::CollapseSubtree(VssTreeNode *sroot, const int time)
	1443	{
	1444	// first count all rays in the subtree
	1445	// use mail 1 for this purpose
	1446	stack<VssTreeNode *> tstack;
	1447	int rayCount = 0;
	1448	int totalRayCount = 0;
	1449	int collapsedNodes = 0;
	1450
	1451	#if DEBUG_COLLAPSE
	1452	cout<<"Collapsing subtree"<<endl;
	1453	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
	1454	cout<<"depth="<<(int)sroot->depth<<endl;
	1455	#endif
	1456
[434]	1457	// tstat.collapsedSubtrees++;
	1458	// tstat.collapseDepths += (int)sroot->depth;
	1459	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
[372]	1460
	1461	tstack.push(sroot);
	1462	VssRay::NewMail();
	1463
	1464	while (!tstack.empty()) {
	1465	collapsedNodes++;
	1466	VssTreeNode *node = tstack.top();
	1467	tstack.pop();
	1468
	1469	if (node->IsLeaf()) {
	1470	VssTreeLeaf leaf = (VssTreeLeaf ) node;
	1471	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	1472	ri != leaf->rays.end();
	1473	ri++) {
[372]	1474
[434]	1475	totalRayCount++;
	1476	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
	1477	(*ri).mRay->Mail();
	1478	rayCount++;
	1479	}
[372]	1480	}
	1481	} else {
	1482	tstack.push(((VssTreeInterior *)node)->back);
	1483	tstack.push(((VssTreeInterior *)node)->front);
	1484	}
	1485	}
	1486
	1487	VssRay::NewMail();
	1488
	1489	// create a new node that will hold the rays
	1490	VssTreeLeaf *newLeaf = new VssTreeLeaf( sroot->parent, rayCount );
	1491	if ( newLeaf->parent )
	1492	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
	1493
	1494
	1495	tstack.push( sroot );
	1496
	1497	while (!tstack.empty()) {
	1498
	1499	VssTreeNode *node = tstack.top();
	1500	tstack.pop();
	1501
	1502	if (node->IsLeaf()) {
	1503	VssTreeLeaf leaf = (VssTreeLeaf ) node;
	1504
	1505	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	1506	ri != leaf->rays.end();
	1507	ri++) {
[372]	1508
[434]	1509	// unref this ray from the old node
[372]	1510
[434]	1511	if ((*ri).mRay->IsActive()) {
	1512	(*ri).mRay->Unref();
	1513	if (!(*ri).mRay->Mailed()) {
	1514	(*ri).mRay->Mail();
	1515	newLeaf->AddRay(*ri);
	1516	}
	1517	} else
	1518	(*ri).mRay->Unref();
[372]	1519
	1520	}
	1521	} else {
	1522	tstack.push(((VssTreeInterior *)node)->back);
	1523	tstack.push(((VssTreeInterior *)node)->front);
	1524	}
	1525	}
	1526
	1527	// delete the node and all its children
	1528	delete sroot;
	1529
[434]	1530	// for(VssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
	1531	// ri != newLeaf->rays.end();
	1532	// ri++)
	1533	// (*ri).ray->UnMail(2);
[372]	1534
	1535
	1536	#if DEBUG_COLLAPSE
	1537	cout<<"Total memory before="<<GetMemUsage()<<endl;
	1538	#endif
	1539
	1540	stat.nodes -= collapsedNodes - 1;
	1541	stat.rayRefs -= totalRayCount - rayCount;
	1542
	1543	#if DEBUG_COLLAPSE
	1544	cout<<"collapsed nodes"<<collapsedNodes<<endl;
	1545	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
	1546	cout<<"Total memory after="<<GetMemUsage()<<endl;
	1547	cout<<"================================"<<endl;
	1548	#endif
	1549
[434]	1550	// tstat.collapsedNodes += collapsedNodes;
	1551	// tstat.collapsedRays += totalRayCount - rayCount;
[372]	1552
	1553	return totalRayCount - rayCount;
	1554	}
[387]	1555
	1556
	1557	int
[434]	1558	VssTree::GetPvsSize(const AxisAlignedBox3 &box) const
[387]	1559	{
[434]	1560	stack<VssTreeNode *> tstack;
[387]	1561	tstack.push(root);
	1562
[434]	1563	Intersectable::NewMail();
	1564	int pvsSize = 0;
[387]	1565
	1566	while (!tstack.empty()) {
	1567	VssTreeNode *node = tstack.top();
	1568	tstack.pop();
	1569
	1570
	1571	if (node->IsLeaf()) {
[434]	1572	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1573	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1574	ri != leaf->rays.end();
	1575	ri++)
	1576	if ((*ri).mRay->IsActive()) {
	1577	Intersectable *object;
[387]	1578	#if BIDIRECTIONAL_RAY
[434]	1579	object = (*ri).mRay->mOriginObject;
	1580	if (object && !object->Mailed()) {
	1581	pvsSize++;
	1582	object->Mail();
	1583	}
[387]	1584	#endif
[434]	1585	object = (*ri).mRay->mTerminationObject;
	1586	if (object && !object->Mailed()) {
	1587	pvsSize++;
	1588	object->Mail();
	1589	}
	1590	}
	1591	} else {
	1592	VssTreeInterior in = (VssTreeInterior )node;
	1593	if (in->axis < 3) {
	1594	if (box.Max(in->axis) >= in->position )
	1595	tstack.push(in->front);
[387]	1596
[434]	1597	if (box.Min(in->axis) <= in->position )
	1598	tstack.push(in->back);
	1599	} else {
	1600	// both nodes for directional splits
	1601	tstack.push(in->front);
	1602	tstack.push(in->back);
	1603	}
[387]	1604	}
[434]	1605	}
	1606	return pvsSize;
[387]	1607	}
[401]	1608
	1609	void
	1610	VssTree::GetRayContributionStatistics(
[434]	1611	float &minRayContribution,
	1612	float &maxRayContribution,
	1613	float &avgRayContribution
	1614	)
[401]	1615	{
[434]	1616	stack<VssTreeNode *> tstack;
[401]	1617	tstack.push(root);
	1618
[434]	1619	minRayContribution = 1.0f;
	1620	maxRayContribution = 0.0f;
	1621	float sumRayContribution = 0.0f;
	1622	int leaves = 0;
[401]	1623
	1624	while (!tstack.empty()) {
	1625	VssTreeNode *node = tstack.top();
	1626	tstack.pop();
	1627
	1628	if (node->IsLeaf()) {
[434]	1629	leaves++;
	1630	VssTreeLeaf leaf = (VssTreeLeaf )node;
[435]	1631	float c = leaf->GetImportance();
[434]	1632	if (c > maxRayContribution)
	1633	maxRayContribution = c;
	1634	if (c < minRayContribution)
	1635	minRayContribution = c;
	1636	sumRayContribution += c;
[401]	1637
[434]	1638	} else {
	1639	VssTreeInterior in = (VssTreeInterior )node;
	1640	// both nodes for directional splits
	1641	tstack.push(in->front);
	1642	tstack.push(in->back);
[401]	1643	}
[434]	1644	}
[401]	1645
[434]	1646	cout<<"sum="<<sumRayContribution<<endl;
	1647	cout<<"leaves="<<leaves<<endl;
	1648	avgRayContribution = sumRayContribution/(float)leaves;
[401]	1649	}
	1650
	1651
	1652	int
	1653	VssTree::GenerateRays(const float ratioPerLeaf,
[434]	1654	SimpleRayContainer &rays)
[401]	1655	{
[434]	1656	stack<VssTreeNode *> tstack;
[401]	1657	tstack.push(root);
	1658
	1659	while (!tstack.empty()) {
	1660	VssTreeNode *node = tstack.top();
	1661	tstack.pop();
	1662
	1663	if (node->IsLeaf()) {
[434]	1664	VssTreeLeaf leaf = (VssTreeLeaf )node;
[435]	1665	float c = leaf->GetImportance();
[1302]	1666	int num = int(c*ratioPerLeaf + 0.5f);
[434]	1667	// cout<<num<<" ";
[401]	1668
[434]	1669	for (int i=0; i < num; i++) {
	1670	Vector3 origin = GetBBox(leaf).GetRandomPoint();
	1671	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
	1672	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
	1673	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	1674	rays.push_back(SimpleRay(origin, direction));
	1675	}
[401]	1676
[434]	1677	} else {
	1678	VssTreeInterior in = (VssTreeInterior )node;
	1679	// both nodes for directional splits
	1680	tstack.push(in->front);
	1681	tstack.push(in->back);
[401]	1682	}
[434]	1683	}
[401]	1684
[1302]	1685	return (int)rays.size();
[401]	1686	}
	1687
[427]	1688	void
	1689	VssTree::CollectLeaves(vector<VssTreeLeaf *> &leaves)
	1690	{
[434]	1691	stack<VssTreeNode *> tstack;
[427]	1692	tstack.push(root);
	1693
	1694	while (!tstack.empty()) {
	1695	VssTreeNode *node = tstack.top();
	1696	tstack.pop();
	1697
	1698	if (node->IsLeaf()) {
[434]	1699	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1700	leaves.push_back(leaf);
	1701	} else {
	1702	VssTreeInterior in = (VssTreeInterior )node;
	1703	// both nodes for directional splits
	1704	tstack.push(in->front);
	1705	tstack.push(in->back);
[427]	1706	}
[434]	1707	}
[427]	1708	}
[401]	1709
[434]	1710	bool
	1711	VssTree::ValidLeaf(VssTreeLeaf *leaf) const
	1712	{
	1713	return leaf->rays.size() > termMinRays/4;
	1714	}
	1715
	1716
	1717	void
	1718	VssTree::GenerateLeafRays(VssTreeLeaf *leaf,
	1719	const int numberOfRays,
	1720	SimpleRayContainer &rays)
	1721	{
[469]	1722	int nrays = (int)leaf->rays.size();
[1297]	1723	for (int i = 0; i < numberOfRays; ++ i) {
[434]	1724	// pickup 3 random rays
[1297]	1725	int r1 = (int)RandomValue(0, (float)(nrays-1));
	1726	int r2 = (int)RandomValue(0, (float)(nrays-1));
	1727	int r3 = (int)RandomValue(0, (float)(nrays-1));
[434]	1728
	1729	Vector3 o1 = leaf->rays[r1].Extrap(RandomValue(leaf->rays[r1].GetMinT(),
	1730	leaf->rays[r1].GetMaxT()));
	1731
	1732	Vector3 o2 = leaf->rays[r2].Extrap(RandomValue(leaf->rays[r2].GetMinT(),
	1733	leaf->rays[r2].GetMaxT()));
	1734
	1735	Vector3 o3 = leaf->rays[r3].Extrap(RandomValue(leaf->rays[r3].GetMinT(),
	1736	leaf->rays[r3].GetMaxT()));
	1737
[469]	1738	const float overlap = 0.1f;
[434]	1739
	1740	Vector3 origin, direction;
	1741	bool useExtendedConvexCombination = true;
	1742	if (useExtendedConvexCombination) {
	1743	float w1, w2, w3;
	1744	GenerateExtendedConvexCombinationWeights(w1, w2, w3, overlap);
	1745	origin = w1o1 + w2o2 + w3*o3;
	1746	direction =
	1747	w1*leaf->rays[r1].mRay->GetDir() +
	1748	w2*leaf->rays[r2].mRay->GetDir() +
	1749	w3*leaf->rays[r3].mRay->GetDir();
	1750	} else {
	1751	origin = GetBBox(leaf).GetRandomPoint();
	1752	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
	1753	direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
	1754	}
	1755	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	1756	rays.push_back(SimpleRay(origin, direction));
	1757	}
	1758	}
	1759
[427]	1760	int
	1761	VssTree::GenerateRays(const int numberOfRays,
[434]	1762	const int numberOfLeaves,
	1763	SimpleRayContainer &rays)
[427]	1764	{
	1765
[434]	1766	vector<VssTreeLeaf *> leaves;
[427]	1767
[434]	1768	CollectLeaves(leaves);
[427]	1769
[434]	1770	sort(leaves.begin(),
	1771	leaves.end(),
	1772	GreaterContribution);
[427]	1773
	1774
[434]	1775	float sumContrib = 0.0;
	1776	int i;
	1777	int k = 0;
	1778	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	1779	if (ValidLeaf(leaves[i])) {
[435]	1780	float c = leaves[i]->GetImportance();
[434]	1781	sumContrib += c;
	1782	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
	1783	k++;
[427]	1784	}
[434]	1785
	1786	float avgContrib = sumContrib/numberOfLeaves;
	1787	float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
	1788	k = 0;
	1789	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	1790	if (ValidLeaf(leaves[i])) {
	1791	k++;
	1792	VssTreeLeaf *leaf = leaves[i];
[435]	1793	float c = leaf->GetImportance();
[469]	1794	int num = (int)(c*ratioPerLeaf + 0.5f);
[434]	1795	GenerateLeafRays(leaf, num, rays);
	1796	}
[427]	1797
[469]	1798	return (int)rays.size();
[427]	1799	}
	1800
	1801
[401]	1802	float
	1803	VssTree::GetAvgPvsSize()
	1804	{
[434]	1805	stack<VssTreeNode *> tstack;
[401]	1806	tstack.push(root);
	1807
[434]	1808	int sumPvs = 0;
	1809	int leaves = 0;
[401]	1810	while (!tstack.empty()) {
	1811	VssTreeNode *node = tstack.top();
	1812	tstack.pop();
	1813
	1814	if (node->IsLeaf()) {
[434]	1815	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1816	// update pvs size
	1817	leaf->UpdatePvsSize();
	1818	sumPvs += leaf->GetPvsSize();
	1819	leaves++;
	1820	} else {
	1821	VssTreeInterior in = (VssTreeInterior )node;
	1822	// both nodes for directional splits
	1823	tstack.push(in->front);
	1824	tstack.push(in->back);
[401]	1825	}
[434]	1826	}
[401]	1827
	1828
[434]	1829	return sumPvs/(float)leaves;
[401]	1830	}
[427]	1831
	1832
[435]	1833	float
	1834	VssTreeLeaf::GetImportance() const
	1835	{
[438]	1836
	1837	if (1) {
	1838	// return GetAvgRayContribution();
[863]	1839	return (float)GetPvsSize();
[438]	1840	} else {
[446]	1841	// return GetAvgRayContribution()*mEntropyImportance;
	1842	//return GetAvgRayContribution();
	1843	return mEntropyImportance;
[438]	1844	}
[435]	1845	}
[438]	1846
	1847
	1848	float
	1849	VssTreeLeaf::ComputePvsEntropy()
	1850	{
	1851	int samples = 0;
	1852	Intersectable::NewMail();
	1853	// set all object as belonging to the fron pvs
	1854	for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
	1855	ri != rays.end();
	1856	ri++)
	1857	if ((*ri).mRay->IsActive()) {
	1858	Intersectable object = (ri).mRay->mTerminationObject;
	1859	if (object) {
	1860	if (!object->Mailed()) {
	1861	object->Mail();
	1862	object->mCounter = 1;
	1863	} else
	1864	object->mCounter++;
	1865	samples++;
	1866	}
	1867	}
	1868
	1869	float entropy = 0.0f;
	1870
	1871	if (samples > 1) {
	1872	Intersectable::NewMail();
	1873	for(RayInfoContainer::const_iterator ri = rays.begin();
	1874	ri != rays.end();
	1875	ri++)
	1876	if ((*ri).mRay->IsActive()) {
	1877	Intersectable object = (ri).mRay->mTerminationObject;
	1878	if (object) {
	1879	if (!object->Mailed()) {
	1880	object->Mail();
	1881	float p = object->mCounter/(float)samples;
	1882	entropy -= p*log(p);
	1883	}
	1884	}
	1885	}
	1886	entropy = entropy/log((float)samples);
	1887	}
	1888	else
	1889	entropy = 1.0f;
	1890
	1891	return entropy;
	1892	}
	1893
	1894	float
	1895	VssTreeLeaf::ComputeRayLengthEntropy()
	1896	{
	1897	// get sum of all ray lengths
	1898	// consider only passing rays or originating rays
	1899	float sum = 0.0f;
	1900	int samples = 0;
	1901
	1902	for(RayInfoContainer::const_iterator ri = rays.begin();
	1903	ri != rays.end();
	1904	ri++) {
	1905	int rayClass = (*ri).GetRayClass();
	1906	if (
	1907	rayClass == RayInfo::PASSING_RAY
	1908	// \|\|
	1909	// rayClass == RayInfo::SOURCE_RAY
	1910	// rayClass == RayInfo::TERMINATION_RAY
	1911	) {
	1912	float c = 1.0f - (*ri).GetMaxT();
	1913	sum += (ri).mRay->GetSize()c;
	1914	samples++;
	1915	}
	1916	}
	1917
	1918	float entropy = 0.0f;
	1919
	1920	if (samples > 1) {
	1921	for(RayInfoContainer::const_iterator ri = rays.begin();
	1922	ri != rays.end();
	1923	ri++) {
	1924	int rayClass = (*ri).GetRayClass();
	1925	if (
	1926	rayClass == RayInfo::PASSING_RAY
	1927	// \|\|
	1928	// rayClass == RayInfo::SOURCE_RAY
	1929	// rayClass == RayInfo::TERMINATION_RAY
	1930	) {
	1931	float c = 1.0f - (*ri).GetMaxT();
	1932	float p = (ri).mRay->GetSize()c/sum;
	1933	entropy -= p*log(p);
	1934	}
	1935	}
	1936	entropy = entropy/log((float)samples);
	1937	} else
	1938	entropy = 1.0f;
	1939
	1940	return entropy;
	1941	}
	1942
	1943	float
	1944	VssTreeLeaf::ComputeRayTerminationEntropy()
	1945	{
	1946
	1947	return 0.0f;
	1948	}
	1949
	1950
	1951	void
	1952	VssTreeLeaf::ComputeEntropyImportance()
	1953	{
	1954	// mEntropy = 1.0f - ComputeRayLengthEntropy();
	1955	mEntropyImportance = 1.0f - ComputePvsEntropy();
	1956
	1957	// cout<<"ei="<<mEntropyImportance<<" ";
	1958	}
[466]	1959
	1960
	1961	int
	1962	VssTree::CollectRays(VssRayContainer &rays,
	1963	const int number)
	1964	{
	1965	VssRayContainer allRays;
	1966	CollectRays(allRays);
	1967
	1968	int desired = min(number, (int)allRays.size());
	1969	float prob = desired/(float)allRays.size();
	1970	while (rays.size() < desired) {
	1971	VssRayContainer::const_iterator it = allRays.begin();
	1972	for (; it != allRays.end() && rays.size() < desired; it++) {
	1973	if (Random(1.0f) < prob)
	1974	rays.push_back(*it);
	1975	}
	1976	}
[1302]	1977	return (int)rays.size();
[466]	1978	}
	1979
	1980
	1981	int
	1982	VssTree::CollectRays(VssRayContainer &rays
	1983	)
	1984	{
	1985	VssRay::NewMail();
	1986
	1987	stack<VssTreeNode *> tstack;
	1988	tstack.push(root);
	1989
	1990	while (!tstack.empty()) {
	1991	VssTreeNode *node = tstack.top();
	1992	tstack.pop();
	1993
	1994	if (node->IsLeaf()) {
	1995	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1996	// update pvs size
	1997	VssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
	1998	for (;it != leaf->rays.end(); ++it)
	1999	if (!(*it).mRay->Mailed()) {
	2000	(*it).mRay->Mail();
	2001	rays.push_back((*it).mRay);
	2002	}
	2003	} else {
	2004	VssTreeInterior in = (VssTreeInterior )node;
	2005	// both nodes for directional splits
	2006	tstack.push(in->front);
	2007	tstack.push(in->back);
	2008	}
	2009	}
	2010
[1302]	2011	return (int)rays.size();
[466]	2012	}
[860]	2013
	2014	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: