Context Navigation

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

Revision 426, 39.5 KB checked in by mattausch, 19 years ago (diff)
fixed ray bug in vspkdtree added visualizations

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

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

Download in other formats: