Context Navigation

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

Revision 1884, 75.2 KB checked in by bittner, 18 years ago (diff)
temporary version, rss preprocessor not functional

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

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

Download in other formats: