Context Navigation

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

Revision 1906, 75.8 KB checked in by mattausch, 18 years ago (diff)
fixed bug in priority computation of object space

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

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

Download in other formats: