Context Navigation

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

Revision 453, 44.0 KB checked in by mattausch, 19 years ago (diff)

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

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

Download in other formats: