Context Navigation

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

Revision 434, 42.9 KB checked in by bittner, 19 years ago (diff)
vssbsp merge, tab change

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

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

Download in other formats: