Context Navigation

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

Revision 1891, 75.6 KB checked in by bittner, 18 years ago (diff)
combined preprocessor

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

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

Download in other formats: