Context Navigation

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

Revision 1233, 47.1 KB checked in by mattausch, 18 years ago (diff)

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
[372]	282	root = new VssTreeLeaf(NULL, rays.size());
	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
[372]	323	stat.rays = 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;
	545	float oldCost = pvsSize;
	546	ratio = newCost/oldCost;
	547	} else {
	548	// importance based cost
	549	#if 0
	550	float newContrib =
	551	((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
	552	(maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
	553
	554	// float newContrib =
	555	// sqr(pvsBack/(raysBack + Limits::Small)) +
	556	// sqr(pvsFront/(raysFront + Limits::Small));
	557	float oldContrib = sqr(leaf->GetAvgRayContribution());
	558	ratio = oldContrib/newContrib;
	559	#else
	560	#if 1
	561	float newCost = raysBackpvsBack + raysFrontpvsFront;
[485]	562	float oldCost = (float)leaf->rays.size()*pvsSize;
[434]	563	ratio = newCost/oldCost;
	564	#else
	565	float newCost = (pvsBack + pvsFront)*0.5f;
	566	float oldCost = pvsSize;
	567	ratio = newCost/oldCost;
	568	#endif
	569	#endif
	570	}
	571
	572	return ratio;
	573	}
[382]	574
	575
	576	float
	577	VssTree::EvalCostRatio(
[434]	578	VssTreeLeaf *leaf,
	579	const int axis,
	580	const float position,
	581	int &raysBack,
	582	int &raysFront,
	583	int &pvsBack,
	584	int &pvsFront
	585	)
[382]	586	{
[434]	587	raysBack = 0;
	588	raysFront = 0;
	589	pvsFront = 0;
	590	pvsBack = 0;
[382]	591
[434]	592
	593	Intersectable::NewMail(3);
	594
	595	if (axis <= VssTreeNode::SPLIT_Z) {
[382]	596	// this is the main ray classification loop!
	597	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	598	ri != leaf->rays.end();
	599	ri++)
[382]	600	if ((*ri).mRay->IsActive()) {
[463]	601	float t;
[434]	602	// determine the side of this ray with respect to the plane
[463]	603	int side = (*ri).ComputeRayIntersection(axis, position, t);
[434]	604	// (*ri).mRay->mSide = side;
[382]	605
[434]	606	if (side <= 0)
	607	raysBack++;
[382]	608
[434]	609	if (side >= 0)
	610	raysFront++;
[382]	611
[434]	612	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
[382]	613	}
	614
	615	} else {
	616
[434]	617	// directional split
[382]	618	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	619	ri != leaf->rays.end();
	620	ri++)
[382]	621	if ((*ri).mRay->IsActive()) {
	622
[434]	623	// determine the side of this ray with respect to the plane
	624	int side;
	625	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
	626	side = 1;
	627	else
	628	side = -1;
[382]	629
[434]	630	if (side <= 0)
	631	raysBack++;
[382]	632
[434]	633	if (side >= 0)
	634	raysFront++;
[382]	635
[434]	636	// (*ri).mRay->mSide = side;
	637	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
[427]	638
[382]	639	}
[434]	640	}
[382]	641
[434]	642	float ratio = GetCostRatio(
	643	leaf,
	644	axis,
	645	position,
	646	raysBack,
	647	raysFront,
	648	pvsBack,
	649	pvsFront);
[427]	650
[434]	651	// cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
	652	// float oldCost = leaf->rays.size();
	653
	654	// cout<<"ratio="<<ratio<<endl;
	655
	656	return ratio;
[382]	657	}
	658
	659	float
[434]	660	VssTree::BestCostRatio(
	661	VssTreeLeaf *leaf,
	662	int &axis,
	663	float &position,
	664	int &raysBack,
	665	int &raysFront,
	666	int &pvsBack,
	667	int &pvsFront
	668	)
	669	{
	670	int nRaysBack[6], nRaysFront[6];
	671	int nPvsBack[6], nPvsFront[6];
	672	float nPosition[6];
	673	float nCostRatio[6];
	674	int bestAxis = -1;
[386]	675
[434]	676	AxisAlignedBox3 sBox = GetBBox(leaf);
	677	AxisAlignedBox3 dBox = GetDirBBox(leaf);
	678	// int sAxis = box.Size().DrivingAxis();
	679	int sAxis = sBox.Size().DrivingAxis();
	680	int dAxis = dBox.Size().DrivingAxis() + 3;
[387]	681
	682
[434]	683	float dirSplitBoxSize = 0.01f;
	684	bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
[427]	685
[438]	686
	687	for (axis = 0; axis < 5; axis++)
	688	if (mInterleaveDirSplits \|\|
	689	(axis < 3 && leaf->depth < mDirSplitDepth) \|\|
	690	(axis >= 3 && leaf->depth >= mDirSplitDepth)
	691	) {
	692	if (!mSplitUseOnlyDrivingAxis \|\| axis == sAxis \|\| axis == dAxis) {
	693
	694
	695	if (splitType == ESplitRegular) {
	696	if (axis < 3)
	697	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
	698	else
	699	nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
	700
	701	nCostRatio[axis] = EvalCostRatio(leaf,
	702	axis,
	703	nPosition[axis],
	704	nRaysBack[axis],
	705	nRaysFront[axis],
	706	nPvsBack[axis],
	707	nPvsFront[axis]
	708	);
	709	} else
	710	if (splitType == ESplitHeuristic) {
	711	nCostRatio[axis] = EvalCostRatioHeuristic(
	712	leaf,
	713	axis,
	714	nPosition[axis],
	715	nRaysBack[axis],
	716	nRaysFront[axis],
	717	nPvsBack[axis],
	718	nPvsFront[axis]);
	719	} else
	720	if (splitType == ESplitHybrid) {
	721	if (leaf->depth > 7)
	722	nCostRatio[axis] = EvalCostRatioHeuristic(
	723	leaf,
	724	axis,
	725	nPosition[axis],
	726	nRaysBack[axis],
	727	nRaysFront[axis],
	728	nPvsBack[axis],
	729	nPvsFront[axis]);
	730	else {
	731	if (axis < 3)
	732	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
	733	else
	734	nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
[434]	735
[438]	736	nCostRatio[axis] = EvalCostRatio(leaf,
	737	axis,
	738	nPosition[axis],
	739	nRaysBack[axis],
	740	nRaysFront[axis],
	741	nPvsBack[axis],
	742	nPvsFront[axis]
	743	);
	744	}
	745	} else {
	746	cerr<<"VssTree: Unknown split heuristics\n";
	747	exit(1);
[434]	748	}
[438]	749
	750
	751	if ( bestAxis == -1)
[434]	752	bestAxis = axis;
[438]	753	else
	754	if ( nCostRatio[axis] < nCostRatio[bestAxis] )
	755	bestAxis = axis;
	756	}
[387]	757	}
[438]	758
[434]	759	axis = bestAxis;
	760	position = nPosition[bestAxis];
[387]	761
[434]	762	raysBack = nRaysBack[bestAxis];
	763	raysFront = nRaysFront[bestAxis];
[387]	764
[434]	765	pvsBack = nPvsBack[bestAxis];
	766	pvsFront = nPvsFront[bestAxis];
[387]	767
[434]	768	return nCostRatio[bestAxis];
[386]	769	}
	770
[434]	771
[386]	772	float
[434]	773	VssTree::EvalCostRatioHeuristic(
	774	VssTreeLeaf *leaf,
	775	const int axis,
	776	float &bestPosition,
	777	int &raysBack,
	778	int &raysFront,
	779	int &pvsBack,
	780	int &pvsFront
	781	)
[386]	782	{
[434]	783	AxisAlignedBox3 box;
	784	float minBox, maxBox;
[387]	785
[434]	786	if (axis < 3) {
	787	box = GetBBox(leaf);
	788	minBox = box.Min(axis);
	789	maxBox = box.Max(axis);
[438]	790	} else {
[434]	791	box = GetDirBBox(leaf);
	792	minBox = box.Min(axis-3);
	793	maxBox = box.Max(axis-3);
	794	}
	795
[1233]	796	SortSubdivisionCandidates(leaf, axis);
[372]	797
[434]	798	// go through the lists, count the number of objects left and right
	799	// and evaluate the following cost funcion:
	800	// C = ct_div_ci + (qlrl + qrrr)/queries
[382]	801
[434]	802	int rl=0, rr = leaf->rays.size();
	803	int pl=0, pr = leaf->GetPvsSize();
	804	float sizeBox = maxBox - minBox;
[387]	805
[434]	806	float minBand = minBox + 0.1*(maxBox - minBox);
	807	float maxBand = minBox + 0.9*(maxBox - minBox);
[387]	808
[434]	809	float minRatio = 1e20;
[387]	810
[434]	811	Intersectable::NewMail();
	812	// set all object as belonging to the fron pvs
	813	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	814	ri != leaf->rays.end();
	815	ri++)
	816	if ((*ri).mRay->IsActive()) {
	817	Intersectable object = (ri).mRay->mTerminationObject;
	818	if (object)
	819	if (!object->Mailed()) {
	820	object->Mail();
	821	object->mCounter = 1;
	822	} else
	823	object->mCounter++;
	824	}
[387]	825
[434]	826	Intersectable::NewMail();
[387]	827
[434]	828	for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
	829	ci < splitCandidates->end();
	830	ci++) {
	831	VssRay *ray;
	832	switch ((*ci).type) {
	833	case SortableEntry::ERayMin: {
	834	rl++;
	835	ray = (VssRay ) (ci).data;
	836	Intersectable *object = ray->mTerminationObject;
	837	if (object && !object->Mailed()) {
	838	object->Mail();
	839	pl++;
	840	}
	841	break;
	842	}
	843	case SortableEntry::ERayMax: {
	844	rr--;
	845	ray = (VssRay ) (ci).data;
	846	Intersectable *object = ray->mTerminationObject;
	847	if (object) {
	848	if (--object->mCounter == 0)
	849	pr--;
	850	}
	851	break;
	852	}
	853	}
	854
	855	float position = (*ci).value;
	856
	857	if (position > minBand && position < maxBand) {
[387]	858
[434]	859	float ratio = GetCostRatio(
	860	leaf,
	861	axis,
	862	position,
	863	rl,
	864	rr,
	865	pl,
	866	pr);
[387]	867
	868
[434]	869	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
	870	// cout<<"cost= "<<sum<<endl;
	871
	872	if (ratio < minRatio) {
	873	minRatio = ratio;
	874	bestPosition = position;
[382]	875
[434]	876	raysBack = rl;
	877	raysFront = rr;
[382]	878
[434]	879	pvsBack = pl;
	880	pvsFront = pr;
[387]	881
[382]	882	}
	883	}
[372]	884	}
[434]	885
[372]	886
[382]	887	// cout<<"===================="<<endl;
	888	// cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
	889	// <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
[434]	890	return minRatio;
[382]	891	}
	892
	893	void
[1233]	894	VssTree::SortSubdivisionCandidates(
[434]	895	VssTreeLeaf *node,
	896	const int axis
	897	)
[382]	898	{
[372]	899
[382]	900	splitCandidates->clear();
	901
	902	int requestedSize = 2*(node->rays.size());
	903	// creates a sorted split candidates array
	904	if (splitCandidates->capacity() > 500000 &&
	905	requestedSize < (int)(splitCandidates->capacity()/10) ) {
	906
	907	delete splitCandidates;
	908	splitCandidates = new vector<SortableEntry>;
[372]	909	}
	910
[382]	911	splitCandidates->reserve(requestedSize);
	912
	913	// insert all queries
	914	for(VssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
	915	ri < node->rays.end();
	916	ri++) {
[434]	917	if ((*ri).mRay->IsActive()) {
	918	if (axis < 3) {
	919	bool positive = (*ri).mRay->HasPosDir(axis);
	920	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
	921	SortableEntry::ERayMax,
	922	(*ri).ExtrapOrigin(axis),
	923	(void )(ri).mRay)
	924	);
	925
	926	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
	927	SortableEntry::ERayMin,
	928	(*ri).ExtrapTermination(axis),
	929	(void )(ri).mRay)
	930	);
	931	} else {
	932	float pos = (*ri).mRay->GetDirParametrization(axis-3);
	933	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
	934	pos - Limits::Small,
	935	(void )(ri).mRay)
	936	);
	937
	938	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMax,
	939	pos + Limits::Small,
	940	(void )(ri).mRay)
	941	);
	942	}
	943	}
[382]	944	}
[434]	945
[382]	946	stable_sort(splitCandidates->begin(), splitCandidates->end());
[372]	947	}
	948
	949
	950	void
	951	VssTree::EvaluateLeafStats(const TraversalData &data)
	952	{
	953
	954	// the node became a leaf -> evaluate stats for leafs
	955	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
	956
[382]	957	if (data.depth >= termMaxDepth)
[372]	958	stat.maxDepthNodes++;
	959
[434]	960	// if ( (int)(leaf->rays.size()) < termMinCost)
	961	// stat.minCostNodes++;
	962	if ( leaf->GetPvsSize() < termMinPvs)
	963	stat.minPvsNodes++;
[403]	964
[434]	965	if ( leaf->GetPvsSize() < termMinRays)
	966	stat.minRaysNodes++;
[403]	967
[434]	968	if (0 && leaf->GetAvgRayContribution() > termMaxRayContribution )
	969	stat.maxRayContribNodes++;
[382]	970
[434]	971	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
	972	stat.minSizeNodes++;
	973	}
[372]	974
[434]	975	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
[469]	976	stat.maxRayRefs = (int)leaf->rays.size();
[382]	977
[372]	978	}
	979
[427]	980	bool
	981	VssTree::TerminationCriteriaSatisfied(VssTreeLeaf *leaf)
	982	{
[434]	983	return ( (leaf->GetPvsSize() < termMinPvs) \|\|
	984	(leaf->rays.size() < termMinRays) \|\|
	985	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) \|\|
	986	(leaf->depth >= termMaxDepth) \|\|
[446]	987	(SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
	988	// \|\|
	989	// (mUseRss && leaf->mPassingRays == leaf->rays.size())
[434]	990	);
[427]	991	}
[372]	992
	993
	994	VssTreeNode *
	995	VssTree::SubdivideNode(
[434]	996	VssTreeLeaf *leaf,
	997	const AxisAlignedBox3 &box,
	998	AxisAlignedBox3 &backBBox,
	999	AxisAlignedBox3 &frontBBox
	1000	)
[372]	1001	{
	1002
[434]	1003	if (TerminationCriteriaSatisfied(leaf)) {
[382]	1004	#if 0
[434]	1005	if (leaf->depth >= termMaxDepth) {
	1006	cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
	1007	cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
	1008	}
[382]	1009	#endif
	1010
[434]	1011	return leaf;
	1012	}
[382]	1013
[434]	1014	float position;
[382]	1015
[434]	1016	// first count ray sides
[382]	1017	int raysBack;
	1018	int raysFront;
[434]	1019	int pvsBack;
	1020	int pvsFront;
[382]	1021
[372]	1022	// select subdivision axis
[382]	1023	int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
[434]	1024	// cout<<axis<<" ";
[427]	1025
[434]	1026	// cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
	1027	// pvsFront<<endl;
[382]	1028
[372]	1029	if (axis == -1) {
	1030	return leaf;
	1031	}
	1032
	1033	stat.nodes+=2;
	1034	stat.splits[axis]++;
	1035
	1036	// add the new nodes to the tree
	1037	VssTreeInterior *node = new VssTreeInterior(leaf->parent);
	1038
	1039	node->axis = axis;
	1040	node->position = position;
	1041	node->bbox = box;
	1042	node->dirBBox = GetDirBBox(leaf);
	1043
	1044	backBBox = box;
	1045	frontBBox = box;
	1046
	1047	VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
	1048	VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);
	1049
	1050	// replace a link from node's parent
	1051	if ( leaf->parent )
	1052	leaf->parent->ReplaceChildLink(leaf, node);
	1053	// and setup child links
	1054	node->SetupChildLinks(back, front);
[382]	1055
[372]	1056	if (axis <= VssTreeNode::SPLIT_Z) {
[434]	1057	backBBox.SetMax(axis, position);
[382]	1058	frontBBox.SetMin(axis, position);
	1059
[434]	1060	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1061	ri != leaf->rays.end();
	1062	ri++) {
[372]	1063	if ((*ri).mRay->IsActive()) {
	1064
[434]	1065	// first unref ray from the former leaf
	1066	(*ri).mRay->Unref();
[382]	1067
[434]	1068	// Debug << "computed t: " << (*ri).mRay->mT << endl;
	1069	// determine the side of this ray with respect to the plane
[463]	1070	float t;
	1071	int side = node->ComputeRayIntersection(*ri, t);
[434]	1072
	1073	if (side == 0) {
	1074	if ((*ri).mRay->HasPosDir(axis)) {
	1075	back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
	1076	(*ri).mMinT,
[463]	1077	t)
[434]	1078	);
	1079	front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
[463]	1080	t,
[434]	1081	(*ri).mMaxT));
	1082	} else {
	1083	back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
[463]	1084	t,
[434]	1085	(*ri).mMaxT));
	1086	front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
	1087	(*ri).mMinT,
[463]	1088	t));
[434]	1089	}
	1090	} else
	1091	if (side == 1)
	1092	front->AddRay(*ri);
	1093	else
	1094	back->AddRay(*ri);
[372]	1095	} else
[434]	1096	(*ri).mRay->Unref();
[372]	1097	}
	1098	} else {
	1099	// rays front/back
	1100
	1101
	1102	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	1103	ri != leaf->rays.end();
	1104	ri++) {
[372]	1105	if ((*ri).mRay->IsActive()) {
[434]	1106	// first unref ray from the former leaf
	1107	(*ri).mRay->Unref();
[382]	1108
[434]	1109	int side;
	1110	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
	1111	side = 1;
	1112	else
	1113	side = -1;
[372]	1114
[434]	1115	if (side == 1)
	1116	front->AddRay(*ri);
	1117	else
	1118	back->AddRay(*ri);
[372]	1119
	1120	} else
[434]	1121	(*ri).mRay->Unref();
[372]	1122	}
	1123	}
[434]	1124
	1125	front->SetPvsSize(pvsFront);
	1126	back->SetPvsSize(pvsBack);
[438]	1127	// compute entropy as well
	1128	front->ComputeEntropyImportance();
	1129	back->ComputeEntropyImportance();
	1130
[372]	1131	// update stats
[469]	1132	stat.rayRefs -= (int)leaf->rays.size();
[372]	1133	stat.rayRefs += raysBack + raysFront;
	1134
	1135
	1136	delete leaf;
	1137	return node;
	1138	}
	1139
	1140
	1141
	1142
	1143
	1144
	1145	int
	1146	VssTree::ReleaseMemory(const int time)
	1147	{
	1148	stack<VssTreeNode *> tstack;
	1149
	1150	// find a node in the tree which subtree will be collapsed
	1151	int maxAccessTime = time - accessTimeThreshold;
	1152	int released;
	1153
	1154	tstack.push(root);
	1155
	1156	while (!tstack.empty()) {
	1157	VssTreeNode *node = tstack.top();
	1158	tstack.pop();
	1159
	1160
	1161	if (!node->IsLeaf()) {
	1162	VssTreeInterior in = (VssTreeInterior )node;
	1163	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
	1164	if (in->depth >= minCollapseDepth &&
[434]	1165	in->lastAccessTime <= maxAccessTime) {
	1166	released = CollapseSubtree(node, time);
	1167	break;
[372]	1168	}
	1169
	1170	if (in->back->GetAccessTime() <
[434]	1171	in->front->GetAccessTime()) {
	1172	tstack.push(in->front);
	1173	tstack.push(in->back);
[372]	1174	} else {
[434]	1175	tstack.push(in->back);
	1176	tstack.push(in->front);
[372]	1177	}
	1178	}
	1179	}
	1180
	1181	while (tstack.empty()) {
	1182	// could find node to collaps...
	1183	// cout<<"Could not find a node to release "<<endl;
	1184	break;
	1185	}
	1186
	1187	return released;
	1188	}
	1189
	1190
	1191
	1192
	1193	VssTreeNode *
	1194	VssTree::SubdivideLeaf(
[434]	1195	VssTreeLeaf *leaf
	1196	)
[372]	1197	{
	1198	VssTreeNode *node = leaf;
[427]	1199
[372]	1200	AxisAlignedBox3 leafBBox = GetBBox(leaf);
	1201
[434]	1202	static int pass = 0;
	1203	pass ++;
[372]	1204
	1205	// check if we should perform a dynamic subdivision of the leaf
[434]	1206	if (!TerminationCriteriaSatisfied(leaf)) {
[372]	1207
[434]	1208	// memory check and realese...
[372]	1209	if (GetMemUsage() > maxTotalMemory) {
	1210	ReleaseMemory( pass );
	1211	}
	1212
	1213	AxisAlignedBox3 backBBox, frontBBox;
	1214
	1215	// subdivide the node
	1216	node =
	1217	SubdivideNode(leaf,
[434]	1218	leafBBox,
	1219	backBBox,
	1220	frontBBox
	1221	);
[372]	1222	}
[427]	1223
[372]	1224	return node;
	1225	}
	1226
	1227
	1228
	1229	void
	1230	VssTree::UpdateRays(VssRayContainer &remove,
[434]	1231	VssRayContainer &add
	1232	)
[372]	1233	{
	1234	VssTreeLeaf::NewMail();
	1235
	1236	// schedule rays for removal
	1237	for(VssRayContainer::const_iterator ri = remove.begin();
	1238	ri != remove.end();
	1239	ri++) {
	1240	(*ri)->ScheduleForRemoval();
	1241	}
	1242
	1243	int inactive=0;
	1244
	1245	for(VssRayContainer::const_iterator ri = remove.begin();
	1246	ri != remove.end();
	1247	ri++) {
	1248	if ((*ri)->ScheduledForRemoval())
[434]	1249	// RemoveRay(*ri, NULL, false);
	1250	// !!! BUG - with true it does not work correctly - aggreated delete
[372]	1251	RemoveRay(*ri, NULL, true);
	1252	else
	1253	inactive++;
	1254	}
	1255
	1256
	1257	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
	1258
	1259	for(VssRayContainer::const_iterator ri = add.begin();
	1260	ri != add.end();
	1261	ri++) {
[434]	1262	VssTreeNode::RayInfo info(*ri);
	1263	if (ClipRay(info, bbox))
	1264	AddRay(info);
[372]	1265	}
	1266	}
	1267
	1268
	1269	void
	1270	VssTree::RemoveRay(VssRay *ray,
[434]	1271	vector<VssTreeLeaf > affectedLeaves,
	1272	const bool removeAllScheduledRays
	1273	)
[372]	1274	{
	1275
	1276	stack<RayTraversalData> tstack;
	1277
	1278	tstack.push(RayTraversalData(root, VssTreeNode::RayInfo(ray)));
	1279
	1280	RayTraversalData data;
	1281
	1282	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
	1283
	1284	while (!tstack.empty()) {
	1285	data = tstack.top();
	1286	tstack.pop();
	1287
	1288	if (!data.node->IsLeaf()) {
	1289	// split the set of rays in two groups intersecting the
	1290	// two subtrees
	1291
	1292	TraverseInternalNode(data, tstack);
	1293
	1294	} else {
	1295	// remove the ray from the leaf
	1296	// find the ray in the leaf and swap it with the last ray...
	1297	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
	1298
	1299	if (!leaf->Mailed()) {
[434]	1300	leaf->Mail();
	1301	if (affectedLeaves)
	1302	affectedLeaves->push_back(leaf);
[372]	1303
[434]	1304	if (removeAllScheduledRays) {
	1305	int tail = leaf->rays.size()-1;
[372]	1306
[434]	1307	for (int i=0; i < (int)(leaf->rays.size()); i++) {
	1308	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
	1309	// find a ray to replace it with
	1310	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
	1311	stat.removedRayRefs++;
	1312	leaf->rays[tail].mRay->Unref();
	1313	leaf->rays.pop_back();
	1314	tail--;
	1315	}
[372]	1316
[434]	1317	if (tail < i)
	1318	break;
[372]	1319
[434]	1320	stat.removedRayRefs++;
	1321	leaf->rays[i].mRay->Unref();
	1322	leaf->rays[i] = leaf->rays[tail];
	1323	leaf->rays.pop_back();
	1324	tail--;
	1325	}
	1326	}
	1327	}
[372]	1328	}
	1329
	1330	if (!removeAllScheduledRays)
[434]	1331	for (int i=0; i < (int)leaf->rays.size(); i++) {
	1332	if (leaf->rays[i].mRay == ray) {
	1333	stat.removedRayRefs++;
	1334	ray->Unref();
	1335	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
	1336	leaf->rays.pop_back();
	1337	// check this ray again
	1338	break;
	1339	}
	1340	}
[372]	1341
	1342	}
	1343	}
	1344
	1345	if (ray->RefCount() != 0) {
	1346	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
	1347	exit(1);
	1348	}
	1349
	1350	}
	1351
	1352
	1353	void
[427]	1354	VssTree::AddRay(VssTreeNode::RayInfo &info)
[372]	1355	{
	1356
	1357	stack<RayTraversalData> tstack;
	1358
[427]	1359	tstack.push(RayTraversalData(root, info));
[372]	1360
	1361	RayTraversalData data;
	1362
	1363	while (!tstack.empty()) {
	1364	data = tstack.top();
	1365	tstack.pop();
	1366
	1367	if (!data.node->IsLeaf()) {
	1368	TraverseInternalNode(data, tstack);
	1369	} else {
	1370	// remove the ray from the leaf
	1371	// find the ray in the leaf and swap it with the last ray...
	1372	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
	1373	leaf->AddRay(data.rayData);
	1374	stat.addedRayRefs++;
	1375	}
	1376	}
	1377	}
	1378
	1379	void
	1380	VssTree::TraverseInternalNode(
[434]	1381	RayTraversalData &data,
	1382	stack<RayTraversalData> &tstack)
[372]	1383	{
	1384	VssTreeInterior in = (VssTreeInterior ) data.node;
	1385
	1386	if (in->axis <= VssTreeNode::SPLIT_Z) {
[463]	1387	float t;
[372]	1388	// determine the side of this ray with respect to the plane
	1389	int side = in->ComputeRayIntersection(data.rayData,
[463]	1390	t);
[372]	1391
	1392
	1393	if (side == 0) {
	1394	if (data.rayData.mRay->HasPosDir(in->axis)) {
[434]	1395	tstack.push(RayTraversalData(in->back,
	1396	VssTreeNode::RayInfo(data.rayData.mRay,
	1397	data.rayData.mMinT,
[463]	1398	t))
[434]	1399	);
[372]	1400
[434]	1401	tstack.push(RayTraversalData(in->front,
	1402	VssTreeNode::RayInfo(data.rayData.mRay,
[463]	1403	t,
[434]	1404	data.rayData.mMaxT
	1405	))
	1406	);
[372]	1407
	1408	} else {
[434]	1409	tstack.push(RayTraversalData(in->back,
	1410	VssTreeNode::RayInfo(data.rayData.mRay,
[463]	1411	t,
[434]	1412	data.rayData.mMaxT
	1413	))
	1414	);
[372]	1415
[434]	1416	tstack.push(RayTraversalData(in->front,
	1417	VssTreeNode::RayInfo(data.rayData.mRay,
	1418	data.rayData.mMinT,
[463]	1419	t))
[434]	1420	);
[372]	1421
	1422
	1423	}
	1424	} else
	1425	if (side == 1)
[434]	1426	tstack.push(RayTraversalData(in->front, data.rayData));
[372]	1427	else
[434]	1428	tstack.push(RayTraversalData(in->back, data.rayData));
[372]	1429	}
	1430	else {
	1431	// directional split
[434]	1432	if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position)
	1433	tstack.push(RayTraversalData(in->front, data.rayData));
	1434	else
	1435	tstack.push(RayTraversalData(in->back, data.rayData));
[372]	1436	}
	1437	}
	1438
	1439
	1440	int
	1441	VssTree::CollapseSubtree(VssTreeNode *sroot, const int time)
	1442	{
	1443	// first count all rays in the subtree
	1444	// use mail 1 for this purpose
	1445	stack<VssTreeNode *> tstack;
	1446	int rayCount = 0;
	1447	int totalRayCount = 0;
	1448	int collapsedNodes = 0;
	1449
	1450	#if DEBUG_COLLAPSE
	1451	cout<<"Collapsing subtree"<<endl;
	1452	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
	1453	cout<<"depth="<<(int)sroot->depth<<endl;
	1454	#endif
	1455
[434]	1456	// tstat.collapsedSubtrees++;
	1457	// tstat.collapseDepths += (int)sroot->depth;
	1458	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
[372]	1459
	1460	tstack.push(sroot);
	1461	VssRay::NewMail();
	1462
	1463	while (!tstack.empty()) {
	1464	collapsedNodes++;
	1465	VssTreeNode *node = tstack.top();
	1466	tstack.pop();
	1467
	1468	if (node->IsLeaf()) {
	1469	VssTreeLeaf leaf = (VssTreeLeaf ) node;
	1470	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	1471	ri != leaf->rays.end();
	1472	ri++) {
[372]	1473
[434]	1474	totalRayCount++;
	1475	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
	1476	(*ri).mRay->Mail();
	1477	rayCount++;
	1478	}
[372]	1479	}
	1480	} else {
	1481	tstack.push(((VssTreeInterior *)node)->back);
	1482	tstack.push(((VssTreeInterior *)node)->front);
	1483	}
	1484	}
	1485
	1486	VssRay::NewMail();
	1487
	1488	// create a new node that will hold the rays
	1489	VssTreeLeaf *newLeaf = new VssTreeLeaf( sroot->parent, rayCount );
	1490	if ( newLeaf->parent )
	1491	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
	1492
	1493
	1494	tstack.push( sroot );
	1495
	1496	while (!tstack.empty()) {
	1497
	1498	VssTreeNode *node = tstack.top();
	1499	tstack.pop();
	1500
	1501	if (node->IsLeaf()) {
	1502	VssTreeLeaf leaf = (VssTreeLeaf ) node;
	1503
	1504	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
[434]	1505	ri != leaf->rays.end();
	1506	ri++) {
[372]	1507
[434]	1508	// unref this ray from the old node
[372]	1509
[434]	1510	if ((*ri).mRay->IsActive()) {
	1511	(*ri).mRay->Unref();
	1512	if (!(*ri).mRay->Mailed()) {
	1513	(*ri).mRay->Mail();
	1514	newLeaf->AddRay(*ri);
	1515	}
	1516	} else
	1517	(*ri).mRay->Unref();
[372]	1518
	1519	}
	1520	} else {
	1521	tstack.push(((VssTreeInterior *)node)->back);
	1522	tstack.push(((VssTreeInterior *)node)->front);
	1523	}
	1524	}
	1525
	1526	// delete the node and all its children
	1527	delete sroot;
	1528
[434]	1529	// for(VssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
	1530	// ri != newLeaf->rays.end();
	1531	// ri++)
	1532	// (*ri).ray->UnMail(2);
[372]	1533
	1534
	1535	#if DEBUG_COLLAPSE
	1536	cout<<"Total memory before="<<GetMemUsage()<<endl;
	1537	#endif
	1538
	1539	stat.nodes -= collapsedNodes - 1;
	1540	stat.rayRefs -= totalRayCount - rayCount;
	1541
	1542	#if DEBUG_COLLAPSE
	1543	cout<<"collapsed nodes"<<collapsedNodes<<endl;
	1544	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
	1545	cout<<"Total memory after="<<GetMemUsage()<<endl;
	1546	cout<<"================================"<<endl;
	1547	#endif
	1548
[434]	1549	// tstat.collapsedNodes += collapsedNodes;
	1550	// tstat.collapsedRays += totalRayCount - rayCount;
[372]	1551
	1552	return totalRayCount - rayCount;
	1553	}
[387]	1554
	1555
	1556	int
[434]	1557	VssTree::GetPvsSize(const AxisAlignedBox3 &box) const
[387]	1558	{
[434]	1559	stack<VssTreeNode *> tstack;
[387]	1560	tstack.push(root);
	1561
[434]	1562	Intersectable::NewMail();
	1563	int pvsSize = 0;
[387]	1564
	1565	while (!tstack.empty()) {
	1566	VssTreeNode *node = tstack.top();
	1567	tstack.pop();
	1568
	1569
	1570	if (node->IsLeaf()) {
[434]	1571	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1572	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1573	ri != leaf->rays.end();
	1574	ri++)
	1575	if ((*ri).mRay->IsActive()) {
	1576	Intersectable *object;
[387]	1577	#if BIDIRECTIONAL_RAY
[434]	1578	object = (*ri).mRay->mOriginObject;
	1579	if (object && !object->Mailed()) {
	1580	pvsSize++;
	1581	object->Mail();
	1582	}
[387]	1583	#endif
[434]	1584	object = (*ri).mRay->mTerminationObject;
	1585	if (object && !object->Mailed()) {
	1586	pvsSize++;
	1587	object->Mail();
	1588	}
	1589	}
	1590	} else {
	1591	VssTreeInterior in = (VssTreeInterior )node;
	1592	if (in->axis < 3) {
	1593	if (box.Max(in->axis) >= in->position )
	1594	tstack.push(in->front);
[387]	1595
[434]	1596	if (box.Min(in->axis) <= in->position )
	1597	tstack.push(in->back);
	1598	} else {
	1599	// both nodes for directional splits
	1600	tstack.push(in->front);
	1601	tstack.push(in->back);
	1602	}
[387]	1603	}
[434]	1604	}
	1605	return pvsSize;
[387]	1606	}
[401]	1607
	1608	void
	1609	VssTree::GetRayContributionStatistics(
[434]	1610	float &minRayContribution,
	1611	float &maxRayContribution,
	1612	float &avgRayContribution
	1613	)
[401]	1614	{
[434]	1615	stack<VssTreeNode *> tstack;
[401]	1616	tstack.push(root);
	1617
[434]	1618	minRayContribution = 1.0f;
	1619	maxRayContribution = 0.0f;
	1620	float sumRayContribution = 0.0f;
	1621	int leaves = 0;
[401]	1622
	1623	while (!tstack.empty()) {
	1624	VssTreeNode *node = tstack.top();
	1625	tstack.pop();
	1626
	1627	if (node->IsLeaf()) {
[434]	1628	leaves++;
	1629	VssTreeLeaf leaf = (VssTreeLeaf )node;
[435]	1630	float c = leaf->GetImportance();
[434]	1631	if (c > maxRayContribution)
	1632	maxRayContribution = c;
	1633	if (c < minRayContribution)
	1634	minRayContribution = c;
	1635	sumRayContribution += c;
[401]	1636
[434]	1637	} else {
	1638	VssTreeInterior in = (VssTreeInterior )node;
	1639	// both nodes for directional splits
	1640	tstack.push(in->front);
	1641	tstack.push(in->back);
[401]	1642	}
[434]	1643	}
[401]	1644
[434]	1645	cout<<"sum="<<sumRayContribution<<endl;
	1646	cout<<"leaves="<<leaves<<endl;
	1647	avgRayContribution = sumRayContribution/(float)leaves;
[401]	1648	}
	1649
	1650
	1651	int
	1652	VssTree::GenerateRays(const float ratioPerLeaf,
[434]	1653	SimpleRayContainer &rays)
[401]	1654	{
[434]	1655	stack<VssTreeNode *> tstack;
[401]	1656	tstack.push(root);
	1657
	1658	while (!tstack.empty()) {
	1659	VssTreeNode *node = tstack.top();
	1660	tstack.pop();
	1661
	1662	if (node->IsLeaf()) {
[434]	1663	VssTreeLeaf leaf = (VssTreeLeaf )node;
[435]	1664	float c = leaf->GetImportance();
[434]	1665	int num = (c*ratioPerLeaf + 0.5);
	1666	// cout<<num<<" ";
[401]	1667
[434]	1668	for (int i=0; i < num; i++) {
	1669	Vector3 origin = GetBBox(leaf).GetRandomPoint();
	1670	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
	1671	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
	1672	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	1673	rays.push_back(SimpleRay(origin, direction));
	1674	}
[401]	1675
[434]	1676	} else {
	1677	VssTreeInterior in = (VssTreeInterior )node;
	1678	// both nodes for directional splits
	1679	tstack.push(in->front);
	1680	tstack.push(in->back);
[401]	1681	}
[434]	1682	}
[401]	1683
[434]	1684	return rays.size();
[401]	1685	}
	1686
[427]	1687	void
	1688	VssTree::CollectLeaves(vector<VssTreeLeaf *> &leaves)
	1689	{
[434]	1690	stack<VssTreeNode *> tstack;
[427]	1691	tstack.push(root);
	1692
	1693	while (!tstack.empty()) {
	1694	VssTreeNode *node = tstack.top();
	1695	tstack.pop();
	1696
	1697	if (node->IsLeaf()) {
[434]	1698	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1699	leaves.push_back(leaf);
	1700	} else {
	1701	VssTreeInterior in = (VssTreeInterior )node;
	1702	// both nodes for directional splits
	1703	tstack.push(in->front);
	1704	tstack.push(in->back);
[427]	1705	}
[434]	1706	}
[427]	1707	}
[401]	1708
[434]	1709	bool
	1710	VssTree::ValidLeaf(VssTreeLeaf *leaf) const
	1711	{
	1712	return leaf->rays.size() > termMinRays/4;
	1713	}
	1714
	1715
	1716	void
	1717	VssTree::GenerateLeafRays(VssTreeLeaf *leaf,
	1718	const int numberOfRays,
	1719	SimpleRayContainer &rays)
	1720	{
[469]	1721	int nrays = (int)leaf->rays.size();
[434]	1722	for (int i=0; i < numberOfRays; i++) {
	1723	// pickup 3 random rays
[473]	1724	int r1 = (int)RandomValue(0, nrays-1);
	1725	int r2 = (int)RandomValue(0, nrays-1);
	1726	int r3 = (int)RandomValue(0, nrays-1);
[434]	1727
	1728	Vector3 o1 = leaf->rays[r1].Extrap(RandomValue(leaf->rays[r1].GetMinT(),
	1729	leaf->rays[r1].GetMaxT()));
	1730
	1731	Vector3 o2 = leaf->rays[r2].Extrap(RandomValue(leaf->rays[r2].GetMinT(),
	1732	leaf->rays[r2].GetMaxT()));
	1733
	1734	Vector3 o3 = leaf->rays[r3].Extrap(RandomValue(leaf->rays[r3].GetMinT(),
	1735	leaf->rays[r3].GetMaxT()));
	1736
[469]	1737	const float overlap = 0.1f;
[434]	1738
	1739	Vector3 origin, direction;
	1740	bool useExtendedConvexCombination = true;
	1741	if (useExtendedConvexCombination) {
	1742	float w1, w2, w3;
	1743	GenerateExtendedConvexCombinationWeights(w1, w2, w3, overlap);
	1744	origin = w1o1 + w2o2 + w3*o3;
	1745	direction =
	1746	w1*leaf->rays[r1].mRay->GetDir() +
	1747	w2*leaf->rays[r2].mRay->GetDir() +
	1748	w3*leaf->rays[r3].mRay->GetDir();
	1749	} else {
	1750	origin = GetBBox(leaf).GetRandomPoint();
	1751	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
	1752	direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
	1753	}
	1754	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	1755	rays.push_back(SimpleRay(origin, direction));
	1756	}
	1757	}
	1758
[427]	1759	int
	1760	VssTree::GenerateRays(const int numberOfRays,
[434]	1761	const int numberOfLeaves,
	1762	SimpleRayContainer &rays)
[427]	1763	{
	1764
[434]	1765	vector<VssTreeLeaf *> leaves;
[427]	1766
[434]	1767	CollectLeaves(leaves);
[427]	1768
[434]	1769	sort(leaves.begin(),
	1770	leaves.end(),
	1771	GreaterContribution);
[427]	1772
	1773
[434]	1774	float sumContrib = 0.0;
	1775	int i;
	1776	int k = 0;
	1777	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	1778	if (ValidLeaf(leaves[i])) {
[435]	1779	float c = leaves[i]->GetImportance();
[434]	1780	sumContrib += c;
	1781	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
	1782	k++;
[427]	1783	}
[434]	1784
	1785	float avgContrib = sumContrib/numberOfLeaves;
	1786	float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
	1787	k = 0;
	1788	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	1789	if (ValidLeaf(leaves[i])) {
	1790	k++;
	1791	VssTreeLeaf *leaf = leaves[i];
[435]	1792	float c = leaf->GetImportance();
[469]	1793	int num = (int)(c*ratioPerLeaf + 0.5f);
[434]	1794	GenerateLeafRays(leaf, num, rays);
	1795	}
[427]	1796
[469]	1797	return (int)rays.size();
[427]	1798	}
	1799
	1800
[401]	1801	float
	1802	VssTree::GetAvgPvsSize()
	1803	{
[434]	1804	stack<VssTreeNode *> tstack;
[401]	1805	tstack.push(root);
	1806
[434]	1807	int sumPvs = 0;
	1808	int leaves = 0;
[401]	1809	while (!tstack.empty()) {
	1810	VssTreeNode *node = tstack.top();
	1811	tstack.pop();
	1812
	1813	if (node->IsLeaf()) {
[434]	1814	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1815	// update pvs size
	1816	leaf->UpdatePvsSize();
	1817	sumPvs += leaf->GetPvsSize();
	1818	leaves++;
	1819	} else {
	1820	VssTreeInterior in = (VssTreeInterior )node;
	1821	// both nodes for directional splits
	1822	tstack.push(in->front);
	1823	tstack.push(in->back);
[401]	1824	}
[434]	1825	}
[401]	1826
	1827
[434]	1828	return sumPvs/(float)leaves;
[401]	1829	}
[427]	1830
	1831
[435]	1832	float
	1833	VssTreeLeaf::GetImportance() const
	1834	{
[438]	1835
	1836	if (1) {
	1837	// return GetAvgRayContribution();
[863]	1838	return (float)GetPvsSize();
[438]	1839	} else {
[446]	1840	// return GetAvgRayContribution()*mEntropyImportance;
	1841	//return GetAvgRayContribution();
	1842	return mEntropyImportance;
[438]	1843	}
[435]	1844	}
[438]	1845
	1846
	1847	float
	1848	VssTreeLeaf::ComputePvsEntropy()
	1849	{
	1850	int samples = 0;
	1851	Intersectable::NewMail();
	1852	// set all object as belonging to the fron pvs
	1853	for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
	1854	ri != rays.end();
	1855	ri++)
	1856	if ((*ri).mRay->IsActive()) {
	1857	Intersectable object = (ri).mRay->mTerminationObject;
	1858	if (object) {
	1859	if (!object->Mailed()) {
	1860	object->Mail();
	1861	object->mCounter = 1;
	1862	} else
	1863	object->mCounter++;
	1864	samples++;
	1865	}
	1866	}
	1867
	1868	float entropy = 0.0f;
	1869
	1870	if (samples > 1) {
	1871	Intersectable::NewMail();
	1872	for(RayInfoContainer::const_iterator ri = rays.begin();
	1873	ri != rays.end();
	1874	ri++)
	1875	if ((*ri).mRay->IsActive()) {
	1876	Intersectable object = (ri).mRay->mTerminationObject;
	1877	if (object) {
	1878	if (!object->Mailed()) {
	1879	object->Mail();
	1880	float p = object->mCounter/(float)samples;
	1881	entropy -= p*log(p);
	1882	}
	1883	}
	1884	}
	1885	entropy = entropy/log((float)samples);
	1886	}
	1887	else
	1888	entropy = 1.0f;
	1889
	1890	return entropy;
	1891	}
	1892
	1893	float
	1894	VssTreeLeaf::ComputeRayLengthEntropy()
	1895	{
	1896	// get sum of all ray lengths
	1897	// consider only passing rays or originating rays
	1898	float sum = 0.0f;
	1899	int samples = 0;
	1900
	1901	for(RayInfoContainer::const_iterator ri = rays.begin();
	1902	ri != rays.end();
	1903	ri++) {
	1904	int rayClass = (*ri).GetRayClass();
	1905	if (
	1906	rayClass == RayInfo::PASSING_RAY
	1907	// \|\|
	1908	// rayClass == RayInfo::SOURCE_RAY
	1909	// rayClass == RayInfo::TERMINATION_RAY
	1910	) {
	1911	float c = 1.0f - (*ri).GetMaxT();
	1912	sum += (ri).mRay->GetSize()c;
	1913	samples++;
	1914	}
	1915	}
	1916
	1917	float entropy = 0.0f;
	1918
	1919	if (samples > 1) {
	1920	for(RayInfoContainer::const_iterator ri = rays.begin();
	1921	ri != rays.end();
	1922	ri++) {
	1923	int rayClass = (*ri).GetRayClass();
	1924	if (
	1925	rayClass == RayInfo::PASSING_RAY
	1926	// \|\|
	1927	// rayClass == RayInfo::SOURCE_RAY
	1928	// rayClass == RayInfo::TERMINATION_RAY
	1929	) {
	1930	float c = 1.0f - (*ri).GetMaxT();
	1931	float p = (ri).mRay->GetSize()c/sum;
	1932	entropy -= p*log(p);
	1933	}
	1934	}
	1935	entropy = entropy/log((float)samples);
	1936	} else
	1937	entropy = 1.0f;
	1938
	1939	return entropy;
	1940	}
	1941
	1942	float
	1943	VssTreeLeaf::ComputeRayTerminationEntropy()
	1944	{
	1945
	1946	return 0.0f;
	1947	}
	1948
	1949
	1950	void
	1951	VssTreeLeaf::ComputeEntropyImportance()
	1952	{
	1953	// mEntropy = 1.0f - ComputeRayLengthEntropy();
	1954	mEntropyImportance = 1.0f - ComputePvsEntropy();
	1955
	1956	// cout<<"ei="<<mEntropyImportance<<" ";
	1957	}
[466]	1958
	1959
	1960	int
	1961	VssTree::CollectRays(VssRayContainer &rays,
	1962	const int number)
	1963	{
	1964	VssRayContainer allRays;
	1965	CollectRays(allRays);
	1966
	1967	int desired = min(number, (int)allRays.size());
	1968	float prob = desired/(float)allRays.size();
	1969	while (rays.size() < desired) {
	1970	VssRayContainer::const_iterator it = allRays.begin();
	1971	for (; it != allRays.end() && rays.size() < desired; it++) {
	1972	if (Random(1.0f) < prob)
	1973	rays.push_back(*it);
	1974	}
	1975	}
	1976	return rays.size();
	1977	}
	1978
	1979
	1980	int
	1981	VssTree::CollectRays(VssRayContainer &rays
	1982	)
	1983	{
	1984	VssRay::NewMail();
	1985
	1986	stack<VssTreeNode *> tstack;
	1987	tstack.push(root);
	1988
	1989	while (!tstack.empty()) {
	1990	VssTreeNode *node = tstack.top();
	1991	tstack.pop();
	1992
	1993	if (node->IsLeaf()) {
	1994	VssTreeLeaf leaf = (VssTreeLeaf )node;
	1995	// update pvs size
	1996	VssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
	1997	for (;it != leaf->rays.end(); ++it)
	1998	if (!(*it).mRay->Mailed()) {
	1999	(*it).mRay->Mail();
	2000	rays.push_back((*it).mRay);
	2001	}
	2002	} else {
	2003	VssTreeInterior in = (VssTreeInterior )node;
	2004	// both nodes for directional splits
	2005	tstack.push(in->front);
	2006	tstack.push(in->back);
	2007	}
	2008	}
	2009
	2010	return rays.size();
	2011	}
[860]	2012
	2013	}

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

Download in other formats: