Context Navigation

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

Revision 2105, 47.3 KB checked in by bittner, 17 years ago (diff)
simple ray separated

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

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

Download in other formats: