Context Navigation

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

Revision 1785, 72.9 KB checked in by bittner, 18 years ago (diff)
merge, filter update, renderebuffer made functional

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	// shift the origin a little bit
2236	direction.Normalize();
2237
2238	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2239	float dist = 0.0f;
2240	float minT, maxT;
2241
2242	// compute interection of the ray with the box
2243	#if 1
2244	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2245	dist = maxT;
2246	#else
2247	dist = 0.5f*boxSize;
2248	#endif
2249
2250	origin += direction*dist;
2251
2252	#if 0
2253	static int counter = 0;
2254	// Debug<<counter++<<flush;
2255
2256	if (counter > 10374968) {
2257	Debug<<"size="<<rays.size()<<endl;
2258	Debug<<"capacity="<<rays.capacity()<<endl;
2259	Debug<<"o="<<origin<<endl;
2260	Debug<<"d="<<direction<<endl;
2261	}
2262	#endif
2263	//Debug<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2264	rays.push_back(SimpleRay(origin, direction));
2265	// Debug<<endl;
2266
2267	}
2268
2269	// Debug<<"D"<<flush;
2270
2271	float density = 1.0f;
2272	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2273	for (int i=startIndex; i < rays.size(); i++) {
2274	rays[i].mPdf = p;
2275	}
2276	}
2277
2278	void
2279	RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
2280	const int numberOfRays,
2281	SimpleRayContainer &rays)
2282	{
2283
2284	if (numberOfRays == 0)
2285	return;
2286
2287
2288	int nrays = (int)leaf->rays.size();
2289	int startIndex = (int)rays.size();
2290
2291	AxisAlignedBox3 box = GetBBox(leaf);
2292	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2293
2294
2295	const float smoothRange = 0.1f;
2296	if (smoothRange != 0.0f) {
2297	box.Scale(1.0f + smoothRange);
2298	dirBox.Scale(1.0f + smoothRange);
2299	}
2300
2301	Exporter *exporter = NULL;
2302	VssRayContainer selectedRays;
2303	static int counter = 0;
2304	if (counter < 0) {
2305	char filename[128];
2306	sprintf(filename, "selected-rays%04d.x3d", counter);
2307	exporter = Exporter::GetExporter(filename);
2308	// exporter->SetWireframe();
2309	// exporter->ExportKdTree(*mKdTree);
2310	exporter->SetWireframe();
2311	// exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
2312	exporter->SetWireframe();
2313	exporter->ExportBox(box);
2314	exporter->SetFilled();
2315	counter++;
2316	}
2317
2318
2319	float extendedConvexCombinationProb = 0.0f; //0.7f
2320	// float silhouetteCheckPercentage = 1.0f; //0.5f
2321	float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
2322	float silhouetteObjectCheckPercentage = 0.8f;
2323
2324	// int numberOfTries = numberOfRays*4;
2325	int numberOfTries = numberOfRays;
2326	int generated = 0;
2327
2328	vector<Intersectable *> pvs(0);
2329	if (silhouetteCheckPercentage != 0.0f) {
2330	Intersectable::NewMail();
2331	pvs.reserve(leaf->GetPvsSize());
2332	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2333	ri != leaf->rays.end();
2334	ri++) {
2335	Intersectable object = (ri).GetObject();
2336	if (object) {
2337	if (!object->Mailed()) {
2338	object->Mail();
2339	object->mCounter = 1;
2340	pvs.push_back(object);
2341	if (exporter)
2342	exporter->ExportIntersectable(object);
2343	} else {
2344	object->mCounter++;
2345	}
2346	}
2347	}
2348	// sort objects based on their mCounter
2349	sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
2350	}
2351
2352	for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
2353	bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
2354	extendedConvexCombinationProb));
2355
2356	Vector3 origin, direction;
2357	// generate half of convex combination and half of random rays
2358	if (useExtendedConvexCombination) {
2359	// pickup 2 random rays
2360	int r1, r2;
2361	PickEdgeRays(leaf, r1, r2);
2362
2363
2364	Vector3 o1 = leaf->rays[r1].GetOrigin();
2365	Vector3 o2 = leaf->rays[r2].GetOrigin();
2366
2367	const float overlap = 0.0f;
2368	float w1, w2;
2369	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2370	origin = w1o1 + w2o2;
2371	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2372	direction =
2373	w1*leaf->rays[r1].GetDir() +
2374	w2*leaf->rays[r2].GetDir();
2375	} else {
2376	Vector3 dirVector;
2377
2378	Vector3 pVector;
2379	Vector3 dVector;
2380
2381	bool useHalton = true;
2382
2383	if (useHalton) {
2384	// generate a random 5D vector
2385	pVector = Vector3(leaf->halton.GetNumber(1),
2386	leaf->halton.GetNumber(2),
2387	leaf->halton.GetNumber(3));
2388
2389	dVector = Vector3(leaf->halton.GetNumber(4),
2390	leaf->halton.GetNumber(5),
2391	0.0f);
2392
2393	leaf->halton.GenerateNext();
2394	} else {
2395	pVector = Vector3(RandomValue(0.0f, 1.0f),
2396	RandomValue(0.0f, 1.0f),
2397	RandomValue(0.0f, 1.0f));
2398
2399	dVector = Vector3(RandomValue(0.0f, 1.0f),
2400	RandomValue(0.0f, 1.0f),
2401	0.0f);
2402	}
2403
2404	origin = box.GetPoint(pVector);
2405	dirVector = dirBox.GetPoint(dVector);
2406
2407	direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
2408	}
2409
2410	// shift the origin a little bit
2411	direction.Normalize();
2412
2413	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2414	float dist = 0.0f;
2415	float minT, maxT;
2416	// compute interection of the ray with the box
2417
2418	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2419	dist = maxT;
2420
2421
2422	origin += direction*dist;
2423
2424	bool intersects = false;
2425
2426	if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
2427	if (exporter) {
2428	VssRay ray = new VssRay(origin, origin + 100direction, NULL, NULL);
2429	selectedRays.push_back(ray);
2430	}
2431
2432	// check whether the ray does not intersect already visible objects
2433	Ray traversalRay;
2434	traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2435	for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2436	oi != pvs.end();
2437	oi++) {
2438	Intersectable object = oi;
2439	// do not test every object
2440	if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
2441	object->CastRay(traversalRay) ) {
2442	intersects = true;
2443	break;
2444	}
2445	}
2446	}
2447
2448
2449	if (!intersects) {
2450	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2451	rays.push_back(SimpleRay(origin, direction));
2452	generated++;
2453	}
2454	}
2455
2456	// cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
2457	// assign pdfs to the generated rays
2458	// float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
2459	float density = 1.0f;
2460
2461	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2462	for (int i=startIndex; i < rays.size(); i++) {
2463	rays[i].mPdf = p;
2464	}
2465
2466
2467	if (exporter) {
2468	exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2469	delete exporter;
2470	CLEAR_CONTAINER(selectedRays);
2471	}
2472
2473	}
2474
2475	int
2476	RssTree::PruneRays(RssTreeLeaf *leaf,
2477	const float contributionThreshold)
2478	{
2479	int i;
2480	int j;
2481
2482	for (j=0, i=0; i < leaf->rays.size(); i++) {
2483
2484	if (leaf->rays[i].mRay->mWeightedPvsContribution > contributionThreshold) {
2485	// copy a valid sample
2486	leaf->rays[j] = leaf->rays[i];
2487	j++;
2488	} else {
2489	// delete the ray
2490	leaf->rays[i].mRay->Unref();
2491	if (leaf->rays[i].mRay->RefCount() != 0) {
2492	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2493	exit(1);
2494	}
2495	delete leaf->rays[i].mRay;
2496	}
2497	}
2498
2499
2500	leaf->rays.resize(j);
2501	int removed = (i-j);
2502	stat.rayRefs -= removed;
2503	return removed;
2504	}
2505
2506	int
2507	RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2508	const float ratio)
2509	{
2510	int i;
2511	int j;
2512
2513	for (j=0, i=0; i < leaf->rays.size(); i++) {
2514
2515	if (Random(1.0f) < ratio) {
2516	// copy a valid sample
2517	if (i!=j)
2518	leaf->rays[j] = leaf->rays[i];
2519	j++;
2520	} else {
2521	// delete the ray
2522	leaf->rays[i].mRay->Unref();
2523	if (leaf->rays[i].mRay->RefCount() != 0) {
2524	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2525	exit(1);
2526	}
2527	delete leaf->rays[i].mRay;
2528	}
2529	}
2530
2531
2532	leaf->rays.resize(j);
2533	int removed = (i-j);
2534	stat.rayRefs -= removed;
2535	return removed;
2536	}
2537
2538	int
2539	RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2540	const float ratio)
2541	{
2542	int i;
2543
2544	if (leaf->rays.size() == 0)
2545	return 0;
2546
2547	sort(leaf->rays.begin(),
2548	leaf->rays.end(),
2549	RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2550
2551	int desired = int(ratio*leaf->rays.size());
2552	int removed = (int)leaf->rays.size() - desired;
2553
2554	for (i=desired; i < (int)leaf->rays.size(); i++) {
2555	// delete the ray
2556	leaf->rays[i].mRay->Unref();
2557	if (leaf->rays[i].mRay->RefCount() != 0) {
2558	cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
2559	cerr<<"desired="<<desired<<"i="<<i<<" size="<<(int)leaf->rays.size()<<endl;
2560	exit(1);
2561	}
2562	delete leaf->rays[i].mRay;
2563	}
2564
2565	leaf->rays.resize(desired);
2566	stat.rayRefs -= removed;
2567	return removed;
2568	}
2569
2570
2571	int
2572	RssTree::PruneRays(
2573	const int maxRays
2574	)
2575	{
2576	bool globalPrunning = false;
2577
2578	stack<RssTreeNode *> tstack;
2579	int prunned = 0;
2580
2581	// prune more rays to amortize the procedure
2582	int desired = (int)(maxRays * 0.8f);
2583
2584	Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;
2585
2586	if (globalPrunning) {
2587	VssRayContainer allRays;
2588	allRays.reserve(stat.rayRefs);
2589	CollectRays(allRays);
2590	sort(allRays.begin(),
2591	allRays.end(),
2592	GreaterWeightedPvsContribution);
2593
2594	if ( maxRays >= allRays.size() )
2595	return 0;
2596
2597	float contributionThreshold = allRays[desired]->mWeightedPvsContribution;
2598
2599	PushRoots(tstack);
2600
2601	while (!tstack.empty()) {
2602	RssTreeNode *node = tstack.top();
2603	tstack.pop();
2604
2605	if (node->IsLeaf()) {
2606	RssTreeLeaf leaf = (RssTreeLeaf )node;
2607	prunned += PruneRays(leaf, contributionThreshold);
2608	} else {
2609	RssTreeInterior in = (RssTreeInterior )node;
2610	// both nodes for directional splits
2611	tstack.push(in->front);
2612	tstack.push(in->back);
2613	}
2614	}
2615	} else {
2616	// prune random rays from each leaf so that the ratio's remain the same
2617	PushRoots(tstack);
2618
2619	if ( maxRays >= stat.rayRefs )
2620	return 0;
2621
2622	float ratio = desired/(float)stat.rayRefs;
2623
2624	while (!tstack.empty()) {
2625	RssTreeNode *node = tstack.top();
2626	tstack.pop();
2627
2628	if (node->IsLeaf()) {
2629	RssTreeLeaf leaf = (RssTreeLeaf )node;
2630	// prunned += PruneRaysRandom(leaf, ratio);
2631	prunned += PruneRaysContribution(leaf, ratio);
2632	} else {
2633	RssTreeInterior in = (RssTreeInterior )node;
2634	// both nodes for directional splits
2635	tstack.push(in->front);
2636	tstack.push(in->back);
2637	}
2638	}
2639
2640
2641
2642
2643	}
2644
2645
2646
2647	Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
2648
2649	return prunned;
2650	}
2651
2652	int
2653	RssTree::GenerateRays(const int numberOfRays,
2654	const int numberOfLeaves,
2655	SimpleRayContainer &rays)
2656	{
2657
2658
2659	vector<RssTreeLeaf *> leaves;
2660
2661	Debug<<"Collecting leaves..."<<endl<<flush;
2662	CollectLeaves(leaves);
2663	Debug<<"done."<<endl<<flush;
2664
2665
2666	if (numberOfLeaves != leaves.size()) {
2667	Debug<<"sorting leaves..."<<endl<<flush;
2668	sort(leaves.begin(),
2669	leaves.end(),
2670	GreaterContribution);
2671	Debug<<"done."<<endl<<flush;
2672	}
2673
2674
2675	Debug<<"Evaluating leaf importance..."<<endl<<flush;
2676
2677	float sumContrib = 0.0;
2678	int i;
2679	int k = 0;
2680	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2681	if (ValidLeaf(leaves[i])) {
2682	float c = leaves[i]->GetImportance();
2683	sumContrib += c;
2684	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
2685	k++;
2686	}
2687
2688	Debug<<"done."<<endl<<flush;
2689
2690	float avgContrib = sumContrib/numberOfLeaves;
2691
2692	// always generate at leat n ray per leaf
2693	int fixedPerLeaf = 1;
2694	// int fixedPerLeaf = 0;
2695	int fixed = fixedPerLeaf*(int)leaves.size();
2696
2697	// check if the number of fixed is not too big
2698
2699	if (fixed >= numberOfRays/2) {
2700	fixedPerLeaf = 0;
2701	fixed = 0;
2702	}
2703
2704	int iGenerated = numberOfRays - fixed;
2705	float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
2706
2707	k = 0;
2708
2709	Debug<<"generating leaf rays..."<<endl<<flush;
2710
2711	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2712	if (ValidLeaf(leaves[i])) {
2713	RssTreeLeaf *leaf = leaves[i];
2714	float c = leaf->GetImportance();
2715
2716	int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
2717	GenerateLeafRaysSimple(leaf, num, rays);
2718
2719	k++;
2720	}
2721
2722	Debug<<"done."<<endl<<flush;
2723
2724	return (int)rays.size();
2725	}
2726
2727
2728	float
2729	RssTree::GetAvgPvsSize()
2730	{
2731	stack<RssTreeNode *> tstack;
2732	PushRoots(tstack);
2733
2734	int sumPvs = 0;
2735	int leaves = 0;
2736	while (!tstack.empty()) {
2737	RssTreeNode *node = tstack.top();
2738	tstack.pop();
2739
2740	if (node->IsLeaf()) {
2741	RssTreeLeaf leaf = (RssTreeLeaf )node;
2742	// update pvs size
2743	UpdatePvsSize(leaf);
2744	sumPvs += leaf->GetPvsSize();
2745	leaves++;
2746	} else {
2747	RssTreeInterior in = (RssTreeInterior )node;
2748	// both nodes for directional splits
2749	tstack.push(in->front);
2750	tstack.push(in->back);
2751	}
2752	}
2753
2754
2755	return sumPvs/(float)leaves;
2756	}
2757
2758	float weightAbsContributions = ABS_CONTRIBUTION_WEIGHT;
2759	// if small very high importance of the last sample
2760	// if 1.0f then weighs = 1 1/2 1/3 1/4
2761	float passSampleWeightDecay = 1.0f;
2762	//float passSampleWeightDecay = 0.0001f;
2763
2764	float
2765	RssTree::GetSampleWeight(const int pass)
2766	{
2767	int passDiff = mCurrentPass - pass;
2768	float weight;
2769	if (1)
2770	weight = 1.0f/sqr(passDiff + passSampleWeightDecay);
2771	else
2772	switch (passDiff) {
2773	case 0:
2774	// weight = 1.0f;
2775	// break;
2776	default:
2777	weight = 1.0f;
2778	break;
2779	// case 1:
2780	// weight = 0.5f;
2781	// break;
2782	// case 2:
2783	// weight = 0.25f;
2784	// break;
2785	// case 3:
2786	// weight = 0.12f;
2787	// break;
2788	// case 4:
2789	// weight = 0.06f;
2790	// break;
2791	// default:
2792	// weight = 0.03f;
2793	// break;
2794	}
2795	return weight;
2796	}
2797
2798	void
2799	RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2800	float &weight,
2801	float &contribution)
2802	{
2803	VssRay *ray = info.mRay;
2804	weight = GetSampleWeight(mCurrentPass);
2805	contribution =
2806	weightAbsContributions*ray->mPvsContribution +
2807	(1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2808	// store the computed value
2809	info.mRay->mWeightedPvsContribution = weight*contribution;
2810	}
2811
2812
2813	void
2814	RssTree::ComputeImportance(RssTreeLeaf *leaf)
2815	{
2816	if (0)
2817	leaf->mImportance = leaf->GetAvgRayContribution();
2818	else {
2819	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2820	it_end = leaf->rays.end();
2821
2822	float sumContributions = 0.0f;
2823	float sumRelContributions = 0.0f;
2824	float sumWeights = 0.0f;
2825	for (; it != it_end; ++it) {
2826	VssRay ray = (it).mRay;
2827	float weight = GetSampleWeight(ray->mPass);
2828
2829	sumContributions += weight*ray->mPvsContribution;
2830	//$$ 20.7. changed to sqr to pronouce higher contribution so that they do not get smoothed!!
2831	//sumRelContributions += weight*ray->mRelativePvsContribution;
2832
2833	sumRelContributions += sqr(weight*ray->mRelativePvsContribution);
2834
2835	//sumWeights += weight;
2836	sumWeights += 1.0f;
2837	}
2838	// $$
2839	// sumWeights = leaf->mTotalRays;
2840
2841	if (sumWeights != 0.0f)
2842	leaf->mImportance =
2843	(weightAbsContributions*sumContributions +
2844	(1.0f - weightAbsContributions)*sumRelContributions)/sumWeights;
2845	else
2846	leaf->mImportance = 0.0f;
2847
2848	}
2849
2850	// return GetAvgRayContribution()*mImportance;
2851	//return GetAvgRayContribution();
2852	}
2853
2854
2855	float
2856	RssTreeLeaf::GetImportance() const
2857	{
2858	//return sqr(mImportance);
2859	return mImportance;
2860	}
2861
2862
2863	float
2864	RssTreeLeaf::ComputePvsEntropy() const
2865	{
2866	int samples = 0;
2867	Intersectable::NewMail();
2868	// set all object as belonging to the fron pvs
2869	for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
2870	ri != rays.end();
2871	ri++)
2872	if ((*ri).mRay->IsActive()) {
2873	Intersectable object = (ri).GetObject();
2874	if (object) {
2875	if (!object->Mailed()) {
2876	object->Mail();
2877	object->mCounter = 1;
2878	} else
2879	object->mCounter++;
2880	samples++;
2881	}
2882	}
2883
2884	float entropy = 0.0f;
2885
2886	if (samples > 1) {
2887	Intersectable::NewMail();
2888	for(RayInfoContainer::const_iterator ri = rays.begin();
2889	ri != rays.end();
2890	ri++)
2891	if ((*ri).mRay->IsActive()) {
2892	Intersectable object = (ri).GetObject();
2893	if (object) {
2894	if (!object->Mailed()) {
2895	object->Mail();
2896	float p = object->mCounter/(float)samples;
2897	entropy -= p*log(p);
2898	}
2899	}
2900	}
2901	entropy = entropy/log((float)samples);
2902	}
2903	else
2904	entropy = 1.0f;
2905
2906	return entropy;
2907	}
2908
2909	float
2910	RssTreeLeaf::ComputeRayLengthEntropy() const
2911	{
2912	// get sum of all ray lengths
2913	// consider only passing rays or originating rays
2914	float sum = 0.0f;
2915	int samples = 0;
2916	int i=0;
2917	for(RayInfoContainer::const_iterator ri = rays.begin();
2918	ri != rays.end();
2919	ri++)
2920	if ((*ri).mRay->IsActive()) {
2921	// float s;
2922	// if (i == 0)
2923	// s = 200;
2924	// else
2925	// s = 100;
2926	// i++;
2927
2928	sum += (*ri).mRay->GetSize();
2929	samples++;
2930	}
2931
2932
2933	float entropy = 0.0f;
2934
2935	if (samples > 1) {
2936	i = 0;
2937	for(RayInfoContainer::const_iterator ri = rays.begin();
2938	ri != rays.end();
2939	ri++)
2940	if ((*ri).mRay->IsActive()) {
2941	// float s;
2942	// if (i==0)
2943	// s = 200;
2944	// else
2945	// s = 100;
2946	// i++;
2947	// float p = s/sum;
2948	float p = (*ri).mRay->GetSize()/sum;
2949	entropy -= p*log(p);
2950	}
2951	entropy = entropy/log((float)samples);
2952	} else
2953	entropy = 1.0f;
2954
2955	return entropy;
2956	}
2957
2958
2959
2960	float
2961	RssTreeLeaf::ComputeEntropyImportance() const
2962	{
2963
2964	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2965	return 1.0f - ComputePvsEntropy();
2966	}
2967
2968	float
2969	RssTreeLeaf::ComputeRayContributionImportance() const
2970	{
2971	// mPvsEntropy = ComputePvsEntropy();
2972	// mRayLengthEntropy = ComputeRayLengthEntropy();
2973
2974	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2975	return 1.0f - ComputePvsEntropy();
2976	}
2977
2978
2979	AxisAlignedBox3
2980	RssTree::GetShrankedBBox(const RssTreeNode *node) const
2981	{
2982	if (node->parent == NULL)
2983	return bbox;
2984
2985	if (!node->IsLeaf())
2986	return ((RssTreeInterior *)node)->bbox;
2987
2988	// evaluate bounding box from the ray origins
2989	AxisAlignedBox3 box;
2990	box.Initialize();
2991	RssTreeLeaf leaf = (RssTreeLeaf )node;
2992	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2993	ri != leaf->rays.end();
2994	ri++)
2995	if ((*ri).mRay->IsActive()) {
2996	box.Include((*ri).GetOrigin());
2997	}
2998	return box;
2999	}
3000
3001
3002	int
3003	RssTree::CollectRays(VssRayContainer &rays,
3004	const int number)
3005	{
3006	VssRayContainer allRays;
3007	CollectRays(allRays);
3008
3009
3010	int desired = min(number, (int)allRays.size());
3011	float prob = desired/(float)allRays.size();
3012	//int skip = allRays.size()/desired;
3013	while (rays.size() < desired) {
3014	VssRayContainer::const_iterator it = allRays.begin();
3015	for (; it != allRays.end() && rays.size() < desired; it++) {
3016	if (Random(1.0f) < prob)
3017	rays.push_back(*it);
3018	}
3019	}
3020	return (int)rays.size();
3021	}
3022
3023
3024	int
3025	RssTree::CollectRays(VssRayContainer &rays
3026	)
3027	{
3028	VssRay::NewMail();
3029
3030	stack<RssTreeNode *> tstack;
3031	PushRoots(tstack);
3032
3033	while (!tstack.empty()) {
3034	RssTreeNode *node = tstack.top();
3035	tstack.pop();
3036
3037	if (node->IsLeaf()) {
3038	RssTreeLeaf leaf = (RssTreeLeaf )node;
3039	// update pvs size
3040	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
3041	for (;it != leaf->rays.end(); ++it)
3042	if (!(*it).mRay->Mailed()) {
3043	(*it).mRay->Mail();
3044	rays.push_back((*it).mRay);
3045	}
3046	} else {
3047	RssTreeInterior in = (RssTreeInterior )node;
3048	// both nodes for directional splits
3049	tstack.push(in->front);
3050	tstack.push(in->back);
3051	}
3052	}
3053
3054	return (int)rays.size();
3055	}
3056
3057
3058	RssTreeNode *
3059	RssTree::GetRoot(Intersectable *object) const
3060	{
3061	if (mPerObjectTree && object) {
3062	int id = object->GetId()-1;
3063	if (id < 0 \|\| id >= mRoots.size()) {
3064	Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<(int)mRoots.size()<<
3065	endl<<flush;
3066	id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
3067	}
3068	return mRoots[id];
3069	} else
3070	return mRoots[0];
3071	}
3072
3073	void
3074	RssTree::PushRoots(priority_queue<TraversalData> &st) const
3075	{
3076	for (int i=0; i < mRoots.size(); i++)
3077	st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
3078	}
3079
3080	void
3081	RssTree::PushRoots(stack<RssTreeNode *> &st) const
3082	{
3083	for (int i=0; i < mRoots.size(); i++)
3084	st.push(mRoots[i]);
3085	}
3086
3087	void
3088	RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info) const
3089	{
3090	for (int i=0; i < mRoots.size(); i++)
3091	st.push(RayTraversalData(mRoots[i], info));
3092	}
3093
3094
3095	}

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

Download in other formats: