Context Navigation

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

Revision 1359, 47.2 KB checked in by mattausch, 18 years ago (diff)
worked on global object sorting

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, (int)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 = (int)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 = (float)pvsSize;
546	ratio = newCost/oldCost;
547
548	} else {
549	// importance based cost
550	#if 0
551	float newContrib =
552	((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
553	(maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
554
555	// float newContrib =
556	// sqr(pvsBack/(raysBack + Limits::Small)) +
557	// sqr(pvsFront/(raysFront + Limits::Small));
558	float oldContrib = sqr(leaf->GetAvgRayContribution());
559	ratio = oldContrib/newContrib;
560	#else
561	#if 1
562	float newCost = float(raysBackpvsBack + raysFrontpvsFront);
563	float oldCost = (float)leaf->rays.size()*pvsSize;
564	ratio = newCost/oldCost;
565	#else
566	float newCost = (pvsBack + pvsFront)*0.5f;
567	float oldCost = pvsSize;
568	ratio = newCost/oldCost;
569	#endif
570	#endif
571	}
572
573	return ratio;
574	}
575
576
577	float
578	VssTree::EvalCostRatio(
579	VssTreeLeaf *leaf,
580	const int axis,
581	const float position,
582	int &raysBack,
583	int &raysFront,
584	int &pvsBack,
585	int &pvsFront
586	)
587	{
588	raysBack = 0;
589	raysFront = 0;
590	pvsFront = 0;
591	pvsBack = 0;
592
593
594	Intersectable::NewMail(3);
595
596	if (axis <= VssTreeNode::SPLIT_Z) {
597	// this is the main ray classification loop!
598	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
599	ri != leaf->rays.end();
600	ri++)
601	if ((*ri).mRay->IsActive()) {
602	float t;
603	// determine the side of this ray with respect to the plane
604	int side = (*ri).ComputeRayIntersection(axis, position, t);
605	// (*ri).mRay->mSide = side;
606
607	if (side <= 0)
608	raysBack++;
609
610	if (side >= 0)
611	raysFront++;
612
613	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
614	}
615
616	} else {
617
618	// directional split
619	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
620	ri != leaf->rays.end();
621	ri++)
622	if ((*ri).mRay->IsActive()) {
623
624	// determine the side of this ray with respect to the plane
625	int side;
626	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
627	side = 1;
628	else
629	side = -1;
630
631	if (side <= 0)
632	raysBack++;
633
634	if (side >= 0)
635	raysFront++;
636
637	// (*ri).mRay->mSide = side;
638	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
639
640	}
641	}
642
643	float ratio = GetCostRatio(
644	leaf,
645	axis,
646	position,
647	raysBack,
648	raysFront,
649	pvsBack,
650	pvsFront);
651
652	// cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
653	// float oldCost = leaf->rays.size();
654
655	// cout<<"ratio="<<ratio<<endl;
656
657	return ratio;
658	}
659
660	float
661	VssTree::BestCostRatio(
662	VssTreeLeaf *leaf,
663	int &axis,
664	float &position,
665	int &raysBack,
666	int &raysFront,
667	int &pvsBack,
668	int &pvsFront
669	)
670	{
671	int nRaysBack[6], nRaysFront[6];
672	int nPvsBack[6], nPvsFront[6];
673	float nPosition[6];
674	float nCostRatio[6];
675	int bestAxis = -1;
676
677	AxisAlignedBox3 sBox = GetBBox(leaf);
678	AxisAlignedBox3 dBox = GetDirBBox(leaf);
679	// int sAxis = box.Size().DrivingAxis();
680	int sAxis = sBox.Size().DrivingAxis();
681	int dAxis = dBox.Size().DrivingAxis() + 3;
682
683
684	float dirSplitBoxSize = 0.01f;
685	bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
686
687
688	for (axis = 0; axis < 5; axis++)
689	if (mInterleaveDirSplits \|\|
690	(axis < 3 && leaf->depth < mDirSplitDepth) \|\|
691	(axis >= 3 && leaf->depth >= mDirSplitDepth)
692	) {
693	if (!mSplitUseOnlyDrivingAxis \|\| axis == sAxis \|\| axis == dAxis) {
694
695
696	if (splitType == ESplitRegular) {
697	if (axis < 3)
698	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
699	else
700	nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
701
702	nCostRatio[axis] = EvalCostRatio(leaf,
703	axis,
704	nPosition[axis],
705	nRaysBack[axis],
706	nRaysFront[axis],
707	nPvsBack[axis],
708	nPvsFront[axis]
709	);
710	} else
711	if (splitType == ESplitHeuristic) {
712	nCostRatio[axis] = EvalCostRatioHeuristic(
713	leaf,
714	axis,
715	nPosition[axis],
716	nRaysBack[axis],
717	nRaysFront[axis],
718	nPvsBack[axis],
719	nPvsFront[axis]);
720	} else
721	if (splitType == ESplitHybrid) {
722	if (leaf->depth > 7)
723	nCostRatio[axis] = EvalCostRatioHeuristic(
724	leaf,
725	axis,
726	nPosition[axis],
727	nRaysBack[axis],
728	nRaysFront[axis],
729	nPvsBack[axis],
730	nPvsFront[axis]);
731	else {
732	if (axis < 3)
733	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
734	else
735	nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
736
737	nCostRatio[axis] = EvalCostRatio(leaf,
738	axis,
739	nPosition[axis],
740	nRaysBack[axis],
741	nRaysFront[axis],
742	nPvsBack[axis],
743	nPvsFront[axis]
744	);
745	}
746	} else {
747	cerr<<"VssTree: Unknown split heuristics\n";
748	exit(1);
749	}
750
751
752	if ( bestAxis == -1)
753	bestAxis = axis;
754	else
755	if ( nCostRatio[axis] < nCostRatio[bestAxis] )
756	bestAxis = axis;
757	}
758	}
759
760	axis = bestAxis;
761	position = nPosition[bestAxis];
762
763	raysBack = nRaysBack[bestAxis];
764	raysFront = nRaysFront[bestAxis];
765
766	pvsBack = nPvsBack[bestAxis];
767	pvsFront = nPvsFront[bestAxis];
768
769	return nCostRatio[bestAxis];
770	}
771
772
773	float
774	VssTree::EvalCostRatioHeuristic(
775	VssTreeLeaf *leaf,
776	const int axis,
777	float &bestPosition,
778	int &raysBack,
779	int &raysFront,
780	int &pvsBack,
781	int &pvsFront
782	)
783	{
784	AxisAlignedBox3 box;
785	float minBox, maxBox;
786
787	if (axis < 3) {
788	box = GetBBox(leaf);
789	minBox = box.Min(axis);
790	maxBox = box.Max(axis);
791	} else {
792	box = GetDirBBox(leaf);
793	minBox = box.Min(axis-3);
794	maxBox = box.Max(axis-3);
795	}
796
797	SortSubdivisionCandidates(leaf, axis);
798
799	// go through the lists, count the number of objects left and right
800	// and evaluate the following cost funcion:
801	// C = ct_div_ci + (qlrl + qrrr)/queries
802
803	int rl=0, rr = (int)leaf->rays.size();
804	int pl=0, pr = (int)leaf->GetPvsSize();
805	float sizeBox = maxBox - minBox;
806
807	float minBand = minBox + 0.1f*(maxBox - minBox);
808	float maxBand = minBox + 0.9f*(maxBox - minBox);
809
810	float minRatio = 1e20f;
811
812	Intersectable::NewMail();
813	// set all object as belonging to the fron pvs
814	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
815	ri != leaf->rays.end();
816	ri++)
817	if ((*ri).mRay->IsActive()) {
818	Intersectable object = (ri).mRay->mTerminationObject;
819	if (object)
820	if (!object->Mailed()) {
821	object->Mail();
822	object->mCounter = 1;
823	} else
824	object->mCounter++;
825	}
826
827	Intersectable::NewMail();
828
829	for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
830	ci < splitCandidates->end();
831	ci++) {
832	VssRay *ray;
833	switch ((*ci).type) {
834	case SortableEntry::ERayMin: {
835	rl++;
836	ray = (VssRay ) (ci).data;
837	Intersectable *object = ray->mTerminationObject;
838	if (object && !object->Mailed()) {
839	object->Mail();
840	pl++;
841	}
842	break;
843	}
844	case SortableEntry::ERayMax: {
845	rr--;
846	ray = (VssRay ) (ci).data;
847	Intersectable *object = ray->mTerminationObject;
848	if (object) {
849	if (--object->mCounter == 0)
850	pr--;
851	}
852	break;
853	}
854	}
855
856	float position = (*ci).value;
857
858	if (position > minBand && position < maxBand) {
859
860	float ratio = GetCostRatio(
861	leaf,
862	axis,
863	position,
864	rl,
865	rr,
866	pl,
867	pr);
868
869
870	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
871	// cout<<"cost= "<<sum<<endl;
872
873	if (ratio < minRatio) {
874	minRatio = ratio;
875	bestPosition = position;
876
877	raysBack = rl;
878	raysFront = rr;
879
880	pvsBack = pl;
881	pvsFront = pr;
882
883	}
884	}
885	}
886
887
888	// cout<<"===================="<<endl;
889	// cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
890	// <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
891	return minRatio;
892	}
893
894	void
895	VssTree::SortSubdivisionCandidates(
896	VssTreeLeaf *node,
897	const int axis
898	)
899	{
900
901	splitCandidates->clear();
902
903	int requestedSize = 2*((int)node->rays.size());
904	// creates a sorted split candidates array
905	if (splitCandidates->capacity() > 500000 &&
906	requestedSize < (int)(splitCandidates->capacity()/10) ) {
907
908	delete splitCandidates;
909	splitCandidates = new vector<SortableEntry>;
910	}
911
912	splitCandidates->reserve(requestedSize);
913
914	// insert all queries
915	for(VssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
916	ri < node->rays.end();
917	ri++) {
918	if ((*ri).mRay->IsActive()) {
919	if (axis < 3) {
920	bool positive = (*ri).mRay->HasPosDir(axis);
921	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
922	SortableEntry::ERayMax,
923	(*ri).ExtrapOrigin(axis),
924	(void )(ri).mRay)
925	);
926
927	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
928	SortableEntry::ERayMin,
929	(*ri).ExtrapTermination(axis),
930	(void )(ri).mRay)
931	);
932	} else {
933	float pos = (*ri).mRay->GetDirParametrization(axis-3);
934	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
935	pos - Limits::Small,
936	(void )(ri).mRay)
937	);
938
939	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMax,
940	pos + Limits::Small,
941	(void )(ri).mRay)
942	);
943	}
944	}
945	}
946
947	stable_sort(splitCandidates->begin(), splitCandidates->end());
948	}
949
950
951	void
952	VssTree::EvaluateLeafStats(const TraversalData &data)
953	{
954
955	// the node became a leaf -> evaluate stats for leafs
956	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
957
958	if (data.depth >= termMaxDepth)
959	stat.maxDepthNodes++;
960
961	// if ( (int)(leaf->rays.size()) < termMinCost)
962	// stat.minCostNodes++;
963	if ( leaf->GetPvsSize() < termMinPvs)
964	stat.minPvsNodes++;
965
966	if ( leaf->GetPvsSize() < termMinRays)
967	stat.minRaysNodes++;
968
969	if (0 && leaf->GetAvgRayContribution() > termMaxRayContribution )
970	stat.maxRayContribNodes++;
971
972	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
973	stat.minSizeNodes++;
974	}
975
976	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
977	stat.maxRayRefs = (int)leaf->rays.size();
978
979	}
980
981	bool
982	VssTree::TerminationCriteriaSatisfied(VssTreeLeaf *leaf)
983	{
984	return ( (leaf->GetPvsSize() < termMinPvs) \|\|
985	(leaf->rays.size() < termMinRays) \|\|
986	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) \|\|
987	(leaf->depth >= termMaxDepth) \|\|
988	(SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
989	// \|\|
990	// (mUseRss && leaf->mPassingRays == leaf->rays.size())
991	);
992	}
993
994
995	VssTreeNode *
996	VssTree::SubdivideNode(
997	VssTreeLeaf *leaf,
998	const AxisAlignedBox3 &box,
999	AxisAlignedBox3 &backBBox,
1000	AxisAlignedBox3 &frontBBox
1001	)
1002	{
1003
1004	if (TerminationCriteriaSatisfied(leaf)) {
1005	#if 0
1006	if (leaf->depth >= termMaxDepth) {
1007	cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
1008	cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
1009	}
1010	#endif
1011
1012	return leaf;
1013	}
1014
1015	float position;
1016
1017	// first count ray sides
1018	int raysBack;
1019	int raysFront;
1020	int pvsBack;
1021	int pvsFront;
1022
1023	// select subdivision axis
1024	int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
1025	// cout<<axis<<" ";
1026
1027	// cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
1028	// pvsFront<<endl;
1029
1030	if (axis == -1) {
1031	return leaf;
1032	}
1033
1034	stat.nodes+=2;
1035	stat.splits[axis]++;
1036
1037	// add the new nodes to the tree
1038	VssTreeInterior *node = new VssTreeInterior(leaf->parent);
1039
1040	node->axis = axis;
1041	node->position = position;
1042	node->bbox = box;
1043	node->dirBBox = GetDirBBox(leaf);
1044
1045	backBBox = box;
1046	frontBBox = box;
1047
1048	VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
1049	VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);
1050
1051	// replace a link from node's parent
1052	if ( leaf->parent )
1053	leaf->parent->ReplaceChildLink(leaf, node);
1054	// and setup child links
1055	node->SetupChildLinks(back, front);
1056
1057	if (axis <= VssTreeNode::SPLIT_Z) {
1058	backBBox.SetMax(axis, position);
1059	frontBBox.SetMin(axis, position);
1060
1061	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1062	ri != leaf->rays.end();
1063	ri++) {
1064	if ((*ri).mRay->IsActive()) {
1065
1066	// first unref ray from the former leaf
1067	(*ri).mRay->Unref();
1068
1069	// Debug << "computed t: " << (*ri).mRay->mT << endl;
1070	// determine the side of this ray with respect to the plane
1071	float t;
1072	int side = node->ComputeRayIntersection(*ri, t);
1073
1074	if (side == 0) {
1075	if ((*ri).mRay->HasPosDir(axis)) {
1076	back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1077	(*ri).mMinT,
1078	t)
1079	);
1080	front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1081	t,
1082	(*ri).mMaxT));
1083	} else {
1084	back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1085	t,
1086	(*ri).mMaxT));
1087	front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1088	(*ri).mMinT,
1089	t));
1090	}
1091	} else
1092	if (side == 1)
1093	front->AddRay(*ri);
1094	else
1095	back->AddRay(*ri);
1096	} else
1097	(*ri).mRay->Unref();
1098	}
1099	} else {
1100	// rays front/back
1101
1102
1103	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1104	ri != leaf->rays.end();
1105	ri++) {
1106	if ((*ri).mRay->IsActive()) {
1107	// first unref ray from the former leaf
1108	(*ri).mRay->Unref();
1109
1110	int side;
1111	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1112	side = 1;
1113	else
1114	side = -1;
1115
1116	if (side == 1)
1117	front->AddRay(*ri);
1118	else
1119	back->AddRay(*ri);
1120
1121	} else
1122	(*ri).mRay->Unref();
1123	}
1124	}
1125
1126	front->SetPvsSize(pvsFront);
1127	back->SetPvsSize(pvsBack);
1128	// compute entropy as well
1129	front->ComputeEntropyImportance();
1130	back->ComputeEntropyImportance();
1131
1132	// update stats
1133	stat.rayRefs -= (int)leaf->rays.size();
1134	stat.rayRefs += raysBack + raysFront;
1135
1136
1137	delete leaf;
1138	return node;
1139	}
1140
1141
1142
1143
1144
1145
1146	int
1147	VssTree::ReleaseMemory(const int time)
1148	{
1149	stack<VssTreeNode *> tstack;
1150
1151	// find a node in the tree which subtree will be collapsed
1152	int maxAccessTime = time - accessTimeThreshold;
1153	int released;
1154
1155	tstack.push(root);
1156
1157	while (!tstack.empty()) {
1158	VssTreeNode *node = tstack.top();
1159	tstack.pop();
1160
1161
1162	if (!node->IsLeaf()) {
1163	VssTreeInterior in = (VssTreeInterior )node;
1164	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1165	if (in->depth >= minCollapseDepth &&
1166	in->lastAccessTime <= maxAccessTime) {
1167	released = CollapseSubtree(node, time);
1168	break;
1169	}
1170
1171	if (in->back->GetAccessTime() <
1172	in->front->GetAccessTime()) {
1173	tstack.push(in->front);
1174	tstack.push(in->back);
1175	} else {
1176	tstack.push(in->back);
1177	tstack.push(in->front);
1178	}
1179	}
1180	}
1181
1182	while (tstack.empty()) {
1183	// could find node to collaps...
1184	// cout<<"Could not find a node to release "<<endl;
1185	break;
1186	}
1187
1188	return released;
1189	}
1190
1191
1192
1193
1194	VssTreeNode *
1195	VssTree::SubdivideLeaf(
1196	VssTreeLeaf *leaf
1197	)
1198	{
1199	VssTreeNode *node = leaf;
1200
1201	AxisAlignedBox3 leafBBox = GetBBox(leaf);
1202
1203	static int pass = 0;
1204	pass ++;
1205
1206	// check if we should perform a dynamic subdivision of the leaf
1207	if (!TerminationCriteriaSatisfied(leaf)) {
1208
1209	// memory check and realese...
1210	if (GetMemUsage() > maxTotalMemory) {
1211	ReleaseMemory( pass );
1212	}
1213
1214	AxisAlignedBox3 backBBox, frontBBox;
1215
1216	// subdivide the node
1217	node =
1218	SubdivideNode(leaf,
1219	leafBBox,
1220	backBBox,
1221	frontBBox
1222	);
1223	}
1224
1225	return node;
1226	}
1227
1228
1229
1230	void
1231	VssTree::UpdateRays(VssRayContainer &remove,
1232	VssRayContainer &add
1233	)
1234	{
1235	VssTreeLeaf::NewMail();
1236
1237	// schedule rays for removal
1238	for(VssRayContainer::const_iterator ri = remove.begin();
1239	ri != remove.end();
1240	ri++) {
1241	(*ri)->ScheduleForRemoval();
1242	}
1243
1244	int inactive=0;
1245
1246	for(VssRayContainer::const_iterator ri = remove.begin();
1247	ri != remove.end();
1248	ri++) {
1249	if ((*ri)->ScheduledForRemoval())
1250	// RemoveRay(*ri, NULL, false);
1251	// !!! BUG - with true it does not work correctly - aggreated delete
1252	RemoveRay(*ri, NULL, true);
1253	else
1254	inactive++;
1255	}
1256
1257
1258	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1259
1260	for(VssRayContainer::const_iterator ri = add.begin();
1261	ri != add.end();
1262	ri++) {
1263	VssTreeNode::RayInfo info(*ri);
1264	if (ClipRay(info, bbox))
1265	AddRay(info);
1266	}
1267	}
1268
1269
1270	void
1271	VssTree::RemoveRay(VssRay *ray,
1272	vector<VssTreeLeaf > affectedLeaves,
1273	const bool removeAllScheduledRays
1274	)
1275	{
1276
1277	stack<RayTraversalData> tstack;
1278
1279	tstack.push(RayTraversalData(root, VssTreeNode::RayInfo(ray)));
1280
1281	RayTraversalData data;
1282
1283	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1284
1285	while (!tstack.empty()) {
1286	data = tstack.top();
1287	tstack.pop();
1288
1289	if (!data.node->IsLeaf()) {
1290	// split the set of rays in two groups intersecting the
1291	// two subtrees
1292
1293	TraverseInternalNode(data, tstack);
1294
1295	} else {
1296	// remove the ray from the leaf
1297	// find the ray in the leaf and swap it with the last ray...
1298	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
1299
1300	if (!leaf->Mailed()) {
1301	leaf->Mail();
1302	if (affectedLeaves)
1303	affectedLeaves->push_back(leaf);
1304
1305	if (removeAllScheduledRays) {
1306	int tail = (int)leaf->rays.size()-1;
1307
1308	for (int i=0; i < (int)(leaf->rays.size()); i++) {
1309	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1310	// find a ray to replace it with
1311	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1312	stat.removedRayRefs++;
1313	leaf->rays[tail].mRay->Unref();
1314	leaf->rays.pop_back();
1315	tail--;
1316	}
1317
1318	if (tail < i)
1319	break;
1320
1321	stat.removedRayRefs++;
1322	leaf->rays[i].mRay->Unref();
1323	leaf->rays[i] = leaf->rays[tail];
1324	leaf->rays.pop_back();
1325	tail--;
1326	}
1327	}
1328	}
1329	}
1330
1331	if (!removeAllScheduledRays)
1332	for (int i=0; i < (int)leaf->rays.size(); i++) {
1333	if (leaf->rays[i].mRay == ray) {
1334	stat.removedRayRefs++;
1335	ray->Unref();
1336	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1337	leaf->rays.pop_back();
1338	// check this ray again
1339	break;
1340	}
1341	}
1342
1343	}
1344	}
1345
1346	if (ray->RefCount() != 0) {
1347	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1348	exit(1);
1349	}
1350
1351	}
1352
1353
1354	void
1355	VssTree::AddRay(VssTreeNode::RayInfo &info)
1356	{
1357
1358	stack<RayTraversalData> tstack;
1359
1360	tstack.push(RayTraversalData(root, info));
1361
1362	RayTraversalData data;
1363
1364	while (!tstack.empty()) {
1365	data = tstack.top();
1366	tstack.pop();
1367
1368	if (!data.node->IsLeaf()) {
1369	TraverseInternalNode(data, tstack);
1370	} else {
1371	// remove the ray from the leaf
1372	// find the ray in the leaf and swap it with the last ray...
1373	VssTreeLeaf leaf = (VssTreeLeaf )data.node;
1374	leaf->AddRay(data.rayData);
1375	stat.addedRayRefs++;
1376	}
1377	}
1378	}
1379
1380	void
1381	VssTree::TraverseInternalNode(
1382	RayTraversalData &data,
1383	stack<RayTraversalData> &tstack)
1384	{
1385	VssTreeInterior in = (VssTreeInterior ) data.node;
1386
1387	if (in->axis <= VssTreeNode::SPLIT_Z) {
1388	float t;
1389	// determine the side of this ray with respect to the plane
1390	int side = in->ComputeRayIntersection(data.rayData,
1391	t);
1392
1393
1394	if (side == 0) {
1395	if (data.rayData.mRay->HasPosDir(in->axis)) {
1396	tstack.push(RayTraversalData(in->back,
1397	VssTreeNode::RayInfo(data.rayData.mRay,
1398	data.rayData.mMinT,
1399	t))
1400	);
1401
1402	tstack.push(RayTraversalData(in->front,
1403	VssTreeNode::RayInfo(data.rayData.mRay,
1404	t,
1405	data.rayData.mMaxT
1406	))
1407	);
1408
1409	} else {
1410	tstack.push(RayTraversalData(in->back,
1411	VssTreeNode::RayInfo(data.rayData.mRay,
1412	t,
1413	data.rayData.mMaxT
1414	))
1415	);
1416
1417	tstack.push(RayTraversalData(in->front,
1418	VssTreeNode::RayInfo(data.rayData.mRay,
1419	data.rayData.mMinT,
1420	t))
1421	);
1422
1423
1424	}
1425	} else
1426	if (side == 1)
1427	tstack.push(RayTraversalData(in->front, data.rayData));
1428	else
1429	tstack.push(RayTraversalData(in->back, data.rayData));
1430	}
1431	else {
1432	// directional split
1433	if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position)
1434	tstack.push(RayTraversalData(in->front, data.rayData));
1435	else
1436	tstack.push(RayTraversalData(in->back, data.rayData));
1437	}
1438	}
1439
1440
1441	int
1442	VssTree::CollapseSubtree(VssTreeNode *sroot, const int time)
1443	{
1444	// first count all rays in the subtree
1445	// use mail 1 for this purpose
1446	stack<VssTreeNode *> tstack;
1447	int rayCount = 0;
1448	int totalRayCount = 0;
1449	int collapsedNodes = 0;
1450
1451	#if DEBUG_COLLAPSE
1452	cout<<"Collapsing subtree"<<endl;
1453	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1454	cout<<"depth="<<(int)sroot->depth<<endl;
1455	#endif
1456
1457	// tstat.collapsedSubtrees++;
1458	// tstat.collapseDepths += (int)sroot->depth;
1459	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1460
1461	tstack.push(sroot);
1462	VssRay::NewMail();
1463
1464	while (!tstack.empty()) {
1465	collapsedNodes++;
1466	VssTreeNode *node = tstack.top();
1467	tstack.pop();
1468
1469	if (node->IsLeaf()) {
1470	VssTreeLeaf leaf = (VssTreeLeaf ) node;
1471	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1472	ri != leaf->rays.end();
1473	ri++) {
1474
1475	totalRayCount++;
1476	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
1477	(*ri).mRay->Mail();
1478	rayCount++;
1479	}
1480	}
1481	} else {
1482	tstack.push(((VssTreeInterior *)node)->back);
1483	tstack.push(((VssTreeInterior *)node)->front);
1484	}
1485	}
1486
1487	VssRay::NewMail();
1488
1489	// create a new node that will hold the rays
1490	VssTreeLeaf *newLeaf = new VssTreeLeaf( sroot->parent, rayCount );
1491	if ( newLeaf->parent )
1492	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1493
1494
1495	tstack.push( sroot );
1496
1497	while (!tstack.empty()) {
1498
1499	VssTreeNode *node = tstack.top();
1500	tstack.pop();
1501
1502	if (node->IsLeaf()) {
1503	VssTreeLeaf leaf = (VssTreeLeaf ) node;
1504
1505	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1506	ri != leaf->rays.end();
1507	ri++) {
1508
1509	// unref this ray from the old node
1510
1511	if ((*ri).mRay->IsActive()) {
1512	(*ri).mRay->Unref();
1513	if (!(*ri).mRay->Mailed()) {
1514	(*ri).mRay->Mail();
1515	newLeaf->AddRay(*ri);
1516	}
1517	} else
1518	(*ri).mRay->Unref();
1519
1520	}
1521	} else {
1522	tstack.push(((VssTreeInterior *)node)->back);
1523	tstack.push(((VssTreeInterior *)node)->front);
1524	}
1525	}
1526
1527	// delete the node and all its children
1528	delete sroot;
1529
1530	// for(VssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1531	// ri != newLeaf->rays.end();
1532	// ri++)
1533	// (*ri).ray->UnMail(2);
1534
1535
1536	#if DEBUG_COLLAPSE
1537	cout<<"Total memory before="<<GetMemUsage()<<endl;
1538	#endif
1539
1540	stat.nodes -= collapsedNodes - 1;
1541	stat.rayRefs -= totalRayCount - rayCount;
1542
1543	#if DEBUG_COLLAPSE
1544	cout<<"collapsed nodes"<<collapsedNodes<<endl;
1545	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1546	cout<<"Total memory after="<<GetMemUsage()<<endl;
1547	cout<<"================================"<<endl;
1548	#endif
1549
1550	// tstat.collapsedNodes += collapsedNodes;
1551	// tstat.collapsedRays += totalRayCount - rayCount;
1552
1553	return totalRayCount - rayCount;
1554	}
1555
1556
1557	int
1558	VssTree::GetPvsSize(const AxisAlignedBox3 &box) const
1559	{
1560	stack<VssTreeNode *> tstack;
1561	tstack.push(root);
1562
1563	Intersectable::NewMail();
1564	int pvsSize = 0;
1565
1566	while (!tstack.empty()) {
1567	VssTreeNode *node = tstack.top();
1568	tstack.pop();
1569
1570
1571	if (node->IsLeaf()) {
1572	VssTreeLeaf leaf = (VssTreeLeaf )node;
1573	for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1574	ri != leaf->rays.end();
1575	ri++)
1576	if ((*ri).mRay->IsActive()) {
1577	Intersectable *object;
1578	#if BIDIRECTIONAL_RAY
1579	object = (*ri).mRay->mOriginObject;
1580	if (object && !object->Mailed()) {
1581	pvsSize++;
1582	object->Mail();
1583	}
1584	#endif
1585	object = (*ri).mRay->mTerminationObject;
1586	if (object && !object->Mailed()) {
1587	pvsSize++;
1588	object->Mail();
1589	}
1590	}
1591	} else {
1592	VssTreeInterior in = (VssTreeInterior )node;
1593	if (in->axis < 3) {
1594	if (box.Max(in->axis) >= in->position )
1595	tstack.push(in->front);
1596
1597	if (box.Min(in->axis) <= in->position )
1598	tstack.push(in->back);
1599	} else {
1600	// both nodes for directional splits
1601	tstack.push(in->front);
1602	tstack.push(in->back);
1603	}
1604	}
1605	}
1606	return pvsSize;
1607	}
1608
1609	void
1610	VssTree::GetRayContributionStatistics(
1611	float &minRayContribution,
1612	float &maxRayContribution,
1613	float &avgRayContribution
1614	)
1615	{
1616	stack<VssTreeNode *> tstack;
1617	tstack.push(root);
1618
1619	minRayContribution = 1.0f;
1620	maxRayContribution = 0.0f;
1621	float sumRayContribution = 0.0f;
1622	int leaves = 0;
1623
1624	while (!tstack.empty()) {
1625	VssTreeNode *node = tstack.top();
1626	tstack.pop();
1627
1628	if (node->IsLeaf()) {
1629	leaves++;
1630	VssTreeLeaf leaf = (VssTreeLeaf )node;
1631	float c = leaf->GetImportance();
1632	if (c > maxRayContribution)
1633	maxRayContribution = c;
1634	if (c < minRayContribution)
1635	minRayContribution = c;
1636	sumRayContribution += c;
1637
1638	} else {
1639	VssTreeInterior in = (VssTreeInterior )node;
1640	// both nodes for directional splits
1641	tstack.push(in->front);
1642	tstack.push(in->back);
1643	}
1644	}
1645
1646	cout<<"sum="<<sumRayContribution<<endl;
1647	cout<<"leaves="<<leaves<<endl;
1648	avgRayContribution = sumRayContribution/(float)leaves;
1649	}
1650
1651
1652	int
1653	VssTree::GenerateRays(const float ratioPerLeaf,
1654	SimpleRayContainer &rays)
1655	{
1656	stack<VssTreeNode *> tstack;
1657	tstack.push(root);
1658
1659	while (!tstack.empty()) {
1660	VssTreeNode *node = tstack.top();
1661	tstack.pop();
1662
1663	if (node->IsLeaf()) {
1664	VssTreeLeaf leaf = (VssTreeLeaf )node;
1665	float c = leaf->GetImportance();
1666	int num = int(c*ratioPerLeaf + 0.5f);
1667	// cout<<num<<" ";
1668
1669	for (int i=0; i < num; i++) {
1670	Vector3 origin = GetBBox(leaf).GetRandomPoint();
1671	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1672	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1673	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1674	rays.push_back(SimpleRay(origin, direction));
1675	}
1676
1677	} else {
1678	VssTreeInterior in = (VssTreeInterior )node;
1679	// both nodes for directional splits
1680	tstack.push(in->front);
1681	tstack.push(in->back);
1682	}
1683	}
1684
1685	return (int)rays.size();
1686	}
1687
1688	void
1689	VssTree::CollectLeaves(vector<VssTreeLeaf *> &leaves)
1690	{
1691	stack<VssTreeNode *> tstack;
1692	tstack.push(root);
1693
1694	while (!tstack.empty()) {
1695	VssTreeNode *node = tstack.top();
1696	tstack.pop();
1697
1698	if (node->IsLeaf()) {
1699	VssTreeLeaf leaf = (VssTreeLeaf )node;
1700	leaves.push_back(leaf);
1701	} else {
1702	VssTreeInterior in = (VssTreeInterior )node;
1703	// both nodes for directional splits
1704	tstack.push(in->front);
1705	tstack.push(in->back);
1706	}
1707	}
1708	}
1709
1710	bool
1711	VssTree::ValidLeaf(VssTreeLeaf *leaf) const
1712	{
1713	return leaf->rays.size() > termMinRays/4;
1714	}
1715
1716
1717	void
1718	VssTree::GenerateLeafRays(VssTreeLeaf *leaf,
1719	const int numberOfRays,
1720	SimpleRayContainer &rays)
1721	{
1722	int nrays = (int)leaf->rays.size();
1723	for (int i = 0; i < numberOfRays; ++ i) {
1724	// pickup 3 random rays
1725	int r1 = (int)RandomValue(0, (float)(nrays-1));
1726	int r2 = (int)RandomValue(0, (float)(nrays-1));
1727	int r3 = (int)RandomValue(0, (float)(nrays-1));
1728
1729	Vector3 o1 = leaf->rays[r1].Extrap(RandomValue(leaf->rays[r1].GetMinT(),
1730	leaf->rays[r1].GetMaxT()));
1731
1732	Vector3 o2 = leaf->rays[r2].Extrap(RandomValue(leaf->rays[r2].GetMinT(),
1733	leaf->rays[r2].GetMaxT()));
1734
1735	Vector3 o3 = leaf->rays[r3].Extrap(RandomValue(leaf->rays[r3].GetMinT(),
1736	leaf->rays[r3].GetMaxT()));
1737
1738	const float overlap = 0.1f;
1739
1740	Vector3 origin, direction;
1741	bool useExtendedConvexCombination = true;
1742	if (useExtendedConvexCombination) {
1743	float w1, w2, w3;
1744	GenerateExtendedConvexCombinationWeights(w1, w2, w3, overlap);
1745	origin = w1o1 + w2o2 + w3*o3;
1746	direction =
1747	w1*leaf->rays[r1].mRay->GetDir() +
1748	w2*leaf->rays[r2].mRay->GetDir() +
1749	w3*leaf->rays[r3].mRay->GetDir();
1750	} else {
1751	origin = GetBBox(leaf).GetRandomPoint();
1752	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1753	direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1754	}
1755	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1756	rays.push_back(SimpleRay(origin, direction));
1757	}
1758	}
1759
1760	int
1761	VssTree::GenerateRays(const int numberOfRays,
1762	const int numberOfLeaves,
1763	SimpleRayContainer &rays)
1764	{
1765
1766	vector<VssTreeLeaf *> leaves;
1767
1768	CollectLeaves(leaves);
1769
1770	sort(leaves.begin(),
1771	leaves.end(),
1772	GreaterContribution);
1773
1774
1775	float sumContrib = 0.0;
1776	int i;
1777	int k = 0;
1778	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
1779	if (ValidLeaf(leaves[i])) {
1780	float c = leaves[i]->GetImportance();
1781	sumContrib += c;
1782	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
1783	k++;
1784	}
1785
1786	float avgContrib = sumContrib/numberOfLeaves;
1787	float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
1788	k = 0;
1789	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
1790	if (ValidLeaf(leaves[i])) {
1791	k++;
1792	VssTreeLeaf *leaf = leaves[i];
1793	float c = leaf->GetImportance();
1794	int num = (int)(c*ratioPerLeaf + 0.5f);
1795	GenerateLeafRays(leaf, num, rays);
1796	}
1797
1798	return (int)rays.size();
1799	}
1800
1801
1802	float
1803	VssTree::GetAvgPvsSize()
1804	{
1805	stack<VssTreeNode *> tstack;
1806	tstack.push(root);
1807
1808	int sumPvs = 0;
1809	int leaves = 0;
1810	while (!tstack.empty()) {
1811	VssTreeNode *node = tstack.top();
1812	tstack.pop();
1813
1814	if (node->IsLeaf()) {
1815	VssTreeLeaf leaf = (VssTreeLeaf )node;
1816	// update pvs size
1817	leaf->UpdatePvsSize();
1818	sumPvs += leaf->GetPvsSize();
1819	leaves++;
1820	} else {
1821	VssTreeInterior in = (VssTreeInterior )node;
1822	// both nodes for directional splits
1823	tstack.push(in->front);
1824	tstack.push(in->back);
1825	}
1826	}
1827
1828
1829	return sumPvs/(float)leaves;
1830	}
1831
1832
1833	float
1834	VssTreeLeaf::GetImportance() const
1835	{
1836
1837	if (1) {
1838	// return GetAvgRayContribution();
1839	return (float)GetPvsSize();
1840	} else {
1841	// return GetAvgRayContribution()*mEntropyImportance;
1842	//return GetAvgRayContribution();
1843	return mEntropyImportance;
1844	}
1845	}
1846
1847
1848	float
1849	VssTreeLeaf::ComputePvsEntropy()
1850	{
1851	int samples = 0;
1852	Intersectable::NewMail();
1853	// set all object as belonging to the fron pvs
1854	for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
1855	ri != rays.end();
1856	ri++)
1857	if ((*ri).mRay->IsActive()) {
1858	Intersectable object = (ri).mRay->mTerminationObject;
1859	if (object) {
1860	if (!object->Mailed()) {
1861	object->Mail();
1862	object->mCounter = 1;
1863	} else
1864	object->mCounter++;
1865	samples++;
1866	}
1867	}
1868
1869	float entropy = 0.0f;
1870
1871	if (samples > 1) {
1872	Intersectable::NewMail();
1873	for(RayInfoContainer::const_iterator ri = rays.begin();
1874	ri != rays.end();
1875	ri++)
1876	if ((*ri).mRay->IsActive()) {
1877	Intersectable object = (ri).mRay->mTerminationObject;
1878	if (object) {
1879	if (!object->Mailed()) {
1880	object->Mail();
1881	float p = object->mCounter/(float)samples;
1882	entropy -= p*log(p);
1883	}
1884	}
1885	}
1886	entropy = entropy/log((float)samples);
1887	}
1888	else
1889	entropy = 1.0f;
1890
1891	return entropy;
1892	}
1893
1894	float
1895	VssTreeLeaf::ComputeRayLengthEntropy()
1896	{
1897	// get sum of all ray lengths
1898	// consider only passing rays or originating rays
1899	float sum = 0.0f;
1900	int samples = 0;
1901
1902	for(RayInfoContainer::const_iterator ri = rays.begin();
1903	ri != rays.end();
1904	ri++) {
1905	int rayClass = (*ri).GetRayClass();
1906	if (
1907	rayClass == RayInfo::PASSING_RAY
1908	// \|\|
1909	// rayClass == RayInfo::SOURCE_RAY
1910	// rayClass == RayInfo::TERMINATION_RAY
1911	) {
1912	float c = 1.0f - (*ri).GetMaxT();
1913	sum += (ri).mRay->GetSize()c;
1914	samples++;
1915	}
1916	}
1917
1918	float entropy = 0.0f;
1919
1920	if (samples > 1) {
1921	for(RayInfoContainer::const_iterator ri = rays.begin();
1922	ri != rays.end();
1923	ri++) {
1924	int rayClass = (*ri).GetRayClass();
1925	if (
1926	rayClass == RayInfo::PASSING_RAY
1927	// \|\|
1928	// rayClass == RayInfo::SOURCE_RAY
1929	// rayClass == RayInfo::TERMINATION_RAY
1930	) {
1931	float c = 1.0f - (*ri).GetMaxT();
1932	float p = (ri).mRay->GetSize()c/sum;
1933	entropy -= p*log(p);
1934	}
1935	}
1936	entropy = entropy/log((float)samples);
1937	} else
1938	entropy = 1.0f;
1939
1940	return entropy;
1941	}
1942
1943	float
1944	VssTreeLeaf::ComputeRayTerminationEntropy()
1945	{
1946
1947	return 0.0f;
1948	}
1949
1950
1951	void
1952	VssTreeLeaf::ComputeEntropyImportance()
1953	{
1954	// mEntropy = 1.0f - ComputeRayLengthEntropy();
1955	mEntropyImportance = 1.0f - ComputePvsEntropy();
1956
1957	// cout<<"ei="<<mEntropyImportance<<" ";
1958	}
1959
1960
1961	int
1962	VssTree::CollectRays(VssRayContainer &rays,
1963	const int number)
1964	{
1965	VssRayContainer allRays;
1966	CollectRays(allRays);
1967
1968	int desired = min(number, (int)allRays.size());
1969	float prob = desired/(float)allRays.size();
1970	while (rays.size() < desired) {
1971	VssRayContainer::const_iterator it = allRays.begin();
1972	for (; it != allRays.end() && rays.size() < desired; it++) {
1973	if (Random(1.0f) < prob)
1974	rays.push_back(*it);
1975	}
1976	}
1977	return (int)rays.size();
1978	}
1979
1980
1981	int
1982	VssTree::CollectRays(VssRayContainer &rays
1983	)
1984	{
1985	VssRay::NewMail();
1986
1987	stack<VssTreeNode *> tstack;
1988	tstack.push(root);
1989
1990	while (!tstack.empty()) {
1991	VssTreeNode *node = tstack.top();
1992	tstack.pop();
1993
1994	if (node->IsLeaf()) {
1995	VssTreeLeaf leaf = (VssTreeLeaf )node;
1996	// update pvs size
1997	VssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
1998	for (;it != leaf->rays.end(); ++it)
1999	if (!(*it).mRay->Mailed()) {
2000	(*it).mRay->Mail();
2001	rays.push_back((*it).mRay);
2002	}
2003	} else {
2004	VssTreeInterior in = (VssTreeInterior )node;
2005	// both nodes for directional splits
2006	tstack.push(in->front);
2007	tstack.push(in->back);
2008	}
2009	}
2010
2011	return (int)rays.size();
2012	}
2013
2014	}

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

Download in other formats: