Context Navigation

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

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

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

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

Download in other formats: