Context Navigation

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

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

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

Download in other formats: