Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/RssTree.cpp @ 1715

Revision 1715, 72.7 KB checked in by bittner, 18 years ago (diff)
new visibility filter support

Line
1	// ================================================================
2	// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
3	// ****************************************************************
4	/**
5	The KD tree based LSDS
6	*/
7	// Initial coding by
8	/**
9	@author Jiri Bittner
10	*/
11
12	// Standard headers
13	#include <stack>
14	#include <queue>
15	#include <algorithm>
16	#include <fstream>
17	#include <string>
18
19	#include "RssTree.h"
20
21	#include "Environment.h"
22	#include "VssRay.h"
23	#include "Intersectable.h"
24	#include "Ray.h"
25	#include "Containers.h"
26	#include "ViewCell.h"
27	#include "Exporter.h"
28	#include "Preprocessor.h"
29	#include "SceneGraph.h"
30
31
32	namespace GtpVisibilityPreprocessor {
33
34
35	#define DEBUG_SPLIT_COST 0
36	#define DEBUG_SPLITS 0
37
38	#define EVAL_VIEWCELLS 0
39
40	// Static variables
41	int
42	RssTreeLeaf::mailID = 0;
43
44	inline void
45	AddObject2Pvs(Intersectable *object,
46	const int side,
47	int &pvsBack,
48	int &pvsFront)
49	{
50
51	if (!object)
52	return;
53
54	if (side <= 0) {
55	if (!object->Mailed() && !object->Mailed(2)) {
56	pvsBack++;
57	if (object->Mailed(1))
58	object->Mail(2);
59	else
60	object->Mail();
61	}
62	}
63
64	if (side >= 0) {
65	if (!object->Mailed(1) && !object->Mailed(2)) {
66	pvsFront++;
67	if (object->Mailed())
68	object->Mail(2);
69	else
70	object->Mail(1);
71	}
72	}
73	}
74
75	inline void
76	AddViewcells2Pvs(const ViewCellContainer &viewcells,
77	const int side,
78	int &viewcellsBack,
79	int &viewcellsFront)
80	{
81	ViewCellContainer::const_iterator it = viewcells.begin();
82
83	for (; it != viewcells.end(); ++it) {
84	ViewCell viewcell = it;
85	if (side <= 0) {
86	if (!viewcell->Mailed() && !viewcell->Mailed(2)) {
87	viewcellsBack++;
88	if (viewcell->Mailed(1))
89	viewcell->Mail(2);
90	else
91	viewcell->Mail();
92	}
93	}
94
95	if (side >= 0) {
96	if (!viewcell->Mailed(1) && !viewcell->Mailed(2)) {
97	viewcellsFront++;
98	if (viewcell->Mailed())
99	viewcell->Mail(2);
100	else
101	viewcell->Mail(1);
102	}
103	}
104	}
105	}
106
107
108	// Constructor
109	RssTree::RssTree()
110	{
111	Environment::GetSingleton()->GetIntValue("RssTree.maxDepth", termMaxDepth);
112	Environment::GetSingleton()->GetIntValue("RssTree.minPvs", termMinPvs);
113	Environment::GetSingleton()->GetIntValue("RssTree.minRays", termMinRays);
114	Environment::GetSingleton()->GetFloatValue("RssTree.maxRayContribution", termMaxRayContribution);
115	Environment::GetSingleton()->GetFloatValue("RssTree.maxCostRatio", termMaxCostRatio);
116
117	Environment::GetSingleton()->GetFloatValue("RssTree.minSize", termMinSize);
118	termMinSize = sqr(termMinSize);
119
120	Environment::GetSingleton()->GetFloatValue("RssTree.refDirBoxMaxSize", refDirBoxMaxSize);
121	refDirBoxMaxSize = sqr(refDirBoxMaxSize);
122
123	Environment::GetSingleton()->GetFloatValue("RssTree.epsilon", epsilon);
124	Environment::GetSingleton()->GetFloatValue("RssTree.ct_div_ci", ct_div_ci);
125
126	Environment::GetSingleton()->GetFloatValue("RssTree.maxTotalMemory", maxTotalMemory);
127	Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
128
129	Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
130
131
132
133	Environment::GetSingleton()->GetIntValue("RssTree.accessTimeThreshold", accessTimeThreshold);
134	//= 1000;
135	Environment::GetSingleton()->GetIntValue("RssTree.minCollapseDepth", minCollapseDepth);
136	// int minCollapseDepth = 4;
137
138	// pRefDirThresh = cos(0.5M_PI - M_PIrefDirAngle/180.0);
139	// cosRefDir = cos(M_PI*refDirAngle/180.0);
140	// sinRefDir = sin(M_PI*refDirAngle/180.0);
141
142
143	// split type
144	char sname[128];
145	Environment::GetSingleton()->GetStringValue("RssTree.splitType", sname);
146	string name(sname);
147
148	if (name.compare("regular") == 0)
149	splitType = ESplitRegular;
150	else
151	if (name.compare("heuristic") == 0)
152	splitType = ESplitHeuristic;
153	else
154	if (name.compare("hybrid") == 0)
155	splitType = ESplitHybrid;
156	else {
157	cerr<<"Invalid RssTree split type "<<name<<endl;
158	exit(1);
159	}
160
161	Environment::GetSingleton()->GetIntValue("RssTree.hybridDepth", mHybridDepth);
162
163	Environment::GetSingleton()->GetBoolValue("RssTree.randomize", randomize);
164	Environment::GetSingleton()->GetBoolValue("RssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
165
166	Environment::GetSingleton()->GetBoolValue("RssTree.interleaveDirSplits", mInterleaveDirSplits);
167	Environment::GetSingleton()->GetIntValue("RssTree.dirSplitDepth", mDirSplitDepth);
168
169	Environment::GetSingleton()->GetBoolValue("RssTree.importanceBasedCost", mImportanceBasedCost);
170
171	Environment::GetSingleton()->GetIntValue("RssTree.maxRays", mMaxRays);
172
173	Environment::GetSingleton()->GetBoolValue("RssTree.perObjectTree", mPerObjectTree);
174
175	// mRoots;
176
177	splitCandidates = new vector<SortableEntry>;
178	}
179
180
181	RssTree::~RssTree()
182	{
183	for (int i=0; i < mRoots.size(); i++)
184	if (mRoots[i])
185	delete mRoots[i];
186	}
187
188
189
190
191	void
192	RssStatistics::Print(ostream &app) const
193	{
194	app << "###### RssTree statistics ######\n";
195
196	app << "#N_RAYS ( Number of rays )\n"
197	<< rays <<endl;
198
199	app << "#N_INITPVS ( Initial PVS size )\n"
200	<< initialPvsSize <<endl;
201
202	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
203
204	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
205
206	app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
207	for (int i=0; i<7; i++)
208	app << splits[i] <<" ";
209	app <<endl;
210
211	app << "#N_RAYREFS ( Number of rayRefs )\n" <<
212	rayRefs << "\n";
213
214	app << "#N_RAYRAYREFS ( Number of rayRefs / ray )\n" <<
215	rayRefs/(double)rays << "\n";
216
217	app << "#N_LEAFRAYREFS ( Number of rayRefs / leaf )\n" <<
218	rayRefs/(double)Leaves() << "\n";
219
220	app << "#N_MAXRAYREFS ( Max number of rayRefs / leaf )\n" <<
221	maxRayRefs << "\n";
222
223
224	// app << setprecision(4);
225
226	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
227	maxDepthNodes*100/(double)Leaves()<<endl;
228
229	app << "#N_PMINPVSLEAVES ( Percentage of leaves with minPvs )\n"<<
230	minPvsNodes*100/(double)Leaves()<<endl;
231
232	app << "#N_PMINRAYSLEAVES ( Percentage of leaves with minRays )\n"<<
233	minRaysNodes*100/(double)Leaves()<<endl;
234
235	app << "#N_PMINSIZELEAVES ( Percentage of leaves with minSize )\n"<<
236	minSizeNodes*100/(double)Leaves()<<endl;
237
238	app << "#N_PMAXRAYCONTRIBLEAVES ( Percentage of leaves with maximal ray contribution )\n"<<
239	maxRayContribNodes*100/(double)Leaves()<<endl;
240
241	app << "#N_PMAXCOSTRATIOLEAVES ( Percentage of leaves with max cost ratio )\n"<<
242	maxCostRatioNodes*100/(double)Leaves()<<endl;
243
244	app << "#N_ADDED_RAYREFS (Number of dynamically added ray references )\n"<<
245	addedRayRefs<<endl;
246
247	app << "#N_REMOVED_RAYREFS (Number of dynamically removed ray references )\n"<<
248	removedRayRefs<<endl;
249
250	// app << setprecision(4);
251
252	app << "#N_CTIME ( Construction time [s] )\n"
253	<< Time() << " \n";
254
255	app << "###### END OF RssTree statistics ######\n";
256
257	}
258
259
260	void
261	RssTree::UpdatePvsSize(RssTreeLeaf *leaf)
262	{
263	if (!leaf->mValidPvs) {
264	Intersectable::NewMail();
265	int pvsSize = 0;
266	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
267	ri != leaf->rays.end();
268	ri++)
269	if ((*ri).mRay->IsActive()) {
270	Intersectable *object;
271	object = (*ri).GetObject();
272	if (object && !object->Mailed()) {
273	pvsSize++;
274	object->Mail();
275	}
276	}
277	leaf->SetPvsSize(pvsSize);
278
279	ComputeImportance(leaf);
280	}
281	}
282
283	bool
284	RssTree::ClipRay(
285	RssTreeNode::RayInfo &rayInfo,
286	const AxisAlignedBox3 &box
287	)
288	{
289	float tmin, tmax;
290	static Ray ray;
291	cerr<<"Clip not reimplented yet...Quiting\n";
292	exit(1);
293	ray.Init(rayInfo.GetOrigin(), rayInfo.GetDir(), Ray::LINE_SEGMENT);
294
295	if (!box.ComputeMinMaxT(ray, &tmin, &tmax) \|\| tmin>=tmax)
296	return false;
297
298	// now check if the ray origin lies inside the box
299	if ( tmax < rayInfo.mRay->GetSize() ) {
300	// this ray does not leave the box
301	rayInfo.mRay->SetupEndPoints(
302	ray.Extrap(tmax),
303	rayInfo.mRay->mTermination
304	);
305	return true;
306	}
307
308	return false;
309	}
310
311
312	void
313	RssTree::Construct(
314	ObjectContainer &objects,
315	VssRayContainer &rays
316	// forced bounding box is only used when computing from-box
317	// visibility
318	// AxisAlignedBox3 *forcedBoundingBox
319	)
320	{
321	cout<<"Constructing rss tree"<<endl<<flush;
322	stat.Start();
323
324	maxMemory = maxStaticMemory;
325
326	// if (root)
327	// delete root;
328
329	if (mPerObjectTree) {
330	// get max id from the rays
331	int i;
332	mRoots.resize(objects.size());
333	for (i = 0; i < objects.size(); i++) {
334	RssTreeLeaf *leaf = new RssTreeLeaf(NULL, 0);
335	// leaf->bbox.Initialize();
336	leaf->dirBBox.Initialize();
337	leaf->dirBBox.SetMin(2, 0.0f);
338	leaf->dirBBox.SetMax(2, 1.0f);
339	mRoots[i] = leaf;
340	}
341
342	// init spatial bounding boxes
343	for (i = 0; i < objects.size(); i++) {
344	GetRoot(objects[i])->bbox = objects[i]->GetBox();
345	}
346	stat.nodes = i;
347	stat.leaves = i;
348	} else {
349	mRoots.resize(1);
350	RssTreeLeaf *leaf = new RssTreeLeaf(NULL, (int)rays.size());
351	leaf->bbox.Initialize();
352	leaf->dirBBox.Initialize();
353	leaf->dirBBox.SetMin(2, 0.0f);
354	leaf->dirBBox.SetMax(2, 1.0f);
355	mRoots[0] = leaf;
356	stat.nodes = 1;
357	stat.leaves = 1;
358	}
359
360	for(VssRayContainer::const_iterator ri = rays.begin();
361	ri != rays.end();
362	ri++) {
363
364	RssTreeNode::RayInfo info(*ri);
365
366	// first construct a leaf that will get subdivide
367	RssTreeLeaf leaf = (RssTreeLeaf ) GetRoot(info.GetSourceObject());
368
369	leaf->AddRay(info);
370
371	// leaf bbox contains bbox of origins only
372
373	// include both origin and terminatin in the global bbox
374	#if USE_ORIGIN
375	leaf->bbox.Include((*ri)->GetOrigin());
376	bbox.Include((*ri)->GetOrigin());
377	#else
378	bbox.Include((*ri)->GetTermination());
379	leaf->bbox.Include((*ri)->GetTermination());
380	#endif
381	Vector3 dVec = Vector3(
382	(*ri)->GetDirParametrization(0),
383	(*ri)->GetDirParametrization(1),
384	0
385	);
386
387	leaf->dirBBox.Include(dVec);
388	dirBBox.Include(dVec);
389	}
390
391	// make the z axis (unused) a unit size
392	// important for volume computation
393
394
395	// if ( forcedBoundingBox )
396	// bbox = *forcedBoundingBox;
397
398	cout<<"Bbox = "<<bbox<<endl;
399	cout<<"Dirr Bbox = "<<dirBBox<<endl;
400
401	stat.rays = rays.size();
402	stat.initialPvsSize = 0;
403	for (int i=0; i < mRoots.size(); i++) {
404	RssTreeLeaf leaf = (RssTreeLeaf )mRoots[i];
405	UpdatePvsSize(leaf);
406	stat.initialPvsSize += leaf->GetPvsSize();
407	mRoots[i] = Subdivide(TraversalData(leaf, GetBBox(leaf), 0));
408	}
409
410	if (splitCandidates) {
411	// force realease of this vector
412	delete splitCandidates;
413	splitCandidates = new vector<SortableEntry>;
414	}
415
416	stat.Stop();
417	stat.Print(cout);
418	cout<<"#Total memory="<<GetMemUsage()<<endl;
419
420	// this rotine also updates importances etc...
421	}
422
423	int
424	RssTree::UpdateSubdivision()
425	{
426	priority_queue<TraversalData> tStack;
427	// stack<TraversalData> tStack;
428
429	// tStack.push(TraversalData(root, bbox, 0));
430	PushRoots(tStack);
431
432	AxisAlignedBox3 backBox;
433	AxisAlignedBox3 frontBox;
434
435	maxMemory = maxTotalMemory;
436	int subdivided = 0;
437	int lastMem = 0;
438	while (!tStack.empty()) {
439
440	float mem = GetMemUsage();
441
442	if ( lastMem/10 != ((int)mem)/10) {
443	cout<<mem<<" MB"<<endl;
444	}
445	lastMem = (int)mem;
446
447	if ( mem > maxMemory ) {
448	// count statistics on unprocessed leafs
449	while (!tStack.empty()) {
450	// EvaluateLeafStats(tStack.top());
451	tStack.pop();
452	}
453	break;
454	}
455
456	TraversalData data = tStack.top();
457	tStack.pop();
458
459	if (data.node->IsLeaf()) {
460	RssTreeNode node = SubdivideNode((RssTreeLeaf ) data.node,
461	data.bbox,
462	backBox,
463	frontBox
464	);
465	if (!node->IsLeaf()) {
466	subdivided++;
467
468
469	RssTreeInterior interior = (RssTreeInterior ) node;
470	// push the children on the stack
471	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
472	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
473	} else {
474	// EvaluateLeafStats(data);
475	}
476	} else {
477	RssTreeInterior interior = (RssTreeInterior ) data.node;
478	tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
479	tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
480	}
481	}
482	return subdivided;
483	}
484
485
486	RssTreeNode *
487	RssTree::Subdivide(const TraversalData &tdata)
488	{
489	RssTreeNode *result = NULL;
490
491	priority_queue<TraversalData> tStack;
492	// stack<TraversalData> tStack;
493
494	tStack.push(tdata);
495
496	AxisAlignedBox3 backBox;
497	AxisAlignedBox3 frontBox;
498
499
500	int lastMem = 0;
501	while (!tStack.empty()) {
502
503	float mem = GetMemUsage();
504
505	if ( lastMem/10 != ((int)mem)/10) {
506	cout<<mem<<" MB"<<endl;
507	}
508	lastMem = (int)mem;
509
510	if ( mem > maxMemory ) {
511	// count statistics on unprocessed leafs
512	while (!tStack.empty()) {
513	EvaluateLeafStats(tStack.top());
514	tStack.pop();
515	}
516	break;
517	}
518
519	TraversalData data = tStack.top();
520	tStack.pop();
521
522	#if DEBUG_SPLITS
523	Debug<<"#Splitting node"<<endl;
524	data.node->Print(Debug);
525	#endif
526	RssTreeNode node = SubdivideNode((RssTreeLeaf ) data.node,
527	data.bbox,
528	backBox,
529	frontBox
530	);
531
532
533	if (result == NULL)
534	result = node;
535
536	if (!node->IsLeaf()) {
537
538	RssTreeInterior interior = (RssTreeInterior ) node;
539	// push the children on the stack
540	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
541	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
542
543	#if DEBUG_SPLITS
544	Debug<<"#New nodes"<<endl;
545	interior->back->Print(Debug);
546	interior->front->Print(Debug);
547	Debug<<"#####################################"<<endl;
548	#endif
549
550	} else {
551	EvaluateLeafStats(data);
552	}
553	}
554
555	return result;
556	}
557
558
559	// returns selected plane for subdivision
560	int
561	RssTree::SelectPlane(
562	RssTreeLeaf *leaf,
563	const AxisAlignedBox3 &box,
564	SplitInfo &info
565	)
566	{
567
568	BestCostRatio(leaf,
569	info);
570
571	#if DEBUG_SPLIT_COST
572	Debug<<"Split Info:"<<endl;
573	Debug<<"axis="<<info.axis<<" ratio="<<info.costRatio<<endl;
574	Debug<<"rays="<<info.rays<<
575	" rays back="<<info.raysBack<<
576	" rays front="<<info.raysFront<<endl;
577	Debug<<"c="<<info.contribution<<
578	" c back="<<info.contributionBack<<
579	" c front="<<info.contributionFront<<endl;
580
581	// Debug<<"viewcells="<<info.viewCells<<
582	// " viewcells back="<<info.viewCellsBack<<
583	// " viewcells back="<<info.viewCellsFront<<endl;
584	#endif
585
586	if (info.costRatio > termMaxCostRatio) {
587	// cout<<"Too big cost ratio "<<costRatio<<endl;
588	stat.maxCostRatioNodes++;
589	return -1;
590	}
591
592	#if 0
593	cout<<
594	"pvs="<<leaf->mPvsSize<<
595	" rays="<<leaf->rays.size()<<
596	" rc="<<leaf->GetAvgRayContribution()<<
597	" axis="<<info.axis<<endl;
598	#endif
599
600	return info.axis;
601	}
602
603
604	void
605	RssTree::GetCostRatio(
606	RssTreeLeaf *leaf,
607	SplitInfo &info
608	)
609	{
610
611	AxisAlignedBox3 box;
612	float minBox, maxBox;
613
614	if (info.axis < 3) {
615	box = GetBBox(leaf);
616	minBox = box.Min(info.axis);
617	maxBox = box.Max(info.axis);
618	} else {
619	box = GetDirBBox(leaf);
620	minBox = box.Min(info.axis-3);
621	maxBox = box.Max(info.axis-3);
622	}
623
624	float sizeBox = maxBox - minBox;
625
626	int pvsSize = leaf->GetPvsSize();
627
628	if (!mImportanceBasedCost) {
629	const int costMethod = 0;
630
631	switch (costMethod) {
632	case 0: {
633	float sum = info.pvsBack(info.position - minBox) + info.pvsFront(maxBox - info.position);
634	float newCost = ct_div_ci + sum/sizeBox;
635	float oldCost = (float)pvsSize;
636	info.costRatio = newCost/oldCost;
637	break;
638	}
639	case 6: {
640	float sum = info.raysBack(info.position - minBox) + info.raysFront(maxBox - info.position);
641	float newCost = ct_div_ci + sum/sizeBox;
642	float oldCost = (float)info.rays;
643	info.costRatio = newCost/oldCost;
644	break;
645	}
646	case 7: {
647	float benefit =
648	(sqr(info.contributionBack) + sqr(info.contributionFront))/
649	(2.0f*sqr(info.contribution));
650	info.costRatio = 1.0f/(benefit + Limits::Small);
651	break;
652	}
653
654
655	case 1: {
656	float newContrib =
657	(info.contributionBack*(info.position - minBox) +
658	+
659	info.contributionFront*(maxBox - info.position))/sizeBox;
660	float oldContrib = info.contribution;
661	info.costRatio = newContrib/oldContrib;
662	// cout<<info.contribution<<" "<<info.contributionBack<<" "<<info.contributionFront<<endl;
663	break;
664	}
665	case 2: {
666	float sum =
667	info.viewCellsBack*(info.position - minBox) +
668	info.viewCellsFront*(maxBox - info.position);
669	float newCost = ct_div_ci + sum/sizeBox;
670	float oldCost = (float)info.viewCells;
671	info.costRatio = newCost/oldCost;
672	break;
673	}
674	case 3: {
675	float newCost = (float)(info.raysBackinfo.pvsBack + info.raysFrontinfo.pvsFront);
676	float oldCost = (float)leaf->rays.size()*pvsSize;
677	info.costRatio = newCost/oldCost;
678	}
679	}
680	} else {
681	const int costMethod = 1;
682	// importance based cost
683	switch (costMethod) {
684	case 0: {
685	break;
686	}
687	case 1: {
688	float newContrib =
689	sqr(info.pvsBack/(info.raysBack + Limits::Small)) +
690	sqr(info.pvsFront/(info.raysFront + Limits::Small));
691	float oldContrib = sqr(leaf->GetAvgRayContribution());
692	info.costRatio = oldContrib/newContrib;
693	break;
694	}
695	case 2: {
696	float newCost = (info.pvsBack + info.pvsFront)*0.5f;
697	float oldCost = (float)pvsSize;
698	info.costRatio = newCost/oldCost;
699	break;
700	}
701	case 3: {
702	float newCost = (float)abs(info.raysBack - info.raysFront);
703	float oldCost = (float)leaf->rays.size();
704	info.costRatio = newCost/oldCost;
705	break;
706	}
707	}
708	}
709	}
710
711
712	void
713	RssTree::EvalCostRatio(
714	RssTreeLeaf *leaf,
715	SplitInfo &info
716	)
717	{
718	info.rays = 0;
719	info.raysBack = 0;
720	info.raysFront = 0;
721	info.pvsFront = 0;
722	info.pvsBack = 0;
723	info.viewCells = 0;
724	info.viewCellsFront = 0;
725	info.viewCellsBack = 0;
726
727
728
729	float sumWeights = Limits::Small;
730	float sumWeightsBack = Limits::Small;
731	float sumWeightsFront = Limits::Small;
732
733	float sumContribution = 0.0f;
734	float sumContributionBack = 0.0f;
735	float sumContributionFront = 0.0f;
736
737	Intersectable::NewMail();
738
739	#if EVAL_VIEWCELLS
740	// count the numebr of viewcells first
741	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
742	ri != leaf->rays.end();
743	ri++)
744	if ((*ri).mRay->IsActive()) {
745
746	ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
747	for (; it != (*ri).mRay->mViewCells.end(); ++it) {
748	if (!(*it)->Mailed()) {
749	(*it)->Mail();
750	info.viewCells++;
751	}
752	}
753	}
754	#endif
755
756	Intersectable::NewMail(3);
757
758	// this is the main ray classification loop!
759	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
760	ri != leaf->rays.end();
761	ri++)
762	if ((*ri).mRay->IsActive()) {
763	int side;
764
765
766	// determine the side of this ray with respect to the plane
767	side = (*ri).ComputeRaySide(info.axis, info.position);
768
769	float weight, contribution;
770
771	GetRayContribution(*ri, weight, contribution);
772	sumWeights += weight;
773	sumContribution += contribution;
774	info.rays++;
775	if (side <= 0) {
776	info.raysBack++;
777	sumWeightsBack += weight;
778	sumContributionBack += contribution;
779	}
780
781	if (side >= 0) {
782	info.raysFront++;
783	sumWeightsFront += weight;
784	sumContributionFront += contribution;
785	}
786
787	AddObject2Pvs((*ri).GetObject(), side, info.pvsBack, info.pvsFront);
788
789	#if EVAL_VIEWCELLS
790	AddViewcells2Pvs((*ri).mRay->mViewCells,
791	side,
792	info.viewCellsBack,
793	info.viewCellsFront);
794	#endif
795	}
796
797	if (sumWeights!=0.0f)
798	info.contribution = sumContribution/sumWeights;
799	else
800	info.contribution = 0.0f;
801
802	if (sumWeightsBack!=0.0f)
803	info.contributionBack = sumContributionBack/sumWeightsBack;
804	else
805	info.contributionBack = 0.0f;
806
807	if (sumWeightsFront!=0.0f)
808	info.contributionFront = sumContributionFront/sumWeightsFront;
809	else
810	info.contributionFront = 0.0f;
811
812	// info.costRatio = 0.1 + Random(0.5f);
813	// return;
814
815	GetCostRatio(
816	leaf,
817	info);
818
819	//cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
820	// float oldCost = leaf->rays.size();
821
822	// cout<<"ratio="<<ratio<<endl;
823	}
824
825	void
826	RssTree::BestCostRatio(
827	RssTreeLeaf *leaf,
828	SplitInfo &info
829	)
830	{
831	SplitInfo nInfo[5];
832	int bestAxis = -1;
833
834	AxisAlignedBox3 sBox = GetBBox(leaf);
835	AxisAlignedBox3 dBox = GetDirBBox(leaf);
836	// int sAxis = box.Size().DrivingAxis();
837	int sAxis = sBox.Size().DrivingAxis();
838	int dAxis = dBox.Size().DrivingAxis() + 3;
839
840	float dirSplitBoxSize = 0.01f;
841	bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
842
843	bool useCyclingAxis = true;
844
845	if (splitType == ESplitHeuristic \|\|
846	(splitType == ESplitHybrid && leaf->depth > mHybridDepth))
847	useCyclingAxis = false;
848
849	int cyclingAxis = 0;
850
851	if (leaf->parent)
852	cyclingAxis = (leaf->parent->axis + 1)%5;
853
854	if (0 && useCyclingAxis) {
855	int axis = cyclingAxis;
856	info.position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
857	info.axis = axis;
858	info.costRatio = 0.5f; // not true but good for speeding up the subdivision at the top levels!
859	info.raysBack = info.raysFront = leaf->rays.size()/2;
860	return;
861	}
862
863
864	for (int axis = 0; axis < 5; axis++)
865	if (
866	(axis < 3 && (leaf->depth < mDirSplitDepth \|\| mInterleaveDirSplits)) \|\|
867	(axis >= 3 && (leaf->depth >= mDirSplitDepth))
868	) {
869	nInfo[axis].axis = axis;
870	if ((!useCyclingAxis && (!mSplitUseOnlyDrivingAxis \|\| axis == sAxis \|\| axis == dAxis))
871	\|\| (useCyclingAxis && (axis == cyclingAxis))
872	) {
873
874	if (splitType == ESplitRegular) {
875	if (axis < 3)
876	nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
877	else
878	nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
879	EvalCostRatio(leaf,
880	nInfo[axis]);
881	} else
882	if (splitType == ESplitHeuristic) {
883	EvalCostRatioHeuristic(
884	leaf,
885	nInfo[axis]
886	);
887	} else
888	if (splitType == ESplitHybrid) {
889	if (leaf->depth > mHybridDepth)
890	EvalCostRatioHeuristic(
891	leaf,
892	nInfo[axis]
893	);
894	else {
895	if (axis < 3)
896	nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
897	else
898	nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
899
900	EvalCostRatio(leaf,
901	nInfo[axis]
902	);
903	}
904	} else {
905	cerr<<"RssTree: Unknown split heuristics\n";
906	exit(1);
907	}
908
909	if ( bestAxis == -1)
910	bestAxis = axis;
911	else
912	if ( nInfo[axis].costRatio < nInfo[bestAxis].costRatio )
913	bestAxis = axis;
914	}
915	}
916
917	info = nInfo[bestAxis];
918	}
919
920
921	void
922	RssTree::EvalCostRatioHeuristic(
923	RssTreeLeaf *leaf,
924	SplitInfo &info
925	)
926	{
927	AxisAlignedBox3 box;
928	float minBox, maxBox;
929
930	if (info.axis < 3) {
931	box = GetBBox(leaf);
932	minBox = box.Min(info.axis);
933	maxBox = box.Max(info.axis);
934	} else {
935	box = GetDirBBox(leaf);
936	minBox = box.Min(info.axis-3);
937	maxBox = box.Max(info.axis-3);
938	}
939
940	SortSubdivisionCandidates(leaf, info.axis);
941
942	// go through the lists, count the number of objects left and right
943	// and evaluate the following cost funcion:
944	// C = ct_div_ci + (qlrl + qrrr)/queries
945
946	SplitInfo currInfo;
947	currInfo.axis = info.axis;
948
949	currInfo.raysBack = 0;
950	currInfo.raysFront = leaf->rays.size();
951
952	currInfo.pvsBack = 0;
953	currInfo.pvsFront = leaf->GetPvsSize();
954
955
956	float sumWeights = Limits::Small;
957	float sumWeightsBack = Limits::Small;
958	float sumWeightsFront = Limits::Small;
959
960	float sumContribution = 0.0f;
961	float sumContributionBack = 0.0f;
962	float sumContributionFront = 0.0f;
963
964
965	float sizeBox = maxBox - minBox;
966
967	float minBand = minBox + 0.1f*(maxBox - minBox);
968	float maxBand = minBox + 0.9f*(maxBox - minBox);
969
970	// best cost ratio
971	info.costRatio = 1e20f;
972
973	currInfo.viewCells = 0;
974	currInfo.rays = 0;
975
976	Intersectable::NewMail();
977	// set all object as belonging to the fron pvs
978	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
979	ri != leaf->rays.end();
980	ri++)
981	if ((*ri).mRay->IsActive()) {
982
983	float weight, contribution;
984	GetRayContribution(*ri, weight, contribution);
985
986	sumWeights += weight;
987	sumContribution += contribution;
988	currInfo.rays++;
989
990	Intersectable object = (ri).GetObject();
991	if (object)
992	if (!object->Mailed()) {
993	object->Mail();
994	object->mCounter = 1;
995	} else
996	object->mCounter++;
997
998	#if EVAL_VIEWCELLS
999	// and do the same for all viewcells
1000	ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
1001
1002	for (; it != (*ri).mRay->mViewCells.end(); ++it) {
1003	ViewCell viewcell = it;
1004	if (!viewcell->Mailed()) {
1005	currInfo.viewCells++;
1006	viewcell->Mail();
1007	viewcell->mCounter = 1;
1008	} else
1009	viewcell->mCounter++;
1010	}
1011	#endif
1012	}
1013
1014	if (sumWeights!=0.0f)
1015	currInfo.contribution = sumContribution/sumWeights;
1016	else
1017	currInfo.contribution = 0.0f;
1018
1019	sumWeightsFront = sumWeights;
1020	sumContributionFront = sumContribution;
1021
1022	sumWeightsBack = 0;
1023	sumContributionBack = 0;
1024
1025	currInfo.viewCellsBack = 0;
1026	currInfo.viewCellsFront = currInfo.viewCells;
1027
1028	Intersectable::NewMail();
1029
1030	for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
1031	ci < splitCandidates->end();
1032	ci++) {
1033	switch ((*ci).type) {
1034	case SortableEntry::ERayMin: {
1035	currInfo.raysFront--;
1036	currInfo.raysBack++;
1037
1038	RssTreeNode::RayInfo rayInfo = (RssTreeNode::RayInfo ) (*ci).data;
1039
1040	float weight, contribution;
1041	GetRayContribution(*rayInfo, weight, contribution);
1042
1043	sumWeightsBack += weight;
1044	sumContributionBack += contribution;
1045
1046	sumWeightsFront -= weight;
1047	sumContributionFront -= contribution;
1048
1049	Intersectable *object = rayInfo->GetObject();
1050	if (object) {
1051	if (!object->Mailed()) {
1052	object->Mail();
1053	currInfo.pvsBack++;
1054	}
1055	if (--object->mCounter == 0)
1056	currInfo.pvsFront--;
1057	}
1058
1059	#if EVAL_VIEWCELLS
1060	ViewCellContainer::const_iterator it = rayInfo->mRay->mViewCells.begin();
1061	for (; it != rayInfo->mRay->mViewCells.end(); ++it) {
1062	ViewCell viewcell = it;
1063	if (!viewcell->Mailed()) {
1064	viewcell->Mail();
1065	currInfo.viewCellsBack++;
1066	}
1067	if (--viewcell->mCounter == 0)
1068	currInfo.viewCellsFront--;
1069	}
1070	#endif
1071
1072	break;
1073	}
1074	}
1075
1076	float position = (*ci).value;
1077
1078	if (position > minBand && position < maxBand) {
1079	currInfo.position = position;
1080
1081
1082	if (sumWeightsBack!=0.0f)
1083	info.contributionBack = sumContributionBack/sumWeightsBack;
1084	else
1085	info.contributionBack = 0.0f;
1086
1087	if (sumWeightsFront!=0.0f)
1088	info.contributionFront = sumContributionFront/sumWeightsFront;
1089	else
1090	info.contributionFront = 0.0f;
1091
1092	GetCostRatio(
1093	leaf,
1094	currInfo);
1095
1096
1097	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
1098	// cout<<"cost= "<<sum<<endl;
1099
1100	if (currInfo.costRatio < info.costRatio) {
1101	info = currInfo;
1102	}
1103	}
1104	}
1105
1106
1107	// cout<<"===================="<<endl;
1108	// cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
1109	// <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
1110	}
1111
1112	void
1113	RssTree::SortSubdivisionCandidates(
1114	RssTreeLeaf *node,
1115	const int axis
1116	)
1117	{
1118
1119	splitCandidates->clear();
1120
1121	int requestedSize = 2*(node->rays.size());
1122	// creates a sorted split candidates array
1123	if (splitCandidates->capacity() > 500000 &&
1124	requestedSize < (int)(splitCandidates->capacity()/10) ) {
1125
1126	delete splitCandidates;
1127	splitCandidates = new vector<SortableEntry>;
1128	}
1129
1130	splitCandidates->reserve(requestedSize);
1131
1132	// insert all queries
1133	for(RssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
1134	ri < node->rays.end();
1135	ri++) {
1136	if ((*ri).mRay->IsActive()) {
1137	if (axis < 3) {
1138	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
1139	(*ri).GetOrigin(axis),
1140	(void )&(ri))
1141	);
1142	} else {
1143	float pos = (*ri).GetDirParametrization(axis-3);
1144	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
1145	pos,
1146	(void )&(ri))
1147	);
1148	}
1149	}
1150	}
1151
1152	stable_sort(splitCandidates->begin(), splitCandidates->end());
1153	}
1154
1155
1156	void
1157	RssTree::EvaluateLeafStats(const TraversalData &data)
1158	{
1159
1160	// the node became a leaf -> evaluate stats for leafs
1161	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1162
1163	if (data.depth >= termMaxDepth)
1164	stat.maxDepthNodes++;
1165
1166	// if ( (int)(leaf->rays.size()) < termMinCost)
1167	// stat.minCostNodes++;
1168	if ( leaf->GetPvsSize() <= termMinPvs)
1169	stat.minPvsNodes++;
1170
1171	if ( leaf->GetPvsSize() <= termMinRays)
1172	stat.minRaysNodes++;
1173
1174	if (leaf->GetAvgRayContribution() > termMaxRayContribution )
1175	stat.maxRayContribNodes++;
1176
1177	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
1178	stat.minSizeNodes++;
1179	}
1180
1181	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
1182	stat.maxRayRefs = leaf->rays.size();
1183
1184	}
1185
1186	bool
1187	RssTree::TerminationCriteriaSatisfied(RssTreeLeaf *leaf)
1188	{
1189	return ( (leaf->GetPvsSize() <= termMinPvs) \|\|
1190	(leaf->rays.size() <= termMinRays) \|\|
1191	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) \|\|
1192	(leaf->depth >= termMaxDepth) \|\|
1193	(SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
1194	);
1195	}
1196
1197
1198	RssTreeNode *
1199	RssTree::SubdivideNode(
1200	RssTreeLeaf *leaf,
1201	const AxisAlignedBox3 &box,
1202	AxisAlignedBox3 &backBBox,
1203	AxisAlignedBox3 &frontBBox
1204	)
1205	{
1206
1207	if (TerminationCriteriaSatisfied(leaf)) {
1208	#if 0
1209	if (leaf->depth >= termMaxDepth) {
1210	cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
1211	cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
1212	}
1213	#endif
1214
1215	return leaf;
1216	}
1217
1218	SplitInfo info;
1219
1220	// select subdivision axis
1221	int axis = SelectPlane( leaf,
1222	box,
1223	info
1224	);
1225	// Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
1226
1227	if (axis == -1) {
1228	return leaf;
1229	}
1230
1231	stat.nodes+=2;
1232	stat.leaves += 1;
1233	stat.splits[axis]++;
1234
1235	// add the new nodes to the tree
1236	RssTreeInterior *node = new RssTreeInterior(leaf->parent);
1237
1238	node->axis = axis;
1239	node->position = info.position;
1240	// node->bbox = box;
1241	// node->dirBBox = GetDirBBox(leaf);
1242
1243
1244	RssTreeLeaf *back = new RssTreeLeaf(node, info.raysBack);
1245	RssTreeLeaf *front = new RssTreeLeaf(node, info.raysFront);
1246
1247
1248	// update halton generator
1249	back->halton.index = leaf->halton.index;
1250	front->halton.index = leaf->halton.index;
1251
1252	// replace a link from node's parent
1253	if ( leaf->parent )
1254	leaf->parent->ReplaceChildLink(leaf, node);
1255	// and setup child links
1256	node->SetupChildLinks(back, front);
1257
1258	back->bbox = leaf->bbox;
1259	front->bbox = leaf->bbox;
1260	back->dirBBox = leaf->dirBBox;
1261	front->dirBBox = leaf->dirBBox;
1262
1263	if (axis <= RssTreeNode::SPLIT_Z) {
1264	back->bbox.SetMax(axis, info.position);
1265	front->bbox.SetMin(axis, info.position);
1266	} else {
1267	back->dirBBox.SetMax(axis-3, info.position);
1268	front->dirBBox.SetMin(axis-3, info.position);
1269	}
1270
1271	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1272	ri != leaf->rays.end();
1273	ri++) {
1274	if ((*ri).mRay->IsActive()) {
1275
1276	// first unref ray from the former leaf
1277	(*ri).mRay->Unref();
1278
1279	// Debug << "computed t: " << (*ri).mRay->mT << endl;
1280	// determine the side of this ray with respect to the plane
1281	int side = node->ComputeRaySide(*ri);
1282
1283	if (side == 1)
1284	front->AddRay(*ri);
1285	else
1286	back->AddRay(*ri);
1287	} else
1288	(*ri).mRay->Unref();
1289	}
1290
1291	// distribute the total number of rays according to the distribution
1292	// of rays which remained
1293	// front->mTotalRays = front->rays.size()*leaf->mTotalRays/leaf->rays.size();
1294	// back->mTotalRays = back->rays.size()*leaf->mTotalRays/leaf->rays.size();
1295
1296	#if 0
1297	front->SetPvsSize(pvsFront);
1298	back->SetPvsSize(pvsBack);
1299	// compute entropy as well
1300	front->ComputeEntropyImportance();
1301	back->ComputeEntropyImportance();
1302	#else
1303	UpdatePvsSize(front);
1304	UpdatePvsSize(back);
1305	#endif
1306
1307
1308	// update stats
1309	stat.rayRefs -= (int)leaf->rays.size();
1310	stat.rayRefs += back->rays.size() + front->rays.size();
1311
1312
1313	delete leaf;
1314	return node;
1315	}
1316
1317
1318
1319
1320
1321
1322	int
1323	RssTree::ReleaseMemory(const int time)
1324	{
1325	stack<RssTreeNode *> tstack;
1326
1327	// find a node in the tree which subtree will be collapsed
1328	int maxAccessTime = time - accessTimeThreshold;
1329	int released;
1330
1331	PushRoots(tstack);
1332
1333	while (!tstack.empty()) {
1334	RssTreeNode *node = tstack.top();
1335	tstack.pop();
1336
1337
1338	if (!node->IsLeaf()) {
1339	RssTreeInterior in = (RssTreeInterior )node;
1340	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1341	if (in->depth >= minCollapseDepth &&
1342	in->lastAccessTime <= maxAccessTime) {
1343	released = CollapseSubtree(node, time);
1344	break;
1345	}
1346
1347	if (in->back->GetAccessTime() <
1348	in->front->GetAccessTime()) {
1349	tstack.push(in->front);
1350	tstack.push(in->back);
1351	} else {
1352	tstack.push(in->back);
1353	tstack.push(in->front);
1354	}
1355	}
1356	}
1357
1358	while (tstack.empty()) {
1359	// could find node to collaps...
1360	// cout<<"Could not find a node to release "<<endl;
1361	break;
1362	}
1363
1364	return released;
1365	}
1366
1367
1368
1369
1370	RssTreeNode *
1371	RssTree::SubdivideLeaf(
1372	RssTreeLeaf *leaf
1373	)
1374	{
1375	RssTreeNode *node = leaf;
1376
1377	AxisAlignedBox3 leafBBox = GetBBox(leaf);
1378
1379	static int pass = 0;
1380	pass ++;
1381
1382	// check if we should perform a dynamic subdivision of the leaf
1383	if (!TerminationCriteriaSatisfied(leaf)) {
1384
1385	// memory check and realese...
1386	if (GetMemUsage() > maxTotalMemory) {
1387	ReleaseMemory( pass );
1388	}
1389
1390	AxisAlignedBox3 backBBox, frontBBox;
1391
1392	// subdivide the node
1393	node =
1394	SubdivideNode(leaf,
1395	leafBBox,
1396	backBBox,
1397	frontBBox
1398	);
1399	}
1400
1401	return node;
1402	}
1403
1404
1405
1406	void
1407	RssTree::UpdateRays(
1408	VssRayContainer &remove,
1409	VssRayContainer &add
1410	)
1411	{
1412	RssTreeLeaf::NewMail();
1413
1414	// schedule rays for removal
1415	for(VssRayContainer::const_iterator ri = remove.begin();
1416	ri != remove.end();
1417	ri++) {
1418	(*ri)->ScheduleForRemoval();
1419	}
1420
1421	int inactive=0;
1422
1423	for(VssRayContainer::const_iterator ri = remove.begin();
1424	ri != remove.end();
1425	ri++) {
1426	if ((*ri)->ScheduledForRemoval())
1427	// RemoveRay(*ri, NULL, false);
1428	// !!! BUG - with true it does not work correctly - aggreated delete
1429	RemoveRay(*ri, NULL, true);
1430	else
1431	inactive++;
1432	}
1433
1434
1435	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1436
1437	for(VssRayContainer::const_iterator ri = add.begin();
1438	ri != add.end();
1439	ri++) {
1440	RssTreeNode::RayInfo info(*ri);
1441	// if (mForcedBoundingBox==NULL \|\| ClipRay(info, bbox))
1442	AddRay(info);
1443	}
1444
1445	stat.rayRefs += add.size() - remove.size();
1446
1447	UpdateTreeStatistics();
1448	// check whether the tree should be prunned
1449	if (stat.rayRefs > mMaxRays) {
1450	PruneRays(mMaxRays);
1451	// UpdateTreeStatistics();
1452	}
1453
1454	}
1455
1456
1457
1458
1459	void
1460	RssTree::RemoveRay(VssRay *ray,
1461	vector<RssTreeLeaf > affectedLeaves,
1462	const bool removeAllScheduledRays
1463	)
1464	{
1465
1466	stack<RayTraversalData> tstack;
1467
1468	PushRoots(tstack, RssTreeLeaf::RayInfo(ray));
1469
1470	RayTraversalData data;
1471
1472	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1473
1474	while (!tstack.empty()) {
1475	data = tstack.top();
1476	tstack.pop();
1477
1478	if (!data.node->IsLeaf()) {
1479	// split the set of rays in two groups intersecting the
1480	// two subtrees
1481
1482	TraverseInternalNode(data, tstack);
1483
1484	} else {
1485	// remove the ray from the leaf
1486	// find the ray in the leaf and swap it with the last ray...
1487	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1488
1489	if (!leaf->Mailed()) {
1490	leaf->Mail();
1491	if (affectedLeaves)
1492	affectedLeaves->push_back(leaf);
1493
1494	if (removeAllScheduledRays) {
1495	int tail = (int)leaf->rays.size()-1;
1496
1497	for (int i=0; i < (int)(leaf->rays.size()); i++) {
1498	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1499	// find a ray to replace it with
1500	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1501	stat.removedRayRefs++;
1502	leaf->rays[tail].mRay->Unref();
1503	leaf->rays.pop_back();
1504	tail--;
1505	}
1506
1507	if (tail < i)
1508	break;
1509
1510	stat.removedRayRefs++;
1511	leaf->rays[i].mRay->Unref();
1512	leaf->rays[i] = leaf->rays[tail];
1513	leaf->rays.pop_back();
1514	tail--;
1515	}
1516	}
1517	}
1518	}
1519
1520	if (!removeAllScheduledRays)
1521	for (int i=0; i < (int)leaf->rays.size(); i++) {
1522	if (leaf->rays[i].mRay == ray) {
1523	stat.removedRayRefs++;
1524	ray->Unref();
1525	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1526	leaf->rays.pop_back();
1527	// check this ray again
1528	break;
1529	}
1530	}
1531
1532	}
1533	}
1534
1535	if (ray->RefCount() != 0) {
1536	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1537	exit(1);
1538	}
1539
1540	}
1541
1542
1543	void
1544	RssTree::AddRay(RssTreeNode::RayInfo &info)
1545	{
1546
1547	stack<RayTraversalData> tstack;
1548
1549	RssTreeNode *root = GetRoot(info.GetSourceObject());
1550	tstack.push(RayTraversalData(root, info));
1551
1552	RayTraversalData data;
1553
1554	while (!tstack.empty()) {
1555	data = tstack.top();
1556	tstack.pop();
1557
1558	if (!data.node->IsLeaf()) {
1559	TraverseInternalNode(data, tstack);
1560	} else {
1561	// remove the ray from the leaf
1562	// find the ray in the leaf and swap it with the last ray...
1563	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1564	leaf->AddRay(data.rayData);
1565	stat.addedRayRefs++;
1566	}
1567	}
1568	}
1569
1570	void
1571	RssTree::TraverseInternalNode(
1572	RayTraversalData &data,
1573	stack<RayTraversalData> &tstack)
1574	{
1575	RssTreeInterior in = (RssTreeInterior ) data.node;
1576
1577	// determine the side of this ray with respect to the plane
1578	int side = in->ComputeRaySide(data.rayData
1579	);
1580
1581	if (side == 1)
1582	tstack.push(RayTraversalData(in->front, data.rayData));
1583	else
1584	tstack.push(RayTraversalData(in->back, data.rayData));
1585
1586	}
1587
1588
1589	int
1590	RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
1591	{
1592	// first count all rays in the subtree
1593	// use mail 1 for this purpose
1594	stack<RssTreeNode *> tstack;
1595	int rayCount = 0;
1596	int totalRayCount = 0;
1597	int collapsedNodes = 0;
1598
1599	#if DEBUG_COLLAPSE
1600	cout<<"Collapsing subtree"<<endl;
1601	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1602	cout<<"depth="<<(int)sroot->depth<<endl;
1603	#endif
1604
1605	// tstat.collapsedSubtrees++;
1606	// tstat.collapseDepths += (int)sroot->depth;
1607	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1608
1609	tstack.push(sroot);
1610	VssRay::NewMail();
1611
1612	while (!tstack.empty()) {
1613	collapsedNodes++;
1614	RssTreeNode *node = tstack.top();
1615	tstack.pop();
1616
1617	if (node->IsLeaf()) {
1618	RssTreeLeaf leaf = (RssTreeLeaf ) node;
1619	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1620	ri != leaf->rays.end();
1621	ri++) {
1622
1623	totalRayCount++;
1624	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
1625	(*ri).mRay->Mail();
1626	rayCount++;
1627	}
1628	}
1629	} else {
1630	tstack.push(((RssTreeInterior *)node)->back);
1631	tstack.push(((RssTreeInterior *)node)->front);
1632	}
1633	}
1634
1635	VssRay::NewMail();
1636
1637	// create a new node that will hold the rays
1638	RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
1639	if ( newLeaf->parent )
1640	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1641
1642
1643	tstack.push( sroot );
1644
1645	while (!tstack.empty()) {
1646
1647	RssTreeNode *node = tstack.top();
1648	tstack.pop();
1649
1650	if (node->IsLeaf()) {
1651	RssTreeLeaf leaf = (RssTreeLeaf ) node;
1652
1653	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1654	ri != leaf->rays.end();
1655	ri++) {
1656
1657	// unref this ray from the old node
1658
1659	if ((*ri).mRay->IsActive()) {
1660	(*ri).mRay->Unref();
1661	if (!(*ri).mRay->Mailed()) {
1662	(*ri).mRay->Mail();
1663	newLeaf->AddRay(*ri);
1664	}
1665	} else
1666	(*ri).mRay->Unref();
1667
1668	}
1669	} else {
1670	tstack.push(((RssTreeInterior *)node)->back);
1671	tstack.push(((RssTreeInterior *)node)->front);
1672	}
1673	}
1674
1675	// delete the node and all its children
1676	delete sroot;
1677
1678	// for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1679	// ri != newLeaf->rays.end();
1680	// ri++)
1681	// (*ri).ray->UnMail(2);
1682
1683
1684	#if DEBUG_COLLAPSE
1685	cout<<"Total memory before="<<GetMemUsage()<<endl;
1686	#endif
1687
1688	stat.nodes -= collapsedNodes - 1;
1689	stat.leaves -= collapsedNodes/2 - 1;
1690	stat.rayRefs -= totalRayCount - rayCount;
1691
1692	#if DEBUG_COLLAPSE
1693	cout<<"collapsed nodes"<<collapsedNodes<<endl;
1694	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1695	cout<<"Total memory after="<<GetMemUsage()<<endl;
1696	cout<<"================================"<<endl;
1697	#endif
1698
1699	// tstat.collapsedNodes += collapsedNodes;
1700	// tstat.collapsedRays += totalRayCount - rayCount;
1701
1702	return totalRayCount - rayCount;
1703	}
1704
1705
1706	int
1707	RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
1708	{
1709	stack<RssTreeNode *> tstack;
1710	PushRoots(tstack);
1711
1712	Intersectable::NewMail();
1713	int pvsSize = 0;
1714
1715	while (!tstack.empty()) {
1716	RssTreeNode *node = tstack.top();
1717	tstack.pop();
1718
1719
1720	if (node->IsLeaf()) {
1721	RssTreeLeaf leaf = (RssTreeLeaf )node;
1722	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1723	ri != leaf->rays.end();
1724	ri++)
1725	if ((*ri).mRay->IsActive()) {
1726	Intersectable *object;
1727	object = (*ri).GetObject();
1728	if (object && !object->Mailed()) {
1729	pvsSize++;
1730	object->Mail();
1731	}
1732	}
1733	} else {
1734	RssTreeInterior in = (RssTreeInterior )node;
1735	if (in->axis < 3) {
1736	if (box.Max(in->axis) >= in->position )
1737	tstack.push(in->front);
1738
1739	if (box.Min(in->axis) <= in->position )
1740	tstack.push(in->back);
1741	} else {
1742	// both nodes for directional splits
1743	tstack.push(in->front);
1744	tstack.push(in->back);
1745	}
1746	}
1747	}
1748	return pvsSize;
1749	}
1750
1751	int
1752	RssTree::CollectPvs(const AxisAlignedBox3 &box,
1753	ObjectContainer &pvs
1754	) const
1755	{
1756	stack<RssTreeNode *> tstack;
1757	PushRoots(tstack);
1758
1759	Intersectable::NewMail();
1760
1761	while (!tstack.empty()) {
1762	RssTreeNode *node = tstack.top();
1763	tstack.pop();
1764
1765
1766	if (node->IsLeaf()) {
1767	RssTreeLeaf leaf = (RssTreeLeaf )node;
1768	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1769	ri != leaf->rays.end();
1770	ri++)
1771	if ((*ri).mRay->IsActive()) {
1772	Intersectable *object;
1773	object = (*ri).GetObject();
1774	if (object && !object->Mailed()) {
1775	pvs.push_back(object);
1776	object->Mail();
1777	}
1778	}
1779	} else {
1780	RssTreeInterior in = (RssTreeInterior )node;
1781	if (in->axis < 3) {
1782	if (box.Max(in->axis) >= in->position )
1783	tstack.push(in->front);
1784
1785	if (box.Min(in->axis) <= in->position )
1786	tstack.push(in->back);
1787	} else {
1788	// both nodes for directional splits
1789	tstack.push(in->front);
1790	tstack.push(in->back);
1791	}
1792	}
1793	}
1794	return pvs.size();
1795	}
1796
1797	void
1798	RssTree::UpdateTreeStatistics()
1799	{
1800	stack<RssTreeNode *> tstack;
1801	PushRoots(tstack);
1802
1803	float sumPvsSize = 0.0f;
1804	float sumRayContribution = 0.0f;
1805	float sumWeightedRayContribution = 0.0f;
1806	float sumImportance = 0.0f;
1807	float sumPvsEntropy = 0.0f;
1808	float sumRayLengthEntropy = 0.0f;
1809	float sumRays = 0.0f;
1810
1811	int leaves = 0;
1812
1813	while (!tstack.empty()) {
1814	RssTreeNode *node = tstack.top();
1815	tstack.pop();
1816
1817	if (node->IsLeaf()) {
1818	leaves++;
1819	RssTreeLeaf leaf = (RssTreeLeaf )node;
1820	UpdatePvsSize(leaf);
1821
1822	sumPvsSize += leaf->GetPvsSize();
1823	sumRayContribution += leaf->GetAvgRayContribution();
1824
1825
1826	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
1827	for (;it != leaf->rays.end(); ++it) {
1828	float weight, contribution;
1829	GetRayContribution(*it, weight, contribution);
1830	sumWeightedRayContribution += weight*contribution;
1831	}
1832
1833	// sumPvsEntropy += leaf->mPvsEntropy;
1834	// sumRayLengthEntropy += leaf->mRayLengthEntropy;
1835	sumRays += leaf->rays.size();
1836
1837	float imp = leaf->GetImportance();
1838
1839	// if (imp > 1.0f)
1840	// cout<<"warning imp > 1.0f:"<<imp<<endl;
1841
1842	sumImportance += imp;
1843
1844
1845	} else {
1846	RssTreeInterior in = (RssTreeInterior )node;
1847	// both nodes for directional splits
1848	tstack.push(in->front);
1849	tstack.push(in->back);
1850	}
1851	}
1852
1853	stat.avgPvsSize = sumPvsSize/(float)leaves;
1854	stat.avgRays = sumRays/(float)leaves;
1855	stat.avgRayContribution = sumRayContribution/(float)leaves;
1856	// avgPvsEntropy = sumPvsEntropy/(float)leaves;
1857	// avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
1858	stat.avgImportance = sumImportance/(float)leaves;
1859	stat.avgWeightedRayContribution = sumWeightedRayContribution/(float)sumRays;
1860	stat.rayRefs = (int)sumRays;
1861	}
1862
1863
1864
1865	int
1866	RssTree::GenerateRays(
1867	const float ratioPerLeaf,
1868	SimpleRayContainer &rays)
1869	{
1870	stack<RssTreeNode *> tstack;
1871	PushRoots(tstack);
1872
1873	while (!tstack.empty()) {
1874	RssTreeNode *node = tstack.top();
1875	tstack.pop();
1876
1877	if (node->IsLeaf()) {
1878	RssTreeLeaf leaf = (RssTreeLeaf )node;
1879	float c = leaf->GetImportance();
1880	int num = int(c*ratioPerLeaf + 0.5f);
1881	// cout<<num<<" ";
1882
1883	for (int i=0; i < num; i++) {
1884	Vector3 origin = GetBBox(leaf).GetRandomPoint();
1885	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1886	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1887	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1888	rays.push_back(SimpleRay(origin, direction));
1889	}
1890	} else {
1891	RssTreeInterior in = (RssTreeInterior )node;
1892	// both nodes for directional splits
1893	tstack.push(in->front);
1894	tstack.push(in->back);
1895	}
1896	}
1897
1898	return (int)rays.size();
1899	}
1900
1901	void
1902	RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
1903	{
1904	stack<RssTreeNode *> tstack;
1905	PushRoots(tstack);
1906
1907	while (!tstack.empty()) {
1908	RssTreeNode *node = tstack.top();
1909	tstack.pop();
1910
1911	if (node->IsLeaf()) {
1912	RssTreeLeaf leaf = (RssTreeLeaf )node;
1913	leaves.push_back(leaf);
1914	} else {
1915	RssTreeInterior in = (RssTreeInterior )node;
1916	// both nodes for directional splits
1917	tstack.push(in->front);
1918	tstack.push(in->back);
1919	}
1920	}
1921	}
1922
1923	bool
1924	RssTree::ValidLeaf(RssTreeLeaf *leaf) const
1925	{
1926	return true;
1927	//return leaf->rays.size() > termMinRays/4;
1928	}
1929
1930
1931	void
1932	RssTree::PickEdgeRays(RssTreeLeaf *leaf,
1933	int &r1,
1934	int &r2
1935	)
1936	{
1937	int nrays = leaf->rays.size();
1938
1939	if (nrays == 2) {
1940	r1 = 0;
1941	r2 = 1;
1942	return;
1943	}
1944
1945	#if 1
1946	int tries = min(20, nrays);
1947
1948	while (--tries) {
1949	r1 = Random(nrays);
1950	r2 = Random(nrays);
1951	if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
1952	break;
1953	}
1954
1955	if (r1 == r2)
1956	r2 = (r1+1)%leaf->rays.size();
1957
1958	#else
1959	// pick a random ray
1960	int base = Random(nrays);
1961	RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
1962	Intersectable *baseObject = baseRay->GetObject();
1963
1964	// and a random 5D derivative which will be used to find the minimal projected distances
1965	Vector3 spatialDerivative;
1966	Vector3 directionalDerivative;
1967
1968	while (1) {
1969	spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
1970	RandomValue(-1.0f, 1.0f),
1971	RandomValue(-1.0f, 1.0f));
1972	float m = Magnitude(spatialDerivative);
1973	if (m != 0) {
1974	spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
1975	break;
1976	}
1977	}
1978
1979	while (1) {
1980	directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
1981	RandomValue(-1.0f, 1.0f),
1982	0.0f);
1983	float m = Magnitude(directionalDerivative);
1984	if (m != 0) {
1985	directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
1986	break;
1987	}
1988	}
1989
1990	// now find the furthest sample from the same object and the closest from a different object
1991	int i = 0;
1992	float minDist = MAX_FLOAT;
1993	float maxDist = -MAX_FLOAT;
1994	r1 = base;
1995	r2 = base;
1996	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1997	ri != leaf->rays.end();
1998	ri++, i++) {
1999	float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
2000	*ri,
2001	spatialDerivative,
2002	directionalDerivative
2003	);
2004
2005	if ((*ri).GetObject() == baseObject) {
2006	if (dist > maxDist) {
2007	maxDist = dist;
2008	r1 = i;
2009	}
2010	} else {
2011	if (dist > 0.0f && dist < minDist) {
2012	minDist = dist;
2013	r2 = i;
2014	}
2015	}
2016	}
2017
2018	#endif
2019	}
2020
2021
2022	struct RayNeighbor {
2023	int rayInfo;
2024	float distance;
2025	RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
2026	};
2027
2028	void
2029	RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
2030	vector<RssTreeLeaf::SilhouetteRays> &rays
2031	)
2032	{
2033	// for every leaf find its neares neighbor from a different object
2034	vector<RayNeighbor> neighbors(leaf->rays.size());
2035	int i, j;
2036	for (i=0; i < leaf->rays.size(); i++)
2037	for (j=0; j < leaf->rays.size(); j++)
2038	if (i!=j) {
2039	float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
2040	if (d < neighbors[i].distance) {
2041	neighbors[i].distance = d;
2042	neighbors[i].rayInfo = j;
2043	}
2044	}
2045
2046	// now check which are the pairs of nearest neighbors
2047	for (i=0; i < leaf->rays.size(); i++) {
2048	int j = neighbors[i].rayInfo;
2049	if (neighbors[j].rayInfo == i) {
2050	// this is a silhouette edge pair
2051	if (i < j) {
2052	// generate an silhoutte ray pair
2053	rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
2054	&leaf->rays[j]));
2055	}
2056	} else {
2057	// this is not a silhouette ray - delete???
2058	}
2059	}
2060
2061	}
2062
2063
2064	bool
2065	GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
2066	const int index,
2067	int &i1,
2068	int &i2,
2069	int &i3)
2070	{
2071	int found = 0;
2072	int indices[3];
2073
2074	int size = rays.size();
2075	// use russian roulete selection for the tripple
2076	// number of free positions for the bullet
2077	int positions = size - 1;
2078	int i;
2079	for (i=0; i < size; i++) {
2080	if (rays[i]->mRay->Mailed())
2081	positions--;
2082	}
2083
2084	if (positions < 3)
2085	return false;
2086
2087	for (i=0; i < size; i++) {
2088	if (i != index && !rays[i]->mRay->Mailed()) {
2089	float p = (3 - found)/(float)positions;
2090	if (Random(1.0f) < p) {
2091	indices[found] = i;
2092	found++;
2093	}
2094	}
2095	positions--;
2096	}
2097	return true; // corr. matt
2098	}
2099
2100	bool
2101	RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
2102	const RssTreeNode::RayInfo &r1,
2103	const RssTreeNode::RayInfo &r2,
2104	const RssTreeNode::RayInfo &r3)
2105	{
2106
2107
2108	return false;
2109	}
2110
2111	int
2112	RssTree::RemoveInteriorRays(
2113	RssTreeLeaf *leaf
2114	)
2115	{
2116	#if 1
2117	// first collect all objects refered in this leaf
2118	map<Intersectable , vector<RssTreeLeaf::RayInfo > > rayMap;
2119
2120	RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
2121	for (; it != leaf->rays.end(); ++it) {
2122	Intersectable object = (it).GetObject();
2123
2124	rayMap[object].push_back(&(*it));
2125	// vector<RayInfo > data = rayMap.Find(object);
2126	// if (data) {
2127	// data->push_back(&(*it));
2128	// } else {
2129	// // rayMap[object] = vector<RayInfo *>;
2130	// rayMap[object].push_back(&(*it));
2131	// }
2132	}
2133
2134	// now go through all objects
2135	map<Intersectable , vector<RssTreeLeaf::RayInfo > >::iterator mi;
2136
2137	// infos of mailed rays are scheduled for removal
2138	VssRay::NewMail();
2139	for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
2140	vector<RssTreeLeaf::RayInfo > &rays = (mi).second;
2141	vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
2142	int rayIndex = 0;
2143
2144	for (; ri != rays.end(); ++ri, ++rayIndex) {
2145	RssTreeNode::RayInfo ray = ri;
2146	int tries = rays.size();
2147	for (int i = 0; i < tries; i++) {
2148	int r1, r2, r3;
2149	GetRandomTripple(rays,
2150	rayIndex,
2151	r1,
2152	r2,
2153	r3);
2154	if (IsRayConvexCombination(*ray,
2155	*rays[r1],
2156	*rays[r2],
2157	*rays[r3])) {
2158	ray->mRay->Mail();
2159	}
2160	}
2161	}
2162	}
2163
2164
2165	#endif
2166	return 0;
2167	}
2168
2169
2170	void
2171	RssTree::GenerateLeafRaysSimple(RssTreeLeaf *leaf,
2172	const int numberOfRays,
2173	SimpleRayContainer &rays)
2174	{
2175
2176	if (numberOfRays == 0)
2177	return;
2178
2179	int startIndex = (int)rays.size();
2180
2181	Debug<<"B"<<flush;
2182
2183	AxisAlignedBox3 box = GetBBox(leaf);
2184	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2185
2186	float boxSize = Magnitude(box.Diagonal());
2187
2188	const float smoothRange = 0.1f;
2189	if (smoothRange != 0.0f) {
2190	box.Scale(1.0f + smoothRange);
2191	dirBox.Scale(1.0f + smoothRange);
2192	}
2193
2194	Debug<<"R"<<flush;
2195
2196	for (int i=0; i < numberOfRays; i++) {
2197	// Debug<<i<<flush;
2198
2199	Vector3 origin, direction;
2200	Vector3 dirVector;
2201
2202	Vector3 pVector;
2203	Vector3 dVector;
2204
2205	bool useHalton = false;
2206
2207
2208	if (useHalton) {
2209	// generate a random 5D vector
2210	pVector = Vector3(leaf->halton.GetNumber(1),
2211	leaf->halton.GetNumber(2),
2212	leaf->halton.GetNumber(3));
2213
2214	dVector = Vector3(leaf->halton.GetNumber(4),
2215	leaf->halton.GetNumber(5),
2216	0.0f);
2217
2218	leaf->halton.GenerateNext();
2219	} else {
2220	pVector = Vector3(RandomValue(0.0f, 1.0f),
2221	RandomValue(0.0f, 1.0f),
2222	RandomValue(0.0f, 1.0f));
2223
2224	dVector = Vector3(RandomValue(0.0f, 1.0f),
2225	RandomValue(0.0f, 1.0f),
2226	0.0f);
2227	}
2228
2229	origin = box.GetPoint(pVector);
2230	dirVector = dirBox.GetPoint(dVector);
2231
2232	direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
2233
2234
2235
2236	// shift the origin a little bit
2237	direction.Normalize();
2238
2239	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2240	float dist = 0.0f;
2241	float minT, maxT;
2242
2243	// compute interection of the ray with the box
2244	#if 1
2245	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2246	dist = maxT;
2247	#else
2248	dist = 0.5f*boxSize;
2249	#endif
2250
2251	origin += direction*dist;
2252
2253	#if 0
2254	static int counter = 0;
2255	Debug<<counter++<<flush;
2256
2257	if (counter > 10374968) {
2258	Debug<<"size="<<rays.size()<<endl;
2259	Debug<<"capacity="<<rays.capacity()<<endl;
2260	Debug<<"o="<<origin<<endl;
2261	Debug<<"d="<<direction<<endl;
2262	}
2263	#endif
2264	//Debug<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2265	rays.push_back(SimpleRay(origin, direction));
2266	// Debug<<endl;
2267
2268	}
2269
2270	Debug<<"D"<<flush;
2271
2272	float density = 1.0f;
2273	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2274	for (i=startIndex; i < rays.size(); i++) {
2275	rays[i].mPdf = p;
2276	}
2277	}
2278
2279	void
2280	RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
2281	const int numberOfRays,
2282	SimpleRayContainer &rays)
2283	{
2284
2285	if (numberOfRays == 0)
2286	return;
2287
2288
2289	int nrays = (int)leaf->rays.size();
2290	int startIndex = (int)rays.size();
2291
2292	AxisAlignedBox3 box = GetBBox(leaf);
2293	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2294
2295
2296	const float smoothRange = 0.1f;
2297	if (smoothRange != 0.0f) {
2298	box.Scale(1.0f + smoothRange);
2299	dirBox.Scale(1.0f + smoothRange);
2300	}
2301
2302	Exporter *exporter = NULL;
2303	VssRayContainer selectedRays;
2304	static int counter = 0;
2305	if (counter < 0) {
2306	char filename[128];
2307	sprintf(filename, "selected-rays%04d.x3d", counter);
2308	exporter = Exporter::GetExporter(filename);
2309	// exporter->SetWireframe();
2310	// exporter->ExportKdTree(*mKdTree);
2311	exporter->SetWireframe();
2312	// exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
2313	exporter->SetWireframe();
2314	exporter->ExportBox(box);
2315	exporter->SetFilled();
2316	counter++;
2317	}
2318
2319
2320	float extendedConvexCombinationProb = 0.0f; //0.7f
2321	// float silhouetteCheckPercentage = 1.0f; //0.5f
2322	float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
2323	float silhouetteObjectCheckPercentage = 0.8f;
2324
2325	// int numberOfTries = numberOfRays*4;
2326	int numberOfTries = numberOfRays;
2327	int generated = 0;
2328
2329	vector<Intersectable *> pvs(0);
2330	if (silhouetteCheckPercentage != 0.0f) {
2331	Intersectable::NewMail();
2332	pvs.reserve(leaf->GetPvsSize());
2333	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2334	ri != leaf->rays.end();
2335	ri++) {
2336	Intersectable object = (ri).GetObject();
2337	if (object) {
2338	if (!object->Mailed()) {
2339	object->Mail();
2340	object->mCounter = 1;
2341	pvs.push_back(object);
2342	if (exporter)
2343	exporter->ExportIntersectable(object);
2344	} else {
2345	object->mCounter++;
2346	}
2347	}
2348	}
2349	// sort objects based on their mCounter
2350	sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
2351	}
2352
2353	for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
2354	bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
2355	extendedConvexCombinationProb));
2356
2357	Vector3 origin, direction;
2358	// generate half of convex combination and half of random rays
2359	if (useExtendedConvexCombination) {
2360	// pickup 2 random rays
2361	int r1, r2;
2362	PickEdgeRays(leaf, r1, r2);
2363
2364
2365	Vector3 o1 = leaf->rays[r1].GetOrigin();
2366	Vector3 o2 = leaf->rays[r2].GetOrigin();
2367
2368	const float overlap = 0.0f;
2369	float w1, w2;
2370	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2371	origin = w1o1 + w2o2;
2372	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2373	direction =
2374	w1*leaf->rays[r1].GetDir() +
2375	w2*leaf->rays[r2].GetDir();
2376	} else {
2377	Vector3 dirVector;
2378
2379	Vector3 pVector;
2380	Vector3 dVector;
2381
2382	bool useHalton = true;
2383
2384	if (useHalton) {
2385	// generate a random 5D vector
2386	pVector = Vector3(leaf->halton.GetNumber(1),
2387	leaf->halton.GetNumber(2),
2388	leaf->halton.GetNumber(3));
2389
2390	dVector = Vector3(leaf->halton.GetNumber(4),
2391	leaf->halton.GetNumber(5),
2392	0.0f);
2393
2394	leaf->halton.GenerateNext();
2395	} else {
2396	pVector = Vector3(RandomValue(0.0f, 1.0f),
2397	RandomValue(0.0f, 1.0f),
2398	RandomValue(0.0f, 1.0f));
2399
2400	dVector = Vector3(RandomValue(0.0f, 1.0f),
2401	RandomValue(0.0f, 1.0f),
2402	0.0f);
2403	}
2404
2405	origin = box.GetPoint(pVector);
2406	dirVector = dirBox.GetPoint(dVector);
2407
2408	direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
2409	}
2410
2411	// shift the origin a little bit
2412	direction.Normalize();
2413
2414	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2415	float dist = 0.0f;
2416	float minT, maxT;
2417	// compute interection of the ray with the box
2418
2419	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2420	dist = maxT;
2421
2422
2423	origin += direction*dist;
2424
2425	bool intersects = false;
2426
2427	if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
2428	if (exporter) {
2429	VssRay ray = new VssRay(origin, origin + 100direction, NULL, NULL);
2430	selectedRays.push_back(ray);
2431	}
2432
2433	// check whether the ray does not intersect already visible objects
2434	Ray traversalRay;
2435	traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2436	for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2437	oi != pvs.end();
2438	oi++) {
2439	Intersectable object = oi;
2440	// do not test every object
2441	if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
2442	object->CastRay(traversalRay) ) {
2443	intersects = true;
2444	break;
2445	}
2446	}
2447	}
2448
2449
2450	if (!intersects) {
2451	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2452	rays.push_back(SimpleRay(origin, direction));
2453	generated++;
2454	}
2455	}
2456
2457	// cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
2458	// assign pdfs to the generated rays
2459	// float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
2460	float density = 1.0f;
2461
2462	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2463	for (i=startIndex; i < rays.size(); i++) {
2464	rays[i].mPdf = p;
2465	}
2466
2467
2468	if (exporter) {
2469	exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2470	delete exporter;
2471	CLEAR_CONTAINER(selectedRays);
2472	}
2473
2474	}
2475
2476	int
2477	RssTree::PruneRays(RssTreeLeaf *leaf,
2478	const float contributionThreshold)
2479	{
2480	int i;
2481	int j;
2482
2483	for (j=0, i=0; i < leaf->rays.size(); i++) {
2484
2485	if (leaf->rays[i].mRay->mWeightedPvsContribution > contributionThreshold) {
2486	// copy a valid sample
2487	leaf->rays[j] = leaf->rays[i];
2488	j++;
2489	} else {
2490	// delete the ray
2491	leaf->rays[i].mRay->Unref();
2492	if (leaf->rays[i].mRay->RefCount() != 0) {
2493	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2494	exit(1);
2495	}
2496	delete leaf->rays[i].mRay;
2497	}
2498	}
2499
2500
2501	leaf->rays.resize(j);
2502	int removed = (i-j);
2503	stat.rayRefs -= removed;
2504	return removed;
2505	}
2506
2507	int
2508	RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2509	const float ratio)
2510	{
2511	int i;
2512	int j;
2513
2514	for (j=0, i=0; i < leaf->rays.size(); i++) {
2515
2516	if (Random(1.0f) < ratio) {
2517	// copy a valid sample
2518	if (i!=j)
2519	leaf->rays[j] = leaf->rays[i];
2520	j++;
2521	} else {
2522	// delete the ray
2523	leaf->rays[i].mRay->Unref();
2524	if (leaf->rays[i].mRay->RefCount() != 0) {
2525	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2526	exit(1);
2527	}
2528	delete leaf->rays[i].mRay;
2529	}
2530	}
2531
2532
2533	leaf->rays.resize(j);
2534	int removed = (i-j);
2535	stat.rayRefs -= removed;
2536	return removed;
2537	}
2538
2539	int
2540	RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2541	const float ratio)
2542	{
2543	int i;
2544
2545	if (leaf->rays.size() == 0)
2546	return 0;
2547
2548	sort(leaf->rays.begin(),
2549	leaf->rays.end(),
2550	RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2551
2552	int desired = int(ratio*leaf->rays.size());
2553	int removed = (int)leaf->rays.size() - desired;
2554
2555	for (i=desired; i < (int)leaf->rays.size(); i++) {
2556	// delete the ray
2557	leaf->rays[i].mRay->Unref();
2558	if (leaf->rays[i].mRay->RefCount() != 0) {
2559	cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
2560	cerr<<"desired="<<desired<<"i="<<i<<" size="<<(int)leaf->rays.size()<<endl;
2561	exit(1);
2562	}
2563	delete leaf->rays[i].mRay;
2564	}
2565
2566	leaf->rays.resize(desired);
2567	stat.rayRefs -= removed;
2568	return removed;
2569	}
2570
2571
2572	int
2573	RssTree::PruneRays(
2574	const int maxRays
2575	)
2576	{
2577	bool globalPrunning = false;
2578
2579	stack<RssTreeNode *> tstack;
2580	int prunned = 0;
2581
2582	// prune more rays to amortize the procedure
2583	int desired = (int)(maxRays * 0.8f);
2584
2585	Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;
2586
2587	if (globalPrunning) {
2588	VssRayContainer allRays;
2589	allRays.reserve(stat.rayRefs);
2590	CollectRays(allRays);
2591	sort(allRays.begin(),
2592	allRays.end(),
2593	GreaterWeightedPvsContribution);
2594
2595	if ( maxRays >= allRays.size() )
2596	return 0;
2597
2598	float contributionThreshold = allRays[desired]->mWeightedPvsContribution;
2599
2600	PushRoots(tstack);
2601
2602	while (!tstack.empty()) {
2603	RssTreeNode *node = tstack.top();
2604	tstack.pop();
2605
2606	if (node->IsLeaf()) {
2607	RssTreeLeaf leaf = (RssTreeLeaf )node;
2608	prunned += PruneRays(leaf, contributionThreshold);
2609	} else {
2610	RssTreeInterior in = (RssTreeInterior )node;
2611	// both nodes for directional splits
2612	tstack.push(in->front);
2613	tstack.push(in->back);
2614	}
2615	}
2616	} else {
2617	// prune random rays from each leaf so that the ratio's remain the same
2618	PushRoots(tstack);
2619
2620	if ( maxRays >= stat.rayRefs )
2621	return 0;
2622
2623	float ratio = desired/(float)stat.rayRefs;
2624
2625	while (!tstack.empty()) {
2626	RssTreeNode *node = tstack.top();
2627	tstack.pop();
2628
2629	if (node->IsLeaf()) {
2630	RssTreeLeaf leaf = (RssTreeLeaf )node;
2631	// prunned += PruneRaysRandom(leaf, ratio);
2632	prunned += PruneRaysContribution(leaf, ratio);
2633	} else {
2634	RssTreeInterior in = (RssTreeInterior )node;
2635	// both nodes for directional splits
2636	tstack.push(in->front);
2637	tstack.push(in->back);
2638	}
2639	}
2640
2641
2642
2643
2644	}
2645
2646
2647
2648	Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
2649
2650	return prunned;
2651	}
2652
2653	int
2654	RssTree::GenerateRays(const int numberOfRays,
2655	const int numberOfLeaves,
2656	SimpleRayContainer &rays)
2657	{
2658
2659
2660	vector<RssTreeLeaf *> leaves;
2661
2662	Debug<<"Collecting leaves..."<<endl<<flush;
2663	CollectLeaves(leaves);
2664	Debug<<"done."<<endl<<flush;
2665
2666
2667	if (numberOfLeaves != leaves.size()) {
2668	Debug<<"sorting leaves..."<<endl<<flush;
2669	sort(leaves.begin(),
2670	leaves.end(),
2671	GreaterContribution);
2672	Debug<<"done."<<endl<<flush;
2673	}
2674
2675
2676	Debug<<"Evaluating leaf importance..."<<endl<<flush;
2677
2678	float sumContrib = 0.0;
2679	int i;
2680	int k = 0;
2681	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2682	if (ValidLeaf(leaves[i])) {
2683	float c = leaves[i]->GetImportance();
2684	sumContrib += c;
2685	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
2686	k++;
2687	}
2688
2689	Debug<<"done."<<endl<<flush;
2690
2691	float avgContrib = sumContrib/numberOfLeaves;
2692
2693	// always generate at leat n ray per leaf
2694	int fixedPerLeaf = 0;
2695	// int fixedPerLeaf = 0;
2696	int fixed = fixedPerLeaf*(int)leaves.size();
2697	int iGenerated = numberOfRays;
2698	float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
2699
2700	k = 0;
2701
2702	Debug<<"generating leaf rays..."<<endl<<flush;
2703
2704	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2705	if (ValidLeaf(leaves[i])) {
2706	RssTreeLeaf *leaf = leaves[i];
2707	float c = leaf->GetImportance();
2708
2709	int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
2710	GenerateLeafRaysSimple(leaf, num, rays);
2711
2712	k++;
2713	}
2714
2715	Debug<<"done."<<endl<<flush;
2716
2717	return (int)rays.size();
2718	}
2719
2720
2721	float
2722	RssTree::GetAvgPvsSize()
2723	{
2724	stack<RssTreeNode *> tstack;
2725	PushRoots(tstack);
2726
2727	int sumPvs = 0;
2728	int leaves = 0;
2729	while (!tstack.empty()) {
2730	RssTreeNode *node = tstack.top();
2731	tstack.pop();
2732
2733	if (node->IsLeaf()) {
2734	RssTreeLeaf leaf = (RssTreeLeaf )node;
2735	// update pvs size
2736	UpdatePvsSize(leaf);
2737	sumPvs += leaf->GetPvsSize();
2738	leaves++;
2739	} else {
2740	RssTreeInterior in = (RssTreeInterior )node;
2741	// both nodes for directional splits
2742	tstack.push(in->front);
2743	tstack.push(in->back);
2744	}
2745	}
2746
2747
2748	return sumPvs/(float)leaves;
2749	}
2750
2751	float weightAbsContributions = 0.0f;
2752	// if small very high importance of the last sample
2753	// if 1.0f then weighs = 1 1/2 1/3 1/4
2754	float passSampleWeightDecay = 1.0f;
2755	//float passSampleWeightDecay = 0.0001f;
2756
2757
2758	float
2759	RssTree::GetSampleWeight(const int pass)
2760	{
2761	int passDiff = mCurrentPass - pass;
2762	float weight;
2763	if (1)
2764	weight = 1.0f/sqr(passDiff + passSampleWeightDecay);
2765	else
2766	switch (passDiff) {
2767	case 0:
2768	// weight = 1.0f;
2769	// break;
2770	default:
2771	weight = 1.0f;
2772	break;
2773	// case 1:
2774	// weight = 0.5f;
2775	// break;
2776	// case 2:
2777	// weight = 0.25f;
2778	// break;
2779	// case 3:
2780	// weight = 0.12f;
2781	// break;
2782	// case 4:
2783	// weight = 0.06f;
2784	// break;
2785	// default:
2786	// weight = 0.03f;
2787	// break;
2788	}
2789	return weight;
2790	}
2791
2792	void
2793	RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2794	float &weight,
2795	float &contribution)
2796	{
2797	VssRay *ray = info.mRay;
2798	weight = GetSampleWeight(mCurrentPass);
2799	contribution =
2800	weightAbsContributions*ray->mPvsContribution +
2801	(1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2802	// store the computed value
2803	info.mRay->mWeightedPvsContribution = weight*contribution;
2804	}
2805
2806
2807	void
2808	RssTree::ComputeImportance(RssTreeLeaf *leaf)
2809	{
2810	if (0)
2811	leaf->mImportance = leaf->GetAvgRayContribution();
2812	else {
2813	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2814	it_end = leaf->rays.end();
2815
2816	float sumContributions = 0.0f;
2817	float sumRelContributions = 0.0f;
2818	float sumWeights = 0.0f;
2819	for (; it != it_end; ++it) {
2820	VssRay ray = (it).mRay;
2821	float weight = GetSampleWeight(ray->mPass);
2822
2823	sumContributions += weight*ray->mPvsContribution;
2824	//$$ 20.7. changed to sqr to pronouce higher contribution so that they do not get smoothed!!
2825	//sumRelContributions += weight*ray->mRelativePvsContribution;
2826
2827	sumRelContributions += sqr(weight*ray->mRelativePvsContribution);
2828
2829	//sumWeights += weight;
2830	sumWeights += 1.0f;
2831	}
2832	// $$
2833	// sumWeights = leaf->mTotalRays;
2834
2835	if (sumWeights != 0.0f)
2836	leaf->mImportance =
2837	(weightAbsContributions*sumContributions +
2838	(1.0f - weightAbsContributions)*sumRelContributions)/sumWeights;
2839	else
2840	leaf->mImportance = 0.0f;
2841
2842	}
2843
2844	// return GetAvgRayContribution()*mImportance;
2845	//return GetAvgRayContribution();
2846	}
2847
2848
2849	float
2850	RssTreeLeaf::GetImportance() const
2851	{
2852	//return sqr(mImportance);
2853	return mImportance;
2854	}
2855
2856
2857	float
2858	RssTreeLeaf::ComputePvsEntropy() const
2859	{
2860	int samples = 0;
2861	Intersectable::NewMail();
2862	// set all object as belonging to the fron pvs
2863	for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
2864	ri != rays.end();
2865	ri++)
2866	if ((*ri).mRay->IsActive()) {
2867	Intersectable object = (ri).GetObject();
2868	if (object) {
2869	if (!object->Mailed()) {
2870	object->Mail();
2871	object->mCounter = 1;
2872	} else
2873	object->mCounter++;
2874	samples++;
2875	}
2876	}
2877
2878	float entropy = 0.0f;
2879
2880	if (samples > 1) {
2881	Intersectable::NewMail();
2882	for(RayInfoContainer::const_iterator ri = rays.begin();
2883	ri != rays.end();
2884	ri++)
2885	if ((*ri).mRay->IsActive()) {
2886	Intersectable object = (ri).GetObject();
2887	if (object) {
2888	if (!object->Mailed()) {
2889	object->Mail();
2890	float p = object->mCounter/(float)samples;
2891	entropy -= p*log(p);
2892	}
2893	}
2894	}
2895	entropy = entropy/log((float)samples);
2896	}
2897	else
2898	entropy = 1.0f;
2899
2900	return entropy;
2901	}
2902
2903	float
2904	RssTreeLeaf::ComputeRayLengthEntropy() const
2905	{
2906	// get sum of all ray lengths
2907	// consider only passing rays or originating rays
2908	float sum = 0.0f;
2909	int samples = 0;
2910	int i=0;
2911	for(RayInfoContainer::const_iterator ri = rays.begin();
2912	ri != rays.end();
2913	ri++)
2914	if ((*ri).mRay->IsActive()) {
2915	// float s;
2916	// if (i == 0)
2917	// s = 200;
2918	// else
2919	// s = 100;
2920	// i++;
2921
2922	sum += (*ri).mRay->GetSize();
2923	samples++;
2924	}
2925
2926
2927	float entropy = 0.0f;
2928
2929	if (samples > 1) {
2930	i = 0;
2931	for(RayInfoContainer::const_iterator ri = rays.begin();
2932	ri != rays.end();
2933	ri++)
2934	if ((*ri).mRay->IsActive()) {
2935	// float s;
2936	// if (i==0)
2937	// s = 200;
2938	// else
2939	// s = 100;
2940	// i++;
2941	// float p = s/sum;
2942	float p = (*ri).mRay->GetSize()/sum;
2943	entropy -= p*log(p);
2944	}
2945	entropy = entropy/log((float)samples);
2946	} else
2947	entropy = 1.0f;
2948
2949	return entropy;
2950	}
2951
2952
2953
2954	float
2955	RssTreeLeaf::ComputeEntropyImportance() const
2956	{
2957
2958	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2959	return 1.0f - ComputePvsEntropy();
2960	}
2961
2962	float
2963	RssTreeLeaf::ComputeRayContributionImportance() const
2964	{
2965	// mPvsEntropy = ComputePvsEntropy();
2966	// mRayLengthEntropy = ComputeRayLengthEntropy();
2967
2968	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2969	return 1.0f - ComputePvsEntropy();
2970	}
2971
2972
2973	AxisAlignedBox3
2974	RssTree::GetShrankedBBox(const RssTreeNode *node) const
2975	{
2976	if (node->parent == NULL)
2977	return bbox;
2978
2979	if (!node->IsLeaf())
2980	return ((RssTreeInterior *)node)->bbox;
2981
2982	// evaluate bounding box from the ray origins
2983	AxisAlignedBox3 box;
2984	box.Initialize();
2985	RssTreeLeaf leaf = (RssTreeLeaf )node;
2986	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2987	ri != leaf->rays.end();
2988	ri++)
2989	if ((*ri).mRay->IsActive()) {
2990	box.Include((*ri).GetOrigin());
2991	}
2992	return box;
2993	}
2994
2995
2996	int
2997	RssTree::CollectRays(VssRayContainer &rays,
2998	const int number)
2999	{
3000	VssRayContainer allRays;
3001	CollectRays(allRays);
3002
3003
3004	int desired = min(number, (int)allRays.size());
3005	float prob = desired/(float)allRays.size();
3006	//int skip = allRays.size()/desired;
3007	while (rays.size() < desired) {
3008	VssRayContainer::const_iterator it = allRays.begin();
3009	for (; it != allRays.end() && rays.size() < desired; it++) {
3010	if (Random(1.0f) < prob)
3011	rays.push_back(*it);
3012	}
3013	}
3014	return (int)rays.size();
3015	}
3016
3017
3018	int
3019	RssTree::CollectRays(VssRayContainer &rays
3020	)
3021	{
3022	VssRay::NewMail();
3023
3024	stack<RssTreeNode *> tstack;
3025	PushRoots(tstack);
3026
3027	while (!tstack.empty()) {
3028	RssTreeNode *node = tstack.top();
3029	tstack.pop();
3030
3031	if (node->IsLeaf()) {
3032	RssTreeLeaf leaf = (RssTreeLeaf )node;
3033	// update pvs size
3034	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
3035	for (;it != leaf->rays.end(); ++it)
3036	if (!(*it).mRay->Mailed()) {
3037	(*it).mRay->Mail();
3038	rays.push_back((*it).mRay);
3039	}
3040	} else {
3041	RssTreeInterior in = (RssTreeInterior )node;
3042	// both nodes for directional splits
3043	tstack.push(in->front);
3044	tstack.push(in->back);
3045	}
3046	}
3047
3048	return (int)rays.size();
3049	}
3050
3051
3052	RssTreeNode *
3053	RssTree::GetRoot(Intersectable *object) const
3054	{
3055	if (mPerObjectTree && object) {
3056	int id = object->GetId()-1;
3057	if (id < 0 \|\| id >= mRoots.size()) {
3058	Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<(int)mRoots.size()<<
3059	endl<<flush;
3060	id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
3061	}
3062	return mRoots[id];
3063	} else
3064	return mRoots[0];
3065	}
3066
3067	void
3068	RssTree::PushRoots(priority_queue<TraversalData> &st) const
3069	{
3070	for (int i=0; i < mRoots.size(); i++)
3071	st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
3072	}
3073
3074	void
3075	RssTree::PushRoots(stack<RssTreeNode *> &st) const
3076	{
3077	for (int i=0; i < mRoots.size(); i++)
3078	st.push(mRoots[i]);
3079	}
3080
3081	void
3082	RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info) const
3083	{
3084	for (int i=0; i < mRoots.size(); i++)
3085	st.push(RayTraversalData(mRoots[i], info));
3086	}
3087
3088
3089	}

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

Download in other formats: