Context Navigation

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

Revision 449, 41.4 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
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)
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	{Debug << "here888" << endl;
450	mStat.Start();
451	Debug << "here1000" << endl;
452	mMaxMemory = mMaxStaticMemory;
453	Debug << "here1" << endl;
454	DEL_PTR(mRoot);
455
456	mRoot = new VspKdTreeLeaf(NULL, (int)rays.size());
457
458	// first construct a leaf that will get subdivided
459	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(mRoot);
460
461	mStat.nodes = 1;
462	mBox.Initialize();
463	Debug << "here2" << endl;
464	//-- compute bounding box
465	if (forcedBoundingBox)
466	mBox = *forcedBoundingBox;
467	else
468	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
469	{
470	mBox.Include((*ri)->GetOrigin());
471	mBox.Include((*ri)->GetTermination());
472	}
473
474	Debug << "bbox = " << mBox << endl;
475
476	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
477	{
478	//leaf->AddRay(RayInfo(*ri));
479	VssRay ray = ri;
480
481	float minT, maxT;
482
483	// TODO: not very efficient to implictly cast between rays types ...
484	if (mBox.GetRaySegment(*ray, minT, maxT))
485	{
486	float len = ray->Length();
487
488	if (!len)
489	len = Limits::Small;
490
491	leaf->AddRay(RayInfo(ray, minT / len, maxT / len));
492	}
493	}
494
495	mStat.rays = (int)leaf->mRays.size();
496	leaf->UpdatePvsSize();
497
498	mStat.initialPvsSize = leaf->GetPvsSize();
499
500	// Subdivide();
501	mRoot = Subdivide(TraversalData(leaf, mBox, 0));
502
503	if (mSplitCandidates)
504	{
505	// force release of this vector
506	delete mSplitCandidates;
507	mSplitCandidates = new vector<SortableEntry>;
508	}
509
510	mStat.Stop();
511
512	mStat.Print(cout);
513	Debug << "#Total memory=" << GetMemUsage() << endl;
514	}
515
516
517
518	VspKdTreeNode *VspKdTree::Subdivide(const TraversalData &tdata)
519	{
520	VspKdTreeNode *result = NULL;
521
522	priority_queue<TraversalData> tStack;
523	//stack<TraversalData> tStack;
524
525	tStack.push(tdata);
526
527	AxisAlignedBox3 backBox;
528	AxisAlignedBox3 frontBox;
529
530	int lastMem = 0;
531
532	while (!tStack.empty())
533	{
534	float mem = GetMemUsage();
535
536	if (lastMem / 10 != ((int)mem) / 10)
537	{
538	Debug << mem << " MB" << endl;
539	}
540	lastMem = (int)mem;
541
542	if (1 && mem > mMaxMemory)
543	{
544	Debug << "memory limit reached: " << mem << endl;
545	// count statistics on unprocessed leafs
546	while (!tStack.empty())
547	{
548	EvaluateLeafStats(tStack.top());
549	tStack.pop();
550	}
551	break;
552	}
553
554	TraversalData data = tStack.top();
555	tStack.pop();
556
557	VspKdTreeNode node = SubdivideNode((VspKdTreeLeaf ) data.mNode,
558	data.mBox, backBox, frontBox);
559	if (result == NULL)
560	result = node;
561
562	if (!node->IsLeaf())
563	{
564	VspKdTreeInterior interior = dynamic_cast<VspKdTreeInterior >(node);
565
566	// push the children on the stack
567	tStack.push(TraversalData(interior->GetBack(), backBox, data.mDepth + 1));
568	tStack.push(TraversalData(interior->GetFront(), frontBox, data.mDepth + 1));
569	}
570	else
571	{
572	EvaluateLeafStats(data);
573	}
574	}
575
576	return result;
577	}
578
579
580	// returns selected plane for subdivision
581	int VspKdTree::SelectPlane(VspKdTreeLeaf *leaf,
582	const AxisAlignedBox3 &box,
583	float &position,
584	int &raysBack,
585	int &raysFront,
586	int &pvsBack,
587	int &pvsFront)
588	{
589	int minDirDepth = 6;
590	int axis = 0;
591	float costRatio = 0;
592
593	if (splitType == ESplitRegular)
594	{
595	costRatio = BestCostRatioRegular(leaf,
596	axis,
597	position,
598	raysBack,
599	raysFront,
600	pvsBack,
601	pvsFront);
602	}
603	else if (splitType == ESplitHeuristic)
604	{
605	costRatio = BestCostRatioHeuristic(leaf,
606	axis,
607	position,
608	raysBack,
609	raysFront,
610	pvsBack,
611	pvsFront);
612	}
613	else
614	{
615	cerr << "VspKdTree: Unknown split heuristics\n";
616	exit(1);
617	}
618
619	if (costRatio > mTermMaxCostRatio)
620	{
621	Debug << "Too big cost ratio " << costRatio << endl;
622	return -1;
623	}
624
625	if (0)
626	Debug << "pvs=" << leaf->mPvsSize
627	<< " rays=" << (int)leaf->mRays.size()
628	<< " rc=" << leaf->GetAvgRayContribution()
629	<< " axis=" << axis << endl;
630
631	return axis;
632	}
633
634
635
636	float VspKdTree::EvalCostRatio(VspKdTreeLeaf *leaf,
637	const int axis,
638	const float position,
639	int &raysBack,
640	int &raysFront,
641	int &pvsBack,
642	int &pvsFront)
643	{
644	raysBack = 0;
645	raysFront = 0;
646	pvsFront = 0;
647	pvsBack = 0;
648
649	float newCost;
650
651	Intersectable::NewMail(3);
652
653	// eval pvs size
654	int pvsSize = leaf->GetPvsSize();
655
656	Intersectable::NewMail(3);
657
658	// this is the main ray classification loop!
659	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
660	ri != leaf->mRays.end(); ++ ri)
661	{
662	if (!(*ri).mRay->IsActive())
663	continue;
664
665	// determine the side of this ray with respect to the plane
666	int side = (ri).ComputeRayIntersection(axis, position, (ri).mRay->mT);
667	// (*ri).mRay->mSide = side;
668
669	if (side <= 0)
670	++ raysBack;
671
672	if (side >= 0)
673	++ raysFront;
674
675	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
676	}
677
678	float ratio = 0;
679
680	if (0)
681	{
682	const float sum = float(pvsBack + pvsFront);
683	const float oldCost = (float)pvsSize * 2;
684
685	return sum / oldCost;
686	}
687	else
688	{
689	AxisAlignedBox3 box = GetBBox(leaf);
690
691	float minBox = box.Min(axis);
692	float maxBox = box.Max(axis);
693	float sizeBox = maxBox - minBox;
694
695	// float sum = raysBack(position - minBox) + raysFront(maxBox - position);
696	const float sum = pvsBack * (position - minBox) + pvsFront * (maxBox - position);
697
698	newCost = mCtDivCi + sum / sizeBox;
699
700	//Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
701	// float oldCost = leaf->mRays.size();
702	const float oldCost = (float)pvsSize;
703
704	float ratio = newCost / oldCost;
705	}
706
707	return ratio;
708	}
709
710	float VspKdTree::BestCostRatioRegular(VspKdTreeLeaf *leaf,
711	int &axis,
712	float &position,
713	int &raysBack,
714	int &raysFront,
715	int &pvsBack,
716	int &pvsFront)
717	{
718	int nRaysBack[3], nRaysFront[3];
719	int nPvsBack[3], nPvsFront[3];
720
721	float nPosition[3];
722	float nCostRatio[3];
723	int bestAxis = -1;
724
725	AxisAlignedBox3 sBox = GetBBox(leaf);
726
727	const int sAxis = sBox.Size().DrivingAxis();
728
729	for (axis = 0; axis < 3; ++ axis)
730	{
731	if (!mOnlyDrivingAxis \|\| axis == sAxis)
732	{
733	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis]) * 0.5f;
734
735	nCostRatio[axis] = EvalCostRatio(leaf,
736	axis,
737	nPosition[axis],
738	nRaysBack[axis],
739	nRaysFront[axis],
740	nPvsBack[axis],
741	nPvsFront[axis]);
742
743	if (bestAxis == -1)
744	bestAxis = axis;
745	else
746	if (nCostRatio[axis] < nCostRatio[bestAxis])
747	{
748	Debug << "pvs front " << nPvsBack[axis]
749	<< " pvs back " << nPvsFront[axis]
750	<< " overall pvs " << leaf->GetPvsSize() << endl;
751	bestAxis = axis;
752	}
753	}
754	}
755
756	//-- assign best axis
757	axis = bestAxis;
758	position = nPosition[bestAxis];
759
760	raysBack = nRaysBack[bestAxis];
761	raysFront = nRaysFront[bestAxis];
762
763	pvsBack = nPvsBack[bestAxis];
764	pvsFront = nPvsFront[bestAxis];
765
766	return nCostRatio[bestAxis];
767	}
768
769	float VspKdTree::BestCostRatioHeuristic(VspKdTreeLeaf *leaf,
770	int &axis,
771	float &position,
772	int &raysBack,
773	int &raysFront,
774	int &pvsBack,
775	int &pvsFront)
776	{
777	AxisAlignedBox3 box = GetBBox(leaf);
778	// AxisAlignedBox3 dirBox = GetDirBBox(node);
779
780	axis = box.Size().DrivingAxis();
781
782	SortSplitCandidates(leaf, axis);
783
784	// go through the lists, count the number of objects left and right
785	// and evaluate the following cost funcion:
786	// C = ct_div_ci + (qlrl + qrrr)/queries
787
788	int rl = 0, rr = (int)leaf->mRays.size();
789	int pl = 0, pr = leaf->GetPvsSize();
790
791	float minBox = box.Min(axis);
792	float maxBox = box.Max(axis);
793	float sizeBox = maxBox - minBox;
794
795	float minBand = minBox + 0.1f*(maxBox - minBox);
796	float maxBand = minBox + 0.9f*(maxBox - minBox);
797
798	float sum = rr*sizeBox;
799	float minSum = 1e20f;
800
801	Intersectable::NewMail();
802
803	// set all object as belonging to the fron pvs
804	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
805	ri != leaf->mRays.end(); ++ ri)
806	{
807	if ((*ri).mRay->IsActive())
808	{
809	Intersectable object = (ri).mRay->mTerminationObject;
810
811	if (object)
812	if (!object->Mailed())
813	{
814	object->Mail();
815	object->mCounter = 1;
816	}
817	else
818	++ object->mCounter;
819	}
820	}
821
822	Intersectable::NewMail();
823
824	for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
825	ci < mSplitCandidates->end(); ++ ci)
826	{
827	VssRay *ray;
828
829	switch ((*ci).type)
830	{
831	case SortableEntry::ERayMin:
832	{
833	++ rl;
834	ray = (VssRay ) (ci).data;
835	Intersectable *object = ray->mTerminationObject;
836	if (object && !object->Mailed())
837	{
838	object->Mail();
839	++ pl;
840	}
841	break;
842	}
843	case SortableEntry::ERayMax:
844	{
845	-- rr;
846
847	ray = (VssRay ) (ci).data;
848	Intersectable *object = ray->mTerminationObject;
849
850	if (object)
851	{
852	if (-- object->mCounter == 0)
853	-- pr;
854	}
855
856	break;
857	}
858	}
859
860	if ((ci).value > minBand && (ci).value < maxBand)
861	{
862	sum = pl((ci).value - minBox) + pr(maxBox - (ci).value);
863
864	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
865	// cout<<"cost= "<<sum<<endl;
866
867	if (sum < minSum)
868	{
869	minSum = sum;
870	position = (*ci).value;
871
872	raysBack = rl;
873	raysFront = rr;
874
875	pvsBack = pl;
876	pvsFront = pr;
877
878	}
879	}
880	}
881
882	float oldCost = (float)leaf->GetPvsSize();
883	float newCost = mCtDivCi + minSum / sizeBox;
884	float ratio = newCost / oldCost;
885
886	//Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
887	// <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
888
889	return ratio;
890	}
891
892	void VspKdTree::SortSplitCandidates(VspKdTreeLeaf *node,
893	const int axis)
894	{
895	mSplitCandidates->clear();
896
897	int requestedSize = 2 * (int)(node->mRays.size());
898	// creates a sorted split candidates array
899	if (mSplitCandidates->capacity() > 500000 &&
900	requestedSize < (int)(mSplitCandidates->capacity()/10) )
901	{
902	DEL_PTR(mSplitCandidates);
903	mSplitCandidates = new vector<SortableEntry>;
904	}
905
906	mSplitCandidates->reserve(requestedSize);
907
908	// insert all queries
909	for(RayInfoContainer::const_iterator ri = node->mRays.begin();
910	ri < node->mRays.end(); ++ ri)
911	{
912	bool positive = (*ri).mRay->HasPosDir(axis);
913	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
914	(ri).ExtrapOrigin(axis), (void )&*ri));
915
916	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
917	(ri).ExtrapTermination(axis), (void )&*ri));
918	}
919
920	stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
921	}
922
923
924	void VspKdTree::EvaluateLeafStats(const TraversalData &data)
925	{
926	// the node became a leaf -> evaluate stats for leafs
927	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(data.mNode);
928
929	if (leaf->GetPvsSize() > mStat.maxPvsSize)
930	mStat.maxPvsSize = leaf->GetPvsSize();
931
932	if ((int)(leaf->mRays.size()) > mStat.maxRayRefs)
933	mStat.maxRayRefs = (int)leaf->mRays.size();
934
935
936	if (data.mDepth >= mTermMaxDepth)
937	++ mStat.maxDepthNodes;
938
939	if (leaf->GetPvsSize() < mTermMinPvs)
940	++ mStat.minPvsNodes;
941
942	if ((int)leaf->GetRays().size() < mTermMinRays)
943	++ mStat.minRaysNodes;
944
945	if (leaf->GetAvgRayContribution() > mTermMaxRayContribution)
946	++ mStat.maxRayContribNodes;
947
948	if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
949	++ mStat.minSizeNodes;
950	}
951
952
953	inline bool VspKdTree::TerminationCriteriaMet(const VspKdTreeLeaf *leaf,
954	const AxisAlignedBox3 &box) const
955	{
956	return ((leaf->GetPvsSize() < mTermMinPvs) \|\|
957	(leaf->mRays.size() < mTermMinRays) \|\|
958	(leaf->GetAvgRayContribution() > mTermMaxRayContribution ) \|\|
959	(leaf->mDepth >= mTermMaxDepth) \|\|
960	(SqrMagnitude(box.Size()) <= mTermMinSize));
961	}
962
963
964	VspKdTreeNode VspKdTree::SubdivideNode(VspKdTreeLeaf leaf,
965	const AxisAlignedBox3 &box,
966	AxisAlignedBox3 &backBBox,
967	AxisAlignedBox3 &frontBBox)
968	{
969	if (TerminationCriteriaMet(leaf, box))
970	{
971	if (1 && leaf->mDepth >= mTermMaxDepth)
972	{
973	Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
974	Debug << "Bbox: " << GetBBox(leaf) << endl;
975	}
976	//Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
977
978	return leaf;
979	}
980
981	float position;
982	// first count ray sides
983	int raysBack;
984	int raysFront;
985	int pvsBack;
986	int pvsFront;
987
988	// select subdivision axis
989	const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
990
991	//Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" << pvsFront << endl;
992
993	if (axis == -1)
994	{
995	return leaf;
996	}
997
998	mStat.nodes += 2;
999	//++ mStat.splits[axis];
1000
1001	// add the new nodes to the tree
1002	VspKdTreeInterior *node = new VspKdTreeInterior(leaf->mParent);
1003
1004	node->mAxis = axis;
1005	node->mPosition = position;
1006	node->mBox = box;
1007
1008	backBBox = box;
1009	frontBBox = box;
1010
1011	VspKdTreeLeaf *back = new VspKdTreeLeaf(node, raysBack);
1012	back->SetPvsSize(pvsBack);
1013	VspKdTreeLeaf *front = new VspKdTreeLeaf(node, raysFront);
1014	front->SetPvsSize(pvsFront);
1015
1016	// replace a link from node's parent
1017	if (leaf->mParent)
1018	leaf->mParent->ReplaceChildLink(leaf, node);
1019	// and setup child links
1020	node->SetupChildLinks(back, front);
1021
1022	if (axis <= VspKdTreeNode::SPLIT_Z)
1023	{
1024	backBBox.SetMax(axis, position);
1025	frontBBox.SetMin(axis, position);
1026
1027	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1028	ri != leaf->mRays.end(); ++ ri)
1029	{
1030	if ((*ri).mRay->IsActive())
1031	{
1032	// first unref ray from the former leaf
1033	(*ri).mRay->Unref();
1034
1035	// determine the side of this ray with respect to the plane
1036	int side = node->ComputeRayIntersection(ri, (ri).mRay->mT);
1037
1038	if (side == 0)
1039	{
1040	if ((*ri).mRay->HasPosDir(axis))
1041	{
1042	back->AddRay(RayInfo((*ri).mRay,
1043	(*ri).mMinT,
1044	(*ri).mRay->mT));
1045	front->AddRay(RayInfo((*ri).mRay,
1046	(*ri).mRay->mT,
1047	(*ri).mMaxT));
1048	}
1049	else
1050	{
1051	back->AddRay(RayInfo((*ri).mRay,
1052	(*ri).mRay->mT,
1053	(*ri).mMaxT));
1054	front->AddRay(RayInfo((*ri).mRay,
1055	(*ri).mMinT,
1056	(*ri).mRay->mT));
1057	}
1058	}
1059	else
1060	if (side == 1)
1061	front->AddRay(*ri);
1062	else
1063	back->AddRay(*ri);
1064	} else
1065	(*ri).mRay->Unref();
1066	}
1067	}
1068	else
1069	{
1070	// rays front/back
1071
1072	for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1073	ri != leaf->mRays.end(); ++ ri)
1074	{
1075	if ((*ri).mRay->IsActive())
1076	{
1077	// first unref ray from the former leaf
1078	(*ri).mRay->Unref();
1079
1080	int side;
1081	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1082	side = 1;
1083	else
1084	side = -1;
1085
1086	if (side == 1)
1087	front->AddRay(*ri);
1088	else
1089	back->AddRay(*ri);
1090	}
1091	else
1092	(*ri).mRay->Unref();
1093	}
1094	}
1095
1096	// update stats
1097	mStat.rayRefs -= (int)leaf->mRays.size();
1098	mStat.rayRefs += raysBack + raysFront;
1099
1100	DEL_PTR(leaf);
1101
1102	return node;
1103	}
1104
1105	int VspKdTree::ReleaseMemory(const int time)
1106	{
1107	stack<VspKdTreeNode *> tstack;
1108
1109	// find a node in the tree which subtree will be collapsed
1110	int maxAccessTime = time - mAccessTimeThreshold;
1111	int released = 0;
1112
1113	tstack.push(mRoot);
1114
1115	while (!tstack.empty())
1116	{
1117	VspKdTreeNode *node = tstack.top();
1118	tstack.pop();
1119
1120	if (!node->IsLeaf())
1121	{
1122	VspKdTreeInterior in = dynamic_cast<VspKdTreeInterior >(node);
1123	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1124	if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
1125	{
1126	released = CollapseSubtree(node, time);
1127	break;
1128	}
1129	if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
1130	{
1131	tstack.push(in->GetFront());
1132	tstack.push(in->GetBack());
1133	}
1134	else
1135	{
1136	tstack.push(in->GetBack());
1137	tstack.push(in->GetFront());
1138	}
1139	}
1140	}
1141
1142	while (tstack.empty())
1143	{
1144	// could find node to collaps...
1145	// cout<<"Could not find a node to release "<<endl;
1146	break;
1147	}
1148
1149	return released;
1150	}
1151
1152
1153	VspKdTreeNode VspKdTree::SubdivideLeaf(VspKdTreeLeaf leaf,
1154	const float sizeThreshold)
1155	{
1156	VspKdTreeNode *node = leaf;
1157
1158	AxisAlignedBox3 leafBBox = GetBBox(leaf);
1159
1160	static int pass = 0;
1161	++ pass;
1162
1163	// check if we should perform a dynamic subdivision of the leaf
1164	if (// leaf->mRays.size() > (unsigned)termMinCost &&
1165	(leaf->GetPvsSize() >= mTermMinPvs) &&
1166	(SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
1167	{
1168	// memory check and release
1169	if (GetMemUsage() > mMaxTotalMemory)
1170	ReleaseMemory(pass);
1171
1172	AxisAlignedBox3 backBBox, frontBBox;
1173
1174	// subdivide the node
1175	node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
1176	}
1177	return node;
1178	}
1179
1180
1181
1182	void VspKdTree::UpdateRays(VssRayContainer &remove,
1183	VssRayContainer &add)
1184	{
1185	VspKdTreeLeaf::NewMail();
1186
1187	// schedule rays for removal
1188	for(VssRayContainer::const_iterator ri = remove.begin();
1189	ri != remove.end(); ++ ri)
1190	{
1191	(*ri)->ScheduleForRemoval();
1192	}
1193
1194	int inactive = 0;
1195
1196	for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
1197	{
1198	if ((*ri)->ScheduledForRemoval())
1199	// RemoveRay(*ri, NULL, false);
1200	// !!! BUG - with true it does not work correctly - aggreated delete
1201	RemoveRay(*ri, NULL, true);
1202	else
1203	++ inactive;
1204	}
1205
1206	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1207	for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
1208	{
1209	AddRay(*ri);
1210	}
1211	}
1212
1213
1214	void VspKdTree::RemoveRay(VssRay *ray,
1215	vector<VspKdTreeLeaf > affectedLeaves,
1216	const bool removeAllScheduledRays)
1217	{
1218	stack<RayTraversalData> tstack;
1219
1220	tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1221
1222	RayTraversalData data;
1223
1224	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1225	while (!tstack.empty())
1226	{
1227	data = tstack.top();
1228	tstack.pop();
1229
1230	if (!data.mNode->IsLeaf())
1231	{
1232	// split the set of rays in two groups intersecting the
1233	// two subtrees
1234	TraverseInternalNode(data, tstack);
1235	}
1236	else
1237	{
1238	// remove the ray from the leaf
1239	// find the ray in the leaf and swap it with the last ray...
1240	VspKdTreeLeaf leaf = (VspKdTreeLeaf )data.mNode;
1241
1242	if (!leaf->Mailed())
1243	{
1244	leaf->Mail();
1245
1246	if (affectedLeaves)
1247	affectedLeaves->push_back(leaf);
1248
1249	if (removeAllScheduledRays)
1250	{
1251	int tail = (int)leaf->mRays.size() - 1;
1252
1253	for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
1254	{
1255	if (leaf->mRays[i].mRay->ScheduledForRemoval())
1256	{
1257	// find a ray to replace it with
1258	while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
1259	{
1260	++ mStat.removedRayRefs;
1261	leaf->mRays[tail].mRay->Unref();
1262	leaf->mRays.pop_back();
1263
1264	-- tail;
1265	}
1266
1267	if (tail < i)
1268	break;
1269
1270	++ mStat.removedRayRefs;
1271
1272	leaf->mRays[i].mRay->Unref();
1273	leaf->mRays[i] = leaf->mRays[tail];
1274	leaf->mRays.pop_back();
1275	-- tail;
1276	}
1277	}
1278	}
1279	}
1280
1281	if (!removeAllScheduledRays)
1282	for (int i=0; i < (int)leaf->mRays.size(); i++)
1283	{
1284	if (leaf->mRays[i].mRay == ray)
1285	{
1286	++ mStat.removedRayRefs;
1287	ray->Unref();
1288	leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
1289	leaf->mRays.pop_back();
1290	// check this ray again
1291	break;
1292	}
1293	}
1294	}
1295	}
1296
1297	if (ray->RefCount() != 0)
1298	{
1299	cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
1300	exit(1);
1301	}
1302
1303	}
1304
1305	void VspKdTree::AddRay(VssRay *ray)
1306	{
1307	stack<RayTraversalData> tstack;
1308
1309	tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
1310
1311	RayTraversalData data;
1312
1313	while (!tstack.empty())
1314	{
1315	data = tstack.top();
1316	tstack.pop();
1317
1318	if (!data.mNode->IsLeaf())
1319	{
1320	TraverseInternalNode(data, tstack);
1321	}
1322	else
1323	{
1324	// remove the ray from the leaf
1325	// find the ray in the leaf and swap it with the last ray
1326	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(data.mNode);
1327
1328	leaf->AddRay(data.mRayData);
1329	++ mStat.addedRayRefs;
1330	}
1331	}
1332	}
1333
1334	void VspKdTree::TraverseInternalNode(RayTraversalData &data,
1335	stack<RayTraversalData> &tstack)
1336	{
1337	VspKdTreeInterior in = dynamic_cast<VspKdTreeInterior >(data.mNode);
1338
1339	if (in->mAxis <= VspKdTreeNode::SPLIT_Z)
1340	{
1341	// determine the side of this ray with respect to the plane
1342	int side = in->ComputeRayIntersection(data.mRayData, data.mRayData.mRay->mT);
1343
1344	if (side == 0)
1345	{
1346	if (data.mRayData.mRay->HasPosDir(in->mAxis))
1347	{
1348	tstack.push(RayTraversalData(in->GetBack(),
1349	RayInfo(data.mRayData.mRay, data.mRayData.mMinT, data.mRayData.mRay->mT)));
1350
1351	tstack.push(RayTraversalData(in->GetFront(),
1352	RayInfo(data.mRayData.mRay, data.mRayData.mRay->mT, data.mRayData.mMaxT)));
1353
1354	}
1355	else
1356	{
1357	tstack.push(RayTraversalData(in->GetBack(),
1358	RayInfo(data.mRayData.mRay,
1359	data.mRayData.mRay->mT,
1360	data.mRayData.mMaxT)));
1361	tstack.push(RayTraversalData(in->GetFront(),
1362	RayInfo(data.mRayData.mRay,
1363	data.mRayData.mMinT,
1364	data.mRayData.mRay->mT)));
1365	}
1366	}
1367	else
1368	if (side == 1)
1369	tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1370	else
1371	tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1372	}
1373	else
1374	{
1375	// directional split
1376	if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
1377	tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
1378	else
1379	tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
1380	}
1381	}
1382
1383
1384	int VspKdTree::CollapseSubtree(VspKdTreeNode *sroot, const int time)
1385	{
1386	// first count all rays in the subtree
1387	// use mail 1 for this purpose
1388	stack<VspKdTreeNode *> tstack;
1389
1390	int rayCount = 0;
1391	int totalRayCount = 0;
1392	int collapsedNodes = 0;
1393
1394	#if DEBUG_COLLAPSE
1395	cout << "Collapsing subtree" << endl;
1396	cout << "accessTime=" << sroot->GetAccessTime() << endl;
1397	cout << "depth=" << (int)sroot->depth << endl;
1398	#endif
1399
1400	// tstat.collapsedSubtrees++;
1401	// tstat.collapseDepths += (int)sroot->depth;
1402	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1403	tstack.push(sroot);
1404	VssRay::NewMail();
1405
1406	while (!tstack.empty())
1407	{
1408	++ collapsedNodes;
1409
1410	VspKdTreeNode *node = tstack.top();
1411	tstack.pop();
1412	if (node->IsLeaf())
1413	{
1414	VspKdTreeLeaf leaf = (VspKdTreeLeaf ) node;
1415
1416	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1417	ri != leaf->mRays.end(); ++ ri)
1418	{
1419	++ totalRayCount;
1420
1421	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed())
1422	{
1423	(*ri).mRay->Mail();
1424	++ rayCount;
1425	}
1426	}
1427	}
1428	else
1429	{
1430	tstack.push(((VspKdTreeInterior *)node)->GetFront());
1431	tstack.push(((VspKdTreeInterior *)node)->GetBack());
1432	}
1433	}
1434
1435	VssRay::NewMail();
1436
1437	// create a new node that will hold the rays
1438	VspKdTreeLeaf *newLeaf = new VspKdTreeLeaf(sroot->mParent, rayCount);
1439
1440	if (newLeaf->mParent)
1441	newLeaf->mParent->ReplaceChildLink(sroot, newLeaf);
1442
1443	tstack.push(sroot);
1444
1445	while (!tstack.empty())
1446	{
1447	VspKdTreeNode *node = tstack.top();
1448	tstack.pop();
1449
1450	if (node->IsLeaf())
1451	{
1452	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(node);
1453
1454	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
1455	ri != leaf->mRays.end(); ++ ri)
1456	{
1457	// unref this ray from the old node
1458	if ((*ri).mRay->IsActive())
1459	{
1460	(*ri).mRay->Unref();
1461	if (!(*ri).mRay->Mailed())
1462	{
1463	(*ri).mRay->Mail();
1464	newLeaf->AddRay(*ri);
1465	}
1466	}
1467	else
1468	(*ri).mRay->Unref();
1469	}
1470	}
1471	else
1472	{
1473	VspKdTreeInterior *interior =
1474	dynamic_cast<VspKdTreeInterior *>(node);
1475	tstack.push(interior->GetBack());
1476	tstack.push(interior->GetFront());
1477	}
1478	}
1479
1480	// delete the node and all its children
1481	DEL_PTR(sroot);
1482
1483	// for(VspKdTreeNode::SRayContainer::iterator ri = newleaf->mRays.begin();
1484	// ri != newleaf->mRays.end(); ++ ri)
1485	// (*ri).ray->UnMail(2);
1486
1487	#if DEBUG_COLLAPSE
1488	cout << "Total memory before=" << GetMemUsage() << endl;
1489	#endif
1490
1491	mStat.nodes -= collapsedNodes - 1;
1492	mStat.rayRefs -= totalRayCount - rayCount;
1493
1494	#if DEBUG_COLLAPSE
1495	cout << "collapsed nodes" << collapsedNodes << endl;
1496	cout << "collapsed rays" << totalRayCount - rayCount << endl;
1497	cout << "Total memory after=" << GetMemUsage() << endl;
1498	cout << "================================" << endl;
1499	#endif
1500
1501	// tstat.collapsedNodes += collapsedNodes;
1502	// tstat.collapsedRays += totalRayCount - rayCount;
1503	return totalRayCount - rayCount;
1504	}
1505
1506
1507	int VspKdTree::GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const
1508	{
1509	stack<VspKdTreeNode *> tstack;
1510	tstack.push(mRoot);
1511
1512	Intersectable::NewMail();
1513	int pvsSize = 0;
1514
1515	while (!tstack.empty())
1516	{
1517	VspKdTreeNode *node = tstack.top();
1518	tstack.pop();
1519
1520	if (node->IsLeaf())
1521	{
1522	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
1523	for (RayInfoContainer::iterator ri = leaf->mRays.begin();
1524	ri != leaf->mRays.end(); ++ ri)
1525	{
1526	if ((*ri).mRay->IsActive())
1527	{
1528	Intersectable *object;
1529	#if BIDIRECTIONAL_RAY
1530	object = (*ri).mRay->mOriginObject;
1531
1532	if (object && !object->Mailed())
1533	{
1534	++ pvsSize;
1535	object->Mail();
1536	}
1537	#endif
1538	object = (*ri).mRay->mTerminationObject;
1539	if (object && !object->Mailed())
1540	{
1541	++ pvsSize;
1542	object->Mail();
1543	}
1544	}
1545	}
1546	}
1547	else
1548	{
1549	VspKdTreeInterior in = dynamic_cast<VspKdTreeInterior >(node);
1550
1551	if (in->mAxis < 3)
1552	{
1553	if (box.Max(in->mAxis) >= in->mPosition)
1554	tstack.push(in->GetFront());
1555
1556	if (box.Min(in->mAxis) <= in->mPosition)
1557	tstack.push(in->GetBack());
1558	}
1559	else
1560	{
1561	// both nodes for directional splits
1562	tstack.push(in->GetFront());
1563	tstack.push(in->GetBack());
1564	}
1565	}
1566	}
1567
1568	return pvsSize;
1569	}
1570
1571	void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
1572	float &maxRayContribution,
1573	float &avgRayContribution)
1574	{
1575	stack<VspKdTreeNode *> tstack;
1576	tstack.push(mRoot);
1577
1578	minRayContribution = 1.0f;
1579	maxRayContribution = 0.0f;
1580
1581	float sumRayContribution = 0.0f;
1582	int leaves = 0;
1583
1584	while (!tstack.empty())
1585	{
1586	VspKdTreeNode *node = tstack.top();
1587	tstack.pop();
1588	if (node->IsLeaf())
1589	{
1590	leaves++;
1591	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
1592	float c = leaf->GetAvgRayContribution();
1593
1594	if (c > maxRayContribution)
1595	maxRayContribution = c;
1596	if (c < minRayContribution)
1597	minRayContribution = c;
1598	sumRayContribution += c;
1599
1600	}
1601	else {
1602	VspKdTreeInterior in = (VspKdTreeInterior )node;
1603	// both nodes for directional splits
1604	tstack.push(in->GetFront());
1605	tstack.push(in->GetBack());
1606	}
1607	}
1608
1609	Debug << "sum=" << sumRayContribution << endl;
1610	Debug << "leaves=" << leaves << endl;
1611	avgRayContribution = sumRayContribution / (float)leaves;
1612	}
1613
1614
1615	int VspKdTree::GenerateRays(const float ratioPerLeaf,
1616	SimpleRayContainer &rays)
1617	{
1618	stack<VspKdTreeNode *> tstack;
1619	tstack.push(mRoot);
1620
1621	while (!tstack.empty())
1622	{
1623	VspKdTreeNode *node = tstack.top();
1624	tstack.pop();
1625
1626	if (node->IsLeaf())
1627	{
1628	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(node);
1629
1630	float c = leaf->GetAvgRayContribution();
1631	int num = (int)(c*ratioPerLeaf + 0.5);
1632	// cout<<num<<" ";
1633
1634	for (int i=0; i < num; i++)
1635	{
1636	Vector3 origin = GetBBox(leaf).GetRandomPoint();
1637	/*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1638	Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1639	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1640	rays.push_back(SimpleRay(origin, direction));*/
1641	}
1642
1643	}
1644	else
1645	{
1646	VspKdTreeInterior *in =
1647	dynamic_cast<VspKdTreeInterior *>(node);
1648	// both nodes for directional splits
1649	tstack.push(in->GetFront());
1650	tstack.push(in->GetBack());
1651	}
1652	}
1653
1654	return (int)rays.size();
1655	}
1656
1657
1658	float VspKdTree::GetAvgPvsSize()
1659	{
1660	stack<VspKdTreeNode *> tstack;
1661	tstack.push(mRoot);
1662
1663	int sumPvs = 0;
1664	int leaves = 0;
1665
1666	while (!tstack.empty())
1667	{
1668	VspKdTreeNode *node = tstack.top();
1669	tstack.pop();
1670
1671	if (node->IsLeaf())
1672	{
1673	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
1674	// update pvs size
1675	leaf->UpdatePvsSize();
1676	sumPvs += leaf->GetPvsSize();
1677	leaves++;
1678	}
1679	else
1680	{
1681	VspKdTreeInterior in = (VspKdTreeInterior )node;
1682	// both nodes for directional splits
1683	tstack.push(in->GetFront());
1684	tstack.push(in->GetBack());
1685	}
1686	}
1687
1688	return sumPvs / (float)leaves;
1689	}
1690
1691	VspKdTreeNode *VspKdTree::GetRoot() const
1692	{
1693	return mRoot;
1694	}
1695
1696	AxisAlignedBox3 VspKdTree::GetBBox(VspKdTreeNode *node) const
1697	{
1698	if (node->mParent == NULL)
1699	return mBox;
1700
1701	if (!node->IsLeaf())
1702	return (dynamic_cast<VspKdTreeInterior *>(node))->mBox;
1703
1704	if (node->mParent->mAxis >= 3)
1705	return node->mParent->mBox;
1706
1707	AxisAlignedBox3 box(node->mParent->mBox);
1708	if (node->mParent->GetFront() == node)
1709	box.SetMin(node->mParent->mAxis, node->mParent->mPosition);
1710	else
1711	box.SetMax(node->mParent->mAxis, node->mParent->mPosition);
1712
1713	return box;
1714	}
1715
1716	int VspKdTree::GetRootPvsSize() const
1717	{
1718	return GetPvsSize(mRoot, mBox);
1719	}
1720
1721	const VspKdStatistics &VspKdTree::GetStatistics() const
1722	{
1723	return mStat;
1724	}
1725
1726	void VspKdTree::AddRays(VssRayContainer &add)
1727	{
1728	VssRayContainer remove;
1729	UpdateRays(remove, add);
1730	}
1731
1732	// return memory usage in MB
1733	float VspKdTree::GetMemUsage() const
1734	{
1735	return
1736	(sizeof(VspKdTree) +
1737	mStat.Leaves() * sizeof(VspKdTreeLeaf) +
1738	mStat.Interior() * sizeof(VspKdTreeInterior) +
1739	mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
1740	}
1741
1742	float VspKdTree::GetRayMemUsage() const
1743	{
1744	return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
1745	}
1746
1747	int VspKdTree::ComputePvsSize(VspKdTreeNode *node,
1748	const RayInfoContainer &globalRays) const
1749	{
1750	int pvsSize = 0;
1751
1752	const AxisAlignedBox3 box = GetBBox(node);
1753
1754	RayInfoContainer::const_iterator it, it_end = globalRays.end();
1755
1756	// warning: implicit conversion from VssRay to Ray
1757	for (it = globalRays.begin(); it != globalRays.end(); ++ it)
1758	pvsSize += box.GetMinMaxT((it).mRay, NULL, NULL);
1759
1760	return pvsSize;
1761	}
1762
1763	void VspKdTree::CollectLeaves(vector<VspKdTreeLeaf *> &leaves) const
1764	{
1765	stack<VspKdTreeNode *> nodeStack;
1766	nodeStack.push(mRoot);
1767
1768	while (!nodeStack.empty())
1769	{
1770	VspKdTreeNode *node = nodeStack.top();
1771
1772	nodeStack.pop();
1773
1774	if (node->IsLeaf())
1775	{
1776	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(node);
1777	leaves.push_back(leaf);
1778	}
1779	else
1780	{
1781	VspKdTreeInterior interior = dynamic_cast<VspKdTreeInterior >(node);
1782
1783	nodeStack.push(interior->GetBack());
1784	nodeStack.push(interior->GetFront());
1785	}
1786	}
1787	}

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

Download in other formats: