Context Navigation

source: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp @ 483

Revision 483, 54.8 KB checked in by mattausch, 19 years ago (diff)

Rev	Line
[408]	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
[406]	14	#include <stack>
[408]	15	#include <queue>
[406]	16	#include <algorithm>
[408]	17	#include <fstream>
	18	#include <string>
	19
	20	#include "VspKdTree.h"
	21
[406]	22	#include "Environment.h"
[408]	23	#include "VssRay.h"
	24	#include "Intersectable.h"
	25	#include "Ray.h"
[420]	26	#include "RayInfo.h"
[453]	27	#include "ViewCellsManager.h"
[462]	28	#include "ViewCellBsp.h"
[406]	29
[418]	30	// Static variables
[462]	31	int VspKdLeaf::sMailId = 0;
	32	int MergeCandidate::sMaxPvsSize = 150;
[472]	33	float MergeCandidate::sOverallCost = Limits::Small;
[428]	34	#define USE_FIXEDPOINT_T 0
	35
[420]	36	/// Adds object to the pvs of the front and back node
[409]	37	inline void AddObject2Pvs(Intersectable *object,
	38	const int side,
	39	int &pvsBack,
	40	int &pvsFront)
[406]	41	{
[408]	42	if (!object)
	43	return;
[465]	44
	45	if (side <= 0)
[409]	46	{
[465]	47	if (!object->Mailed() && !object->Mailed(2))
[409]	48	{
	49	++ pvsBack;
	50
[408]	51	if (object->Mailed(1))
	52	object->Mail(2);
	53	else
	54	object->Mail();
	55	}
	56	}
[465]	57
[428]	58	if (side >= 0)
[409]	59	{
[465]	60	if (!object->Mailed(1) && !object->Mailed(2))
[409]	61	{
	62	++ pvsFront;
[428]	63
[408]	64	if (object->Mailed())
	65	object->Mail(2);
	66	else
	67	object->Mail(1);
	68	}
	69	}
[406]	70	}
	71
[418]	72	/**************************************************************/
[462]	73	/* class VspKdNode implementation */
[420]	74	/**************************************************************/
	75
	76	// Inline constructor
[462]	77	VspKdNode::VspKdNode(VspKdInterior *p):
[465]	78	mParent(p), mAxis(-1), mDepth(p ? p->mDepth + 1 : 0)
[420]	79	{}
	80
[465]	81	VspKdNode::~VspKdNode()
[480]	82	{
	83	};
[420]	84
[462]	85	inline VspKdInterior *VspKdNode::GetParent() const
[420]	86	{
	87	return mParent;
	88	}
	89
[462]	90	inline void VspKdNode::SetParent(VspKdInterior *p)
[420]	91	{
	92	mParent = p;
	93	}
	94
[465]	95	bool VspKdNode::IsLeaf() const
	96	{
	97	return mAxis == -1;
[420]	98	}
[465]	99
	100	int VspKdNode::GetAccessTime()
[420]	101	{
	102	return 0x7FFFFFF;
	103	}
	104
	105	/**************************************************************/
[462]	106	/* VspKdInterior implementation */
[418]	107	/**************************************************************/
[408]	108
[462]	109	VspKdInterior::VspKdInterior(VspKdInterior *p):
	110	VspKdNode(p), mBack(NULL), mFront(NULL), mAccesses(0), mLastAccessTime(-1)
[418]	111	{
	112	}
	113
[462]	114	int VspKdInterior::GetAccessTime()
[418]	115	{
[419]	116	return mLastAccessTime;
[418]	117	}
	118
[465]	119	void VspKdInterior::SetupChildLinks(VspKdNode b, VspKdNode f)
[418]	120	{
	121	mBack = b;
	122	mFront = f;
[420]	123	b->SetParent(this);
	124	f->SetParent(this);
[418]	125	}
	126
[465]	127	void VspKdInterior::ReplaceChildLink(VspKdNode *oldChild,
[462]	128	VspKdNode *newChild)
[418]	129	{
	130	if (mBack == oldChild)
	131	mBack = newChild;
	132	else
	133	mFront = newChild;
	134	}
	135
[465]	136	int VspKdInterior::Type() const
	137	{
	138	return EInterior;
[418]	139	}
[465]	140
[462]	141	VspKdInterior::~VspKdInterior()
[418]	142	{
	143	DEL_PTR(mBack);
	144	DEL_PTR(mFront);
	145	}
[465]	146
	147	void VspKdInterior::Print(ostream &s) const
[418]	148	{
[423]	149	switch (mAxis)
	150	{
	151	case 0:
	152	s << "x "; break;
	153	case 1:
	154	s << "y "; break;
	155	case 2:
	156	s << "z "; break;
	157	}
	158
[419]	159	s << mPosition << " ";
[418]	160
	161	mBack->Print(s);
	162	mFront->Print(s);
	163	}
[465]	164
	165	int VspKdInterior::ComputeRayIntersection(const RayInfo &rayData, float &t)
[418]	166	{
[419]	167	return rayData.ComputeRayIntersection(mAxis, mPosition, t);
[418]	168	}
	169
[462]	170	VspKdNode *VspKdInterior::GetBack() const
[419]	171	{
	172	return mBack;
	173	}
[418]	174
[462]	175	VspKdNode *VspKdInterior::GetFront() const
[406]	176	{
[419]	177	return mFront;
	178	}
[406]	179
[420]	180
[425]	181	/**************************************************************/
[462]	182	/* class VspKdLeaf implementation */
[425]	183	/**************************************************************/
[462]	184
	185
[472]	186	VspKdLeaf::VspKdLeaf(VspKdInterior *p,
	187	const int nRays,
	188	const int maxCostMisses):
	189	VspKdNode(p), mRays(), mPvsSize(0), mValidPvs(false), mViewCell(NULL),
	190	mMaxCostMisses(maxCostMisses)
[420]	191	{
	192	mRays.reserve(nRays);
	193	}
	194
[465]	195	VspKdLeaf::~VspKdLeaf()
[462]	196	{
	197	}
[420]	198
[472]	199
	200	int VspKdLeaf::GetMaxCostMisses()
	201	{
	202	return mMaxCostMisses;
	203	}
	204
	205
[465]	206	int VspKdLeaf::Type() const
	207	{
	208	return ELeaf;
[420]	209	}
	210
[465]	211	void VspKdLeaf::Print(ostream &s) const
[420]	212	{
	213	s << endl << "L: r = " << (int)mRays.size() << endl;
	214	};
	215
[465]	216	void VspKdLeaf::AddRay(const RayInfo &data)
[420]	217	{
	218	mValidPvs = false;
	219	mRays.push_back(data);
	220	data.mRay->Ref();
	221	}
	222
[465]	223	int VspKdLeaf::GetPvsSize() const
[420]	224	{
	225	return mPvsSize;
	226	}
	227
[465]	228	void VspKdLeaf::SetPvsSize(const int s)
[420]	229	{
	230	mPvsSize = s;
	231	}
	232
[462]	233	void VspKdLeaf::Mail()
[465]	234	{
	235	mMailbox = sMailId;
[420]	236	}
	237
[465]	238	void VspKdLeaf::NewMail()
	239	{
	240	++ sMailId;
[420]	241	}
	242
[465]	243	bool VspKdLeaf::Mailed() const
	244	{
	245	return mMailbox == sMailId;
[420]	246	}
	247
[462]	248	bool VspKdLeaf::Mailed(const int mail) const
[420]	249	{
[462]	250	return mMailbox == mail;
[420]	251	}
	252
[462]	253	int VspKdLeaf::GetMailbox() const
[420]	254	{
[462]	255	return mMailbox;
	256	}
	257
[465]	258	float VspKdLeaf::GetAvgRayContribution() const
[462]	259	{
[420]	260	return GetPvsSize() / ((float)mRays.size() + Limits::Small);
	261	}
	262
[465]	263	float VspKdLeaf::GetSqrRayContribution() const
[420]	264	{
	265	return sqr(GetAvgRayContribution());
	266	}
	267
[462]	268	RayInfoContainer &VspKdLeaf::GetRays()
[426]	269	{
	270	return mRays;
	271	}
[428]	272
[462]	273	void VspKdLeaf::ExtractPvs(ObjectContainer &objects) const
[428]	274	{
	275	RayInfoContainer::const_iterator it, it_end = mRays.end();
	276
	277	for (it = mRays.begin(); it != it_end; ++ it)
	278	{
	279	if ((*it).mRay->mTerminationObject)
	280	objects.push_back((*it).mRay->mTerminationObject);
	281	if ((*it).mRay->mOriginObject)
	282	objects.push_back((*it).mRay->mOriginObject);
	283	}
	284	}
	285
[462]	286	void VspKdLeaf::GetRays(VssRayContainer &rays)
[428]	287	{
	288	RayInfoContainer::const_iterator it, it_end = mRays.end();
	289
	290	for (it = mRays.begin(); it != mRays.end(); ++ it)
	291	rays.push_back((*it).mRay);
	292	}
	293
[462]	294	VspKdViewCell *VspKdLeaf::GetViewCell()
[453]	295	{
	296	return mViewCell;
	297	}
	298
[462]	299	void VspKdLeaf::SetViewCell(VspKdViewCell *viewCell)
	300	{
	301	mViewCell = viewCell;
	302	}
	303
	304
	305	void VspKdLeaf::UpdatePvsSize()
	306	{
[465]	307	if (!mValidPvs)
[462]	308	{
	309	Intersectable::NewMail();
	310	int pvsSize = 0;
[465]	311	for(RayInfoContainer::iterator ri = mRays.begin();
[462]	312	ri != mRays.end(); ++ ri)
	313	{
[465]	314	if ((*ri).mRay->IsActive())
[462]	315	{
	316	Intersectable *object;
	317	#if BIDIRECTIONAL_RAY
	318	object = (*ri).mRay->mOriginObject;
[465]	319
	320	if (object && !object->Mailed())
[462]	321	{
	322	++ pvsSize;
	323	object->Mail();
	324	}
	325	#endif
	326	object = (*ri).mRay->mTerminationObject;
[465]	327	if (object && !object->Mailed())
[462]	328	{
	329	++ pvsSize;
	330	object->Mail();
	331	}
	332	}
	333	}
	334
	335	mPvsSize = pvsSize;
	336	mValidPvs = true;
	337	}
	338	}
	339
[420]	340	/*********************************************************/
	341	/* class VspKdTree implementation */
	342	/*********************************************************/
	343
	344
[482]	345	VspKdTree::VspKdTree(): mOnlyDrivingAxis(false)
[465]	346	{
[421]	347	environment->GetIntValue("VspKdTree.Termination.maxDepth", mTermMaxDepth);
	348	environment->GetIntValue("VspKdTree.Termination.minPvs", mTermMinPvs);
	349	environment->GetIntValue("VspKdTree.Termination.minRays", mTermMinRays);
	350	environment->GetFloatValue("VspKdTree.Termination.maxRayContribution", mTermMaxRayContribution);
	351	environment->GetFloatValue("VspKdTree.Termination.maxCostRatio", mTermMaxCostRatio);
[422]	352	environment->GetFloatValue("VspKdTree.Termination.minSize", mTermMinSize);
[421]	353
[472]	354	environment->GetIntValue("VspKdTree.Termination.missTolerance", mTermMissTolerance);
	355
[421]	356	mTermMinSize = sqr(mTermMinSize);
	357
[422]	358	environment->GetFloatValue("VspKdTree.epsilon", mEpsilon);
	359	environment->GetFloatValue("VspKdTree.ct_div_ci", mCtDivCi);
[465]	360
[422]	361	environment->GetFloatValue("VspKdTree.maxTotalMemory", mMaxTotalMemory);
	362	environment->GetFloatValue("VspKdTree.maxStaticMemory", mMaxStaticMemory);
[465]	363
[422]	364	environment->GetIntValue("VspKdTree.accessTimeThreshold", mAccessTimeThreshold);
	365	environment->GetIntValue("VspKdTree.minCollapseDepth", mMinCollapseDepth);
[408]	366
[465]	367	//environment->GetIntValue("ViewCells.maxViewCells", mMaxViewCells);
[462]	368
[472]	369	environment->GetIntValue("VspKdTree.PostProcess.minViewCells", mMergeMinViewCells);
	370	environment->GetFloatValue("VspKdTree.PostProcess.maxCostRatio", mMergeMaxCostRatio);
[465]	371	environment->GetIntValue("VspKdTree.PostProcess.maxPvsSize",
	372	MergeCandidate::sMaxPvsSize);
[462]	373
[409]	374	// split type
[408]	375	char sname[128];
	376	environment->GetStringValue("VspKdTree.splitType", sname);
[409]	377	string name(sname);
[465]	378
[423]	379	Debug << "======= vsp kd tree options ========" << endl;
	380	Debug << "max depth: "<< mTermMaxDepth << endl;
	381	Debug << "min pvs: "<< mTermMinPvs << endl;
	382	Debug << "min rays: "<< mTermMinRays << endl;
	383	Debug << "max ray contribution: "<< mTermMaxRayContribution << endl;
[473]	384	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
	385	Debug << "min size: " << mTermMinSize << endl;
[423]	386
[409]	387	if (name.compare("regular") == 0)
[423]	388	{
	389	Debug << "using regular split" << endl;
[409]	390	splitType = ESplitRegular;
[423]	391	}
[409]	392	else
	393	{
	394	if (name.compare("heuristic") == 0)
[423]	395	{
	396	Debug << "using heuristic split" << endl;
[409]	397	splitType = ESplitHeuristic;
[423]	398	}
[465]	399	else
[409]	400	{
	401	cerr << "Invalid VspKdTree split type " << name << endl;
[408]	402	exit(1);
	403	}
[409]	404	}
[408]	405
[418]	406	mRoot = NULL;
[409]	407
[422]	408	mSplitCandidates = new vector<SortableEntry>;
[406]	409	}
	410
[408]	411
	412	VspKdTree::~VspKdTree()
[406]	413	{
[419]	414	DEL_PTR(mRoot);
[449]	415	DEL_PTR(mSplitCandidates);
[408]	416	}
[406]	417
[418]	418	void VspKdStatistics::Print(ostream &app) const
[408]	419	{
[426]	420	app << "===== VspKdTree statistics ===============\n";
[408]	421
[426]	422	app << "#N_RAYS ( Number of rays )\n"
	423	<< rays << endl;
[465]	424
[426]	425	app << "#N_INITPVS ( Initial PVS size )\n"
	426	<< initialPvsSize << endl;
[465]	427
[426]	428	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
[406]	429
[426]	430	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
[406]	431
[482]	432	app << "#N_SPLITS ( Number of splits in axes x y z\n";
[465]	433
[482]	434	for (int i = 0; i < 3; ++ i)
	435	app << splits[i] << " ";
[426]	436	app << endl;
[406]	437
[426]	438	app << "#N_RAYREFS ( Number of rayRefs )\n" <<
	439	rayRefs << "\n";
[408]	440
[426]	441	app << "#N_RAYRAYREFS ( Number of rayRefs / ray )\n" <<
	442	rayRefs / (double)rays << "\n";
[408]	443
[426]	444	app << "#N_LEAFRAYREFS ( Number of rayRefs / leaf )\n" <<
	445	rayRefs / (double)Leaves() << "\n";
[408]	446
[426]	447	app << "#N_MAXRAYREFS ( Max number of rayRefs / leaf )\n" <<
	448	maxRayRefs << "\n";
[408]	449
	450	// app << setprecision(4);
	451
[426]	452	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n" <<
	453	maxDepthNodes * 100 / (double)Leaves() << endl;
[408]	454
[426]	455	app << "#N_PMINPVSLEAVES ( Percentage of leaves with mininimal PVS )\n" <<
	456	minPvsNodes * 100 / (double)Leaves() << endl;
[408]	457
[426]	458	app << "#N_PMINRAYSLEAVES ( Percentage of leaves with minimal number of rays)\n" <<
	459	minRaysNodes * 100 / (double)Leaves() << endl;
	460
[408]	461	app << "#N_PMINSIZELEAVES ( Percentage of leaves with minSize )\n"<<
[426]	462	minSizeNodes * 100 / (double)Leaves() << endl;
[408]	463
[426]	464	app << "#N_PMAXRAYCONTRIBLEAVES ( Percentage of leaves with maximal ray contribution )\n" <<
	465	maxRayContribNodes * 100 / (double)Leaves() << endl;
[408]	466
[482]	467	app << "#N_PMAXRCOSTLEAVES ( Percentage of leaves with maximal cost ratio )\n" <<
	468	maxCostNodes * 100 / (double)Leaves() << endl;
	469
[426]	470	app << "#N_ADDED_RAYREFS ( Number of dynamically added ray references )\n"<<
	471	addedRayRefs << endl;
[408]	472
[426]	473	app << "#N_REMOVED_RAYREFS ( Number of dynamically removed ray references )\n"<<
	474	removedRayRefs << endl;
[408]	475
[426]	476	// app << setprecision(4);
[408]	477
[462]	478	app << "#N_( Maximal PVS size / leaf)\n"
[426]	479	<< maxPvsSize << endl;
[408]	480
[482]	481	app << "#N_( Average PVS size / leaf)\n"
	482	<< (double) accPvsSize / (double)Leaves() << endl;
	483
[426]	484	app << "#N_CTIME ( Construction time [s] )\n"
	485	<< Time() << " \n";
[408]	486
[426]	487	app << "===== END OF VspKdTree statistics ==========\n";
[406]	488	}
	489
[408]	490
[463]	491	void VspKdTree::CollectViewCells(ViewCellContainer &viewCells) const
	492	{
	493	stack<VspKdNode *> nodeStack;
	494	nodeStack.push(mRoot);
	495
	496	ViewCell::NewMail();
	497
[465]	498	while (!nodeStack.empty())
[463]	499	{
	500	VspKdNode *node = nodeStack.top();
	501	nodeStack.pop();
	502
[465]	503	if (node->IsLeaf())
[463]	504	{
	505	ViewCell viewCell = dynamic_cast<VspKdLeaf >(node)->mViewCell;
	506
[465]	507	if (!viewCell->Mailed())
[463]	508	{
	509	viewCell->Mail();
	510	viewCells.push_back(viewCell);
	511	}
	512	}
[465]	513	else
[463]	514	{
	515	VspKdInterior interior = dynamic_cast<VspKdInterior >(node);
	516
	517	nodeStack.push(interior->GetFront());
	518	nodeStack.push(interior->GetBack());
	519	}
	520	}
[408]	521	}
[406]	522
[462]	523	void VspKdTree::Construct(const VssRayContainer &rays,
	524	AxisAlignedBox3 *forcedBoundingBox)
[452]	525	{
[414]	526	mStat.Start();
[465]	527
[422]	528	mMaxMemory = mMaxStaticMemory;
[419]	529	DEL_PTR(mRoot);
[408]	530
[462]	531	mRoot = new VspKdLeaf(NULL, (int)rays.size());
[408]	532
[423]	533	// first construct a leaf that will get subdivided
[462]	534	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(mRoot);
[465]	535
[414]	536	mStat.nodes = 1;
[448]	537	mBox.Initialize();
[452]	538
[448]	539	//-- compute bounding box
[409]	540	if (forcedBoundingBox)
[419]	541	mBox = *forcedBoundingBox;
[448]	542	else
	543	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
	544	{
[452]	545	if ((*ri)->mOriginObject)
	546	mBox.Include((*ri)->GetOrigin());
	547	if ((*ri)->mTerminationObject)
[465]	548	mBox.Include((*ri)->GetTermination());
[448]	549	}
[408]	550
[422]	551	Debug << "bbox = " << mBox << endl;
[448]	552
	553	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
	554	{
	555	//leaf->AddRay(RayInfo(*ri));
[463]	556	VssRay ray = ri;
	557	float minT, maxT;
[482]	558
	559	// TODO: not very efficient to implictly cast between rays types!
[463]	560	if (mBox.GetRaySegment(*ray, minT, maxT))
	561	{
	562	float len = ray->Length();
[465]	563
	564	if (!len)
[463]	565	len = Limits::Small;
[465]	566
	567	leaf->AddRay(RayInfo(ray, minT / len, maxT / len));
[448]	568	}
	569	}
[465]	570
[420]	571	mStat.rays = (int)leaf->mRays.size();
[408]	572	leaf->UpdatePvsSize();
[465]	573
[414]	574	mStat.initialPvsSize = leaf->GetPvsSize();
[465]	575
[409]	576	// Subdivide();
[419]	577	mRoot = Subdivide(TraversalData(leaf, mBox, 0));
[408]	578
[465]	579	if (mSplitCandidates)
[409]	580	{
[422]	581	// force release of this vector
	582	delete mSplitCandidates;
	583	mSplitCandidates = new vector<SortableEntry>;
[409]	584	}
[465]	585
[414]	586	mStat.Stop();
[408]	587
[414]	588	mStat.Print(cout);
[423]	589	Debug << "#Total memory=" << GetMemUsage() << endl;
[408]	590	}
	591
	592
	593
[462]	594	VspKdNode *VspKdTree::Subdivide(const TraversalData &tdata)
[408]	595	{
[462]	596	VspKdNode *result = NULL;
[408]	597
[409]	598	priority_queue<TraversalData> tStack;
[423]	599	//stack<TraversalData> tStack;
[465]	600
[409]	601	tStack.push(tdata);
[406]	602
[409]	603	AxisAlignedBox3 backBox;
	604	AxisAlignedBox3 frontBox;
[465]	605
[408]	606	int lastMem = 0;
[423]	607
[465]	608	while (!tStack.empty())
[409]	609	{
[408]	610	float mem = GetMemUsage();
[465]	611
	612	if (lastMem / 10 != ((int)mem) / 10)
[409]	613	{
[425]	614	Debug << mem << " MB" << endl;
[408]	615	}
	616	lastMem = (int)mem;
[465]	617
	618	if (1 && mem > mMaxMemory)
[409]	619	{
[423]	620	Debug << "memory limit reached: " << mem << endl;
[409]	621	// count statistics on unprocessed leafs
[465]	622	while (!tStack.empty())
[409]	623	{
[408]	624	EvaluateLeafStats(tStack.top());
	625	tStack.pop();
[409]	626	}
	627	break;
	628	}
[465]	629
[409]	630	TraversalData data = tStack.top();
[465]	631	tStack.pop();
	632
[472]	633	VspKdNode node = SubdivideNode(dynamic_cast<VspKdLeaf >(data.mNode),
	634	data.mBox, backBox, frontBox);
[409]	635	if (result == NULL)
	636	result = node;
[465]	637
	638	if (!node->IsLeaf())
[409]	639	{
[462]	640	VspKdInterior interior = dynamic_cast<VspKdInterior >(node);
[408]	641
[409]	642	// push the children on the stack
[419]	643	tStack.push(TraversalData(interior->GetBack(), backBox, data.mDepth + 1));
	644	tStack.push(TraversalData(interior->GetFront(), frontBox, data.mDepth + 1));
[465]	645	}
	646	else
[409]	647	{
	648	EvaluateLeafStats(data);
	649	}
	650	}
	651
	652	return result;
[408]	653	}
[406]	654
[408]	655
	656	// returns selected plane for subdivision
[462]	657	int VspKdTree::SelectPlane(VspKdLeaf *leaf,
[424]	658	const AxisAlignedBox3 &box,
	659	float &position,
	660	int &raysBack,
	661	int &raysFront,
	662	int &pvsBack,
	663	int &pvsFront)
[408]	664	{
[428]	665	int axis = 0;
	666	float costRatio = 0;
[465]	667
	668	if (splitType == ESplitRegular)
[422]	669	{
[408]	670	costRatio = BestCostRatioRegular(leaf,
[480]	671	box,
[410]	672	axis,
	673	position,
	674	raysBack,
	675	raysFront,
	676	pvsBack,
	677	pvsFront);
[465]	678	}
[423]	679	else if (splitType == ESplitHeuristic)
[410]	680	{
[472]	681	costRatio = BestCostRatioHeuristic(leaf,
[480]	682	box,
[472]	683	axis,
	684	position,
	685	raysBack,
	686	raysFront,
	687	pvsBack,
	688	pvsFront);
[423]	689	}
[465]	690	else
[423]	691	{
	692	cerr << "VspKdTree: Unknown split heuristics\n";
	693	exit(1);
	694	}
[408]	695
[465]	696	if (costRatio > mTermMaxCostRatio)
[410]	697	{
[482]	698	//Debug << "Too big cost ratio " << costRatio << endl;
	699	++ leaf->mMaxCostMisses;
	700	if (leaf->mMaxCostMisses > mTermMissTolerance)
	701	{
	702	++ mStat.maxCostNodes;
	703	return -1;
	704	}
[408]	705	}
	706
[428]	707	if (0)
[465]	708	Debug << "pvs=" << leaf->mPvsSize
	709	<< " rays=" << (int)leaf->mRays.size()
	710	<< " rc=" << leaf->GetAvgRayContribution()
[428]	711	<< " axis=" << axis << endl;
[465]	712
[408]	713	return axis;
[406]	714	}
	715
	716
[465]	717
[462]	718	float VspKdTree::EvalCostRatio(VspKdLeaf *leaf,
[480]	719	const AxisAlignedBox3 &box,
[425]	720	const int axis,
	721	const float position,
	722	int &raysBack,
	723	int &raysFront,
	724	int &pvsBack,
	725	int &pvsFront)
[406]	726	{
[408]	727	raysBack = 0;
	728	raysFront = 0;
	729	pvsFront = 0;
	730	pvsBack = 0;
	731
	732	Intersectable::NewMail(3);
[465]	733
[408]	734	// eval pvs size
[482]	735	const int pvsSize = leaf->GetPvsSize();
[408]	736
[419]	737	// this is the main ray classification loop!
[480]	738	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[463]	739	ri != leaf->mRays.end(); ++ ri)
[480]	740	{
	741	if (!(*ri).mRay->IsActive())
	742	continue;
[465]	743
[480]	744	// determine the side of this ray with respect to the plane
	745	float t;
	746	int side = (*ri).ComputeRayIntersection(axis, position, t);
	747
[419]	748	if (side <= 0)
	749	++ raysBack;
[465]	750
[419]	751	if (side >= 0)
	752	++ raysFront;
[465]	753
[419]	754	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
[465]	755	}
[408]	756
[482]	757	//-- only one of these options should be one
	758
	759	if (0) //-- only pvs
[425]	760	{
[428]	761	const float sum = float(pvsBack + pvsFront);
[482]	762	const float oldCost = (float)pvsSize;
[425]	763
	764	return sum / oldCost;
	765	}
[482]	766
	767	//-- pvs + probability heuristics
	768	float pBack, pFront, pOverall;
	769
	770	if (0)
[428]	771	{
[482]	772	// box length substitute for probability
	773	const float minBox = box.Min(axis);
	774	const float maxBox = box.Max(axis);
[480]	775
[482]	776	pBack = position - minBox;
	777	pFront = maxBox - position;
	778	pOverall = maxBox - minBox;
	779	}
[465]	780
[482]	781	if (1) //-- area substitute for probability
	782	{
	783
	784	const bool useMidSplit = true;
	785	//const bool useMidSplit = false;
	786
	787	pOverall = box.SurfaceArea();
	788
	789	if (!useMidSplit)
	790	{
	791	Vector3 pos = box.Max(); pos[axis] = position;
	792	pBack = AxisAlignedBox3(box.Min(), pos).SurfaceArea();
[480]	793
[482]	794	pos = box.Min(); pos[axis] = position;
	795	pFront = AxisAlignedBox3(pos, box.Max()).SurfaceArea();
	796	}
	797	else
	798	{
	799	//-- simplified computation for mid split
	800	const int axis2 = (axis + 1) % 3;
	801	const int axis3 = (axis + 2) % 3;
[480]	802
[482]	803	const float faceArea =
	804	(box.Max(axis2) - box.Min(axis2)) *
	805	(box.Max(axis3) - box.Min(axis3));
[408]	806
[482]	807	pBack = pFront = pOverall * 0.5f + faceArea;
	808	}
	809	}
[465]	810
[482]	811	//-- ray density substitute for probability
[483]	812	if (0)
[482]	813	{
	814	pBack = (float)raysBack;
	815	pFront = (float)raysFront;
	816	pOverall = (float)leaf->mRays.size();
	817	}
[465]	818
[482]	819	//Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
	820	//Debug << pFront << " " << pBack << " " << pOverall << endl;
	821
	822	// float sum = raysBack(position - minBox) + raysFront(maxBox - position);
	823	const float newCost = pvsBack * pBack + pvsFront * pFront;
	824	// float oldCost = leaf->mRays.size();
	825	const float oldCost = (float)pvsSize * pOverall;
	826
	827	return (mCtDivCi + newCost) / oldCost;
[406]	828	}
	829
[462]	830	float VspKdTree::BestCostRatioRegular(VspKdLeaf *leaf,
[480]	831	const AxisAlignedBox3 &box,
[422]	832	int &axis,
	833	float &position,
	834	int &raysBack,
	835	int &raysFront,
	836	int &pvsBack,
	837	int &pvsFront)
[408]	838	{
[422]	839	int nRaysBack[3], nRaysFront[3];
	840	int nPvsBack[3], nPvsFront[3];
	841
	842	float nPosition[3];
	843	float nCostRatio[3];
[408]	844	int bestAxis = -1;
[465]	845
[480]	846	//AxisAlignedBox3 sBox = GetBBox(leaf);
	847	const int sAxis = box.Size().DrivingAxis();
[408]	848
[465]	849	for (axis = 0; axis < 3; ++ axis)
[410]	850	{
[465]	851	if (!mOnlyDrivingAxis \|\| axis == sAxis)
[410]	852	{
[465]	853
[480]	854	nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
[465]	855
[408]	856	nCostRatio[axis] = EvalCostRatio(leaf,
[480]	857	box,
[410]	858	axis,
	859	nPosition[axis],
	860	nRaysBack[axis],
	861	nRaysFront[axis],
	862	nPvsBack[axis],
	863	nPvsFront[axis]);
	864	if (bestAxis == -1)
[452]	865	{
[408]	866	bestAxis = axis;
[452]	867	}
	868	else if (nCostRatio[axis] < nCostRatio[bestAxis])
	869	{
[465]	870	/*Debug << "pvs front " << nPvsBack[axis]
[452]	871	<< " pvs back " << nPvsFront[axis]
	872	<< " overall pvs " << leaf->GetPvsSize() << endl;*/
	873	bestAxis = axis;
	874	}
[465]	875
[408]	876	}
	877	}
	878
[428]	879	//-- assign best axis
[408]	880	axis = bestAxis;
	881	position = nPosition[bestAxis];
	882
	883	raysBack = nRaysBack[bestAxis];
	884	raysFront = nRaysFront[bestAxis];
	885
	886	pvsBack = nPvsBack[bestAxis];
	887	pvsFront = nPvsFront[bestAxis];
[465]	888
[408]	889	return nCostRatio[bestAxis];
	890	}
	891
[462]	892	float VspKdTree::BestCostRatioHeuristic(VspKdLeaf *leaf,
[480]	893	const AxisAlignedBox3 &box,
[410]	894	int &axis,
	895	float &position,
	896	int &raysBack,
	897	int &raysFront,
	898	int &pvsBack,
	899	int &pvsFront)
[406]	900	{
[480]	901	//AxisAlignedBox3 box = GetBBox(leaf);
[408]	902	// AxisAlignedBox3 dirBox = GetDirBBox(node);
[465]	903
[408]	904	axis = box.Size().DrivingAxis();
[465]	905
[408]	906	SortSplitCandidates(leaf, axis);
[465]	907
[408]	908	// go through the lists, count the number of objects left and right
	909	// and evaluate the following cost funcion:
	910	// C = ct_div_ci + (qlrl + qrrr)/queries
[465]	911
[420]	912	int rl = 0, rr = (int)leaf->mRays.size();
	913	int pl = 0, pr = leaf->GetPvsSize();
[424]	914
[408]	915	float minBox = box.Min(axis);
	916	float maxBox = box.Max(axis);
	917	float sizeBox = maxBox - minBox;
[465]	918
[419]	919	float minBand = minBox + 0.1f*(maxBox - minBox);
	920	float maxBand = minBox + 0.9f*(maxBox - minBox);
[465]	921
[408]	922	float sum = rr*sizeBox;
[411]	923	float minSum = 1e20f;
[465]	924
[408]	925	Intersectable::NewMail();
[424]	926
[408]	927	// set all object as belonging to the fron pvs
[420]	928	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[424]	929	ri != leaf->mRays.end(); ++ ri)
[410]	930	{
	931	if ((*ri).mRay->IsActive())
	932	{
[408]	933	Intersectable object = (ri).mRay->mTerminationObject;
[465]	934
[408]	935	if (object)
[465]	936	if (!object->Mailed())
[410]	937	{
[408]	938	object->Mail();
	939	object->mCounter = 1;
[465]	940	}
[410]	941	else
	942	++ object->mCounter;
[408]	943	}
[410]	944	}
[465]	945
[408]	946	Intersectable::NewMail();
[465]	947
[422]	948	for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
[465]	949	ci < mSplitCandidates->end(); ++ ci)
[410]	950	{
[408]	951	VssRay *ray;
[465]	952
	953	switch ((*ci).type)
[410]	954	{
[465]	955	case SortableEntry::ERayMin:
[410]	956	{
	957	++ rl;
[482]	958	ray = (*ci).ray;
[410]	959	Intersectable *object = ray->mTerminationObject;
[465]	960	if (object && !object->Mailed())
[410]	961	{
	962	object->Mail();
	963	++ pl;
	964	}
	965	break;
	966	}
[465]	967	case SortableEntry::ERayMax:
[410]	968	{
	969	-- rr;
[465]	970
[482]	971	ray = (*ci).ray;
[410]	972	Intersectable *object = ray->mTerminationObject;
[465]	973
[410]	974	if (object)
	975	{
	976	if (-- object->mCounter == 0)
	977	-- pr;
	978	}
[465]	979
[410]	980	break;
	981	}
[408]	982	}
[465]	983
	984	if ((ci).value > minBand && (ci).value < maxBand)
[410]	985	{
[408]	986	sum = pl((ci).value - minBox) + pr(maxBox - (ci).value);
[465]	987
[410]	988	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
	989	// cout<<"cost= "<<sum<<endl;
[465]	990
	991	if (sum < minSum)
[410]	992	{
[408]	993	minSum = sum;
	994	position = (*ci).value;
[465]	995
[408]	996	raysBack = rl;
	997	raysFront = rr;
[465]	998
[408]	999	pvsBack = pl;
	1000	pvsFront = pr;
[465]	1001
[410]	1002	}
	1003	}
	1004	}
	1005
[419]	1006	float oldCost = (float)leaf->GetPvsSize();
[422]	1007	float newCost = mCtDivCi + minSum / sizeBox;
[410]	1008	float ratio = newCost / oldCost;
[465]	1009
[422]	1010	//Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
	1011	// <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
[465]	1012
[408]	1013	return ratio;
	1014	}
	1015
[462]	1016	void VspKdTree::SortSplitCandidates(VspKdLeaf *node,
[410]	1017	const int axis)
[408]	1018	{
[422]	1019	mSplitCandidates->clear();
[465]	1020
[420]	1021	int requestedSize = 2 * (int)(node->mRays.size());
[410]	1022	// creates a sorted split candidates array
[422]	1023	if (mSplitCandidates->capacity() > 500000 &&
[465]	1024	requestedSize < (int)(mSplitCandidates->capacity()/10) )
[410]	1025	{
[449]	1026	DEL_PTR(mSplitCandidates);
[422]	1027	mSplitCandidates = new vector<SortableEntry>;
[410]	1028	}
[465]	1029
[422]	1030	mSplitCandidates->reserve(requestedSize);
[408]	1031
[465]	1032	// insert all queries
[420]	1033	for(RayInfoContainer::const_iterator ri = node->mRays.begin();
[465]	1034	ri < node->mRays.end(); ++ ri)
[410]	1035	{
	1036	bool positive = (*ri).mRay->HasPosDir(axis);
[465]	1037	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
[482]	1038	(ri).ExtrapOrigin(axis), (ri).mRay));
[465]	1039
[422]	1040	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
[482]	1041	(ri).ExtrapTermination(axis), (ri).mRay));
[410]	1042	}
[465]	1043
[422]	1044	stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
[406]	1045	}
	1046
[408]	1047
[410]	1048	void VspKdTree::EvaluateLeafStats(const TraversalData &data)
[406]	1049	{
[410]	1050	// the node became a leaf -> evaluate stats for leafs
[462]	1051	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(data.mNode);
[408]	1052
[426]	1053	if (leaf->GetPvsSize() > mStat.maxPvsSize)
	1054	mStat.maxPvsSize = leaf->GetPvsSize();
	1055
	1056	if ((int)(leaf->mRays.size()) > mStat.maxRayRefs)
	1057	mStat.maxRayRefs = (int)leaf->mRays.size();
	1058
[482]	1059	mStat.rays += (int)leaf->mRays.size();
[426]	1060
[421]	1061	if (data.mDepth >= mTermMaxDepth)
[413]	1062	++ mStat.maxDepthNodes;
[465]	1063
[421]	1064	if (leaf->GetPvsSize() < mTermMinPvs)
[413]	1065	++ mStat.minPvsNodes;
[408]	1066
[482]	1067	mStat.accPvsSize += leaf->GetPvsSize();
	1068
[426]	1069	if ((int)leaf->GetRays().size() < mTermMinRays)
[413]	1070	++ mStat.minRaysNodes;
[408]	1071
[426]	1072	if (leaf->GetAvgRayContribution() > mTermMaxRayContribution)
[413]	1073	++ mStat.maxRayContribNodes;
[465]	1074
	1075	if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
[413]	1076	++ mStat.minSizeNodes;
[408]	1077	}
	1078
	1079
[465]	1080	inline bool VspKdTree::TerminationCriteriaMet(const VspKdLeaf *leaf,
[421]	1081	const AxisAlignedBox3 &box) const
	1082	{
[465]	1083	return ((leaf->GetPvsSize() < mTermMinPvs) \|\|
[426]	1084	(leaf->mRays.size() < mTermMinRays) \|\|
[482]	1085	(leaf->GetAvgRayContribution() > mTermMaxRayContribution ) \|\|
[465]	1086	(leaf->mDepth >= mTermMaxDepth) \|\|
[423]	1087	(SqrMagnitude(box.Size()) <= mTermMinSize));
[421]	1088	}
[408]	1089
[421]	1090
[462]	1091	VspKdNode VspKdTree::SubdivideNode(VspKdLeaf leaf,
[472]	1092	const AxisAlignedBox3 &box,
	1093	AxisAlignedBox3 &backBBox,
	1094	AxisAlignedBox3 &frontBBox)
[408]	1095	{
[421]	1096	if (TerminationCriteriaMet(leaf, box))
[410]	1097	{
[465]	1098	if (1 && leaf->mDepth >= mTermMaxDepth)
[422]	1099	{
[428]	1100	Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
	1101	Debug << "Bbox: " << GetBBox(leaf) << endl;
[408]	1102	}
[428]	1103	//Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
[465]	1104
[408]	1105	return leaf;
	1106	}
[465]	1107
[408]	1108	float position;
	1109	// first count ray sides
[410]	1110	int raysBack;
	1111	int raysFront;
[408]	1112	int pvsBack;
	1113	int pvsFront;
[465]	1114
[410]	1115	// select subdivision axis
[425]	1116	const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
[426]	1117	//Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" << pvsFront << endl;
[408]	1118
[465]	1119	if (axis == -1)
[482]	1120	return leaf;
	1121
[413]	1122	mStat.nodes += 2;
[482]	1123	++ mStat.splits[axis];
[408]	1124
[410]	1125	// add the new nodes to the tree
[462]	1126	VspKdInterior *node = new VspKdInterior(leaf->mParent);
[406]	1127
[419]	1128	node->mAxis = axis;
	1129	node->mPosition = position;
	1130	node->mBox = box;
[465]	1131
[410]	1132	backBBox = box;
	1133	frontBBox = box;
[406]	1134
[482]	1135	VspKdLeaf *back = new VspKdLeaf(node, raysBack, leaf->GetMaxCostMisses());
[408]	1136	back->SetPvsSize(pvsBack);
[482]	1137	VspKdLeaf *front = new VspKdLeaf(node, raysFront, leaf->GetMaxCostMisses());
[408]	1138	front->SetPvsSize(pvsFront);
[406]	1139
[410]	1140	// replace a link from node's parent
[419]	1141	if (leaf->mParent)
	1142	leaf->mParent->ReplaceChildLink(leaf, node);
[410]	1143	// and setup child links
	1144	node->SetupChildLinks(back, front);
[465]	1145
	1146	if (axis <= VspKdNode::SPLIT_Z)
[410]	1147	{
[408]	1148	backBBox.SetMax(axis, position);
[410]	1149	frontBBox.SetMin(axis, position);
[465]	1150
[420]	1151	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]	1152	ri != leaf->mRays.end(); ++ ri)
[410]	1153	{
[465]	1154	if ((*ri).mRay->IsActive())
[410]	1155	{
[408]	1156	// first unref ray from the former leaf
	1157	(*ri).mRay->Unref();
[463]	1158	float t;
[408]	1159	// determine the side of this ray with respect to the plane
[463]	1160	int side = node->ComputeRayIntersection(*ri, t);
[408]	1161
[465]	1162	if (side == 0)
[410]	1163	{
[465]	1164	if ((*ri).mRay->HasPosDir(axis))
[410]	1165	{
[420]	1166	back->AddRay(RayInfo((*ri).mRay,
	1167	(*ri).mMinT,
[463]	1168	t));
[420]	1169	front->AddRay(RayInfo((*ri).mRay,
[463]	1170	t,
[420]	1171	(*ri).mMaxT));
[465]	1172	}
	1173	else
[410]	1174	{
[420]	1175	back->AddRay(RayInfo((*ri).mRay,
[463]	1176	t,
[420]	1177	(*ri).mMaxT));
	1178	front->AddRay(RayInfo((*ri).mRay,
	1179	(*ri).mMinT,
[463]	1180	t));
[408]	1181	}
[465]	1182	}
[410]	1183	else
[465]	1184	if (side == 1)
[408]	1185	front->AddRay(*ri);
	1186	else
	1187	back->AddRay(*ri);
[410]	1188	} else
[408]	1189	(*ri).mRay->Unref();
[410]	1190	}
[465]	1191	}
	1192	else
[410]	1193	{
	1194	// rays front/back
[465]	1195
[420]	1196	for (RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]	1197	ri != leaf->mRays.end(); ++ ri)
[410]	1198	{
[465]	1199	if ((*ri).mRay->IsActive())
[410]	1200	{
[408]	1201	// first unref ray from the former leaf
	1202	(*ri).mRay->Unref();
	1203
	1204	int side;
	1205	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
	1206	side = 1;
	1207	else
	1208	side = -1;
[465]	1209
[408]	1210	if (side == 1)
	1211	front->AddRay(*ri);
	1212	else
[465]	1213	back->AddRay(*ri);
	1214	}
[410]	1215	else
[408]	1216	(*ri).mRay->Unref();
[410]	1217	}
	1218	}
[465]	1219
[410]	1220	// update stats
[420]	1221	mStat.rayRefs -= (int)leaf->mRays.size();
[414]	1222	mStat.rayRefs += raysBack + raysFront;
[408]	1223
[420]	1224	DEL_PTR(leaf);
	1225
[410]	1226	return node;
[406]	1227	}
	1228
[410]	1229	int VspKdTree::ReleaseMemory(const int time)
	1230	{
[462]	1231	stack<VspKdNode *> tstack;
[406]	1232
[410]	1233	// find a node in the tree which subtree will be collapsed
[422]	1234	int maxAccessTime = time - mAccessTimeThreshold;
[420]	1235	int released = 0;
[406]	1236
[419]	1237	tstack.push(mRoot);
[406]	1238
[423]	1239	while (!tstack.empty())
[410]	1240	{
[462]	1241	VspKdNode *node = tstack.top();
[410]	1242	tstack.pop();
[465]	1243
	1244	if (!node->IsLeaf())
[410]	1245	{
[462]	1246	VspKdInterior in = dynamic_cast<VspKdInterior >(node);
[410]	1247	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
[465]	1248	if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
[410]	1249	{
	1250	released = CollapseSubtree(node, time);
	1251	break;
	1252	}
[465]	1253	if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
[410]	1254	{
[419]	1255	tstack.push(in->GetFront());
	1256	tstack.push(in->GetBack());
[465]	1257	}
	1258	else
[410]	1259	{
[419]	1260	tstack.push(in->GetBack());
	1261	tstack.push(in->GetFront());
[410]	1262	}
	1263	}
	1264	}
[406]	1265
[465]	1266	while (tstack.empty())
[410]	1267	{
	1268	// could find node to collaps...
	1269	// cout<<"Could not find a node to release "<<endl;
	1270	break;
	1271	}
[465]	1272
[410]	1273	return released;
[408]	1274	}
[406]	1275
	1276
[462]	1277	VspKdNode VspKdTree::SubdivideLeaf(VspKdLeaf leaf,
[410]	1278	const float sizeThreshold)
[408]	1279	{
[462]	1280	VspKdNode *node = leaf;
[406]	1281
[410]	1282	AxisAlignedBox3 leafBBox = GetBBox(leaf);
[406]	1283
[408]	1284	static int pass = 0;
[410]	1285	++ pass;
[465]	1286
[410]	1287	// check if we should perform a dynamic subdivision of the leaf
[420]	1288	if (// leaf->mRays.size() > (unsigned)termMinCost &&
[465]	1289	(leaf->GetPvsSize() >= mTermMinPvs) &&
[410]	1290	(SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
	1291	{
[423]	1292	// memory check and release
[465]	1293	if (GetMemUsage() > mMaxTotalMemory)
[410]	1294	ReleaseMemory(pass);
[465]	1295
[410]	1296	AxisAlignedBox3 backBBox, frontBBox;
[406]	1297
[410]	1298	// subdivide the node
	1299	node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
	1300	}
	1301	return node;
[406]	1302	}
	1303
	1304
	1305
[423]	1306	void VspKdTree::UpdateRays(VssRayContainer &remove,
	1307	VssRayContainer &add)
[406]	1308	{
[462]	1309	VspKdLeaf::NewMail();
[406]	1310
[410]	1311	// schedule rays for removal
	1312	for(VssRayContainer::const_iterator ri = remove.begin();
[465]	1313	ri != remove.end(); ++ ri)
[410]	1314	{
[465]	1315	(*ri)->ScheduleForRemoval();
[410]	1316	}
[406]	1317
[410]	1318	int inactive = 0;
[406]	1319
[465]	1320	for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
[410]	1321	{
	1322	if ((*ri)->ScheduledForRemoval())
	1323	// RemoveRay(*ri, NULL, false);
	1324	// !!! BUG - with true it does not work correctly - aggreated delete
	1325	RemoveRay(*ri, NULL, true);
	1326	else
	1327	++ inactive;
	1328	}
[406]	1329
[410]	1330	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
[465]	1331	for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
[410]	1332	{
	1333	AddRay(*ri);
	1334	}
[406]	1335	}
	1336
	1337
[410]	1338	void VspKdTree::RemoveRay(VssRay *ray,
[462]	1339	vector<VspKdLeaf > affectedLeaves,
[410]	1340	const bool removeAllScheduledRays)
[406]	1341	{
[410]	1342	stack<RayTraversalData> tstack;
[406]	1343
[420]	1344	tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[465]	1345
[410]	1346	RayTraversalData data;
[406]	1347
[410]	1348	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
[465]	1349	while (!tstack.empty())
[410]	1350	{
	1351	data = tstack.top();
	1352	tstack.pop();
[406]	1353
[465]	1354	if (!data.mNode->IsLeaf())
[410]	1355	{
	1356	// split the set of rays in two groups intersecting the
	1357	// two subtrees
	1358	TraverseInternalNode(data, tstack);
[465]	1359	}
	1360	else
[410]	1361	{
	1362	// remove the ray from the leaf
	1363	// find the ray in the leaf and swap it with the last ray...
[462]	1364	VspKdLeaf leaf = (VspKdLeaf )data.mNode;
[465]	1365
	1366	if (!leaf->Mailed())
[410]	1367	{
	1368	leaf->Mail();
[465]	1369
[410]	1370	if (affectedLeaves)
	1371	affectedLeaves->push_back(leaf);
[465]	1372
	1373	if (removeAllScheduledRays)
[410]	1374	{
[420]	1375	int tail = (int)leaf->mRays.size() - 1;
[406]	1376
[465]	1377	for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
[410]	1378	{
[465]	1379	if (leaf->mRays[i].mRay->ScheduledForRemoval())
[410]	1380	{
	1381	// find a ray to replace it with
[465]	1382	while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
[410]	1383	{
[414]	1384	++ mStat.removedRayRefs;
[420]	1385	leaf->mRays[tail].mRay->Unref();
	1386	leaf->mRays.pop_back();
[465]	1387
[410]	1388	-- tail;
	1389	}
[465]	1390
[410]	1391	if (tail < i)
	1392	break;
[465]	1393
[414]	1394	++ mStat.removedRayRefs;
[465]	1395
[420]	1396	leaf->mRays[i].mRay->Unref();
	1397	leaf->mRays[i] = leaf->mRays[tail];
	1398	leaf->mRays.pop_back();
[410]	1399	-- tail;
	1400	}
	1401	}
	1402	}
	1403	}
[465]	1404
[410]	1405	if (!removeAllScheduledRays)
[465]	1406	for (int i=0; i < (int)leaf->mRays.size(); i++)
[410]	1407	{
[465]	1408	if (leaf->mRays[i].mRay == ray)
[410]	1409	{
[414]	1410	++ mStat.removedRayRefs;
[410]	1411	ray->Unref();
[420]	1412	leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
	1413	leaf->mRays.pop_back();
[410]	1414	// check this ray again
	1415	break;
	1416	}
	1417	}
	1418	}
[408]	1419	}
	1420
[465]	1421	if (ray->RefCount() != 0)
[410]	1422	{
	1423	cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
	1424	exit(1);
	1425	}
	1426
[406]	1427	}
	1428
[410]	1429	void VspKdTree::AddRay(VssRay *ray)
[406]	1430	{
[410]	1431	stack<RayTraversalData> tstack;
[465]	1432
[420]	1433	tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[465]	1434
[410]	1435	RayTraversalData data;
[465]	1436
	1437	while (!tstack.empty())
[410]	1438	{
	1439	data = tstack.top();
	1440	tstack.pop();
[406]	1441
[465]	1442	if (!data.mNode->IsLeaf())
[410]	1443	{
	1444	TraverseInternalNode(data, tstack);
[465]	1445	}
	1446	else
[410]	1447	{
	1448	// remove the ray from the leaf
[420]	1449	// find the ray in the leaf and swap it with the last ray
[462]	1450	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(data.mNode);
[423]	1451
[419]	1452	leaf->AddRay(data.mRayData);
[413]	1453	++ mStat.addedRayRefs;
[410]	1454	}
	1455	}
[406]	1456	}
	1457
[465]	1458	void VspKdTree::TraverseInternalNode(RayTraversalData &data,
[410]	1459	stack<RayTraversalData> &tstack)
[406]	1460	{
[462]	1461	VspKdInterior in = dynamic_cast<VspKdInterior >(data.mNode);
[465]	1462
	1463	if (in->mAxis <= VspKdNode::SPLIT_Z)
[410]	1464	{
	1465	// determine the side of this ray with respect to the plane
[465]	1466	float t;
[463]	1467	int side = in->ComputeRayIntersection(data.mRayData, t);
[465]	1468
	1469	if (side == 0)
[410]	1470	{
[465]	1471	if (data.mRayData.mRay->HasPosDir(in->mAxis))
[433]	1472	{
[465]	1473	tstack.push(RayTraversalData(in->GetBack(),
[463]	1474	RayInfo(data.mRayData.mRay, data.mRayData.mMinT, t)));
[465]	1475
[433]	1476	tstack.push(RayTraversalData(in->GetFront(),
[463]	1477	RayInfo(data.mRayData.mRay, t, data.mRayData.mMaxT)));
[465]	1478
	1479	}
	1480	else
[433]	1481	{
	1482	tstack.push(RayTraversalData(in->GetBack(),
	1483	RayInfo(data.mRayData.mRay,
[463]	1484	t,
[433]	1485	data.mRayData.mMaxT)));
	1486	tstack.push(RayTraversalData(in->GetFront(),
	1487	RayInfo(data.mRayData.mRay,
	1488	data.mRayData.mMinT,
[463]	1489	t)));
[433]	1490	}
[410]	1491	}
[433]	1492	else
	1493	if (side == 1)
	1494	tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
	1495	else
	1496	tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[465]	1497	}
	1498	else
[410]	1499	{
	1500	// directional split
[419]	1501	if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
	1502	tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
[408]	1503	else
[419]	1504	tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[410]	1505	}
[406]	1506	}
	1507
[408]	1508
[462]	1509	int VspKdTree::CollapseSubtree(VspKdNode *sroot, const int time)
[406]	1510	{
[410]	1511	// first count all rays in the subtree
	1512	// use mail 1 for this purpose
[462]	1513	stack<VspKdNode *> tstack;
[465]	1514
[410]	1515	int rayCount = 0;
	1516	int totalRayCount = 0;
	1517	int collapsedNodes = 0;
[408]	1518
	1519	#if DEBUG_COLLAPSE
[420]	1520	cout << "Collapsing subtree" << endl;
[428]	1521	cout << "accessTime=" << sroot->GetAccessTime() << endl;
[420]	1522	cout << "depth=" << (int)sroot->depth << endl;
[408]	1523	#endif
	1524
[410]	1525	// tstat.collapsedSubtrees++;
	1526	// tstat.collapseDepths += (int)sroot->depth;
	1527	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
	1528	tstack.push(sroot);
	1529	VssRay::NewMail();
[465]	1530
	1531	while (!tstack.empty())
[410]	1532	{
	1533	++ collapsedNodes;
[465]	1534
[462]	1535	VspKdNode *node = tstack.top();
[410]	1536	tstack.pop();
[465]	1537	if (node->IsLeaf())
[410]	1538	{
[462]	1539	VspKdLeaf leaf = (VspKdLeaf ) node;
[465]	1540
[420]	1541	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]	1542	ri != leaf->mRays.end(); ++ ri)
[410]	1543	{
	1544	++ totalRayCount;
[465]	1545
	1546	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed())
[410]	1547	{
[408]	1548	(*ri).mRay->Mail();
[410]	1549	++ rayCount;
[408]	1550	}
[410]	1551	}
[465]	1552	}
	1553	else
[410]	1554	{
[472]	1555	tstack.push(dynamic_cast<VspKdInterior *>(node)->GetFront());
	1556	tstack.push(dynamic_cast<VspKdInterior *>(node)->GetBack());
[410]	1557	}
	1558	}
[465]	1559
[410]	1560	VssRay::NewMail();
[406]	1561
[410]	1562	// create a new node that will hold the rays
[462]	1563	VspKdLeaf *newLeaf = new VspKdLeaf(sroot->mParent, rayCount);
[465]	1564
[419]	1565	if (newLeaf->mParent)
[472]	1566	{
[419]	1567	newLeaf->mParent->ReplaceChildLink(sroot, newLeaf);
[472]	1568	}
[410]	1569	tstack.push(sroot);
[406]	1570
[465]	1571	while (!tstack.empty())
[410]	1572	{
[462]	1573	VspKdNode *node = tstack.top();
[410]	1574	tstack.pop();
[408]	1575
[465]	1576	if (node->IsLeaf())
[410]	1577	{
[462]	1578	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(node);
[408]	1579
[420]	1580	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[465]	1581	ri != leaf->mRays.end(); ++ ri)
[410]	1582	{
[465]	1583	// unref this ray from the old node
	1584	if ((*ri).mRay->IsActive())
[410]	1585	{
[408]	1586	(*ri).mRay->Unref();
[465]	1587	if (!(*ri).mRay->Mailed())
[410]	1588	{
[408]	1589	(*ri).mRay->Mail();
	1590	newLeaf->AddRay(*ri);
	1591	}
[465]	1592	}
[410]	1593	else
[408]	1594	(*ri).mRay->Unref();
[410]	1595	}
[465]	1596	}
	1597	else
[410]	1598	{
[465]	1599	VspKdInterior *interior =
[462]	1600	dynamic_cast<VspKdInterior *>(node);
[419]	1601	tstack.push(interior->GetBack());
	1602	tstack.push(interior->GetFront());
[410]	1603	}
	1604	}
[465]	1605
[410]	1606	// delete the node and all its children
[420]	1607	DEL_PTR(sroot);
[465]	1608
[462]	1609	// for(VspKdNode::SRayContainer::iterator ri = newleaf->mRays.begin();
[420]	1610	// ri != newleaf->mRays.end(); ++ ri)
[410]	1611	// (*ri).ray->UnMail(2);
[408]	1612
	1613	#if DEBUG_COLLAPSE
[423]	1614	cout << "Total memory before=" << GetMemUsage() << endl;
[408]	1615	#endif
	1616
[414]	1617	mStat.nodes -= collapsedNodes - 1;
	1618	mStat.rayRefs -= totalRayCount - rayCount;
[465]	1619
[408]	1620	#if DEBUG_COLLAPSE
[410]	1621	cout << "collapsed nodes" << collapsedNodes << endl;
	1622	cout << "collapsed rays" << totalRayCount - rayCount << endl;
	1623	cout << "Total memory after=" << GetMemUsage() << endl;
	1624	cout << "================================" << endl;
[408]	1625	#endif
	1626
	1627	// tstat.collapsedNodes += collapsedNodes;
	1628	// tstat.collapsedRays += totalRayCount - rayCount;
[410]	1629	return totalRayCount - rayCount;
[406]	1630	}
	1631
[408]	1632
[462]	1633	int VspKdTree::GetPvsSize(VspKdNode *node, const AxisAlignedBox3 &box) const
[406]	1634	{
[462]	1635	stack<VspKdNode *> tstack;
[419]	1636	tstack.push(mRoot);
[408]	1637
	1638	Intersectable::NewMail();
	1639	int pvsSize = 0;
[465]	1640
	1641	while (!tstack.empty())
[410]	1642	{
[462]	1643	VspKdNode *node = tstack.top();
[410]	1644	tstack.pop();
[465]	1645
	1646	if (node->IsLeaf())
[410]	1647	{
[462]	1648	VspKdLeaf leaf = (VspKdLeaf )node;
[420]	1649	for (RayInfoContainer::iterator ri = leaf->mRays.begin();
	1650	ri != leaf->mRays.end(); ++ ri)
[410]	1651	{
[465]	1652	if ((*ri).mRay->IsActive())
[410]	1653	{
[408]	1654	Intersectable *object;
	1655	#if BIDIRECTIONAL_RAY
	1656	object = (*ri).mRay->mOriginObject;
[465]	1657
	1658	if (object && !object->Mailed())
[410]	1659	{
	1660	++ pvsSize;
[408]	1661	object->Mail();
	1662	}
	1663	#endif
	1664	object = (*ri).mRay->mTerminationObject;
[465]	1665	if (object && !object->Mailed())
[410]	1666	{
	1667	++ pvsSize;
[408]	1668	object->Mail();
	1669	}
	1670	}
[465]	1671	}
[410]	1672	}
[465]	1673	else
[410]	1674	{
[462]	1675	VspKdInterior in = dynamic_cast<VspKdInterior >(node);
[465]	1676
	1677	if (in->mAxis < 3)
[410]	1678	{
[419]	1679	if (box.Max(in->mAxis) >= in->mPosition)
	1680	tstack.push(in->GetFront());
[465]	1681
[419]	1682	if (box.Min(in->mAxis) <= in->mPosition)
	1683	tstack.push(in->GetBack());
[465]	1684	}
	1685	else
[410]	1686	{
[408]	1687	// both nodes for directional splits
[419]	1688	tstack.push(in->GetFront());
	1689	tstack.push(in->GetBack());
[408]	1690	}
	1691	}
	1692	}
[410]	1693
[408]	1694	return pvsSize;
[406]	1695	}
	1696
[410]	1697	void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
	1698	float &maxRayContribution,
	1699	float &avgRayContribution)
[406]	1700	{
[462]	1701	stack<VspKdNode *> tstack;
[419]	1702	tstack.push(mRoot);
[465]	1703
[408]	1704	minRayContribution = 1.0f;
	1705	maxRayContribution = 0.0f;
[410]	1706
[408]	1707	float sumRayContribution = 0.0f;
	1708	int leaves = 0;
[406]	1709
[465]	1710	while (!tstack.empty())
[410]	1711	{
[462]	1712	VspKdNode *node = tstack.top();
[410]	1713	tstack.pop();
[465]	1714	if (node->IsLeaf())
[410]	1715	{
[408]	1716	leaves++;
[462]	1717	VspKdLeaf leaf = (VspKdLeaf )node;
[408]	1718	float c = leaf->GetAvgRayContribution();
[410]	1719
[408]	1720	if (c > maxRayContribution)
	1721	maxRayContribution = c;
	1722	if (c < minRayContribution)
	1723	minRayContribution = c;
	1724	sumRayContribution += c;
[465]	1725
	1726	}
[482]	1727	else
	1728	{
[462]	1729	VspKdInterior in = (VspKdInterior )node;
[482]	1730
[419]	1731	tstack.push(in->GetFront());
	1732	tstack.push(in->GetBack());
[408]	1733	}
	1734	}
[465]	1735
[423]	1736	Debug << "sum=" << sumRayContribution << endl;
	1737	Debug << "leaves=" << leaves << endl;
	1738	avgRayContribution = sumRayContribution / (float)leaves;
[406]	1739	}
	1740
	1741
[410]	1742	int VspKdTree::GenerateRays(const float ratioPerLeaf,
	1743	SimpleRayContainer &rays)
[406]	1744	{
[462]	1745	stack<VspKdNode *> tstack;
[419]	1746	tstack.push(mRoot);
[465]	1747
	1748	while (!tstack.empty())
[410]	1749	{
[462]	1750	VspKdNode *node = tstack.top();
[410]	1751	tstack.pop();
[465]	1752
	1753	if (node->IsLeaf())
[410]	1754	{
[462]	1755	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(node);
[423]	1756
[408]	1757	float c = leaf->GetAvgRayContribution();
[411]	1758	int num = (int)(c*ratioPerLeaf + 0.5);
[408]	1759	// cout<<num<<" ";
	1760
[465]	1761	for (int i=0; i < num; i++)
[410]	1762	{
[408]	1763	Vector3 origin = GetBBox(leaf).GetRandomPoint();
[419]	1764	/*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
[408]	1765	Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
	1766	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
[419]	1767	rays.push_back(SimpleRay(origin, direction));*/
[408]	1768	}
[465]	1769
	1770	}
	1771	else
[410]	1772	{
[465]	1773	VspKdInterior *in =
[462]	1774	dynamic_cast<VspKdInterior *>(node);
[482]	1775
[419]	1776	tstack.push(in->GetFront());
	1777	tstack.push(in->GetBack());
[408]	1778	}
[406]	1779	}
	1780
[411]	1781	return (int)rays.size();
[406]	1782	}
	1783
	1784
[410]	1785	float VspKdTree::GetAvgPvsSize()
[406]	1786	{
[462]	1787	stack<VspKdNode *> tstack;
[419]	1788	tstack.push(mRoot);
[406]	1789
[408]	1790	int sumPvs = 0;
	1791	int leaves = 0;
[465]	1792
	1793	while (!tstack.empty())
[410]	1794	{
[462]	1795	VspKdNode *node = tstack.top();
[410]	1796	tstack.pop();
[465]	1797
	1798	if (node->IsLeaf())
[410]	1799	{
[462]	1800	VspKdLeaf leaf = (VspKdLeaf )node;
[408]	1801	// update pvs size
	1802	leaf->UpdatePvsSize();
	1803	sumPvs += leaf->GetPvsSize();
	1804	leaves++;
[465]	1805	}
	1806	else
[410]	1807	{
[462]	1808	VspKdInterior in = (VspKdInterior )node;
[482]	1809
[419]	1810	tstack.push(in->GetFront());
	1811	tstack.push(in->GetBack());
[406]	1812	}
	1813	}
[408]	1814
[410]	1815	return sumPvs / (float)leaves;
[418]	1816	}
	1817
[462]	1818	VspKdNode *VspKdTree::GetRoot() const
[418]	1819	{
	1820	return mRoot;
	1821	}
	1822
[462]	1823	AxisAlignedBox3 VspKdTree::GetBBox(VspKdNode *node) const
[418]	1824	{
[419]	1825	if (node->mParent == NULL)
	1826	return mBox;
[418]	1827
	1828	if (!node->IsLeaf())
[462]	1829	return (dynamic_cast<VspKdInterior *>(node))->mBox;
[418]	1830
[419]	1831	if (node->mParent->mAxis >= 3)
	1832	return node->mParent->mBox;
[465]	1833
[419]	1834	AxisAlignedBox3 box(node->mParent->mBox);
[448]	1835	if (node->mParent->GetFront() == node)
[419]	1836	box.SetMin(node->mParent->mAxis, node->mParent->mPosition);
[418]	1837	else
[419]	1838	box.SetMax(node->mParent->mAxis, node->mParent->mPosition);
	1839
	1840	return box;
[418]	1841	}
	1842
[465]	1843	int VspKdTree::GetRootPvsSize() const
[418]	1844	{
	1845	return GetPvsSize(mRoot, mBox);
	1846	}
[419]	1847
[465]	1848	const VspKdStatistics &VspKdTree::GetStatistics() const
[419]	1849	{
	1850	return mStat;
[420]	1851	}
	1852
	1853	void VspKdTree::AddRays(VssRayContainer &add)
	1854	{
	1855	VssRayContainer remove;
	1856	UpdateRays(remove, add);
	1857	}
[465]	1858
[420]	1859	// return memory usage in MB
[465]	1860	float VspKdTree::GetMemUsage() const
[420]	1861	{
[465]	1862	return
	1863	(sizeof(VspKdTree) +
[462]	1864	mStat.Leaves() * sizeof(VspKdLeaf) +
	1865	mStat.Interior() * sizeof(VspKdInterior) +
[420]	1866	mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
	1867	}
[465]	1868
	1869	float VspKdTree::GetRayMemUsage() const
[420]	1870	{
	1871	return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
[425]	1872	}
	1873
[482]	1874
[465]	1875	int VspKdTree::ComputePvsSize(VspKdNode *node,
[426]	1876	const RayInfoContainer &globalRays) const
[425]	1877	{
	1878	int pvsSize = 0;
	1879
	1880	const AxisAlignedBox3 box = GetBBox(node);
	1881
	1882	RayInfoContainer::const_iterator it, it_end = globalRays.end();
	1883
[482]	1884	// TODO
	1885
[426]	1886	// warning: implicit conversion from VssRay to Ray
[425]	1887	for (it = globalRays.begin(); it != globalRays.end(); ++ it)
[426]	1888	pvsSize += box.GetMinMaxT((it).mRay, NULL, NULL);
[465]	1889
[425]	1890	return pvsSize;
[428]	1891	}
	1892
[482]	1893
[462]	1894	void VspKdTree::CollectLeaves(vector<VspKdLeaf *> &leaves) const
[428]	1895	{
[462]	1896	stack<VspKdNode *> nodeStack;
[428]	1897	nodeStack.push(mRoot);
	1898
[465]	1899	while (!nodeStack.empty())
[428]	1900	{
[462]	1901	VspKdNode *node = nodeStack.top();
[465]	1902
[428]	1903	nodeStack.pop();
[465]	1904
	1905	if (node->IsLeaf())
[428]	1906	{
[462]	1907	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(node);
[428]	1908	leaves.push_back(leaf);
[465]	1909	}
	1910	else
[428]	1911	{
[462]	1912	VspKdInterior interior = dynamic_cast<VspKdInterior >(node);
[428]	1913
[448]	1914	nodeStack.push(interior->GetBack());
	1915	nodeStack.push(interior->GetFront());
[428]	1916	}
	1917	}
[452]	1918	}
	1919
[462]	1920	int VspKdTree::FindNeighbors(VspKdLeaf *n,
	1921	vector<VspKdLeaf *> &neighbors,
[452]	1922	bool onlyUnmailed)
	1923	{
[462]	1924	stack<VspKdNode *> nodeStack;
[465]	1925
[452]	1926	nodeStack.push(mRoot);
	1927
	1928	AxisAlignedBox3 box = GetBBox(n);
	1929
[465]	1930	while (!nodeStack.empty())
[452]	1931	{
[462]	1932	VspKdNode *node = nodeStack.top();
[452]	1933	nodeStack.pop();
	1934
[465]	1935	if (node->IsLeaf())
[452]	1936	{
[462]	1937	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(node);
[452]	1938
[465]	1939	if (leaf != n && (!onlyUnmailed \|\| !leaf->Mailed()))
[452]	1940	neighbors.push_back(leaf);
[465]	1941	}
	1942	else
[452]	1943	{
[462]	1944	VspKdInterior interior = dynamic_cast<VspKdInterior >(node);
[452]	1945
	1946	if (interior->mPosition > box.Max(interior->mAxis))
	1947	nodeStack.push(interior->mBack);
	1948	else
	1949	{
	1950	if (interior->mPosition < box.Min(interior->mAxis))
	1951	nodeStack.push(interior->mFront);
[465]	1952	else
[452]	1953	{
	1954	// random decision
	1955	nodeStack.push(interior->mBack);
	1956	nodeStack.push(interior->mFront);
	1957	}
	1958	}
	1959	}
	1960	}
	1961
	1962	return (int)neighbors.size();
	1963	}
	1964
[462]	1965
	1966	void VspKdTree::SetViewCellsManager(ViewCellsManager *vcm)
[452]	1967	{
[462]	1968	mViewCellsManager = vcm;
	1969	}
[452]	1970
[462]	1971
[468]	1972	int VspKdTree::CastLineSegment(const Vector3 &origin,
	1973	const Vector3 &termination,
	1974	vector<ViewCell *> &viewcells)
[462]	1975	{
	1976	int hits = 0;
[465]	1977
[468]	1978	float mint = 0.0f, maxt = 1.0f;
	1979	const Vector3 dir = termination - origin;
[465]	1980
[483]	1981	stack<LineTraversalData> tStack;
[462]	1982
	1983	Intersectable::NewMail();
	1984
[468]	1985	Vector3 entp = origin;
	1986	Vector3 extp = termination;
[465]	1987
[462]	1988	VspKdNode *node = mRoot;
	1989	VspKdNode *farChild;
	1990
	1991	float position;
	1992	int axis;
[465]	1993
	1994	while (1)
[462]	1995	{
[465]	1996	if (!node->IsLeaf())
[462]	1997	{
	1998	VspKdInterior in = dynamic_cast<VspKdInterior >(node);
	1999	position = in->mPosition;
	2000	axis = in->mAxis;
	2001
[465]	2002	if (entp[axis] <= position)
[462]	2003	{
[465]	2004	if (extp[axis] <= position)
[462]	2005	{
	2006	node = in->mBack;
	2007	// cases N1,N2,N3,P5,Z2,Z3
	2008	continue;
[465]	2009	} else
[462]	2010	{
	2011	// case N4
	2012	node = in->mBack;
	2013	farChild = in->mFront;
	2014	}
	2015	}
[465]	2016	else
[462]	2017	{
[465]	2018	if (position <= extp[axis])
[462]	2019	{
	2020	node = in->mFront;
	2021	// cases P1,P2,P3,N5,Z1
	2022	continue;
[465]	2023	}
	2024	else
[462]	2025	{
	2026	node = in->mFront;
	2027	farChild = in->mBack;
	2028	// case P4
	2029	}
	2030	}
	2031
	2032	// $$ modification 3.5.2004 - hints from Kamil Ghais
	2033	// case N4 or P4
[468]	2034	float tdist = (position - origin[axis]) / dir[axis];
[483]	2035	tStack.push(LineTraversalData(farChild, extp, maxt)); //TODO
[468]	2036	extp = origin + dir * tdist;
[462]	2037	maxt = tdist;
[465]	2038	}
	2039	else
[462]	2040	{
	2041	// compute intersection with all objects in this leaf
[468]	2042	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(node);
	2043	ViewCell *vc = leaf->GetViewCell();
[465]	2044
[468]	2045	if (!vc->Mailed())
[462]	2046	{
[468]	2047	vc->Mail();
	2048	viewcells.push_back(vc);
[469]	2049	++ hits;
[462]	2050	}
[465]	2051
[483]	2052	if (1) //matt TODO: REMOVE LATER
	2053	leaf->mRays.push_back(RayInfo(new VssRay(origin, termination, NULL, NULL, 0)));
	2054
[462]	2055	// get the next node from the stack
	2056	if (tStack.empty())
	2057	break;
[465]	2058
[462]	2059	entp = extp;
	2060	mint = maxt;
[468]	2061
[483]	2062	LineTraversalData &s = tStack.top();
[462]	2063	node = s.mNode;
	2064	extp = s.mExitPoint;
	2065	maxt = s.mMaxT;
	2066	tStack.pop();
	2067	}
	2068	}
[465]	2069
[462]	2070	return hits;
	2071	}
	2072
	2073
	2074	void VspKdTree::GenerateViewCell(VspKdLeaf *leaf)
	2075	{
	2076	RayInfoContainer::const_iterator it, it_end = leaf->GetRays().end();
	2077
	2078	VspKdViewCell vc = dynamic_cast<VspKdViewCell >(mViewCellsManager->GenerateViewCell());
	2079	leaf->SetViewCell(vc);
	2080
[479]	2081	vc->SetVolume(GetBBox(leaf).GetVolume());
	2082	vc->SetArea(GetBBox(leaf).SurfaceArea());
	2083
[469]	2084	vc->mLeaves.push_back(leaf);
[462]	2085
	2086	for (it = leaf->GetRays().begin(); it != it_end; ++ it)
	2087	{
	2088	VssRay ray = (it).mRay;
	2089
[463]	2090	if (ray->mTerminationObject)
[466]	2091	vc->GetPvs().AddSample(ray->mTerminationObject);
[465]	2092
[463]	2093	if (ray->mOriginObject)
[466]	2094	vc->GetPvs().AddSample(ray->mOriginObject);
[462]	2095	}
	2096	}
	2097
	2098
[465]	2099	bool VspKdTree::MergeViewCells(VspKdLeaf l1, VspKdLeaf l2)
	2100	{
[472]	2101	//-- change pointer to view cells of all leaves associated
	2102	//-- with the previous view cells
[465]	2103	VspKdViewCell *fVc = l1->GetViewCell();
	2104	VspKdViewCell *bVc = l2->GetViewCell();
	2105
	2106	VspKdViewCell vc = dynamic_cast<VspKdViewCell >(
	2107	mViewCellsManager->MergeViewCells(fVc, bVc));
	2108
	2109	// if merge was unsuccessful
	2110	if (!vc) return false;
	2111
	2112	// set new size of view cell
[479]	2113	vc->SetVolume(fVc->GetVolume() + bVc->GetVolume());
[482]	2114	// new area
	2115	vc->SetArea(fVc->GetArea() + bVc->GetArea());
[465]	2116
[469]	2117	vector<VspKdLeaf *> fLeaves = fVc->mLeaves;
	2118	vector<VspKdLeaf *> bLeaves = bVc->mLeaves;
[465]	2119
	2120	vector<VspKdLeaf *>::const_iterator it;
	2121
	2122	//-- change view cells of all the other leaves the view cell belongs to
	2123	for (it = fLeaves.begin(); it != fLeaves.end(); ++ it)
	2124	{
	2125	(*it)->SetViewCell(vc);
[469]	2126	vc->mLeaves.push_back(*it);
[465]	2127	}
	2128
	2129	for (it = bLeaves.begin(); it != bLeaves.end(); ++ it)
	2130	{
	2131	(*it)->SetViewCell(vc);
[469]	2132	vc->mLeaves.push_back(*it);
[465]	2133	}
	2134
	2135	vc->Mail();
	2136
	2137	// clean up old view cells
	2138	DEL_PTR(fVc);
	2139	DEL_PTR(bVc);
	2140
[463]	2141	return true;
[462]	2142	}
	2143
[483]	2144	void VspKdTree::CollectMergeCandidates()
[462]	2145	{
[483]	2146	MergeCandidate::sOverallCost = 0;
[462]	2147	vector<VspKdLeaf *> leaves;
[465]	2148
[462]	2149	// collect the leaves, e.g., the "voxels" that will build the view cells
[452]	2150	CollectLeaves(leaves);
	2151
[462]	2152	VspKdLeaf::NewMail();
[452]	2153
[462]	2154	vector<VspKdLeaf *>::const_iterator it, it_end = leaves.end();
[452]	2155
[462]	2156	// generate view cells
[452]	2157	for (it = leaves.begin(); it != it_end; ++ it)
[462]	2158	GenerateViewCell(*it);
	2159
	2160	// find merge candidates and push them into queue
	2161	for (it = leaves.begin(); it != it_end; ++ it)
[452]	2162	{
[462]	2163	VspKdLeaf leaf = it;
[465]	2164
[462]	2165	// no leaf is part of two merge candidates
[452]	2166	if (!leaf->Mailed())
	2167	{
[462]	2168	leaf->Mail();
	2169
	2170	vector<VspKdLeaf *> neighbors;
[452]	2171	FindNeighbors(leaf, neighbors, true);
	2172
[465]	2173	vector<VspKdLeaf *>::const_iterator nit,
[462]	2174	nit_end = neighbors.end();
[452]	2175
[462]	2176	for (nit = neighbors.begin(); nit != nit_end; ++ nit)
[472]	2177	{
[483]	2178	// TODO: test if at least one ray goes from one leaf to the other
[472]	2179	MergeCandidate mc = MergeCandidate(leaf, *nit);
[483]	2180	mMergeQueue.push(mc);
[472]	2181
	2182	MergeCandidate::sOverallCost += mc.GetLeaf1Cost();
	2183	MergeCandidate::sOverallCost += mc.GetLeaf2Cost();
	2184	}
[462]	2185	}
	2186	}
[483]	2187	}
[462]	2188
[483]	2189	/*
	2190	void VspKdTree::CollectMergeCandidates(const vector<VspKdRay *> &rays)
	2191	{
	2192	MergeCandidate::sOverallCost = 0;
	2193
	2194	vector<VspKdIntersection>::const_iterator iit;
	2195
	2196	for (int i = 0; i < (int)rays.size(); ++ i)
	2197	{
	2198	//VspKdLeaf::NewMail();
	2199
	2200	VspKdRay *ray = rays[i];
	2201
	2202	// traverse leaves stored in the rays and compare and merge consecutive
	2203	// leaves (i.e., the neighbors in the tree)
	2204	if (ray->intersections.size() < 2)
	2205	continue;
	2206
	2207	iit = ray->intersections.begin();
	2208
	2209	BspLeaf previousLeaf = (iit).mLeaf;
	2210
	2211	++ iit;
	2212
	2213	for (; iit != ray->intersections.end(); ++ iit)
	2214	{
	2215	BspLeaf leaf = (iit).mLeaf;
	2216
	2217	// TODO: how to sort out doubles?
	2218	MergeCandidate mc = MergeCandidate(leaf, previousLeaf);
	2219	mMergeQueue.push(mc);
	2220
	2221	MergeCandidate::sOverallCost += mc.GetLeaf1Cost();
	2222	MergeCandidate::sOverallCost += mc.GetLeaf2Cost();
	2223
	2224	previousLeaf = leaf;
	2225	}
	2226	}
	2227	}
	2228	*/
	2229
	2230	int VspKdTree::MergeLeaves()
	2231	{
[462]	2232	int merged = 0;
	2233
[483]	2234	int vcSize = mStat.nodes / 2 + 1;
[462]	2235	// use priority queue to merge leaves
[483]	2236	while (!mMergeQueue.empty() && (vcSize > mMergeMinViewCells) &&
	2237	(mMergeQueue.top().GetMergeCost() <
[482]	2238	mMergeMaxCostRatio * MergeCandidate::sOverallCost))
[462]	2239	{
[482]	2240	//Debug << "mergecost: " << mergeQueue.top().GetMergeCost() / MergeCandidate::sOverallCost << " " << mMergeMaxCostRatio << endl;
[483]	2241	MergeCandidate mc = mMergeQueue.top();
	2242	mMergeQueue.pop();
[465]	2243
	2244	// both view cells equal!
[462]	2245	if (mc.GetLeaf1()->GetViewCell() == mc.GetLeaf2()->GetViewCell())
	2246	continue;
	2247
	2248	if (mc.Valid())
	2249	{
[465]	2250	ViewCell::NewMail();
	2251	MergeViewCells(mc.GetLeaf1(), mc.GetLeaf2());
[462]	2252
	2253	++ merged;
	2254	-- vcSize;
[472]	2255	// increase absolute merge cost
	2256	MergeCandidate::sOverallCost += mc.GetMergeCost();
[452]	2257	}
[462]	2258	// merge candidate not valid, because one of the leaves was already
	2259	// merged with another one
[465]	2260	else
[462]	2261	{
	2262	// validate and reinsert into queue
	2263	mc.SetValid();
[483]	2264	mMergeQueue.push(mc);
[465]	2265	//Debug << "validate " << mc.GetMergeCost() << endl;
[462]	2266	}
[452]	2267	}
[462]	2268
[479]	2269	// collapse siblings belonging to the same view cell
[482]	2270	CollapseTree(mRoot);
[465]	2271	// revalidate leaves
[482]	2272	RepairVcLeafLists();
[465]	2273
[483]	2274	Debug << "merged " << merged << " of " << mStat.nodes / 2 + 1 << " leaves" << endl;
[465]	2275
[469]	2276	//TODO: should return sample contributions
[462]	2277	return merged;
[452]	2278	}
	2279
	2280
[479]	2281	void VspKdTree::RepairVcLeafLists()
[465]	2282	{
	2283	// list not valid anymore => clear
	2284	stack<VspKdNode *> nodeStack;
	2285	nodeStack.push(mRoot);
	2286
	2287	ViewCell::NewMail();
	2288
	2289	while (!nodeStack.empty())
	2290	{
	2291	VspKdNode *node = nodeStack.top();
	2292	nodeStack.pop();
	2293
	2294	if (node->IsLeaf())
	2295	{
	2296	VspKdLeaf leaf = dynamic_cast<VspKdLeaf >(node);
	2297
	2298	VspKdViewCell *viewCell = leaf->GetViewCell();
	2299
	2300	if (!viewCell->Mailed())
	2301	{
[469]	2302	viewCell->mLeaves.clear();
[465]	2303	viewCell->Mail();
	2304	}
	2305
[469]	2306	viewCell->mLeaves.push_back(leaf);
[465]	2307	}
	2308	else
	2309	{
	2310	VspKdInterior interior = dynamic_cast<VspKdInterior >(node);
	2311
	2312	nodeStack.push(interior->GetFront());
	2313	nodeStack.push(interior->GetBack());
	2314	}
	2315	}
	2316	}
	2317
[479]	2318
[465]	2319	VspKdNode VspKdTree::CollapseTree(VspKdNode node)
	2320	{
	2321	if (node->IsLeaf())
	2322	return node;
	2323
	2324	VspKdInterior interior = dynamic_cast<VspKdInterior >(node);
	2325
	2326	VspKdNode *front = CollapseTree(interior->GetFront());
	2327	VspKdNode *back = CollapseTree(interior->GetBack());
	2328
	2329	if (front->IsLeaf() && back->IsLeaf())
	2330	{
	2331	VspKdLeaf frontLeaf = dynamic_cast<VspKdLeaf >(front);
	2332	VspKdLeaf backLeaf = dynamic_cast<VspKdLeaf >(back);
	2333
	2334	//-- collapse tree
	2335	if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
	2336	{
	2337	VspKdViewCell *vc = frontLeaf->GetViewCell();
	2338
	2339	VspKdLeaf *leaf = new VspKdLeaf(interior->GetParent(), 0);
	2340	leaf->SetViewCell(vc);
	2341
	2342	// replace a link from node's parent
	2343	if (leaf->mParent)
	2344	leaf->mParent->ReplaceChildLink(node, leaf);
	2345
	2346	delete interior;
	2347
	2348	return leaf;
	2349	}
	2350	}
	2351
	2352	return node;
	2353	}
	2354
	2355
[473]	2356	void VspKdTree::RefineViewCells()
	2357	{
	2358	}
	2359
	2360
[452]	2361	/*********************************************************************/
	2362	/* MergeCandidate implementation */
	2363	/*********************************************************************/
	2364
	2365
[462]	2366	MergeCandidate::MergeCandidate(VspKdLeaf l1, VspKdLeaf l2):
	2367	mMergeCost(0),
	2368	mLeaf1(l1),
	2369	mLeaf2(l2),
	2370	mLeaf1Id(l1->GetViewCell()->mMailbox),
	2371	mLeaf2Id(l2->GetViewCell()->mMailbox)
	2372	{
	2373	EvalMergeCost();
	2374	}
[452]	2375
[472]	2376	float MergeCandidate::GetLeaf1Cost() const
	2377	{
	2378	VspKdViewCell *vc = mLeaf1->GetViewCell();
	2379	return vc->GetPvs().GetSize() * vc->GetVolume();
	2380	}
[462]	2381
[472]	2382	float MergeCandidate::GetLeaf2Cost() const
	2383	{
	2384	VspKdViewCell *vc = mLeaf2->GetViewCell();
	2385	return vc->GetPvs().GetSize() * vc->GetVolume();
	2386	}
	2387
[462]	2388	void MergeCandidate::EvalMergeCost()
[452]	2389	{
[465]	2390	//-- compute pvs difference
[463]	2391	VspKdViewCell *vc1 = mLeaf1->GetViewCell();
[462]	2392	VspKdViewCell *vc2 = mLeaf2->GetViewCell();
[452]	2393
[463]	2394	const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
[462]	2395	const int vcPvs = diff1 + vc1->GetPvs().GetSize();
[453]	2396
[465]	2397	//-- compute ratio of old cost
[462]	2398	//-- (added size of left and right view cell times pvs size)
	2399	//-- to new rendering cost (size of merged view cell times pvs size)
[472]	2400	const float oldCost = GetLeaf1Cost() + GetLeaf2Cost();
	2401
[465]	2402	const float newCost =
[469]	2403	(float)vcPvs * (vc1->GetVolume() + vc2->GetVolume());
[453]	2404
[465]	2405	mMergeCost = newCost - oldCost;
[462]	2406
[465]	2407	// if (vcPvs > sMaxPvsSize) // strong penalty if pvs size too large
	2408	// mMergeCost += 1.0;
[462]	2409	}
	2410
	2411	void MergeCandidate::SetLeaf1(VspKdLeaf *l)
	2412	{
	2413	mLeaf1 = l;
	2414	}
	2415
	2416	void MergeCandidate::SetLeaf2(VspKdLeaf *l)
	2417	{
	2418	mLeaf2 = l;
	2419	}
	2420
[482]	2421
[462]	2422	VspKdLeaf *MergeCandidate::GetLeaf1()
	2423	{
	2424	return mLeaf1;
	2425	}
	2426
[482]	2427
[462]	2428	VspKdLeaf *MergeCandidate::GetLeaf2()
	2429	{
	2430	return mLeaf2;
	2431	}
	2432
[482]	2433
[462]	2434	bool MergeCandidate::Valid() const
[465]	2435	{
	2436	return
	2437	(mLeaf1->GetViewCell()->mMailbox == mLeaf1Id) &&
[462]	2438	(mLeaf2->GetViewCell()->mMailbox == mLeaf2Id);
	2439	}
	2440
[482]	2441
[462]	2442	float MergeCandidate::GetMergeCost() const
	2443	{
	2444	return mMergeCost;
	2445	}
	2446
[482]	2447
[462]	2448	void MergeCandidate::SetValid()
	2449	{
	2450	mLeaf1Id = mLeaf1->GetViewCell()->mMailbox;
	2451	mLeaf2Id = mLeaf2->GetViewCell()->mMailbox;
[465]	2452
[462]	2453	EvalMergeCost();
[463]	2454	}

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

Download in other formats: