Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/VssTree.cpp @ 1233

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

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

Download in other formats: