Context Navigation

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

Revision 1874, 74.5 KB checked in by bittner, 18 years ago (diff)

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	}
678	}
679	} else {
680	const int costMethod = 1;
681	// importance based cost
682	switch (costMethod) {
683	case 0: {
684	break;
685	}
686	case 1: {
687	float newContrib =
688	sqr(info.pvsBack/(info.raysBack + Limits::Small)) +
689	sqr(info.pvsFront/(info.raysFront + Limits::Small));
690	float oldContrib = sqr(leaf->GetAvgRayContribution());
691	info.costRatio = oldContrib/newContrib;
692	break;
693	}
694	case 2: {
695	float newCost = (info.pvsBack + info.pvsFront)*0.5f;
696	float oldCost = (float)pvsSize;
697	info.costRatio = newCost/oldCost;
698	break;
699	}
700	case 3: {
701	float newCost = (float)abs(info.raysBack - info.raysFront);
702	float oldCost = (float)leaf->rays.size();
703	info.costRatio = newCost/oldCost;
704	break;
705	}
706	}
707	}
708	}
709
710
711	void
712	RssTree::EvalCostRatio(
713	RssTreeLeaf *leaf,
714	SplitInfo &info
715	)
716	{
717	info.rays = 0;
718	info.raysBack = 0;
719	info.raysFront = 0;
720	info.pvsFront = 0;
721	info.pvsBack = 0;
722	info.viewCells = 0;
723	info.viewCellsFront = 0;
724	info.viewCellsBack = 0;
725
726
727
728	float sumWeights = Limits::Small;
729	float sumWeightsBack = Limits::Small;
730	float sumWeightsFront = Limits::Small;
731
732	float sumContribution = 0.0f;
733	float sumContributionBack = 0.0f;
734	float sumContributionFront = 0.0f;
735
736	Intersectable::NewMail();
737
738	#if EVAL_VIEWCELLS
739	// count the numebr of viewcells first
740	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
741	ri != leaf->rays.end();
742	ri++)
743	if ((*ri).mRay->IsActive()) {
744
745	ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
746	for (; it != (*ri).mRay->mViewCells.end(); ++it) {
747	if (!(*it)->Mailed()) {
748	(*it)->Mail();
749	info.viewCells++;
750	}
751	}
752	}
753	#endif
754
755	Intersectable::NewMail(3);
756
757	// this is the main ray classification loop!
758	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
759	ri != leaf->rays.end();
760	ri++)
761	if ((*ri).mRay->IsActive()) {
762	int side;
763
764
765	// determine the side of this ray with respect to the plane
766	side = (*ri).ComputeRaySide(info.axis, info.position);
767
768	float weight, contribution;
769
770	GetRayContribution(*ri, weight, contribution);
771	sumWeights += weight;
772	sumContribution += contribution;
773	info.rays++;
774	if (side <= 0) {
775	info.raysBack++;
776	sumWeightsBack += weight;
777	sumContributionBack += contribution;
778	}
779
780	if (side >= 0) {
781	info.raysFront++;
782	sumWeightsFront += weight;
783	sumContributionFront += contribution;
784	}
785
786	AddObject2Pvs((*ri).GetObject(), side, info.pvsBack, info.pvsFront);
787
788	#if EVAL_VIEWCELLS
789	AddViewcells2Pvs((*ri).mRay->mViewCells,
790	side,
791	info.viewCellsBack,
792	info.viewCellsFront);
793	#endif
794	}
795
796	if (sumWeights!=0.0f)
797	info.contribution = sumContribution/sumWeights;
798	else
799	info.contribution = 0.0f;
800
801	if (sumWeightsBack!=0.0f)
802	info.contributionBack = sumContributionBack/sumWeightsBack;
803	else
804	info.contributionBack = 0.0f;
805
806	if (sumWeightsFront!=0.0f)
807	info.contributionFront = sumContributionFront/sumWeightsFront;
808	else
809	info.contributionFront = 0.0f;
810
811	// info.costRatio = 0.1 + Random(0.5f);
812	// return;
813
814	GetCostRatio(
815	leaf,
816	info);
817
818	//cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
819	// float oldCost = leaf->rays.size();
820
821	// cout<<"ratio="<<ratio<<endl;
822	}
823
824	void
825	RssTree::BestCostRatio(
826	RssTreeLeaf *leaf,
827	SplitInfo &info
828	)
829	{
830	SplitInfo nInfo[5];
831	int bestAxis = -1;
832
833	AxisAlignedBox3 sBox = GetBBox(leaf);
834	AxisAlignedBox3 dBox = GetDirBBox(leaf);
835	// int sAxis = box.Size().DrivingAxis();
836	int sAxis = sBox.Size().DrivingAxis();
837	int dAxis = dBox.Size().DrivingAxis() + 3;
838
839	float dirSplitBoxSize = 0.01f;
840	bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
841
842	bool useCyclingAxis = true;
843
844	if (splitType == ESplitHeuristic \|\|
845	(splitType == ESplitHybrid && leaf->depth > mHybridDepth))
846	useCyclingAxis = false;
847
848	int cyclingAxis = 0;
849
850	if (leaf->parent)
851	cyclingAxis = (leaf->parent->axis + 1)%5;
852
853	if (0 && useCyclingAxis) {
854	int axis = cyclingAxis;
855	info.position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
856	info.axis = axis;
857	info.costRatio = 0.5f; // not true but good for speeding up the subdivision at the top levels!
858	info.raysBack = info.raysFront = leaf->rays.size()/2;
859	return;
860	}
861
862
863	for (int axis = 0; axis < 5; axis++)
864	if (
865	(axis < 3 && (leaf->depth < mDirSplitDepth \|\| mInterleaveDirSplits)) \|\|
866	(axis >= 3 && (leaf->depth >= mDirSplitDepth))
867	) {
868	nInfo[axis].axis = axis;
869	if ((!useCyclingAxis && (!mSplitUseOnlyDrivingAxis \|\| axis == sAxis \|\| axis == dAxis))
870	\|\| (useCyclingAxis && (axis == cyclingAxis))
871	) {
872
873	if (splitType == ESplitRegular) {
874	if (axis < 3)
875	nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
876	else
877	nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
878	EvalCostRatio(leaf,
879	nInfo[axis]);
880	} else
881	if (splitType == ESplitHeuristic) {
882	EvalCostRatioHeuristic(
883	leaf,
884	nInfo[axis]
885	);
886	} else
887	if (splitType == ESplitHybrid) {
888	if (leaf->depth > mHybridDepth)
889	EvalCostRatioHeuristic(
890	leaf,
891	nInfo[axis]
892	);
893	else {
894	if (axis < 3)
895	nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
896	else
897	nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
898
899	EvalCostRatio(leaf,
900	nInfo[axis]
901	);
902	}
903	} else {
904	cerr<<"RssTree: Unknown split heuristics\n";
905	exit(1);
906	}
907
908	if ( bestAxis == -1)
909	bestAxis = axis;
910	else
911	if ( nInfo[axis].costRatio < nInfo[bestAxis].costRatio )
912	bestAxis = axis;
913	}
914	}
915
916	info = nInfo[bestAxis];
917	}
918
919
920	void
921	RssTree::EvalCostRatioHeuristic(
922	RssTreeLeaf *leaf,
923	SplitInfo &info
924	)
925	{
926	AxisAlignedBox3 box;
927	float minBox, maxBox;
928
929	if (info.axis < 3) {
930	box = GetBBox(leaf);
931	minBox = box.Min(info.axis);
932	maxBox = box.Max(info.axis);
933	} else {
934	box = GetDirBBox(leaf);
935	minBox = box.Min(info.axis-3);
936	maxBox = box.Max(info.axis-3);
937	}
938
939	SortSubdivisionCandidates(leaf, info.axis);
940
941	// go through the lists, count the number of objects left and right
942	// and evaluate the following cost funcion:
943	// C = ct_div_ci + (qlrl + qrrr)/queries
944
945	SplitInfo currInfo;
946	currInfo.axis = info.axis;
947
948	currInfo.raysBack = 0;
949	currInfo.raysFront = leaf->rays.size();
950
951	currInfo.pvsBack = 0;
952	currInfo.pvsFront = leaf->GetPvsSize();
953
954
955	float sumWeights = Limits::Small;
956	float sumWeightsBack = Limits::Small;
957	float sumWeightsFront = Limits::Small;
958
959	float sumContribution = 0.0f;
960	float sumContributionBack = 0.0f;
961	float sumContributionFront = 0.0f;
962
963
964	float sizeBox = maxBox - minBox;
965
966	float minBand = minBox + 0.1f*(maxBox - minBox);
967	float maxBand = minBox + 0.9f*(maxBox - minBox);
968
969	// best cost ratio
970	info.costRatio = 1e20f;
971
972	currInfo.viewCells = 0;
973	currInfo.rays = 0;
974
975	Intersectable::NewMail();
976	// set all object as belonging to the fron pvs
977	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
978	ri != leaf->rays.end();
979	ri++)
980	if ((*ri).mRay->IsActive()) {
981
982	float weight, contribution;
983	GetRayContribution(*ri, weight, contribution);
984
985	sumWeights += weight;
986	sumContribution += contribution;
987	currInfo.rays++;
988
989	Intersectable object = (ri).GetObject();
990	if (object)
991	if (!object->Mailed()) {
992	object->Mail();
993	object->mCounter = 1;
994	} else
995	object->mCounter++;
996
997	#if EVAL_VIEWCELLS
998	// and do the same for all viewcells
999	ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
1000
1001	for (; it != (*ri).mRay->mViewCells.end(); ++it) {
1002	ViewCell viewcell = it;
1003	if (!viewcell->Mailed()) {
1004	currInfo.viewCells++;
1005	viewcell->Mail();
1006	viewcell->mCounter = 1;
1007	} else
1008	viewcell->mCounter++;
1009	}
1010	#endif
1011	}
1012
1013	if (sumWeights!=0.0f)
1014	currInfo.contribution = sumContribution/sumWeights;
1015	else
1016	currInfo.contribution = 0.0f;
1017
1018	sumWeightsFront = sumWeights;
1019	sumContributionFront = sumContribution;
1020
1021	sumWeightsBack = 0;
1022	sumContributionBack = 0;
1023
1024	currInfo.viewCellsBack = 0;
1025	currInfo.viewCellsFront = currInfo.viewCells;
1026
1027	Intersectable::NewMail();
1028
1029	for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
1030	ci < splitCandidates->end();
1031	ci++) {
1032	switch ((*ci).type) {
1033	case SortableEntry::ERayMin: {
1034	currInfo.raysFront--;
1035	currInfo.raysBack++;
1036
1037	RssTreeNode::RayInfo rayInfo = (RssTreeNode::RayInfo ) (*ci).data;
1038
1039	float weight, contribution;
1040	GetRayContribution(*rayInfo, weight, contribution);
1041
1042	sumWeightsBack += weight;
1043	sumContributionBack += contribution;
1044
1045	sumWeightsFront -= weight;
1046	sumContributionFront -= contribution;
1047
1048	Intersectable *object = rayInfo->GetObject();
1049	if (object) {
1050	if (!object->Mailed()) {
1051	object->Mail();
1052	currInfo.pvsBack++;
1053	}
1054	if (--object->mCounter == 0)
1055	currInfo.pvsFront--;
1056	}
1057
1058	#if EVAL_VIEWCELLS
1059	ViewCellContainer::const_iterator it = rayInfo->mRay->mViewCells.begin();
1060	for (; it != rayInfo->mRay->mViewCells.end(); ++it) {
1061	ViewCell viewcell = it;
1062	if (!viewcell->Mailed()) {
1063	viewcell->Mail();
1064	currInfo.viewCellsBack++;
1065	}
1066	if (--viewcell->mCounter == 0)
1067	currInfo.viewCellsFront--;
1068	}
1069	#endif
1070
1071	break;
1072	}
1073	}
1074
1075	float position = (*ci).value;
1076
1077	if (position > minBand && position < maxBand) {
1078	currInfo.position = position;
1079
1080
1081	if (sumWeightsBack!=0.0f)
1082	info.contributionBack = sumContributionBack/sumWeightsBack;
1083	else
1084	info.contributionBack = 0.0f;
1085
1086	if (sumWeightsFront!=0.0f)
1087	info.contributionFront = sumContributionFront/sumWeightsFront;
1088	else
1089	info.contributionFront = 0.0f;
1090
1091	GetCostRatio(
1092	leaf,
1093	currInfo);
1094
1095
1096	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
1097	// cout<<"cost= "<<sum<<endl;
1098
1099	if (currInfo.costRatio < info.costRatio) {
1100	info = currInfo;
1101	}
1102	}
1103	}
1104
1105
1106	// cout<<"===================="<<endl;
1107	// cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
1108	// <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
1109	}
1110
1111	void
1112	RssTree::SortSubdivisionCandidates(
1113	RssTreeLeaf *node,
1114	const int axis
1115	)
1116	{
1117
1118	splitCandidates->clear();
1119
1120	int requestedSize = 2*(node->rays.size());
1121	// creates a sorted split candidates array
1122	if (splitCandidates->capacity() > 500000 &&
1123	requestedSize < (int)(splitCandidates->capacity()/10) ) {
1124
1125	delete splitCandidates;
1126	splitCandidates = new vector<SortableEntry>;
1127	}
1128
1129	splitCandidates->reserve(requestedSize);
1130
1131	// insert all queries
1132	for(RssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
1133	ri < node->rays.end();
1134	ri++) {
1135	if ((*ri).mRay->IsActive()) {
1136	if (axis < 3) {
1137	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
1138	(*ri).GetOrigin(axis),
1139	(void )&(ri))
1140	);
1141	} else {
1142	float pos = (*ri).GetDirParametrization(axis-3);
1143	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
1144	pos,
1145	(void )&(ri))
1146	);
1147	}
1148	}
1149	}
1150
1151	stable_sort(splitCandidates->begin(), splitCandidates->end());
1152	}
1153
1154
1155	void
1156	RssTree::EvaluateLeafStats(const TraversalData &data)
1157	{
1158
1159	// the node became a leaf -> evaluate stats for leafs
1160	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1161
1162	if (data.depth >= termMaxDepth)
1163	stat.maxDepthNodes++;
1164
1165	// if ( (int)(leaf->rays.size()) < termMinCost)
1166	// stat.minCostNodes++;
1167	if ( leaf->GetPvsSize() <= termMinPvs)
1168	stat.minPvsNodes++;
1169
1170	if ( leaf->GetPvsSize() <= termMinRays)
1171	stat.minRaysNodes++;
1172
1173	if (leaf->GetAvgRayContribution() > termMaxRayContribution )
1174	stat.maxRayContribNodes++;
1175
1176	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
1177	stat.minSizeNodes++;
1178	}
1179
1180	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
1181	stat.maxRayRefs = leaf->rays.size();
1182
1183	}
1184
1185	bool
1186	RssTree::TerminationCriteriaSatisfied(RssTreeLeaf *leaf)
1187	{
1188	return ( (leaf->GetPvsSize() <= termMinPvs) \|\|
1189	(leaf->rays.size() <= termMinRays) \|\|
1190	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) \|\|
1191	(leaf->depth >= termMaxDepth) \|\|
1192	(SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
1193	);
1194	}
1195
1196
1197	RssTreeNode *
1198	RssTree::SubdivideNode(
1199	RssTreeLeaf *leaf,
1200	const AxisAlignedBox3 &box,
1201	AxisAlignedBox3 &backBBox,
1202	AxisAlignedBox3 &frontBBox
1203	)
1204	{
1205
1206	if (TerminationCriteriaSatisfied(leaf)) {
1207	#if 0
1208	if (leaf->depth >= termMaxDepth) {
1209	cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
1210	cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
1211	}
1212	#endif
1213
1214	return leaf;
1215	}
1216
1217	SplitInfo info;
1218
1219	// select subdivision axis
1220	int axis = SelectPlane( leaf,
1221	box,
1222	info
1223	);
1224	// Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
1225
1226	if (axis == -1) {
1227	return leaf;
1228	}
1229
1230	stat.nodes+=2;
1231	stat.leaves += 1;
1232	stat.splits[axis]++;
1233
1234	// add the new nodes to the tree
1235	RssTreeInterior *node = new RssTreeInterior(leaf->parent);
1236
1237	node->axis = axis;
1238	node->position = info.position;
1239
1240	node->bbox = GetBBox(leaf);
1241	node->dirBBox = GetDirBBox(leaf);
1242
1243	RssTreeLeaf *back = new RssTreeLeaf(node, info.raysBack);
1244	RssTreeLeaf *front = new RssTreeLeaf(node, info.raysFront);
1245
1246
1247
1248	// replace a link from node's parent
1249	if ( leaf->parent )
1250	leaf->parent->ReplaceChildLink(leaf, node);
1251	// and setup child links
1252	node->SetupChildLinks(back, front);
1253
1254	back->bbox = leaf->bbox;
1255	front->bbox = leaf->bbox;
1256	back->dirBBox = leaf->dirBBox;
1257	front->dirBBox = leaf->dirBBox;
1258
1259	if (axis <= RssTreeNode::SPLIT_Z) {
1260	back->bbox.SetMax(axis, info.position);
1261	front->bbox.SetMin(axis, info.position);
1262	} else {
1263	back->dirBBox.SetMax(axis-3, info.position);
1264	front->dirBBox.SetMin(axis-3, info.position);
1265	}
1266
1267	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1268	ri != leaf->rays.end();
1269	ri++) {
1270	if ((*ri).mRay->IsActive()) {
1271
1272	// first unref ray from the former leaf
1273	(*ri).mRay->Unref();
1274
1275	// Debug << "computed t: " << (*ri).mRay->mT << endl;
1276	// determine the side of this ray with respect to the plane
1277	int side = node->ComputeRaySide(*ri);
1278
1279	if (side == 1)
1280	front->AddRay(*ri);
1281	else
1282	back->AddRay(*ri);
1283	} else
1284	(*ri).mRay->Unref();
1285	}
1286
1287	// distribute the total number of rays according to the distribution
1288	// of rays which remained
1289	// front->mTotalRays = front->rays.size()*leaf->mTotalRays/leaf->rays.size();
1290	// back->mTotalRays = back->rays.size()*leaf->mTotalRays/leaf->rays.size();
1291
1292	#if 0
1293	front->SetPvsSize(pvsFront);
1294	back->SetPvsSize(pvsBack);
1295	// compute entropy as well
1296	front->ComputeEntropyImportance();
1297	back->ComputeEntropyImportance();
1298	#else
1299	UpdatePvsSize(front);
1300	UpdatePvsSize(back);
1301	#endif
1302
1303
1304	// update stats
1305	stat.rayRefs -= (int)leaf->rays.size();
1306	stat.rayRefs += back->rays.size() + front->rays.size();
1307
1308
1309	delete leaf;
1310	return node;
1311	}
1312
1313
1314
1315
1316
1317
1318	int
1319	RssTree::ReleaseMemory(const int time)
1320	{
1321	stack<RssTreeNode *> tstack;
1322
1323	// find a node in the tree which subtree will be collapsed
1324	int maxAccessTime = time - accessTimeThreshold;
1325	int released;
1326
1327	PushRoots(tstack);
1328
1329	while (!tstack.empty()) {
1330	RssTreeNode *node = tstack.top();
1331	tstack.pop();
1332
1333
1334	if (!node->IsLeaf()) {
1335	RssTreeInterior in = (RssTreeInterior )node;
1336	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1337	if (in->depth >= minCollapseDepth &&
1338	in->lastAccessTime <= maxAccessTime) {
1339	released = CollapseSubtree(node, time);
1340	break;
1341	}
1342
1343	if (in->back->GetAccessTime() <
1344	in->front->GetAccessTime()) {
1345	tstack.push(in->front);
1346	tstack.push(in->back);
1347	} else {
1348	tstack.push(in->back);
1349	tstack.push(in->front);
1350	}
1351	}
1352	}
1353
1354	while (tstack.empty()) {
1355	// could find node to collaps...
1356	// cout<<"Could not find a node to release "<<endl;
1357	break;
1358	}
1359
1360	return released;
1361	}
1362
1363
1364
1365
1366	RssTreeNode *
1367	RssTree::SubdivideLeaf(
1368	RssTreeLeaf *leaf
1369	)
1370	{
1371	RssTreeNode *node = leaf;
1372
1373	AxisAlignedBox3 leafBBox = GetBBox(leaf);
1374
1375	static int pass = 0;
1376	pass ++;
1377
1378	// check if we should perform a dynamic subdivision of the leaf
1379	if (!TerminationCriteriaSatisfied(leaf)) {
1380
1381	// memory check and realese...
1382	if (GetMemUsage() > maxTotalMemory) {
1383	ReleaseMemory( pass );
1384	}
1385
1386	AxisAlignedBox3 backBBox, frontBBox;
1387
1388	// subdivide the node
1389	node =
1390	SubdivideNode(leaf,
1391	leafBBox,
1392	backBBox,
1393	frontBBox
1394	);
1395	}
1396
1397	return node;
1398	}
1399
1400
1401
1402	void
1403	RssTree::UpdateRays(
1404	VssRayContainer &remove,
1405	VssRayContainer &add
1406	)
1407	{
1408	RssTreeLeaf::NewMail();
1409
1410	// schedule rays for removal
1411	for(VssRayContainer::const_iterator ri = remove.begin();
1412	ri != remove.end();
1413	ri++) {
1414	(*ri)->ScheduleForRemoval();
1415	}
1416
1417	int inactive=0;
1418
1419	for(VssRayContainer::const_iterator ri = remove.begin();
1420	ri != remove.end();
1421	ri++) {
1422	if ((*ri)->ScheduledForRemoval())
1423	// RemoveRay(*ri, NULL, false);
1424	// !!! BUG - with true it does not work correctly - aggreated delete
1425	RemoveRay(*ri, NULL, true);
1426	else
1427	inactive++;
1428	}
1429
1430
1431	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1432
1433	for(VssRayContainer::const_iterator ri = add.begin();
1434	ri != add.end();
1435	ri++) {
1436	RssTreeNode::RayInfo info(*ri);
1437	// if (mForcedBoundingBox==NULL \|\| ClipRay(info, bbox))
1438	AddRay(info);
1439	}
1440
1441	stat.rayRefs += add.size() - remove.size();
1442
1443	UpdateTreeStatistics();
1444	// check whether the tree should be prunned
1445	if (stat.rayRefs > mMaxRays) {
1446	PruneRays(mMaxRays);
1447	// UpdateTreeStatistics();
1448	}
1449
1450	}
1451
1452
1453
1454
1455	void
1456	RssTree::RemoveRay(VssRay *ray,
1457	vector<RssTreeLeaf > affectedLeaves,
1458	const bool removeAllScheduledRays
1459	)
1460	{
1461
1462	stack<RayTraversalData> tstack;
1463
1464	PushRoots(tstack, RssTreeLeaf::RayInfo(ray));
1465
1466	RayTraversalData data;
1467
1468	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1469
1470	while (!tstack.empty()) {
1471	data = tstack.top();
1472	tstack.pop();
1473
1474	if (!data.node->IsLeaf()) {
1475	// split the set of rays in two groups intersecting the
1476	// two subtrees
1477
1478	TraverseInternalNode(data, tstack);
1479
1480	} else {
1481	// remove the ray from the leaf
1482	// find the ray in the leaf and swap it with the last ray...
1483	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1484
1485	if (!leaf->Mailed()) {
1486	leaf->Mail();
1487	if (affectedLeaves)
1488	affectedLeaves->push_back(leaf);
1489
1490	if (removeAllScheduledRays) {
1491	int tail = (int)leaf->rays.size()-1;
1492
1493	for (int i=0; i < (int)(leaf->rays.size()); i++) {
1494	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1495	// find a ray to replace it with
1496	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1497	stat.removedRayRefs++;
1498	leaf->rays[tail].mRay->Unref();
1499	leaf->rays.pop_back();
1500	tail--;
1501	}
1502
1503	if (tail < i)
1504	break;
1505
1506	stat.removedRayRefs++;
1507	leaf->rays[i].mRay->Unref();
1508	leaf->rays[i] = leaf->rays[tail];
1509	leaf->rays.pop_back();
1510	tail--;
1511	}
1512	}
1513	}
1514	}
1515
1516	if (!removeAllScheduledRays)
1517	for (int i=0; i < (int)leaf->rays.size(); i++) {
1518	if (leaf->rays[i].mRay == ray) {
1519	stat.removedRayRefs++;
1520	ray->Unref();
1521	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1522	leaf->rays.pop_back();
1523	// check this ray again
1524	break;
1525	}
1526	}
1527
1528	}
1529	}
1530
1531	if (ray->RefCount() != 0) {
1532	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1533	exit(1);
1534	}
1535
1536	}
1537
1538
1539	void
1540	RssTree::AddRay(RssTreeNode::RayInfo &info)
1541	{
1542
1543	stack<RayTraversalData> tstack;
1544
1545	RssTreeNode *root = GetRoot(info.GetSourceObject());
1546	tstack.push(RayTraversalData(root, info));
1547
1548	RayTraversalData data;
1549
1550	while (!tstack.empty()) {
1551	data = tstack.top();
1552	tstack.pop();
1553
1554	if (!data.node->IsLeaf()) {
1555	TraverseInternalNode(data, tstack);
1556	} else {
1557	// remove the ray from the leaf
1558	// find the ray in the leaf and swap it with the last ray...
1559	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
1560	leaf->AddRay(data.rayData);
1561	stat.addedRayRefs++;
1562	}
1563	}
1564	}
1565
1566	void
1567	RssTree::TraverseInternalNode(
1568	RayTraversalData &data,
1569	stack<RayTraversalData> &tstack)
1570	{
1571	RssTreeInterior in = (RssTreeInterior ) data.node;
1572
1573	// determine the side of this ray with respect to the plane
1574	int side = in->ComputeRaySide(data.rayData
1575	);
1576
1577	if (side == 1)
1578	tstack.push(RayTraversalData(in->front, data.rayData));
1579	else
1580	tstack.push(RayTraversalData(in->back, data.rayData));
1581
1582	}
1583
1584
1585	int
1586	RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
1587	{
1588	// first count all rays in the subtree
1589	// use mail 1 for this purpose
1590	stack<RssTreeNode *> tstack;
1591	int rayCount = 0;
1592	int totalRayCount = 0;
1593	int collapsedNodes = 0;
1594
1595	#if DEBUG_COLLAPSE
1596	cout<<"Collapsing subtree"<<endl;
1597	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1598	cout<<"depth="<<(int)sroot->depth<<endl;
1599	#endif
1600
1601	// tstat.collapsedSubtrees++;
1602	// tstat.collapseDepths += (int)sroot->depth;
1603	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1604
1605	tstack.push(sroot);
1606	VssRay::NewMail();
1607
1608	while (!tstack.empty()) {
1609	collapsedNodes++;
1610	RssTreeNode *node = tstack.top();
1611	tstack.pop();
1612
1613	if (node->IsLeaf()) {
1614	RssTreeLeaf leaf = (RssTreeLeaf ) node;
1615	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1616	ri != leaf->rays.end();
1617	ri++) {
1618
1619	totalRayCount++;
1620	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
1621	(*ri).mRay->Mail();
1622	rayCount++;
1623	}
1624	}
1625	} else {
1626	tstack.push(((RssTreeInterior *)node)->back);
1627	tstack.push(((RssTreeInterior *)node)->front);
1628	}
1629	}
1630
1631	VssRay::NewMail();
1632
1633	// create a new node that will hold the rays
1634	RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
1635	if ( newLeaf->parent )
1636	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1637
1638
1639	tstack.push( sroot );
1640
1641	while (!tstack.empty()) {
1642
1643	RssTreeNode *node = tstack.top();
1644	tstack.pop();
1645
1646	if (node->IsLeaf()) {
1647	RssTreeLeaf leaf = (RssTreeLeaf ) node;
1648
1649	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1650	ri != leaf->rays.end();
1651	ri++) {
1652
1653	// unref this ray from the old node
1654
1655	if ((*ri).mRay->IsActive()) {
1656	(*ri).mRay->Unref();
1657	if (!(*ri).mRay->Mailed()) {
1658	(*ri).mRay->Mail();
1659	newLeaf->AddRay(*ri);
1660	}
1661	} else
1662	(*ri).mRay->Unref();
1663
1664	}
1665	} else {
1666	tstack.push(((RssTreeInterior *)node)->back);
1667	tstack.push(((RssTreeInterior *)node)->front);
1668	}
1669	}
1670
1671	// delete the node and all its children
1672	delete sroot;
1673
1674	// for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1675	// ri != newLeaf->rays.end();
1676	// ri++)
1677	// (*ri).ray->UnMail(2);
1678
1679
1680	#if DEBUG_COLLAPSE
1681	cout<<"Total memory before="<<GetMemUsage()<<endl;
1682	#endif
1683
1684	stat.nodes -= collapsedNodes - 1;
1685	stat.leaves -= collapsedNodes/2 - 1;
1686	stat.rayRefs -= totalRayCount - rayCount;
1687
1688	#if DEBUG_COLLAPSE
1689	cout<<"collapsed nodes"<<collapsedNodes<<endl;
1690	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1691	cout<<"Total memory after="<<GetMemUsage()<<endl;
1692	cout<<"================================"<<endl;
1693	#endif
1694
1695	// tstat.collapsedNodes += collapsedNodes;
1696	// tstat.collapsedRays += totalRayCount - rayCount;
1697
1698	return totalRayCount - rayCount;
1699	}
1700
1701
1702	int
1703	RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
1704	{
1705	stack<RssTreeNode *> tstack;
1706	PushRoots(tstack);
1707
1708	Intersectable::NewMail();
1709	int pvsSize = 0;
1710
1711	while (!tstack.empty()) {
1712	RssTreeNode *node = tstack.top();
1713	tstack.pop();
1714
1715
1716	if (node->IsLeaf()) {
1717	RssTreeLeaf leaf = (RssTreeLeaf )node;
1718	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1719	ri != leaf->rays.end();
1720	ri++)
1721	if ((*ri).mRay->IsActive()) {
1722	Intersectable *object;
1723	object = (*ri).GetObject();
1724	if (object && !object->Mailed()) {
1725	pvsSize++;
1726	object->Mail();
1727	}
1728	}
1729	} else {
1730	RssTreeInterior in = (RssTreeInterior )node;
1731	if (in->axis < 3) {
1732	if (box.Max(in->axis) >= in->position )
1733	tstack.push(in->front);
1734
1735	if (box.Min(in->axis) <= in->position )
1736	tstack.push(in->back);
1737	} else {
1738	// both nodes for directional splits
1739	tstack.push(in->front);
1740	tstack.push(in->back);
1741	}
1742	}
1743	}
1744	return pvsSize;
1745	}
1746
1747	int
1748	RssTree::CollectPvs(const AxisAlignedBox3 &box,
1749	ObjectContainer &pvs
1750	) const
1751	{
1752	stack<RssTreeNode *> tstack;
1753	PushRoots(tstack);
1754
1755	Intersectable::NewMail();
1756
1757	while (!tstack.empty()) {
1758	RssTreeNode *node = tstack.top();
1759	tstack.pop();
1760
1761
1762	if (node->IsLeaf()) {
1763	RssTreeLeaf leaf = (RssTreeLeaf )node;
1764	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1765	ri != leaf->rays.end();
1766	ri++)
1767	if ((*ri).mRay->IsActive()) {
1768	Intersectable *object;
1769	object = (*ri).GetObject();
1770	if (object && !object->Mailed()) {
1771	pvs.push_back(object);
1772	object->Mail();
1773	}
1774	}
1775	} else {
1776	RssTreeInterior in = (RssTreeInterior )node;
1777	if (in->axis < 3) {
1778	if (box.Max(in->axis) >= in->position )
1779	tstack.push(in->front);
1780
1781	if (box.Min(in->axis) <= in->position )
1782	tstack.push(in->back);
1783	} else {
1784	// both nodes for directional splits
1785	tstack.push(in->front);
1786	tstack.push(in->back);
1787	}
1788	}
1789	}
1790	return pvs.size();
1791	}
1792
1793	void
1794	RssTree::UpdateTreeStatistics()
1795	{
1796	stack<RssTreeNode *> tstack;
1797	PushRoots(tstack);
1798
1799	float sumPvsSize = 0.0f;
1800	float sumRayContribution = 0.0f;
1801	float sumWeightedRayContribution = 0.0f;
1802	float sumImportance = 0.0f;
1803	float sumPvsEntropy = 0.0f;
1804	float sumRayLengthEntropy = 0.0f;
1805	float sumRays = 0.0f;
1806
1807	int leaves = 0;
1808
1809	for (int i=0; i < mRoots.size(); i++)
1810	UpdateImportance(mRoots[i]);
1811
1812	while (!tstack.empty()) {
1813	RssTreeNode *node = tstack.top();
1814	tstack.pop();
1815
1816	if (node->IsLeaf()) {
1817	leaves++;
1818	RssTreeLeaf leaf = (RssTreeLeaf )node;
1819	// updated above already
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	} else {
1845	RssTreeInterior in = (RssTreeInterior )node;
1846	tstack.push(in->front);
1847	tstack.push(in->back);
1848	}
1849	}
1850
1851	stat.avgPvsSize = sumPvsSize/(float)leaves;
1852	stat.avgRays = sumRays/(float)leaves;
1853	stat.avgRayContribution = sumRayContribution/(float)leaves;
1854	// avgPvsEntropy = sumPvsEntropy/(float)leaves;
1855	// avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
1856	stat.avgImportance = sumImportance/(float)leaves;
1857	stat.avgWeightedRayContribution = sumWeightedRayContribution/(float)sumRays;
1858	stat.rayRefs = (int)sumRays;
1859
1860	}
1861
1862	void
1863	RssTree::UpdateImportance(
1864	RssTreeNode *node)
1865	{
1866	if (node->IsLeaf()) {
1867	UpdatePvsSize((RssTreeLeaf *)node);
1868	} else {
1869	RssTreeInterior in = (RssTreeInterior )node;
1870	UpdateImportance(in->front);
1871	UpdateImportance(in->back);
1872	node->mImportance = in->front->mImportance + in->back->mImportance;
1873	}
1874	}
1875
1876	// generate a single ray described by parameters in a unit cube
1877	void
1878	RssTree::GenerateRay(float *params, SimpleRayContainer &rays)
1879	{
1880	RssTreeNode *node;
1881
1882	// works only with single root now!
1883	node = mRoots[0];
1884
1885	while (!node->IsLeaf()) {
1886	RssTreeInterior in = (RssTreeInterior )node;
1887	// decide whether to go left or right
1888	int axis = in->axis;
1889	// float p = in->GetRelativePosition();
1890	float p = in->back->GetImportance()/(in->back->GetImportance() + in->front->GetImportance());
1891
1892	// cout<<"axis="<<axis<<" pos="<<pos<<endl;
1893	// cout<<GetBBox(in)<<endl;
1894	// cout<<GetDirBBox(in)<<endl;
1895	// cout<<"********"<<endl;
1896
1897	if (params[axis] < p) {
1898	node = in->back;
1899	params[axis] = params[axis] / p;
1900	} else {
1901	node = in->front;
1902	params[axis] = (params[axis] - p) / ( 1.0f - p);
1903	}
1904	}
1905
1906	// cout<<params[0]<<" "<<params[1]<<" "<<params[2]<<" "<<params[3]<<" "<<params[4]<<endl;
1907
1908	GenerateLeafRay(params, rays, node->bbox, node->dirBBox);
1909	}
1910
1911	int
1912	RssTree::GenerateRays(
1913	const float ratioPerLeaf,
1914	SimpleRayContainer &rays)
1915	{
1916
1917
1918
1919
1920	stack<RssTreeNode *> tstack;
1921	PushRoots(tstack);
1922
1923	while (!tstack.empty()) {
1924	RssTreeNode *node = tstack.top();
1925	tstack.pop();
1926
1927	if (node->IsLeaf()) {
1928	RssTreeLeaf leaf = (RssTreeLeaf )node;
1929	float c = leaf->GetImportance();
1930	int num = int(c*ratioPerLeaf + 0.5f);
1931	// cout<<num<<" ";
1932
1933	for (int i=0; i < num; i++) {
1934	Vector3 origin = GetBBox(leaf).GetRandomPoint();
1935	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1936	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1937	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1938	rays.push_back(SimpleRay(origin, direction));
1939	}
1940	} else {
1941	RssTreeInterior in = (RssTreeInterior )node;
1942	// both nodes for directional splits
1943	tstack.push(in->front);
1944	tstack.push(in->back);
1945	}
1946	}
1947
1948	return (int)rays.size();
1949	}
1950
1951	void
1952	RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
1953	{
1954	stack<RssTreeNode *> tstack;
1955	PushRoots(tstack);
1956
1957	while (!tstack.empty()) {
1958	RssTreeNode *node = tstack.top();
1959	tstack.pop();
1960
1961	if (node->IsLeaf()) {
1962	RssTreeLeaf leaf = (RssTreeLeaf )node;
1963	leaves.push_back(leaf);
1964	} else {
1965	RssTreeInterior in = (RssTreeInterior )node;
1966	// both nodes for directional splits
1967	tstack.push(in->front);
1968	tstack.push(in->back);
1969	}
1970	}
1971	}
1972
1973	bool
1974	RssTree::ValidLeaf(RssTreeLeaf *leaf) const
1975	{
1976	return true;
1977	//return leaf->rays.size() > termMinRays/4;
1978	}
1979
1980
1981	void
1982	RssTree::PickEdgeRays(RssTreeLeaf *leaf,
1983	int &r1,
1984	int &r2
1985	)
1986	{
1987	int nrays = leaf->rays.size();
1988
1989	if (nrays == 2) {
1990	r1 = 0;
1991	r2 = 1;
1992	return;
1993	}
1994
1995	#if 1
1996	int tries = min(20, nrays);
1997
1998	while (--tries) {
1999	r1 = Random(nrays);
2000	r2 = Random(nrays);
2001	if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
2002	break;
2003	}
2004
2005	if (r1 == r2)
2006	r2 = (r1+1)%leaf->rays.size();
2007
2008	#else
2009	// pick a random ray
2010	int base = Random(nrays);
2011	RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
2012	Intersectable *baseObject = baseRay->GetObject();
2013
2014	// and a random 5D derivative which will be used to find the minimal projected distances
2015	Vector3 spatialDerivative;
2016	Vector3 directionalDerivative;
2017
2018	while (1) {
2019	spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2020	RandomValue(-1.0f, 1.0f),
2021	RandomValue(-1.0f, 1.0f));
2022	float m = Magnitude(spatialDerivative);
2023	if (m != 0) {
2024	spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
2025	break;
2026	}
2027	}
2028
2029	while (1) {
2030	directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2031	RandomValue(-1.0f, 1.0f),
2032	0.0f);
2033	float m = Magnitude(directionalDerivative);
2034	if (m != 0) {
2035	directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
2036	break;
2037	}
2038	}
2039
2040	// now find the furthest sample from the same object and the closest from a different object
2041	int i = 0;
2042	float minDist = MAX_FLOAT;
2043	float maxDist = -MAX_FLOAT;
2044	r1 = base;
2045	r2 = base;
2046	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2047	ri != leaf->rays.end();
2048	ri++, i++) {
2049	float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
2050	*ri,
2051	spatialDerivative,
2052	directionalDerivative
2053	);
2054
2055	if ((*ri).GetObject() == baseObject) {
2056	if (dist > maxDist) {
2057	maxDist = dist;
2058	r1 = i;
2059	}
2060	} else {
2061	if (dist > 0.0f && dist < minDist) {
2062	minDist = dist;
2063	r2 = i;
2064	}
2065	}
2066	}
2067
2068	#endif
2069	}
2070
2071
2072	struct RayNeighbor {
2073	int rayInfo;
2074	float distance;
2075	RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
2076	};
2077
2078	void
2079	RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
2080	vector<RssTreeLeaf::SilhouetteRays> &rays
2081	)
2082	{
2083	// for every leaf find its neares neighbor from a different object
2084	vector<RayNeighbor> neighbors(leaf->rays.size());
2085	int i, j;
2086	for (i=0; i < leaf->rays.size(); i++)
2087	for (j=0; j < leaf->rays.size(); j++)
2088	if (i!=j) {
2089	float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
2090	if (d < neighbors[i].distance) {
2091	neighbors[i].distance = d;
2092	neighbors[i].rayInfo = j;
2093	}
2094	}
2095
2096	// now check which are the pairs of nearest neighbors
2097	for (i=0; i < leaf->rays.size(); i++) {
2098	int j = neighbors[i].rayInfo;
2099	if (neighbors[j].rayInfo == i) {
2100	// this is a silhouette edge pair
2101	if (i < j) {
2102	// generate an silhoutte ray pair
2103	rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
2104	&leaf->rays[j]));
2105	}
2106	} else {
2107	// this is not a silhouette ray - delete???
2108	}
2109	}
2110
2111	}
2112
2113
2114	bool
2115	GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
2116	const int index,
2117	int &i1,
2118	int &i2,
2119	int &i3)
2120	{
2121	int found = 0;
2122	int indices[3];
2123
2124	int size = rays.size();
2125	// use russian roulete selection for the tripple
2126	// number of free positions for the bullet
2127	int positions = size - 1;
2128	int i;
2129	for (i=0; i < size; i++) {
2130	if (rays[i]->mRay->Mailed())
2131	positions--;
2132	}
2133
2134	if (positions < 3)
2135	return false;
2136
2137	for (i=0; i < size; i++) {
2138	if (i != index && !rays[i]->mRay->Mailed()) {
2139	float p = (3 - found)/(float)positions;
2140	if (Random(1.0f) < p) {
2141	indices[found] = i;
2142	found++;
2143	}
2144	}
2145	positions--;
2146	}
2147	return true; // corr. matt
2148	}
2149
2150	bool
2151	RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
2152	const RssTreeNode::RayInfo &r1,
2153	const RssTreeNode::RayInfo &r2,
2154	const RssTreeNode::RayInfo &r3)
2155	{
2156
2157
2158	return false;
2159	}
2160
2161	int
2162	RssTree::RemoveInteriorRays(
2163	RssTreeLeaf *leaf
2164	)
2165	{
2166	#if 1
2167	// first collect all objects refered in this leaf
2168	map<Intersectable , vector<RssTreeLeaf::RayInfo > > rayMap;
2169
2170	RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
2171	for (; it != leaf->rays.end(); ++it) {
2172	Intersectable object = (it).GetObject();
2173
2174	rayMap[object].push_back(&(*it));
2175	// vector<RayInfo > data = rayMap.Find(object);
2176	// if (data) {
2177	// data->push_back(&(*it));
2178	// } else {
2179	// // rayMap[object] = vector<RayInfo *>;
2180	// rayMap[object].push_back(&(*it));
2181	// }
2182	}
2183
2184	// now go through all objects
2185	map<Intersectable , vector<RssTreeLeaf::RayInfo > >::iterator mi;
2186
2187	// infos of mailed rays are scheduled for removal
2188	VssRay::NewMail();
2189	for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
2190	vector<RssTreeLeaf::RayInfo > &rays = (mi).second;
2191	vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
2192	int rayIndex = 0;
2193
2194	for (; ri != rays.end(); ++ri, ++rayIndex) {
2195	RssTreeNode::RayInfo ray = ri;
2196	int tries = rays.size();
2197	for (int i = 0; i < tries; i++) {
2198	int r1, r2, r3;
2199	GetRandomTripple(rays,
2200	rayIndex,
2201	r1,
2202	r2,
2203	r3);
2204	if (IsRayConvexCombination(*ray,
2205	*rays[r1],
2206	*rays[r2],
2207	*rays[r3])) {
2208	ray->mRay->Mail();
2209	}
2210	}
2211	}
2212	}
2213
2214
2215	#endif
2216	return 0;
2217	}
2218
2219
2220	void
2221	RssTree::GenerateLeafRaysSimple(RssTreeLeaf *leaf,
2222	const int numberOfRays,
2223	SimpleRayContainer &rays)
2224	{
2225
2226	if (numberOfRays == 0)
2227	return;
2228
2229	int startIndex = (int)rays.size();
2230
2231	// Debug<<"B"<<flush;
2232
2233	AxisAlignedBox3 box = GetBBox(leaf);
2234	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2235
2236	float boxSize = Magnitude(box.Diagonal());
2237
2238	const float smoothRange = 0.1f;
2239	if (smoothRange != 0.0f) {
2240	box.Scale(1.0f + smoothRange);
2241	dirBox.Scale(1.0f + smoothRange);
2242	}
2243
2244	// Debug<<"R"<<flush;
2245
2246	for (int i=0; i < numberOfRays; i++) {
2247	// Debug<<i<<flush;
2248	float r[5];
2249
2250	r[0] = RandomValue(0.0f, 1.0f);
2251	r[1] = RandomValue(0.0f, 1.0f);
2252	r[2] = RandomValue(0.0f, 1.0f);
2253	r[3] = RandomValue(0.0f, 1.0f);
2254	r[4] = RandomValue(0.0f, 1.0f);
2255
2256	GenerateLeafRay(r, rays, box, dirBox);
2257	}
2258
2259	// Debug<<"D"<<flush;
2260
2261	float density = 1.0f;
2262	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2263	for (int i=startIndex; i < rays.size(); i++) {
2264	rays[i].mPdf = p;
2265	}
2266	}
2267
2268	void
2269	RssTree::GenerateLeafRay(
2270	float *params,
2271	SimpleRayContainer &rays,
2272	const AxisAlignedBox3 &box,
2273	const AxisAlignedBox3 &dirBBox
2274	)
2275	{
2276	Vector3 origin = bbox.GetPoint(Vector3(params[0], params[1], params[2]));
2277	Vector3 dirVector = dirBBox.GetPoint(Vector3(params[3], params[4], 0.0f));
2278	Vector3 direction = VssRay::GetInvDirParam(dirVector.x, dirVector.y);
2279
2280	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2281	float dist = 0.0f;
2282	float minT, maxT;
2283
2284	float boxSize = Magnitude(bbox.Size());
2285	// compute interection of the ray with the box
2286	#if 1
2287	if (bbox.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2288	dist = (maxT + 1e-2*boxSize);
2289	#else
2290	dist = 0.5f*boxSize;
2291	#endif
2292
2293	origin += direction*dist;
2294
2295
2296	//Debug<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2297	rays.push_back(SimpleRay(origin, direction));
2298	}
2299
2300	void
2301	RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
2302	const int numberOfRays,
2303	SimpleRayContainer &rays)
2304	{
2305
2306	if (numberOfRays == 0)
2307	return;
2308
2309
2310	int nrays = (int)leaf->rays.size();
2311	int startIndex = (int)rays.size();
2312
2313	AxisAlignedBox3 box = GetBBox(leaf);
2314	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2315
2316
2317	const float smoothRange = 0.1f;
2318	if (smoothRange != 0.0f) {
2319	box.Scale(1.0f + smoothRange);
2320	dirBox.Scale(1.0f + smoothRange);
2321	}
2322
2323	Exporter *exporter = NULL;
2324	VssRayContainer selectedRays;
2325	static int counter = 0;
2326	if (counter < 0) {
2327	char filename[128];
2328	sprintf(filename, "selected-rays%04d.x3d", counter);
2329	exporter = Exporter::GetExporter(filename);
2330	// exporter->SetWireframe();
2331	// exporter->ExportKdTree(*mKdTree);
2332	exporter->SetWireframe();
2333	// exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
2334	exporter->SetWireframe();
2335	exporter->ExportBox(box);
2336	exporter->SetFilled();
2337	counter++;
2338	}
2339
2340
2341	float extendedConvexCombinationProb = 0.0f; //0.7f
2342	// float silhouetteCheckPercentage = 1.0f; //0.5f
2343	float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
2344	float silhouetteObjectCheckPercentage = 0.8f;
2345
2346	// int numberOfTries = numberOfRays*4;
2347	int numberOfTries = numberOfRays;
2348	int generated = 0;
2349
2350	vector<Intersectable *> pvs(0);
2351	if (silhouetteCheckPercentage != 0.0f) {
2352	Intersectable::NewMail();
2353	pvs.reserve(leaf->GetPvsSize());
2354	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2355	ri != leaf->rays.end();
2356	ri++) {
2357	Intersectable object = (ri).GetObject();
2358	if (object) {
2359	if (!object->Mailed()) {
2360	object->Mail();
2361	object->mCounter = 1;
2362	pvs.push_back(object);
2363	if (exporter)
2364	exporter->ExportIntersectable(object);
2365	} else {
2366	object->mCounter++;
2367	}
2368	}
2369	}
2370	// sort objects based on their mCounter
2371	sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
2372	}
2373
2374	for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
2375	bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
2376	extendedConvexCombinationProb));
2377
2378	Vector3 origin, direction;
2379	// generate half of convex combination and half of random rays
2380	if (useExtendedConvexCombination) {
2381	// pickup 2 random rays
2382	int r1, r2;
2383	PickEdgeRays(leaf, r1, r2);
2384
2385
2386	Vector3 o1 = leaf->rays[r1].GetOrigin();
2387	Vector3 o2 = leaf->rays[r2].GetOrigin();
2388
2389	const float overlap = 0.0f;
2390	float w1, w2;
2391	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2392	origin = w1o1 + w2o2;
2393	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2394	direction =
2395	w1*leaf->rays[r1].GetDir() +
2396	w2*leaf->rays[r2].GetDir();
2397
2398	// shift the origin a little bit
2399	direction.Normalize();
2400
2401	} else {
2402	Vector3 dirVector;
2403
2404	Vector3 pVector;
2405	Vector3 dVector;
2406
2407	pVector = Vector3(RandomValue(0.0f, 1.0f),
2408	RandomValue(0.0f, 1.0f),
2409	RandomValue(0.0f, 1.0f));
2410
2411	dVector = Vector3(RandomValue(0.0f, 1.0f),
2412	RandomValue(0.0f, 1.0f),
2413	0.0f);
2414
2415
2416	origin = box.GetPoint(pVector);
2417	dirVector = dirBox.GetPoint(dVector);
2418
2419	direction = VssRay::GetInvDirParam(dirVector.x,dirVector.y);
2420	}
2421
2422
2423	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2424	float dist = 0.0f;
2425	float minT, maxT;
2426	// compute interection of the ray with the box
2427
2428	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2429	dist = maxT;
2430
2431
2432	origin += direction*dist;
2433
2434	bool intersects = false;
2435
2436	if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
2437	if (exporter) {
2438	VssRay ray = new VssRay(origin, origin + 100direction, NULL, NULL);
2439	selectedRays.push_back(ray);
2440	}
2441
2442	// check whether the ray does not intersect already visible objects
2443	Ray traversalRay;
2444	traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2445	for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2446	oi != pvs.end();
2447	oi++) {
2448	Intersectable object = oi;
2449	// do not test every object
2450	if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
2451	object->CastRay(traversalRay) ) {
2452	intersects = true;
2453	break;
2454	}
2455	}
2456	}
2457
2458
2459	if (!intersects) {
2460	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2461	rays.push_back(SimpleRay(origin, direction));
2462	generated++;
2463	}
2464	}
2465
2466	// cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
2467	// assign pdfs to the generated rays
2468	// float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
2469	float density = 1.0f;
2470
2471	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2472	for (int i=startIndex; i < rays.size(); i++) {
2473	rays[i].mPdf = p;
2474	}
2475
2476
2477	if (exporter) {
2478	exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2479	delete exporter;
2480	CLEAR_CONTAINER(selectedRays);
2481	}
2482
2483	}
2484
2485	int
2486	RssTree::PruneRays(RssTreeLeaf *leaf,
2487	const float contributionThreshold)
2488	{
2489	int i;
2490	int j;
2491
2492	for (j=0, i=0; i < leaf->rays.size(); i++) {
2493
2494	if (leaf->rays[i].mRay->mWeightedPvsContribution > contributionThreshold) {
2495	// copy a valid sample
2496	leaf->rays[j] = leaf->rays[i];
2497	j++;
2498	} else {
2499	// delete the ray
2500	leaf->rays[i].mRay->Unref();
2501	if (leaf->rays[i].mRay->RefCount() != 0) {
2502	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2503	exit(1);
2504	}
2505	delete leaf->rays[i].mRay;
2506	}
2507	}
2508
2509
2510	leaf->rays.resize(j);
2511	int removed = (i-j);
2512	stat.rayRefs -= removed;
2513	return removed;
2514	}
2515
2516	int
2517	RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2518	const float ratio)
2519	{
2520	int i;
2521	int j;
2522
2523	for (j=0, i=0; i < leaf->rays.size(); i++) {
2524
2525	if (Random(1.0f) < ratio) {
2526	// copy a valid sample
2527	if (i!=j)
2528	leaf->rays[j] = leaf->rays[i];
2529	j++;
2530	} else {
2531	// delete the ray
2532	leaf->rays[i].mRay->Unref();
2533	if (leaf->rays[i].mRay->RefCount() != 0) {
2534	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2535	exit(1);
2536	}
2537	delete leaf->rays[i].mRay;
2538	}
2539	}
2540
2541
2542	leaf->rays.resize(j);
2543	int removed = (i-j);
2544	stat.rayRefs -= removed;
2545	return removed;
2546	}
2547
2548	int
2549	RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2550	const float ratio)
2551	{
2552	int i;
2553
2554	if (leaf->rays.size() == 0)
2555	return 0;
2556
2557	sort(leaf->rays.begin(),
2558	leaf->rays.end(),
2559	RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2560
2561	int desired = int(ratio*leaf->rays.size());
2562	int removed = (int)leaf->rays.size() - desired;
2563
2564	for (i=desired; i < (int)leaf->rays.size(); i++) {
2565	// delete the ray
2566	leaf->rays[i].mRay->Unref();
2567	if (leaf->rays[i].mRay->RefCount() != 0) {
2568	cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
2569	cerr<<"desired="<<desired<<"i="<<i<<" size="<<(int)leaf->rays.size()<<endl;
2570	exit(1);
2571	}
2572	delete leaf->rays[i].mRay;
2573	}
2574
2575	leaf->rays.resize(desired);
2576	stat.rayRefs -= removed;
2577	return removed;
2578	}
2579
2580
2581	int
2582	RssTree::PruneRays(
2583	const int maxRays
2584	)
2585	{
2586	bool globalPrunning = false;
2587
2588	stack<RssTreeNode *> tstack;
2589	int prunned = 0;
2590
2591	// prune more rays to amortize the procedure
2592	int desired = (int)(maxRays * 0.8f);
2593
2594	Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;
2595
2596	if (globalPrunning) {
2597	VssRayContainer allRays;
2598	allRays.reserve(stat.rayRefs);
2599	CollectRays(allRays);
2600	sort(allRays.begin(),
2601	allRays.end(),
2602	GreaterWeightedPvsContribution);
2603
2604	if ( maxRays >= allRays.size() )
2605	return 0;
2606
2607	float contributionThreshold = allRays[desired]->mWeightedPvsContribution;
2608
2609	PushRoots(tstack);
2610
2611	while (!tstack.empty()) {
2612	RssTreeNode *node = tstack.top();
2613	tstack.pop();
2614
2615	if (node->IsLeaf()) {
2616	RssTreeLeaf leaf = (RssTreeLeaf )node;
2617	prunned += PruneRays(leaf, contributionThreshold);
2618	} else {
2619	RssTreeInterior in = (RssTreeInterior )node;
2620	// both nodes for directional splits
2621	tstack.push(in->front);
2622	tstack.push(in->back);
2623	}
2624	}
2625	} else {
2626	// prune random rays from each leaf so that the ratio's remain the same
2627	PushRoots(tstack);
2628
2629	if ( maxRays >= stat.rayRefs )
2630	return 0;
2631
2632	float ratio = desired/(float)stat.rayRefs;
2633
2634	while (!tstack.empty()) {
2635	RssTreeNode *node = tstack.top();
2636	tstack.pop();
2637
2638	if (node->IsLeaf()) {
2639	RssTreeLeaf leaf = (RssTreeLeaf )node;
2640	// prunned += PruneRaysRandom(leaf, ratio);
2641	prunned += PruneRaysContribution(leaf, ratio);
2642	} else {
2643	RssTreeInterior in = (RssTreeInterior )node;
2644	// both nodes for directional splits
2645	tstack.push(in->front);
2646	tstack.push(in->back);
2647	}
2648	}
2649
2650
2651
2652
2653	}
2654
2655
2656
2657	Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
2658
2659	return prunned;
2660	}
2661
2662	int
2663	RssTree::GenerateRays(const int numberOfRays,
2664	const int numberOfLeaves,
2665	SimpleRayContainer &rays)
2666	{
2667
2668	float param[5];
2669
2670	for (int i=0; i < numberOfRays; i++) {
2671	halton.GetNext(param);
2672	GenerateRay(param, rays);
2673	}
2674
2675	return rays.size();
2676
2677	vector<RssTreeLeaf *> leaves;
2678
2679	Debug<<"Collecting leaves..."<<endl<<flush;
2680	CollectLeaves(leaves);
2681	Debug<<"done."<<endl<<flush;
2682
2683
2684	if (numberOfLeaves != leaves.size()) {
2685	Debug<<"sorting leaves..."<<endl<<flush;
2686	sort(leaves.begin(),
2687	leaves.end(),
2688	GreaterContribution);
2689	Debug<<"done."<<endl<<flush;
2690	}
2691
2692
2693	Debug<<"Evaluating leaf importance..."<<endl<<flush;
2694
2695	float sumContrib = 0.0;
2696	int i;
2697	int k = 0;
2698	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2699	if (ValidLeaf(leaves[i])) {
2700	float c = leaves[i]->GetImportance();
2701	sumContrib += c;
2702	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
2703	k++;
2704	}
2705
2706	Debug<<"done."<<endl<<flush;
2707
2708	float avgContrib = sumContrib/numberOfLeaves;
2709
2710	// always generate at leat n ray per leaf
2711	int fixedPerLeaf = 1;
2712	// int fixedPerLeaf = 0;
2713	int fixed = fixedPerLeaf*(int)leaves.size();
2714
2715	// check if the number of fixed is not too big
2716
2717	if (fixed >= numberOfRays/2) {
2718	fixedPerLeaf = 0;
2719	fixed = 0;
2720	}
2721
2722	int iGenerated = numberOfRays - fixed;
2723	float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
2724
2725	k = 0;
2726
2727	Debug<<"generating leaf rays..."<<endl<<flush;
2728
2729	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2730	if (ValidLeaf(leaves[i])) {
2731	RssTreeLeaf *leaf = leaves[i];
2732	float c = leaf->GetImportance();
2733
2734	int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
2735	GenerateLeafRaysSimple(leaf, num, rays);
2736
2737	k++;
2738	}
2739
2740	Debug<<"done."<<endl<<flush;
2741
2742	return (int)rays.size();
2743	}
2744
2745
2746	float
2747	RssTree::GetAvgPvsSize()
2748	{
2749	stack<RssTreeNode *> tstack;
2750	PushRoots(tstack);
2751
2752	int sumPvs = 0;
2753	int leaves = 0;
2754	while (!tstack.empty()) {
2755	RssTreeNode *node = tstack.top();
2756	tstack.pop();
2757
2758	if (node->IsLeaf()) {
2759	RssTreeLeaf leaf = (RssTreeLeaf )node;
2760	// update pvs size
2761	UpdatePvsSize(leaf);
2762	sumPvs += leaf->GetPvsSize();
2763	leaves++;
2764	} else {
2765	RssTreeInterior in = (RssTreeInterior )node;
2766	// both nodes for directional splits
2767	tstack.push(in->front);
2768	tstack.push(in->back);
2769	}
2770	}
2771
2772
2773	return sumPvs/(float)leaves;
2774	}
2775
2776	float weightAbsContributions = ABS_CONTRIBUTION_WEIGHT;
2777	// if small very high importance of the last sample
2778	// if 1.0f then weighs = 1 1/2 1/3 1/4
2779	//float passSampleWeightDecay = 1.0f;
2780	float passSampleWeightDecay = 1.0f;
2781
2782	float
2783	RssTree::GetSampleWeight(const int pass)
2784	{
2785	int passDiff = mCurrentPass - pass;
2786	float weight;
2787	if (1)
2788	weight = 1.0f/(passDiff + passSampleWeightDecay);
2789	else
2790	switch (passDiff) {
2791	case 0:
2792	// weight = 1.0f;
2793	// break;
2794	default:
2795	weight = 1.0f;
2796	break;
2797	// case 1:
2798	// weight = 0.5f;
2799	// break;
2800	// case 2:
2801	// weight = 0.25f;
2802	// break;
2803	// case 3:
2804	// weight = 0.12f;
2805	// break;
2806	// case 4:
2807	// weight = 0.06f;
2808	// break;
2809	// default:
2810	// weight = 0.03f;
2811	// break;
2812	}
2813	return weight;
2814	}
2815
2816	void
2817	RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2818	float &weight,
2819	float &contribution)
2820	{
2821	VssRay *ray = info.mRay;
2822	weight = GetSampleWeight(mCurrentPass);
2823	contribution =
2824	weightAbsContributions*ray->mPvsContribution +
2825	(1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2826	// store the computed value
2827	info.mRay->mWeightedPvsContribution = weight*contribution;
2828	}
2829
2830
2831	void
2832	RssTree::ComputeImportance(RssTreeLeaf *leaf)
2833	{
2834	if (0)
2835	leaf->mImportance = leaf->GetAvgRayContribution();
2836	else {
2837	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2838	it_end = leaf->rays.end();
2839
2840	float sumContributions = 0.0f;
2841	float sumRelContributions = 0.0f;
2842	float sumWeights = 0.0f;
2843	float maxContribution = 0.0f;
2844	float maxRelContribution = 0.0f;
2845	float sumLength = 0;
2846
2847	for (; it != it_end; ++it) {
2848	VssRay ray = (it).mRay;
2849	float sumLength = ray->Length();
2850
2851	float weight = GetSampleWeight(ray->mPass);
2852
2853	float c1 = weight*ray->mPvsContribution;
2854	float c2 = weight*ray->mRelativePvsContribution;
2855	sumContributions += c1;
2856	if (c1 > maxContribution)
2857	maxContribution = c1;
2858
2859	//$$ 20.7. changed to sqr to pronouce higher contribution so that
2860	// they do not get smoothed!!
2861	//sumRelContributions += weight*ray->mRelativePvsContribution;
2862
2863	sumRelContributions += c2;
2864	if (c2 > maxRelContribution)
2865	maxRelContribution = c2;
2866
2867	//sumWeights += weight;
2868	sumWeights += 1.0f;
2869	}
2870
2871	float distImportance = 0.0f;
2872
2873	float spatialArea = GetBBox(leaf).SurfaceArea();
2874	float dirArea = GetDirBBox(leaf).SurfaceArea();
2875	if (leaf->rays.size()) {
2876	// compute average length of the ray
2877	float avgLength = sumLength/leaf->rays.size();
2878	// compute estimated area of the visible surface corresponding to these rays
2879	float ss = spatialArea/2.0f;
2880	float sd = dirArea;
2881	float s = ss + avgLength(2sqrt(sssd) + avgLengthsd);
2882
2883	distImportance = s/leaf->rays.size();
2884	}
2885
2886	// $$
2887	// sumWeights = leaf->mTotalRays;
2888
2889	//$$ test 30.11. 2006
2890	// normalize per volume of the cell not the number of rays!
2891	// sumWeights = spatialArea*dirArea;
2892
2893	if (sumWeights != 0.0f) {
2894	leaf->mImportance =
2895	sqr(((weightAbsContributions*sumContributions +
2896	(1.0f - weightAbsContributions)*sumRelContributions)/sumWeights));
2897
2898	// leaf->mImportance =
2899	// (weightAbsContributions*maxContribution +
2900	// (1.0f - weightAbsContributions)*maxRelContribution)/sumWeights;
2901
2902	} else
2903	leaf->mImportance = 0.0f;
2904
2905	}
2906
2907	// return GetAvgRayContribution()*mImportance;
2908	//return GetAvgRayContribution();
2909	}
2910
2911
2912	float
2913	RssTreeNode::GetImportance() const
2914	{
2915	// return sqr(mImportance);
2916	return mImportance;
2917	}
2918
2919
2920	float
2921	RssTreeLeaf::ComputePvsEntropy() const
2922	{
2923	int samples = 0;
2924	Intersectable::NewMail();
2925	// set all object as belonging to the fron pvs
2926	for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
2927	ri != rays.end();
2928	ri++)
2929	if ((*ri).mRay->IsActive()) {
2930	Intersectable object = (ri).GetObject();
2931	if (object) {
2932	if (!object->Mailed()) {
2933	object->Mail();
2934	object->mCounter = 1;
2935	} else
2936	object->mCounter++;
2937	samples++;
2938	}
2939	}
2940
2941	float entropy = 0.0f;
2942
2943	if (samples > 1) {
2944	Intersectable::NewMail();
2945	for(RayInfoContainer::const_iterator ri = rays.begin();
2946	ri != rays.end();
2947	ri++)
2948	if ((*ri).mRay->IsActive()) {
2949	Intersectable object = (ri).GetObject();
2950	if (object) {
2951	if (!object->Mailed()) {
2952	object->Mail();
2953	float p = object->mCounter/(float)samples;
2954	entropy -= p*log(p);
2955	}
2956	}
2957	}
2958	entropy = entropy/log((float)samples);
2959	}
2960	else
2961	entropy = 1.0f;
2962
2963	return entropy;
2964	}
2965
2966	float
2967	RssTreeLeaf::ComputeRayLengthEntropy() const
2968	{
2969	// get sum of all ray lengths
2970	// consider only passing rays or originating rays
2971	float sum = 0.0f;
2972	int samples = 0;
2973	int i=0;
2974	for(RayInfoContainer::const_iterator ri = rays.begin();
2975	ri != rays.end();
2976	ri++)
2977	if ((*ri).mRay->IsActive()) {
2978	// float s;
2979	// if (i == 0)
2980	// s = 200;
2981	// else
2982	// s = 100;
2983	// i++;
2984
2985	sum += (*ri).mRay->GetSize();
2986	samples++;
2987	}
2988
2989
2990	float entropy = 0.0f;
2991
2992	if (samples > 1) {
2993	i = 0;
2994	for(RayInfoContainer::const_iterator ri = rays.begin();
2995	ri != rays.end();
2996	ri++)
2997	if ((*ri).mRay->IsActive()) {
2998	// float s;
2999	// if (i==0)
3000	// s = 200;
3001	// else
3002	// s = 100;
3003	// i++;
3004	// float p = s/sum;
3005	float p = (*ri).mRay->GetSize()/sum;
3006	entropy -= p*log(p);
3007	}
3008	entropy = entropy/log((float)samples);
3009	} else
3010	entropy = 1.0f;
3011
3012	return entropy;
3013	}
3014
3015
3016
3017	float
3018	RssTreeLeaf::ComputeEntropyImportance() const
3019	{
3020
3021	// mEntropy = 1.0f - ComputeRayLengthEntropy();
3022	return 1.0f - ComputePvsEntropy();
3023	}
3024
3025	float
3026	RssTreeLeaf::ComputeRayContributionImportance() const
3027	{
3028	// mPvsEntropy = ComputePvsEntropy();
3029	// mRayLengthEntropy = ComputeRayLengthEntropy();
3030
3031	// mEntropy = 1.0f - ComputeRayLengthEntropy();
3032	return 1.0f - ComputePvsEntropy();
3033	}
3034
3035
3036	AxisAlignedBox3
3037	RssTree::GetShrankedBBox(const RssTreeNode *node) const
3038	{
3039	if (node->parent == NULL)
3040	return bbox;
3041
3042	if (!node->IsLeaf())
3043	return ((RssTreeInterior *)node)->bbox;
3044
3045	// evaluate bounding box from the ray origins
3046	AxisAlignedBox3 box;
3047	box.Initialize();
3048	RssTreeLeaf leaf = (RssTreeLeaf )node;
3049	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
3050	ri != leaf->rays.end();
3051	ri++)
3052	if ((*ri).mRay->IsActive()) {
3053	box.Include((*ri).GetOrigin());
3054	}
3055	return box;
3056	}
3057
3058
3059	int
3060	RssTree::CollectRays(VssRayContainer &rays,
3061	const int number)
3062	{
3063	VssRayContainer allRays;
3064	CollectRays(allRays);
3065
3066
3067	int desired = min(number, (int)allRays.size());
3068	float prob = desired/(float)allRays.size();
3069	//int skip = allRays.size()/desired;
3070	while (rays.size() < desired) {
3071	VssRayContainer::const_iterator it = allRays.begin();
3072	for (; it != allRays.end() && rays.size() < desired; it++) {
3073	if (Random(1.0f) < prob)
3074	rays.push_back(*it);
3075	}
3076	}
3077	return (int)rays.size();
3078	}
3079
3080
3081	int
3082	RssTree::CollectRays(VssRayContainer &rays
3083	)
3084	{
3085	VssRay::NewMail();
3086
3087	stack<RssTreeNode *> tstack;
3088	PushRoots(tstack);
3089
3090	while (!tstack.empty()) {
3091	RssTreeNode *node = tstack.top();
3092	tstack.pop();
3093
3094	if (node->IsLeaf()) {
3095	RssTreeLeaf leaf = (RssTreeLeaf )node;
3096	// update pvs size
3097	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
3098	for (;it != leaf->rays.end(); ++it)
3099	if (!(*it).mRay->Mailed()) {
3100	(*it).mRay->Mail();
3101	rays.push_back((*it).mRay);
3102	}
3103	} else {
3104	RssTreeInterior in = (RssTreeInterior )node;
3105	// both nodes for directional splits
3106	tstack.push(in->front);
3107	tstack.push(in->back);
3108	}
3109	}
3110
3111	return (int)rays.size();
3112	}
3113
3114
3115	RssTreeNode *
3116	RssTree::GetRoot(Intersectable *object) const
3117	{
3118	if (mPerObjectTree && object) {
3119	int id = object->GetId()-1;
3120	if (id < 0 \|\| id >= mRoots.size()) {
3121	Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<(int)mRoots.size()<<
3122	endl<<flush;
3123	id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
3124	}
3125	return mRoots[id];
3126	} else
3127	return mRoots[0];
3128	}
3129
3130	void
3131	RssTree::PushRoots(priority_queue<TraversalData> &st) const
3132	{
3133	for (int i=0; i < mRoots.size(); i++)
3134	st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
3135	}
3136
3137	void
3138	RssTree::PushRoots(stack<RssTreeNode *> &st) const
3139	{
3140	for (int i=0; i < mRoots.size(); i++)
3141	st.push(mRoots[i]);
3142	}
3143
3144	void
3145	RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info) const
3146	{
3147	for (int i=0; i < mRoots.size(); i++)
3148	st.push(RayTraversalData(mRoots[i], info));
3149	}
3150
3151
3152
3153
3154
3155	}

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

Download in other formats: