Context Navigation

source: trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp @ 403

Revision 403, 35.8 KB checked in by bittner, 19 years ago (diff)
vvs preprocessor update

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

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

Download in other formats: