Context Navigation

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

Revision 1877, 73.7 KB checked in by bittner, 18 years ago (diff)
sampling updates

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, 0.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, 0.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	// if ( forcedBoundingBox )
395	// bbox = *forcedBoundingBox;
396
397	cout<<"Bbox = "<<bbox<<endl;
398	cout<<"Dirr Bbox = "<<dirBBox<<endl;
399
400	stat.rays = rays.size();
401	stat.initialPvsSize = 0;
402	for (int i=0; i < mRoots.size(); i++) {
403	RssTreeLeaf leaf = (RssTreeLeaf )mRoots[i];
404	UpdatePvsSize(leaf);
405	stat.initialPvsSize += leaf->GetPvsSize();
406	mRoots[i] = Subdivide(TraversalData(leaf, GetBBox(leaf), 0));
407	}
408
409	if (splitCandidates) {
410	// force realease of this vector
411	delete splitCandidates;
412	splitCandidates = new vector<SortableEntry>;
413	}
414
415	stat.Stop();
416	stat.Print(cout);
417	cout<<"#Total memory="<<GetMemUsage()<<endl;
418
419	// this rotine also updates importances etc...
420	}
421
422	int
423	RssTree::UpdateSubdivision()
424	{
425	priority_queue<TraversalData> tStack;
426	// stack<TraversalData> tStack;
427
428	// tStack.push(TraversalData(root, bbox, 0));
429	PushRoots(tStack);
430
431	AxisAlignedBox3 backBox;
432	AxisAlignedBox3 frontBox;
433
434	maxMemory = maxTotalMemory;
435	int subdivided = 0;
436	int lastMem = 0;
437	while (!tStack.empty()) {
438
439	float mem = GetMemUsage();
440
441	if ( lastMem/10 != ((int)mem)/10) {
442	cout<<mem<<" MB"<<endl;
443	}
444	lastMem = (int)mem;
445
446	if ( mem > maxMemory ) {
447	// count statistics on unprocessed leafs
448	while (!tStack.empty()) {
449	// EvaluateLeafStats(tStack.top());
450	tStack.pop();
451	}
452	break;
453	}
454
455	TraversalData data = tStack.top();
456	tStack.pop();
457
458	if (data.node->IsLeaf()) {
459	RssTreeNode node = SubdivideNode((RssTreeLeaf ) data.node,
460	data.bbox,
461	backBox,
462	frontBox
463	);
464	if (!node->IsLeaf()) {
465	subdivided++;
466
467
468	RssTreeInterior interior = (RssTreeInterior ) node;
469	// push the children on the stack
470	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
471	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
472	} else {
473	// EvaluateLeafStats(data);
474	}
475	} else {
476	RssTreeInterior interior = (RssTreeInterior ) data.node;
477	tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
478	tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
479	}
480	}
481	return subdivided;
482	}
483
484
485	RssTreeNode *
486	RssTree::Subdivide(const TraversalData &tdata)
487	{
488	RssTreeNode *result = NULL;
489
490	priority_queue<TraversalData> tStack;
491	// stack<TraversalData> tStack;
492
493	tStack.push(tdata);
494
495	AxisAlignedBox3 backBox;
496	AxisAlignedBox3 frontBox;
497
498
499	int lastMem = 0;
500	while (!tStack.empty()) {
501
502	float mem = GetMemUsage();
503
504	if ( lastMem/10 != ((int)mem)/10) {
505	cout<<mem<<" MB"<<endl;
506	}
507	lastMem = (int)mem;
508
509	if ( mem > maxMemory ) {
510	// count statistics on unprocessed leafs
511	while (!tStack.empty()) {
512	EvaluateLeafStats(tStack.top());
513	tStack.pop();
514	}
515	break;
516	}
517
518	TraversalData data = tStack.top();
519	tStack.pop();
520
521	#if DEBUG_SPLITS
522	Debug<<"#Splitting node"<<endl;
523	data.node->Print(Debug);
524	#endif
525	RssTreeNode node = SubdivideNode((RssTreeLeaf ) data.node,
526	data.bbox,
527	backBox,
528	frontBox
529	);
530
531
532	if (result == NULL)
533	result = node;
534
535	if (!node->IsLeaf()) {
536
537	RssTreeInterior interior = (RssTreeInterior ) node;
538	// push the children on the stack
539	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
540	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
541
542	#if DEBUG_SPLITS
543	Debug<<"#New nodes"<<endl;
544	interior->back->Print(Debug);
545	interior->front->Print(Debug);
546	Debug<<"#####################################"<<endl;
547	#endif
548
549	} else {
550	EvaluateLeafStats(data);
551	}
552	}
553
554	return result;
555	}
556
557
558	// returns selected plane for subdivision
559	int
560	RssTree::SelectPlane(
561	RssTreeLeaf *leaf,
562	const AxisAlignedBox3 &box,
563	SplitInfo &info
564	)
565	{
566
567	BestCostRatio(leaf,
568	info);
569
570	#if DEBUG_SPLIT_COST
571	Debug<<"Split Info:"<<endl;
572	Debug<<"axis="<<info.axis<<" ratio="<<info.costRatio<<endl;
573	Debug<<"rays="<<info.rays<<
574	" rays back="<<info.raysBack<<
575	" rays front="<<info.raysFront<<endl;
576	Debug<<"c="<<info.contribution<<
577	" c back="<<info.contributionBack<<
578	" c front="<<info.contributionFront<<endl;
579
580	// Debug<<"viewcells="<<info.viewCells<<
581	// " viewcells back="<<info.viewCellsBack<<
582	// " viewcells back="<<info.viewCellsFront<<endl;
583	#endif
584
585	if (info.costRatio > termMaxCostRatio) {
586	// cout<<"Too big cost ratio "<<costRatio<<endl;
587	stat.maxCostRatioNodes++;
588	return -1;
589	}
590
591	#if 0
592	cout<<
593	"pvs="<<leaf->mPvsSize<<
594	" rays="<<leaf->rays.size()<<
595	" rc="<<leaf->GetAvgRayContribution()<<
596	" axis="<<info.axis<<endl;
597	#endif
598
599	return info.axis;
600	}
601
602
603	void
604	RssTree::GetCostRatio(
605	RssTreeLeaf *leaf,
606	SplitInfo &info
607	)
608	{
609
610	AxisAlignedBox3 box;
611	float minBox, maxBox;
612
613	if (info.axis < 3) {
614	box = GetBBox(leaf);
615	minBox = box.Min(info.axis);
616	maxBox = box.Max(info.axis);
617	} else {
618	box = GetDirBBox(leaf);
619	minBox = box.Min(info.axis-3);
620	maxBox = box.Max(info.axis-3);
621	}
622
623	float sizeBox = maxBox - minBox;
624
625	int pvsSize = leaf->GetPvsSize();
626
627	if (!mImportanceBasedCost) {
628	const int costMethod = 0;
629
630	switch (costMethod) {
631	case 0: {
632	float sum = info.pvsBack(info.position - minBox) + info.pvsFront(maxBox - info.position);
633	float newCost = ct_div_ci + sum/sizeBox;
634	float oldCost = (float)pvsSize;
635	info.costRatio = newCost/oldCost;
636	break;
637	}
638	case 6: {
639	float sum = info.raysBack(info.position - minBox) + info.raysFront(maxBox - info.position);
640	float newCost = ct_div_ci + sum/sizeBox;
641	float oldCost = (float)info.rays;
642	info.costRatio = newCost/oldCost;
643	break;
644	}
645	case 7: {
646	float benefit =
647	(sqr(info.contributionBack) + sqr(info.contributionFront))/
648	(2.0f*sqr(info.contribution));
649	info.costRatio = 1.0f/(benefit + Limits::Small);
650	break;
651	}
652
653
654	case 1: {
655	float newContrib =
656	(info.contributionBack*(info.position - minBox) +
657	+
658	info.contributionFront*(maxBox - info.position))/sizeBox;
659	float oldContrib = info.contribution;
660	info.costRatio = newContrib/oldContrib;
661	// cout<<info.contribution<<" "<<info.contributionBack<<" "<<info.contributionFront<<endl;
662	break;
663	}
664	case 2: {
665	float sum =
666	info.viewCellsBack*(info.position - minBox) +
667	info.viewCellsFront*(maxBox - info.position);
668	float newCost = ct_div_ci + sum/sizeBox;
669	float oldCost = (float)info.viewCells;
670	info.costRatio = newCost/oldCost;
671	break;
672	}
673	case 3: {
674	float newCost = (float)(info.raysBackinfo.pvsBack + info.raysFrontinfo.pvsFront);
675	float oldCost = (float)leaf->rays.size()*pvsSize;
676	info.costRatio = newCost/oldCost;
677	break;
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
1241	node->bbox = GetBBox(leaf);
1242	node->dirBBox = GetDirBBox(leaf);
1243
1244	RssTreeLeaf *back = new RssTreeLeaf(node, info.raysBack);
1245	RssTreeLeaf *front = new RssTreeLeaf(node, info.raysFront);
1246
1247
1248
1249	// replace a link from node's parent
1250	if ( leaf->parent )
1251	leaf->parent->ReplaceChildLink(leaf, node);
1252	// and setup child links
1253	node->SetupChildLinks(back, front);
1254
1255	back->bbox = leaf->bbox;
1256	front->bbox = leaf->bbox;
1257	back->dirBBox = leaf->dirBBox;
1258	front->dirBBox = leaf->dirBBox;
1259
1260	if (axis <= RssTreeNode::SPLIT_Z) {
1261	back->bbox.SetMax(axis, info.position);
1262	front->bbox.SetMin(axis, info.position);
1263	} else {
1264	back->dirBBox.SetMax(axis-3, info.position);
1265	front->dirBBox.SetMin(axis-3, info.position);
1266	}
1267
1268	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1269	ri != leaf->rays.end();
1270	ri++) {
1271	if ((*ri).mRay->IsActive()) {
1272
1273	// first unref ray from the former leaf
1274	(*ri).mRay->Unref();
1275
1276	// Debug << "computed t: " << (*ri).mRay->mT << endl;
1277	// determine the side of this ray with respect to the plane
1278	int side = node->ComputeRaySide(*ri);
1279
1280	if (side == 1)
1281	front->AddRay(*ri);
1282	else
1283	back->AddRay(*ri);
1284	} else
1285	(*ri).mRay->Unref();
1286	}
1287
1288	// distribute the total number of rays according to the distribution
1289	// of rays which remained
1290	// front->mTotalRays = front->rays.size()*leaf->mTotalRays/leaf->rays.size();
1291	// back->mTotalRays = back->rays.size()*leaf->mTotalRays/leaf->rays.size();
1292
1293	#if 0
1294	front->SetPvsSize(pvsFront);
1295	back->SetPvsSize(pvsBack);
1296	// compute entropy as well
1297	front->ComputeEntropyImportance();
1298	back->ComputeEntropyImportance();
1299	#else
1300	UpdatePvsSize(front);
1301	UpdatePvsSize(back);
1302	#endif
1303
1304
1305	// update stats
1306	stat.rayRefs -= (int)leaf->rays.size();
1307	stat.rayRefs += back->rays.size() + front->rays.size();
1308
1309
1310	delete leaf;
1311	return node;
1312	}
1313
1314
1315
1316
1317
1318
1319	int
1320	RssTree::ReleaseMemory(const int time)
1321	{
1322	stack<RssTreeNode *> tstack;
1323
1324	// find a node in the tree which subtree will be collapsed
1325	int maxAccessTime = time - accessTimeThreshold;
1326	int released;
1327
1328	PushRoots(tstack);
1329
1330	while (!tstack.empty()) {
1331	RssTreeNode *node = tstack.top();
1332	tstack.pop();
1333
1334
1335	if (!node->IsLeaf()) {
1336	RssTreeInterior in = (RssTreeInterior )node;
1337	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1338	if (in->depth >= minCollapseDepth &&
1339	in->lastAccessTime <= maxAccessTime) {
1340	released = CollapseSubtree(node, time);
1341	break;
1342	}
1343
1344	if (in->back->GetAccessTime() <
1345	in->front->GetAccessTime()) {
1346	tstack.push(in->front);
1347	tstack.push(in->back);
1348	} else {
1349	tstack.push(in->back);
1350	tstack.push(in->front);
1351	}
1352	}
1353	}
1354
1355	while (tstack.empty()) {
1356	// could find node to collaps...
1357	// cout<<"Could not find a node to release "<<endl;
1358	break;
1359	}
1360
1361	return released;
1362	}
1363
1364
1365
1366
1367	RssTreeNode *
1368	RssTree::SubdivideLeaf(
1369	RssTreeLeaf *leaf
1370	)
1371	{
1372	RssTreeNode *node = leaf;
1373
1374	AxisAlignedBox3 leafBBox = GetBBox(leaf);
1375
1376	static int pass = 0;
1377	pass ++;
1378
1379	// check if we should perform a dynamic subdivision of the leaf
1380	if (!TerminationCriteriaSatisfied(leaf)) {
1381
1382	// memory check and realese...
1383	if (GetMemUsage() > maxTotalMemory) {
1384	ReleaseMemory( pass );
1385	}
1386
1387	AxisAlignedBox3 backBBox, frontBBox;
1388
1389	// subdivide the node
1390	node =
1391	SubdivideNode(leaf,
1392	leafBBox,
1393	backBBox,
1394	frontBBox
1395	);
1396	}
1397
1398	return node;
1399	}
1400
1401
1402
1403	void
1404	RssTree::UpdateRays(
1405	VssRayContainer &remove,
1406	VssRayContainer &add
1407	)
1408	{
1409	RssTreeLeaf::NewMail();
1410
1411	// schedule rays for removal
1412	for(VssRayContainer::const_iterator ri = remove.begin();
1413	ri != remove.end();
1414	ri++) {
1415	(*ri)->ScheduleForRemoval();
1416	}
1417
1418	int inactive=0;
1419	Debug<<"Deleting rays..."<<endl<<flush;
1420
1421	for(VssRayContainer::const_iterator ri = remove.begin();
1422	ri != remove.end();
1423	ri++) {
1424	if ((*ri)->ScheduledForRemoval())
1425	// RemoveRay(*ri, NULL, false);
1426	// !!! BUG - with true it does not work correctly - aggreated delete
1427	RemoveRay(*ri, NULL, true);
1428	else
1429	inactive++;
1430	}
1431
1432	Debug<<"done."<<endl<<flush;
1433
1434	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1435	Debug<<"Adding rays..."<<endl<<flush;
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	Debug<<"done."<<endl<<flush;
1446
1447	stat.rayRefs += add.size() - remove.size();
1448
1449	Debug<<"Updating statistics..."<<endl<<flush;
1450	UpdateTreeStatistics();
1451	Debug<<"done."<<endl<<flush;
1452
1453	// check whether the tree should be prunned
1454	if (stat.rayRefs > mMaxRays) {
1455	Debug<<"Prunning rays..."<<endl<<flush;
1456	PruneRays(mMaxRays);
1457	Debug<<"done."<<endl<<flush;
1458	// UpdateTreeStatistics();
1459	}
1460
1461	}
1462
1463
1464
1465
1466	void
1467	RssTree::RemoveRay(VssRay *ray,
1468	vector<RssTreeLeaf > affectedLeaves,
1469	const bool removeAllScheduledRays
1470	)
1471	{
1472
1473	stack<RayTraversalData> tstack;
1474
1475	PushRoots(tstack, RssTreeLeaf::RayInfo(ray));
1476
1477	RayTraversalData data;
1478
1479	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1480
1481	while (!tstack.empty()) {
1482	data = tstack.top();
1483	tstack.pop();
1484
1485	if (!data.node->IsLeaf()) {
1486	// split the set of rays in two groups intersecting the
1487	// two subtrees
1488
1489	TraverseInternalNode(data, tstack);
1490
1491	} else {
1492	// remove the ray from the leaf
1493	// find the ray in the leaf and swap it with the last ray...
1494	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1495
1496	if (!leaf->Mailed()) {
1497	leaf->Mail();
1498	if (affectedLeaves)
1499	affectedLeaves->push_back(leaf);
1500
1501	if (removeAllScheduledRays) {
1502	int tail = (int)leaf->rays.size()-1;
1503
1504	for (int i=0; i < (int)(leaf->rays.size()); i++) {
1505	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1506	// find a ray to replace it with
1507	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1508	stat.removedRayRefs++;
1509	leaf->rays[tail].mRay->Unref();
1510	leaf->rays.pop_back();
1511	tail--;
1512	}
1513
1514	if (tail < i)
1515	break;
1516
1517	stat.removedRayRefs++;
1518	leaf->rays[i].mRay->Unref();
1519	leaf->rays[i] = leaf->rays[tail];
1520	leaf->rays.pop_back();
1521	tail--;
1522	}
1523	}
1524	}
1525	}
1526
1527	if (!removeAllScheduledRays)
1528	for (int i=0; i < (int)leaf->rays.size(); i++) {
1529	if (leaf->rays[i].mRay == ray) {
1530	stat.removedRayRefs++;
1531	ray->Unref();
1532	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1533	leaf->rays.pop_back();
1534	// check this ray again
1535	break;
1536	}
1537	}
1538
1539	}
1540	}
1541
1542	if (ray->RefCount() != 0) {
1543	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1544	exit(1);
1545	}
1546
1547	}
1548
1549
1550	void
1551	RssTree::AddRay(RssTreeNode::RayInfo &info)
1552	{
1553
1554	stack<RayTraversalData> tstack;
1555
1556	RssTreeNode *root = GetRoot(info.GetSourceObject());
1557	tstack.push(RayTraversalData(root, info));
1558
1559	RayTraversalData data;
1560
1561	while (!tstack.empty()) {
1562	data = tstack.top();
1563	tstack.pop();
1564
1565	if (!data.node->IsLeaf()) {
1566	TraverseInternalNode(data, tstack);
1567	} else {
1568	// remove the ray from the leaf
1569	// find the ray in the leaf and swap it with the last ray...
1570	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1571	leaf->AddRay(data.rayData);
1572	stat.addedRayRefs++;
1573	}
1574	}
1575	}
1576
1577	void
1578	RssTree::TraverseInternalNode(
1579	RayTraversalData &data,
1580	stack<RayTraversalData> &tstack)
1581	{
1582	RssTreeInterior in = (RssTreeInterior ) data.node;
1583
1584	// determine the side of this ray with respect to the plane
1585	int side = in->ComputeRaySide(data.rayData
1586	);
1587
1588	if (side == 1)
1589	tstack.push(RayTraversalData(in->front, data.rayData));
1590	else
1591	tstack.push(RayTraversalData(in->back, data.rayData));
1592
1593	}
1594
1595
1596	int
1597	RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
1598	{
1599	// first count all rays in the subtree
1600	// use mail 1 for this purpose
1601	stack<RssTreeNode *> tstack;
1602	int rayCount = 0;
1603	int totalRayCount = 0;
1604	int collapsedNodes = 0;
1605
1606	#if DEBUG_COLLAPSE
1607	cout<<"Collapsing subtree"<<endl;
1608	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1609	cout<<"depth="<<(int)sroot->depth<<endl;
1610	#endif
1611
1612	// tstat.collapsedSubtrees++;
1613	// tstat.collapseDepths += (int)sroot->depth;
1614	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1615
1616	tstack.push(sroot);
1617	VssRay::NewMail();
1618
1619	while (!tstack.empty()) {
1620	collapsedNodes++;
1621	RssTreeNode *node = tstack.top();
1622	tstack.pop();
1623
1624	if (node->IsLeaf()) {
1625	RssTreeLeaf leaf = (RssTreeLeaf ) node;
1626	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1627	ri != leaf->rays.end();
1628	ri++) {
1629
1630	totalRayCount++;
1631	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
1632	(*ri).mRay->Mail();
1633	rayCount++;
1634	}
1635	}
1636	} else {
1637	tstack.push(((RssTreeInterior *)node)->back);
1638	tstack.push(((RssTreeInterior *)node)->front);
1639	}
1640	}
1641
1642	VssRay::NewMail();
1643
1644	// create a new node that will hold the rays
1645	RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
1646	if ( newLeaf->parent )
1647	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1648
1649
1650	tstack.push( sroot );
1651
1652	while (!tstack.empty()) {
1653
1654	RssTreeNode *node = tstack.top();
1655	tstack.pop();
1656
1657	if (node->IsLeaf()) {
1658	RssTreeLeaf leaf = (RssTreeLeaf ) node;
1659
1660	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1661	ri != leaf->rays.end();
1662	ri++) {
1663
1664	// unref this ray from the old node
1665
1666	if ((*ri).mRay->IsActive()) {
1667	(*ri).mRay->Unref();
1668	if (!(*ri).mRay->Mailed()) {
1669	(*ri).mRay->Mail();
1670	newLeaf->AddRay(*ri);
1671	}
1672	} else
1673	(*ri).mRay->Unref();
1674
1675	}
1676	} else {
1677	tstack.push(((RssTreeInterior *)node)->back);
1678	tstack.push(((RssTreeInterior *)node)->front);
1679	}
1680	}
1681
1682	// delete the node and all its children
1683	delete sroot;
1684
1685	// for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1686	// ri != newLeaf->rays.end();
1687	// ri++)
1688	// (*ri).ray->UnMail(2);
1689
1690
1691	#if DEBUG_COLLAPSE
1692	cout<<"Total memory before="<<GetMemUsage()<<endl;
1693	#endif
1694
1695	stat.nodes -= collapsedNodes - 1;
1696	stat.leaves -= collapsedNodes/2 - 1;
1697	stat.rayRefs -= totalRayCount - rayCount;
1698
1699	#if DEBUG_COLLAPSE
1700	cout<<"collapsed nodes"<<collapsedNodes<<endl;
1701	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1702	cout<<"Total memory after="<<GetMemUsage()<<endl;
1703	cout<<"================================"<<endl;
1704	#endif
1705
1706	// tstat.collapsedNodes += collapsedNodes;
1707	// tstat.collapsedRays += totalRayCount - rayCount;
1708
1709	return totalRayCount - rayCount;
1710	}
1711
1712
1713	int
1714	RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
1715	{
1716	stack<RssTreeNode *> tstack;
1717	PushRoots(tstack);
1718
1719	Intersectable::NewMail();
1720	int pvsSize = 0;
1721
1722	while (!tstack.empty()) {
1723	RssTreeNode *node = tstack.top();
1724	tstack.pop();
1725
1726
1727	if (node->IsLeaf()) {
1728	RssTreeLeaf leaf = (RssTreeLeaf )node;
1729	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1730	ri != leaf->rays.end();
1731	ri++)
1732	if ((*ri).mRay->IsActive()) {
1733	Intersectable *object;
1734	object = (*ri).GetObject();
1735	if (object && !object->Mailed()) {
1736	pvsSize++;
1737	object->Mail();
1738	}
1739	}
1740	} else {
1741	RssTreeInterior in = (RssTreeInterior )node;
1742	if (in->axis < 3) {
1743	if (box.Max(in->axis) >= in->position )
1744	tstack.push(in->front);
1745
1746	if (box.Min(in->axis) <= in->position )
1747	tstack.push(in->back);
1748	} else {
1749	// both nodes for directional splits
1750	tstack.push(in->front);
1751	tstack.push(in->back);
1752	}
1753	}
1754	}
1755	return pvsSize;
1756	}
1757
1758	int
1759	RssTree::CollectPvs(const AxisAlignedBox3 &box,
1760	ObjectContainer &pvs
1761	) const
1762	{
1763	stack<RssTreeNode *> tstack;
1764	PushRoots(tstack);
1765
1766	Intersectable::NewMail();
1767
1768	while (!tstack.empty()) {
1769	RssTreeNode *node = tstack.top();
1770	tstack.pop();
1771
1772
1773	if (node->IsLeaf()) {
1774	RssTreeLeaf leaf = (RssTreeLeaf )node;
1775	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1776	ri != leaf->rays.end();
1777	ri++)
1778	if ((*ri).mRay->IsActive()) {
1779	Intersectable *object;
1780	object = (*ri).GetObject();
1781	if (object && !object->Mailed()) {
1782	pvs.push_back(object);
1783	object->Mail();
1784	}
1785	}
1786	} else {
1787	RssTreeInterior in = (RssTreeInterior )node;
1788	if (in->axis < 3) {
1789	if (box.Max(in->axis) >= in->position )
1790	tstack.push(in->front);
1791
1792	if (box.Min(in->axis) <= in->position )
1793	tstack.push(in->back);
1794	} else {
1795	// both nodes for directional splits
1796	tstack.push(in->front);
1797	tstack.push(in->back);
1798	}
1799	}
1800	}
1801	return pvs.size();
1802	}
1803
1804	void
1805	RssTree::UpdateTreeStatistics()
1806	{
1807	stack<RssTreeNode *> tstack;
1808
1809	float sumPvsSize = 0.0f;
1810	float sumRayContribution = 0.0f;
1811	float sumWeightedRayContribution = 0.0f;
1812	float sumImportance = 0.0f;
1813	float sumPvsEntropy = 0.0f;
1814	float sumRayLengthEntropy = 0.0f;
1815	float sumRays = 0.0f;
1816
1817	int leaves = 0;
1818
1819	for (int i=0; i < mRoots.size(); i++)
1820	UpdateImportance(mRoots[i]);
1821
1822	PushRoots(tstack);
1823
1824	while (!tstack.empty()) {
1825	RssTreeNode *node = tstack.top();
1826	tstack.pop();
1827
1828	if (node->IsLeaf()) {
1829	leaves++;
1830	RssTreeLeaf leaf = (RssTreeLeaf )node;
1831	// updated above already
1832	// UpdatePvsSize(leaf);
1833
1834	sumPvsSize += leaf->GetPvsSize();
1835	sumRayContribution += leaf->GetAvgRayContribution();
1836
1837
1838	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
1839	for (;it != leaf->rays.end(); ++it) {
1840	float weight, contribution;
1841	GetRayContribution(*it, weight, contribution);
1842	sumWeightedRayContribution += weight*contribution;
1843	}
1844
1845	// sumPvsEntropy += leaf->mPvsEntropy;
1846	// sumRayLengthEntropy += leaf->mRayLengthEntropy;
1847	sumRays += leaf->rays.size();
1848
1849	float imp = leaf->GetImportance();
1850
1851	// if (imp > 1.0f)
1852	// cout<<"warning imp > 1.0f:"<<imp<<endl;
1853
1854	sumImportance += imp;
1855
1856	} else {
1857	RssTreeInterior in = (RssTreeInterior )node;
1858	tstack.push(in->front);
1859	tstack.push(in->back);
1860	}
1861	}
1862
1863	stat.avgPvsSize = sumPvsSize/(float)leaves;
1864	stat.avgRays = sumRays/(float)leaves;
1865	stat.avgRayContribution = sumRayContribution/(float)leaves;
1866	// avgPvsEntropy = sumPvsEntropy/(float)leaves;
1867	// avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
1868	stat.avgImportance = sumImportance/(float)leaves;
1869	stat.avgWeightedRayContribution = sumWeightedRayContribution/(float)sumRays;
1870	stat.rayRefs = (int)sumRays;
1871
1872	}
1873
1874	void
1875	RssTree::UpdateImportance(
1876	RssTreeNode *node)
1877	{
1878	if (node->IsLeaf()) {
1879	UpdatePvsSize((RssTreeLeaf *)node);
1880	} else {
1881	RssTreeInterior in = (RssTreeInterior )node;
1882	UpdateImportance(in->front);
1883	UpdateImportance(in->back);
1884	node->mImportance = in->front->mImportance + in->back->mImportance;
1885	}
1886	}
1887
1888	// generate a single ray described by parameters in a unit cube
1889	void
1890	RssTree::GenerateRay(float *params, SimpleRayContainer &rays)
1891	{
1892	RssTreeNode *node;
1893
1894	// works only with single root now!
1895	node = mRoots[0];
1896	while (!node->IsLeaf()) {
1897	RssTreeInterior in = (RssTreeInterior )node;
1898	// decide whether to go left or right
1899	int axis = in->axis;
1900	// float p = in->GetRelativePosition();
1901	float p = in->back->GetImportance()/(in->back->GetImportance() + in->front->GetImportance());
1902
1903	// cout<<"p = "<<p<<endl;
1904	// cout<<"axis="<<axis<<" pos="<<pos<<endl;
1905	// cout<<GetBBox(in)<<endl;
1906	// cout<<GetDirBBox(in)<<endl;
1907	// cout<<"********"<<endl;
1908
1909	if (params[axis] < p) {
1910	node = in->back;
1911	if (p > Limits::Small)
1912	params[axis] = params[axis] / p;
1913	} else {
1914	node = in->front;
1915	if (p < 1.0f - Limits::Small)
1916	params[axis] = (params[axis] - p) / (1.0f - p);
1917	}
1918	}
1919
1920	// cout<<params[0]<<" "<<params[1]<<" "<<params[2]<<" "<<params[3]<<" "<<params[4]<<endl;
1921
1922	GenerateLeafRay(params, rays, node->bbox, node->dirBBox);
1923
1924
1925	}
1926
1927	int
1928	RssTree::GenerateRays(
1929	const float ratioPerLeaf,
1930	SimpleRayContainer &rays)
1931	{
1932
1933	stack<RssTreeNode *> tstack;
1934	PushRoots(tstack);
1935
1936	while (!tstack.empty()) {
1937	RssTreeNode *node = tstack.top();
1938	tstack.pop();
1939
1940	if (node->IsLeaf()) {
1941	RssTreeLeaf leaf = (RssTreeLeaf )node;
1942	float c = leaf->GetImportance();
1943	int num = int(c*ratioPerLeaf + 0.5f);
1944	// cout<<num<<" ";
1945
1946	for (int i=0; i < num; i++) {
1947	Vector3 origin = GetBBox(leaf).GetRandomPoint();
1948	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1949	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1950	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1951	rays.push_back(SimpleRay(origin, direction));
1952	}
1953	} else {
1954	RssTreeInterior in = (RssTreeInterior )node;
1955	// both nodes for directional splits
1956	tstack.push(in->front);
1957	tstack.push(in->back);
1958	}
1959	}
1960
1961	return (int)rays.size();
1962	}
1963
1964	void
1965	RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
1966	{
1967	stack<RssTreeNode *> tstack;
1968	PushRoots(tstack);
1969
1970	while (!tstack.empty()) {
1971	RssTreeNode *node = tstack.top();
1972	tstack.pop();
1973
1974	if (node->IsLeaf()) {
1975	RssTreeLeaf leaf = (RssTreeLeaf )node;
1976	leaves.push_back(leaf);
1977	} else {
1978	RssTreeInterior in = (RssTreeInterior )node;
1979	// both nodes for directional splits
1980	tstack.push(in->front);
1981	tstack.push(in->back);
1982	}
1983	}
1984	}
1985
1986	bool
1987	RssTree::ValidLeaf(RssTreeLeaf *leaf) const
1988	{
1989	return true;
1990	//return leaf->rays.size() > termMinRays/4;
1991	}
1992
1993
1994	void
1995	RssTree::PickEdgeRays(RssTreeLeaf *leaf,
1996	int &r1,
1997	int &r2
1998	)
1999	{
2000	int nrays = leaf->rays.size();
2001
2002	if (nrays == 2) {
2003	r1 = 0;
2004	r2 = 1;
2005	return;
2006	}
2007
2008	#if 1
2009	int tries = min(20, nrays);
2010
2011	while (--tries) {
2012	r1 = Random(nrays);
2013	r2 = Random(nrays);
2014	if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
2015	break;
2016	}
2017
2018	if (r1 == r2)
2019	r2 = (r1+1)%leaf->rays.size();
2020
2021	#else
2022	// pick a random ray
2023	int base = Random(nrays);
2024	RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
2025	Intersectable *baseObject = baseRay->GetObject();
2026
2027	// and a random 5D derivative which will be used to find the minimal projected distances
2028	Vector3 spatialDerivative;
2029	Vector3 directionalDerivative;
2030
2031	while (1) {
2032	spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2033	RandomValue(-1.0f, 1.0f),
2034	RandomValue(-1.0f, 1.0f));
2035	float m = Magnitude(spatialDerivative);
2036	if (m != 0) {
2037	spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
2038	break;
2039	}
2040	}
2041
2042	while (1) {
2043	directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2044	RandomValue(-1.0f, 1.0f),
2045	0.0f);
2046	float m = Magnitude(directionalDerivative);
2047	if (m != 0) {
2048	directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
2049	break;
2050	}
2051	}
2052
2053	// now find the furthest sample from the same object and the closest from a different object
2054	int i = 0;
2055	float minDist = MAX_FLOAT;
2056	float maxDist = -MAX_FLOAT;
2057	r1 = base;
2058	r2 = base;
2059	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2060	ri != leaf->rays.end();
2061	ri++, i++) {
2062	float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
2063	*ri,
2064	spatialDerivative,
2065	directionalDerivative
2066	);
2067
2068	if ((*ri).GetObject() == baseObject) {
2069	if (dist > maxDist) {
2070	maxDist = dist;
2071	r1 = i;
2072	}
2073	} else {
2074	if (dist > 0.0f && dist < minDist) {
2075	minDist = dist;
2076	r2 = i;
2077	}
2078	}
2079	}
2080
2081	#endif
2082	}
2083
2084
2085	struct RayNeighbor {
2086	int rayInfo;
2087	float distance;
2088	RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
2089	};
2090
2091	void
2092	RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
2093	vector<RssTreeLeaf::SilhouetteRays> &rays
2094	)
2095	{
2096	// for every leaf find its neares neighbor from a different object
2097	vector<RayNeighbor> neighbors(leaf->rays.size());
2098	int i, j;
2099	for (i=0; i < leaf->rays.size(); i++)
2100	for (j=0; j < leaf->rays.size(); j++)
2101	if (i!=j) {
2102	float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
2103	if (d < neighbors[i].distance) {
2104	neighbors[i].distance = d;
2105	neighbors[i].rayInfo = j;
2106	}
2107	}
2108
2109	// now check which are the pairs of nearest neighbors
2110	for (i=0; i < leaf->rays.size(); i++) {
2111	int j = neighbors[i].rayInfo;
2112	if (neighbors[j].rayInfo == i) {
2113	// this is a silhouette edge pair
2114	if (i < j) {
2115	// generate an silhoutte ray pair
2116	rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
2117	&leaf->rays[j]));
2118	}
2119	} else {
2120	// this is not a silhouette ray - delete???
2121	}
2122	}
2123
2124	}
2125
2126
2127	bool
2128	GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
2129	const int index,
2130	int &i1,
2131	int &i2,
2132	int &i3)
2133	{
2134	int found = 0;
2135	int indices[3];
2136
2137	int size = rays.size();
2138	// use russian roulete selection for the tripple
2139	// number of free positions for the bullet
2140	int positions = size - 1;
2141	int i;
2142	for (i=0; i < size; i++) {
2143	if (rays[i]->mRay->Mailed())
2144	positions--;
2145	}
2146
2147	if (positions < 3)
2148	return false;
2149
2150	for (i=0; i < size; i++) {
2151	if (i != index && !rays[i]->mRay->Mailed()) {
2152	float p = (3 - found)/(float)positions;
2153	if (Random(1.0f) < p) {
2154	indices[found] = i;
2155	found++;
2156	}
2157	}
2158	positions--;
2159	}
2160	return true; // corr. matt
2161	}
2162
2163	bool
2164	RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
2165	const RssTreeNode::RayInfo &r1,
2166	const RssTreeNode::RayInfo &r2,
2167	const RssTreeNode::RayInfo &r3)
2168	{
2169
2170
2171	return false;
2172	}
2173
2174	int
2175	RssTree::RemoveInteriorRays(
2176	RssTreeLeaf *leaf
2177	)
2178	{
2179	#if 1
2180	// first collect all objects refered in this leaf
2181	map<Intersectable , vector<RssTreeLeaf::RayInfo > > rayMap;
2182
2183	RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
2184	for (; it != leaf->rays.end(); ++it) {
2185	Intersectable object = (it).GetObject();
2186
2187	rayMap[object].push_back(&(*it));
2188	// vector<RayInfo > data = rayMap.Find(object);
2189	// if (data) {
2190	// data->push_back(&(*it));
2191	// } else {
2192	// // rayMap[object] = vector<RayInfo *>;
2193	// rayMap[object].push_back(&(*it));
2194	// }
2195	}
2196
2197	// now go through all objects
2198	map<Intersectable , vector<RssTreeLeaf::RayInfo > >::iterator mi;
2199
2200	// infos of mailed rays are scheduled for removal
2201	VssRay::NewMail();
2202	for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
2203	vector<RssTreeLeaf::RayInfo > &rays = (mi).second;
2204	vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
2205	int rayIndex = 0;
2206
2207	for (; ri != rays.end(); ++ri, ++rayIndex) {
2208	RssTreeNode::RayInfo ray = ri;
2209	int tries = rays.size();
2210	for (int i = 0; i < tries; i++) {
2211	int r1, r2, r3;
2212	GetRandomTripple(rays,
2213	rayIndex,
2214	r1,
2215	r2,
2216	r3);
2217	if (IsRayConvexCombination(*ray,
2218	*rays[r1],
2219	*rays[r2],
2220	*rays[r3])) {
2221	ray->mRay->Mail();
2222	}
2223	}
2224	}
2225	}
2226
2227
2228	#endif
2229	return 0;
2230	}
2231
2232
2233	void
2234	RssTree::GenerateLeafRaysSimple(RssTreeLeaf *leaf,
2235	const int numberOfRays,
2236	SimpleRayContainer &rays)
2237	{
2238
2239	if (numberOfRays == 0)
2240	return;
2241
2242	int startIndex = (int)rays.size();
2243	// Debug<<"B"<<flush;
2244
2245	AxisAlignedBox3 box = GetBBox(leaf);
2246	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2247
2248	float boxSize = Magnitude(box.Diagonal());
2249
2250	const float smoothRange = 0.1f;
2251	if (smoothRange != 0.0f) {
2252	box.Scale(1.0f + smoothRange);
2253	dirBox.Scale(1.0f + smoothRange);
2254	}
2255
2256	// Debug<<"R"<<flush;
2257
2258	for (int i=0; i < numberOfRays; i++) {
2259	// Debug<<i<<flush;
2260	float r[5];
2261
2262	r[0] = RandomValue(0.0f, 1.0f);
2263	r[1] = RandomValue(0.0f, 1.0f);
2264	r[2] = RandomValue(0.0f, 1.0f);
2265	r[3] = RandomValue(0.0f, 1.0f);
2266	r[4] = RandomValue(0.0f, 1.0f);
2267
2268	GenerateLeafRay(r, rays, box, dirBox);
2269	}
2270
2271	// Debug<<"D"<<flush;
2272
2273	float density = 1.0f;
2274	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2275	for (int i=startIndex; i < rays.size(); i++) {
2276	rays[i].mPdf = p;
2277	}
2278	}
2279
2280	void
2281	RssTree::GenerateLeafRay(
2282	float *params,
2283	SimpleRayContainer &rays,
2284	const AxisAlignedBox3 &sbox,
2285	const AxisAlignedBox3 &sDirBBox
2286	)
2287	{
2288	Vector3 origin = sbox.GetPoint(Vector3(params[0], params[1], params[2]));
2289
2290	Vector3 dirVector = sDirBBox.GetPoint(Vector3(params[3], params[4], 0.0f));
2291
2292	Vector3 direction = VssRay::GetInvDirParam(dirVector.x, dirVector.y);
2293
2294	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2295	float dist = 0.0f;
2296	float minT, maxT;
2297
2298	float boxSize = Magnitude(sbox.Size());
2299	// compute interection of the ray with the box
2300	#if 1
2301	if (sbox.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2302	dist = (maxT + 0.0f*boxSize);
2303	#else
2304	dist = 0.5f*boxSize;
2305	#endif
2306
2307	origin += direction*dist;
2308
2309	//Debug<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2310	rays.push_back(SimpleRay(origin, direction));
2311	}
2312
2313	void
2314	RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
2315	const int numberOfRays,
2316	SimpleRayContainer &rays)
2317	{
2318
2319	if (numberOfRays == 0)
2320	return;
2321
2322
2323	int nrays = (int)leaf->rays.size();
2324	int startIndex = (int)rays.size();
2325
2326	AxisAlignedBox3 box = GetBBox(leaf);
2327	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2328
2329
2330	const float smoothRange = 0.1f;
2331	if (smoothRange != 0.0f) {
2332	box.Scale(1.0f + smoothRange);
2333	dirBox.Scale(1.0f + smoothRange);
2334	}
2335
2336	Exporter *exporter = NULL;
2337	VssRayContainer selectedRays;
2338	static int counter = 0;
2339	if (counter < 0) {
2340	char filename[128];
2341	sprintf(filename, "selected-rays%04d.x3d", counter);
2342	exporter = Exporter::GetExporter(filename);
2343	// exporter->SetWireframe();
2344	// exporter->ExportKdTree(*mKdTree);
2345	exporter->SetWireframe();
2346	// exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
2347	exporter->SetWireframe();
2348	exporter->ExportBox(box);
2349	exporter->SetFilled();
2350	counter++;
2351	}
2352
2353
2354	float extendedConvexCombinationProb = 0.0f; //0.7f
2355	// float silhouetteCheckPercentage = 1.0f; //0.5f
2356	float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
2357	float silhouetteObjectCheckPercentage = 0.8f;
2358
2359	// int numberOfTries = numberOfRays*4;
2360	int numberOfTries = numberOfRays;
2361	int generated = 0;
2362
2363	vector<Intersectable *> pvs(0);
2364	if (silhouetteCheckPercentage != 0.0f) {
2365	Intersectable::NewMail();
2366	pvs.reserve(leaf->GetPvsSize());
2367	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2368	ri != leaf->rays.end();
2369	ri++) {
2370	Intersectable object = (ri).GetObject();
2371	if (object) {
2372	if (!object->Mailed()) {
2373	object->Mail();
2374	object->mCounter = 1;
2375	pvs.push_back(object);
2376	if (exporter)
2377	exporter->ExportIntersectable(object);
2378	} else {
2379	object->mCounter++;
2380	}
2381	}
2382	}
2383	// sort objects based on their mCounter
2384	sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
2385	}
2386
2387	for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
2388	bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
2389	extendedConvexCombinationProb));
2390
2391	Vector3 origin, direction;
2392	// generate half of convex combination and half of random rays
2393	if (useExtendedConvexCombination) {
2394	// pickup 2 random rays
2395	int r1, r2;
2396	PickEdgeRays(leaf, r1, r2);
2397
2398
2399	Vector3 o1 = leaf->rays[r1].GetOrigin();
2400	Vector3 o2 = leaf->rays[r2].GetOrigin();
2401
2402	const float overlap = 0.0f;
2403	float w1, w2;
2404	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2405	origin = w1o1 + w2o2;
2406	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2407	direction =
2408	w1*leaf->rays[r1].GetDir() +
2409	w2*leaf->rays[r2].GetDir();
2410
2411	// shift the origin a little bit
2412	direction.Normalize();
2413
2414	} else {
2415	Vector3 dirVector;
2416
2417	Vector3 pVector;
2418	Vector3 dVector;
2419
2420	pVector = Vector3(RandomValue(0.0f, 1.0f),
2421	RandomValue(0.0f, 1.0f),
2422	RandomValue(0.0f, 1.0f));
2423
2424	dVector = Vector3(RandomValue(0.0f, 1.0f),
2425	RandomValue(0.0f, 1.0f),
2426	0.0f);
2427
2428
2429	origin = box.GetPoint(pVector);
2430	dirVector = dirBox.GetPoint(dVector);
2431
2432	direction = VssRay::GetInvDirParam(dirVector.x,dirVector.y);
2433	}
2434
2435
2436	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2437	float dist = 0.0f;
2438	float minT, maxT;
2439	// compute interection of the ray with the box
2440
2441	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2442	dist = maxT;
2443
2444
2445	origin += direction*dist;
2446
2447	bool intersects = false;
2448
2449	if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
2450	if (exporter) {
2451	VssRay ray = new VssRay(origin, origin + 100direction, NULL, NULL);
2452	selectedRays.push_back(ray);
2453	}
2454
2455	// check whether the ray does not intersect already visible objects
2456	Ray traversalRay;
2457	traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2458	for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2459	oi != pvs.end();
2460	oi++) {
2461	Intersectable object = oi;
2462	// do not test every object
2463	if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
2464	object->CastRay(traversalRay) ) {
2465	intersects = true;
2466	break;
2467	}
2468	}
2469	}
2470
2471
2472	if (!intersects) {
2473	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2474	rays.push_back(SimpleRay(origin, direction));
2475	generated++;
2476	}
2477	}
2478
2479	// cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
2480	// assign pdfs to the generated rays
2481	// float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
2482	float density = 1.0f;
2483
2484	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2485	for (int i=startIndex; i < rays.size(); i++) {
2486	rays[i].mPdf = p;
2487	}
2488
2489
2490	if (exporter) {
2491	exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2492	delete exporter;
2493	CLEAR_CONTAINER(selectedRays);
2494	}
2495
2496	}
2497
2498	int
2499	RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2500	const float ratio)
2501	{
2502	int i;
2503	int j;
2504
2505	for (j=0, i=0; i < leaf->rays.size(); i++) {
2506	if (Random(1.0f) < ratio) {
2507	// copy a valid sample
2508	if (i!=j)
2509	leaf->rays[j] = leaf->rays[i];
2510	j++;
2511	} else {
2512	// delete the ray
2513	leaf->rays[i].mRay->Unref();
2514	if (leaf->rays[i].mRay->RefCount() != 0) {
2515	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2516	exit(1);
2517	}
2518	delete leaf->rays[i].mRay;
2519	}
2520	}
2521
2522
2523	leaf->rays.resize(j);
2524	int removed = (i-j);
2525	stat.rayRefs -= removed;
2526	return removed;
2527	}
2528
2529	int
2530	RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2531	const float ratio)
2532	{
2533	int i;
2534
2535	if (leaf->rays.size() == 0)
2536	return 0;
2537
2538	sort(leaf->rays.begin(),
2539	leaf->rays.end(),
2540	RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2541
2542	int desired = int(ratio*leaf->rays.size());
2543	int removed = (int)leaf->rays.size() - desired;
2544
2545	for (i=desired; i < (int)leaf->rays.size(); i++) {
2546	// delete the ray
2547	leaf->rays[i].mRay->Unref();
2548	if (leaf->rays[i].mRay->RefCount() != 0) {
2549	cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
2550	cerr<<"desired="<<desired<<"i="<<i<<" size="<<(int)leaf->rays.size()<<endl;
2551	exit(1);
2552	}
2553	delete leaf->rays[i].mRay;
2554	}
2555
2556	leaf->rays.resize(desired);
2557	stat.rayRefs -= removed;
2558	return removed;
2559	}
2560
2561
2562	int
2563	RssTree::PruneRays(
2564	const int maxRays
2565	)
2566	{
2567
2568	stack<RssTreeNode *> tstack;
2569	int prunned = 0;
2570
2571	// prune more rays to amortize the procedure
2572	int desired = (int)(maxRays * 0.8f);
2573
2574	Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;
2575
2576	// prune random rays from each leaf so that the ratio's remain the same
2577	PushRoots(tstack);
2578
2579	if ( maxRays >= stat.rayRefs )
2580	return 0;
2581
2582	float ratio = desired/(float)stat.rayRefs;
2583
2584	while (!tstack.empty()) {
2585	RssTreeNode *node = tstack.top();
2586	tstack.pop();
2587
2588	if (node->IsLeaf()) {
2589	RssTreeLeaf leaf = (RssTreeLeaf )node;
2590	prunned += PruneRaysRandom(leaf, ratio);
2591	//prunned += PruneRaysContribution(leaf, ratio);
2592	} else {
2593	RssTreeInterior in = (RssTreeInterior )node;
2594	// both nodes for directional splits
2595	tstack.push(in->front);
2596	tstack.push(in->back);
2597	}
2598	}
2599
2600	Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
2601
2602	return prunned;
2603	}
2604
2605	int
2606	RssTree::GenerateRays(const int numberOfRays,
2607	const int numberOfLeaves,
2608	SimpleRayContainer &rays)
2609	{
2610
2611	if (1) {
2612	float param[5];
2613
2614	for (int i=0; i < numberOfRays; i++) {
2615	mHalton.GetNext(5, param);
2616	GenerateRay(param, rays);
2617	}
2618
2619	return rays.size();
2620	}
2621
2622	vector<RssTreeLeaf *> leaves;
2623
2624	Debug<<"Collecting leaves..."<<endl<<flush;
2625	CollectLeaves(leaves);
2626	Debug<<"done."<<endl<<flush;
2627
2628
2629	if (numberOfLeaves != leaves.size()) {
2630	Debug<<"sorting leaves..."<<endl<<flush;
2631	sort(leaves.begin(),
2632	leaves.end(),
2633	GreaterContribution);
2634	Debug<<"done."<<endl<<flush;
2635	}
2636
2637
2638	Debug<<"Evaluating leaf importance..."<<endl<<flush;
2639
2640	float sumContrib = 0.0;
2641	int i;
2642	int k = 0;
2643	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2644	if (ValidLeaf(leaves[i])) {
2645	float c = leaves[i]->GetImportance();
2646	sumContrib += c;
2647	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
2648	k++;
2649	}
2650
2651	Debug<<"done."<<endl<<flush;
2652
2653	float avgContrib = sumContrib/numberOfLeaves;
2654
2655	// always generate at leat n ray per leaf
2656	int fixedPerLeaf = 1;
2657	// int fixedPerLeaf = 0;
2658	int fixed = fixedPerLeaf*(int)leaves.size();
2659
2660	// check if the number of fixed is not too big
2661
2662	if (fixed >= numberOfRays/2) {
2663	fixedPerLeaf = 0;
2664	fixed = 0;
2665	}
2666
2667	int iGenerated = numberOfRays - fixed;
2668	float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
2669
2670	k = 0;
2671
2672	Debug<<"generating leaf rays..."<<endl<<flush;
2673
2674	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2675	if (ValidLeaf(leaves[i])) {
2676	RssTreeLeaf *leaf = leaves[i];
2677	float c = leaf->GetImportance();
2678
2679	int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
2680	GenerateLeafRaysSimple(leaf, num, rays);
2681
2682	k++;
2683	}
2684
2685	Debug<<"done."<<endl<<flush;
2686
2687	return (int)rays.size();
2688	}
2689
2690
2691	float
2692	RssTree::GetAvgPvsSize()
2693	{
2694	stack<RssTreeNode *> tstack;
2695	PushRoots(tstack);
2696
2697	int sumPvs = 0;
2698	int leaves = 0;
2699	while (!tstack.empty()) {
2700	RssTreeNode *node = tstack.top();
2701	tstack.pop();
2702
2703	if (node->IsLeaf()) {
2704	RssTreeLeaf leaf = (RssTreeLeaf )node;
2705	// update pvs size
2706	UpdatePvsSize(leaf);
2707	sumPvs += leaf->GetPvsSize();
2708	leaves++;
2709	} else {
2710	RssTreeInterior in = (RssTreeInterior )node;
2711	// both nodes for directional splits
2712	tstack.push(in->front);
2713	tstack.push(in->back);
2714	}
2715	}
2716
2717
2718	return sumPvs/(float)leaves;
2719	}
2720
2721	float weightAbsContributions = ABS_CONTRIBUTION_WEIGHT;
2722	// if small very high importance of the last sample
2723	// if 1.0f then weighs = 1 1/2 1/3 1/4
2724	//float passSampleWeightDecay = 1.0f;
2725	float passSampleWeightDecay = 1000.0f;
2726
2727	float
2728	RssTree::GetSampleWeight(const int pass)
2729	{
2730	int passDiff = mCurrentPass - pass;
2731	float weight;
2732	if (1)
2733	weight = 1.0f/(passDiff + passSampleWeightDecay);
2734	else
2735	switch (passDiff) {
2736	case 0:
2737	// weight = 1.0f;
2738	// break;
2739	default:
2740	weight = 1.0f;
2741	break;
2742	// case 1:
2743	// weight = 0.5f;
2744	// break;
2745	// case 2:
2746	// weight = 0.25f;
2747	// break;
2748	// case 3:
2749	// weight = 0.12f;
2750	// break;
2751	// case 4:
2752	// weight = 0.06f;
2753	// break;
2754	// default:
2755	// weight = 0.03f;
2756	// break;
2757	}
2758	return weight;
2759	}
2760
2761	void
2762	RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2763	float &weight,
2764	float &contribution)
2765	{
2766	VssRay *ray = info.mRay;
2767	weight = GetSampleWeight(mCurrentPass);
2768	contribution =
2769	weightAbsContributions*ray->mPvsContribution +
2770	(1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2771	// store the computed value
2772	info.mRay->mWeightedPvsContribution = weight*contribution;
2773	}
2774
2775
2776	void
2777	RssTree::ComputeImportance(RssTreeLeaf *leaf)
2778	{
2779	if (leaf->mTotalRays)
2780	leaf->mImportance = leaf->mTotalContribution / leaf->mTotalRays;
2781	else
2782	leaf->mImportance = 0.0f;
2783	return;
2784
2785	if (0)
2786	leaf->mImportance = leaf->GetAvgRayContribution();
2787	else {
2788	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2789	it_end = leaf->rays.end();
2790
2791	float sumContributions = 0.0f;
2792	float sumRelContributions = 0.0f;
2793	float sumWeights = 0.0f;
2794	float maxContribution = 0.0f;
2795	float maxRelContribution = 0.0f;
2796	float sumLength = 0;
2797
2798	for (; it != it_end; ++it) {
2799	VssRay ray = (it).mRay;
2800	float sumLength = ray->Length();
2801
2802	float weight = GetSampleWeight(ray->mPass);
2803
2804	float c1 = weight*ray->mPvsContribution;
2805	float c2 = weight*ray->mRelativePvsContribution;
2806	sumContributions += c1;
2807	if (c1 > maxContribution)
2808	maxContribution = c1;
2809
2810	//$$ 20.7. changed to sqr to pronouce higher contribution so that
2811	// they do not get smoothed!!
2812	//sumRelContributions += weight*ray->mRelativePvsContribution;
2813
2814	sumRelContributions += c2;
2815	if (c2 > maxRelContribution)
2816	maxRelContribution = c2;
2817
2818	//sumWeights += weight;
2819	sumWeights += 1.0f;
2820	}
2821
2822	float distImportance = 0.0f;
2823
2824	float spatialArea = GetBBox(leaf).SurfaceArea();
2825	float dirArea = GetDirBBox(leaf).SurfaceArea();
2826	if (leaf->rays.size()) {
2827	// compute average length of the ray
2828	float avgLength = sumLength/leaf->rays.size();
2829	// compute estimated area of the visible surface corresponding to these rays
2830	float ss = spatialArea/2.0f;
2831	float sd = dirArea;
2832	float s = ss + avgLength(2sqrt(sssd) + avgLengthsd);
2833
2834	distImportance = s/leaf->rays.size();
2835	}
2836
2837	// $$
2838	// sumWeights = leaf->mTotalRays;
2839
2840	//$$ test 30.11. 2006
2841	// normalize per volume of the cell not the number of rays!
2842	// sumWeights = spatialArea*dirArea;
2843
2844	if (sumWeights != 0.0f) {
2845	leaf->mImportance =
2846	sqr(((weightAbsContributions*sumContributions +
2847	(1.0f - weightAbsContributions)*sumRelContributions)/sumWeights));
2848
2849	// leaf->mImportance =
2850	// (weightAbsContributions*maxContribution +
2851	// (1.0f - weightAbsContributions)*maxRelContribution)/sumWeights;
2852
2853	} else
2854	leaf->mImportance = 0.0f;
2855
2856	}
2857
2858	// return GetAvgRayContribution()*mImportance;
2859	//return GetAvgRayContribution();
2860	}
2861
2862
2863	float
2864	RssTreeNode::GetImportance() const
2865	{
2866	// return sqr(mImportance);
2867	return mImportance;
2868	}
2869
2870
2871	float
2872	RssTreeLeaf::ComputePvsEntropy() const
2873	{
2874	int samples = 0;
2875	Intersectable::NewMail();
2876	// set all object as belonging to the fron pvs
2877	for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
2878	ri != rays.end();
2879	ri++)
2880	if ((*ri).mRay->IsActive()) {
2881	Intersectable object = (ri).GetObject();
2882	if (object) {
2883	if (!object->Mailed()) {
2884	object->Mail();
2885	object->mCounter = 1;
2886	} else
2887	object->mCounter++;
2888	samples++;
2889	}
2890	}
2891
2892	float entropy = 0.0f;
2893
2894	if (samples > 1) {
2895	Intersectable::NewMail();
2896	for(RayInfoContainer::const_iterator ri = rays.begin();
2897	ri != rays.end();
2898	ri++)
2899	if ((*ri).mRay->IsActive()) {
2900	Intersectable object = (ri).GetObject();
2901	if (object) {
2902	if (!object->Mailed()) {
2903	object->Mail();
2904	float p = object->mCounter/(float)samples;
2905	entropy -= p*log(p);
2906	}
2907	}
2908	}
2909	entropy = entropy/log((float)samples);
2910	}
2911	else
2912	entropy = 1.0f;
2913
2914	return entropy;
2915	}
2916
2917	float
2918	RssTreeLeaf::ComputeRayLengthEntropy() const
2919	{
2920	// get sum of all ray lengths
2921	// consider only passing rays or originating rays
2922	float sum = 0.0f;
2923	int samples = 0;
2924	int i=0;
2925	for(RayInfoContainer::const_iterator ri = rays.begin();
2926	ri != rays.end();
2927	ri++)
2928	if ((*ri).mRay->IsActive()) {
2929	// float s;
2930	// if (i == 0)
2931	// s = 200;
2932	// else
2933	// s = 100;
2934	// i++;
2935
2936	sum += (*ri).mRay->GetSize();
2937	samples++;
2938	}
2939
2940
2941	float entropy = 0.0f;
2942
2943	if (samples > 1) {
2944	i = 0;
2945	for(RayInfoContainer::const_iterator ri = rays.begin();
2946	ri != rays.end();
2947	ri++)
2948	if ((*ri).mRay->IsActive()) {
2949	// float s;
2950	// if (i==0)
2951	// s = 200;
2952	// else
2953	// s = 100;
2954	// i++;
2955	// float p = s/sum;
2956	float p = (*ri).mRay->GetSize()/sum;
2957	entropy -= p*log(p);
2958	}
2959	entropy = entropy/log((float)samples);
2960	} else
2961	entropy = 1.0f;
2962
2963	return entropy;
2964	}
2965
2966
2967
2968	float
2969	RssTreeLeaf::ComputeEntropyImportance() const
2970	{
2971
2972	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2973	return 1.0f - ComputePvsEntropy();
2974	}
2975
2976	float
2977	RssTreeLeaf::ComputeRayContributionImportance() const
2978	{
2979	// mPvsEntropy = ComputePvsEntropy();
2980	// mRayLengthEntropy = ComputeRayLengthEntropy();
2981
2982	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2983	return 1.0f - ComputePvsEntropy();
2984	}
2985
2986
2987	AxisAlignedBox3
2988	RssTree::GetShrankedBBox(const RssTreeNode *node) const
2989	{
2990	if (node->parent == NULL)
2991	return bbox;
2992
2993	if (!node->IsLeaf())
2994	return ((RssTreeInterior *)node)->bbox;
2995
2996	// evaluate bounding box from the ray origins
2997	AxisAlignedBox3 box;
2998	box.Initialize();
2999	RssTreeLeaf leaf = (RssTreeLeaf )node;
3000	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
3001	ri != leaf->rays.end();
3002	ri++)
3003	if ((*ri).mRay->IsActive()) {
3004	box.Include((*ri).GetOrigin());
3005	}
3006	return box;
3007	}
3008
3009
3010	int
3011	RssTree::CollectRays(VssRayContainer &rays,
3012	const int number)
3013	{
3014	VssRayContainer allRays;
3015	CollectRays(allRays);
3016
3017
3018	int desired = min(number, (int)allRays.size());
3019	float prob = desired/(float)allRays.size();
3020	//int skip = allRays.size()/desired;
3021	while (rays.size() < desired) {
3022	VssRayContainer::const_iterator it = allRays.begin();
3023	for (; it != allRays.end() && rays.size() < desired; it++) {
3024	if (Random(1.0f) < prob)
3025	rays.push_back(*it);
3026	}
3027	}
3028	return (int)rays.size();
3029	}
3030
3031
3032	int
3033	RssTree::CollectRays(VssRayContainer &rays
3034	)
3035	{
3036	VssRay::NewMail();
3037
3038	stack<RssTreeNode *> tstack;
3039	PushRoots(tstack);
3040
3041	while (!tstack.empty()) {
3042	RssTreeNode *node = tstack.top();
3043	tstack.pop();
3044
3045	if (node->IsLeaf()) {
3046	RssTreeLeaf leaf = (RssTreeLeaf )node;
3047	// update pvs size
3048	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
3049	for (;it != leaf->rays.end(); ++it)
3050	if (!(*it).mRay->Mailed()) {
3051	(*it).mRay->Mail();
3052	rays.push_back((*it).mRay);
3053	}
3054	} else {
3055	RssTreeInterior in = (RssTreeInterior )node;
3056	// both nodes for directional splits
3057	tstack.push(in->front);
3058	tstack.push(in->back);
3059	}
3060	}
3061
3062	return (int)rays.size();
3063	}
3064
3065
3066	RssTreeNode *
3067	RssTree::GetRoot(Intersectable *object) const
3068	{
3069	if (mPerObjectTree && object) {
3070	int id = object->GetId()-1;
3071	if (id < 0 \|\| id >= mRoots.size()) {
3072	Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<(int)mRoots.size()<<
3073	endl<<flush;
3074	id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
3075	}
3076	return mRoots[id];
3077	} else
3078	return mRoots[0];
3079	}
3080
3081	void
3082	RssTree::PushRoots(priority_queue<TraversalData> &st) const
3083	{
3084	for (int i=0; i < mRoots.size(); i++)
3085	st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
3086	}
3087
3088	void
3089	RssTree::PushRoots(stack<RssTreeNode *> &st) const
3090	{
3091	for (int i=0; i < mRoots.size(); i++)
3092	st.push(mRoots[i]);
3093	}
3094
3095	void
3096	RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info) const
3097	{
3098	for (int i=0; i < mRoots.size(); i++)
3099	st.push(RayTraversalData(mRoots[i], info));
3100	}
3101
3102
3103
3104
3105
3106	}

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

Download in other formats: