Context Navigation

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

Revision 452, 43.4 KB checked in by mattausch, 19 years ago (diff)
started to add view cell merging to vsp kd tree

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

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

Download in other formats: