Context Navigation

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

Revision 438, 45.9 KB checked in by bittner, 19 years ago (diff)
vss updates

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

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

Download in other formats: