Context Navigation

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

Revision 1883, 74.0 KB checked in by bittner, 18 years ago (diff)
mixture distribution initial coding

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	}
1927
1928	int
1929	RssTree::GenerateRays(
1930	const float ratioPerLeaf,
1931	SimpleRayContainer &rays)
1932	{
1933
1934	stack<RssTreeNode *> tstack;
1935	PushRoots(tstack);
1936
1937	while (!tstack.empty()) {
1938	RssTreeNode *node = tstack.top();
1939	tstack.pop();
1940
1941	if (node->IsLeaf()) {
1942	RssTreeLeaf leaf = (RssTreeLeaf )node;
1943	float c = leaf->GetImportance();
1944	int num = int(c*ratioPerLeaf + 0.5f);
1945	// cout<<num<<" ";
1946
1947	for (int i=0; i < num; i++) {
1948	Vector3 origin = GetBBox(leaf).GetRandomPoint();
1949	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1950	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1951	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1952	rays.push_back(SimpleRay(origin,
1953	direction,
1954	SamplingStrategy::RSS_BASED_DISTRIBUTION,
1955	1.0f
1956	));
1957	}
1958	} else {
1959	RssTreeInterior in = (RssTreeInterior )node;
1960	// both nodes for directional splits
1961	tstack.push(in->front);
1962	tstack.push(in->back);
1963	}
1964	}
1965
1966	return (int)rays.size();
1967	}
1968
1969	void
1970	RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
1971	{
1972	stack<RssTreeNode *> tstack;
1973	PushRoots(tstack);
1974
1975	while (!tstack.empty()) {
1976	RssTreeNode *node = tstack.top();
1977	tstack.pop();
1978
1979	if (node->IsLeaf()) {
1980	RssTreeLeaf leaf = (RssTreeLeaf )node;
1981	leaves.push_back(leaf);
1982	} else {
1983	RssTreeInterior in = (RssTreeInterior )node;
1984	// both nodes for directional splits
1985	tstack.push(in->front);
1986	tstack.push(in->back);
1987	}
1988	}
1989	}
1990
1991	bool
1992	RssTree::ValidLeaf(RssTreeLeaf *leaf) const
1993	{
1994	return true;
1995	//return leaf->rays.size() > termMinRays/4;
1996	}
1997
1998
1999	void
2000	RssTree::PickEdgeRays(RssTreeLeaf *leaf,
2001	int &r1,
2002	int &r2
2003	)
2004	{
2005	int nrays = leaf->rays.size();
2006
2007	if (nrays == 2) {
2008	r1 = 0;
2009	r2 = 1;
2010	return;
2011	}
2012
2013	#if 1
2014	int tries = min(20, nrays);
2015
2016	while (--tries) {
2017	r1 = Random(nrays);
2018	r2 = Random(nrays);
2019	if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
2020	break;
2021	}
2022
2023	if (r1 == r2)
2024	r2 = (r1+1)%leaf->rays.size();
2025
2026	#else
2027	// pick a random ray
2028	int base = Random(nrays);
2029	RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
2030	Intersectable *baseObject = baseRay->GetObject();
2031
2032	// and a random 5D derivative which will be used to find the minimal projected distances
2033	Vector3 spatialDerivative;
2034	Vector3 directionalDerivative;
2035
2036	while (1) {
2037	spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2038	RandomValue(-1.0f, 1.0f),
2039	RandomValue(-1.0f, 1.0f));
2040	float m = Magnitude(spatialDerivative);
2041	if (m != 0) {
2042	spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
2043	break;
2044	}
2045	}
2046
2047	while (1) {
2048	directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
2049	RandomValue(-1.0f, 1.0f),
2050	0.0f);
2051	float m = Magnitude(directionalDerivative);
2052	if (m != 0) {
2053	directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
2054	break;
2055	}
2056	}
2057
2058	// now find the furthest sample from the same object and the closest from a different object
2059	int i = 0;
2060	float minDist = MAX_FLOAT;
2061	float maxDist = -MAX_FLOAT;
2062	r1 = base;
2063	r2 = base;
2064	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2065	ri != leaf->rays.end();
2066	ri++, i++) {
2067	float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
2068	*ri,
2069	spatialDerivative,
2070	directionalDerivative
2071	);
2072
2073	if ((*ri).GetObject() == baseObject) {
2074	if (dist > maxDist) {
2075	maxDist = dist;
2076	r1 = i;
2077	}
2078	} else {
2079	if (dist > 0.0f && dist < minDist) {
2080	minDist = dist;
2081	r2 = i;
2082	}
2083	}
2084	}
2085
2086	#endif
2087	}
2088
2089
2090	struct RayNeighbor {
2091	int rayInfo;
2092	float distance;
2093	RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
2094	};
2095
2096	void
2097	RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
2098	vector<RssTreeLeaf::SilhouetteRays> &rays
2099	)
2100	{
2101	// for every leaf find its neares neighbor from a different object
2102	vector<RayNeighbor> neighbors(leaf->rays.size());
2103	int i, j;
2104	for (i=0; i < leaf->rays.size(); i++)
2105	for (j=0; j < leaf->rays.size(); j++)
2106	if (i!=j) {
2107	float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
2108	if (d < neighbors[i].distance) {
2109	neighbors[i].distance = d;
2110	neighbors[i].rayInfo = j;
2111	}
2112	}
2113
2114	// now check which are the pairs of nearest neighbors
2115	for (i=0; i < leaf->rays.size(); i++) {
2116	int j = neighbors[i].rayInfo;
2117	if (neighbors[j].rayInfo == i) {
2118	// this is a silhouette edge pair
2119	if (i < j) {
2120	// generate an silhoutte ray pair
2121	rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
2122	&leaf->rays[j]));
2123	}
2124	} else {
2125	// this is not a silhouette ray - delete???
2126	}
2127	}
2128
2129	}
2130
2131
2132	bool
2133	GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
2134	const int index,
2135	int &i1,
2136	int &i2,
2137	int &i3)
2138	{
2139	int found = 0;
2140	int indices[3];
2141
2142	int size = rays.size();
2143	// use russian roulete selection for the tripple
2144	// number of free positions for the bullet
2145	int positions = size - 1;
2146	int i;
2147	for (i=0; i < size; i++) {
2148	if (rays[i]->mRay->Mailed())
2149	positions--;
2150	}
2151
2152	if (positions < 3)
2153	return false;
2154
2155	for (i=0; i < size; i++) {
2156	if (i != index && !rays[i]->mRay->Mailed()) {
2157	float p = (3 - found)/(float)positions;
2158	if (Random(1.0f) < p) {
2159	indices[found] = i;
2160	found++;
2161	}
2162	}
2163	positions--;
2164	}
2165	return true; // corr. matt
2166	}
2167
2168	bool
2169	RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
2170	const RssTreeNode::RayInfo &r1,
2171	const RssTreeNode::RayInfo &r2,
2172	const RssTreeNode::RayInfo &r3)
2173	{
2174
2175
2176	return false;
2177	}
2178
2179	int
2180	RssTree::RemoveInteriorRays(
2181	RssTreeLeaf *leaf
2182	)
2183	{
2184	#if 1
2185	// first collect all objects refered in this leaf
2186	map<Intersectable , vector<RssTreeLeaf::RayInfo > > rayMap;
2187
2188	RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
2189	for (; it != leaf->rays.end(); ++it) {
2190	Intersectable object = (it).GetObject();
2191
2192	rayMap[object].push_back(&(*it));
2193	// vector<RayInfo > data = rayMap.Find(object);
2194	// if (data) {
2195	// data->push_back(&(*it));
2196	// } else {
2197	// // rayMap[object] = vector<RayInfo *>;
2198	// rayMap[object].push_back(&(*it));
2199	// }
2200	}
2201
2202	// now go through all objects
2203	map<Intersectable , vector<RssTreeLeaf::RayInfo > >::iterator mi;
2204
2205	// infos of mailed rays are scheduled for removal
2206	VssRay::NewMail();
2207	for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
2208	vector<RssTreeLeaf::RayInfo > &rays = (mi).second;
2209	vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
2210	int rayIndex = 0;
2211
2212	for (; ri != rays.end(); ++ri, ++rayIndex) {
2213	RssTreeNode::RayInfo ray = ri;
2214	int tries = rays.size();
2215	for (int i = 0; i < tries; i++) {
2216	int r1, r2, r3;
2217	GetRandomTripple(rays,
2218	rayIndex,
2219	r1,
2220	r2,
2221	r3);
2222	if (IsRayConvexCombination(*ray,
2223	*rays[r1],
2224	*rays[r2],
2225	*rays[r3])) {
2226	ray->mRay->Mail();
2227	}
2228	}
2229	}
2230	}
2231
2232
2233	#endif
2234	return 0;
2235	}
2236
2237
2238	void
2239	RssTree::GenerateLeafRaysSimple(RssTreeLeaf *leaf,
2240	const int numberOfRays,
2241	SimpleRayContainer &rays)
2242	{
2243
2244	if (numberOfRays == 0)
2245	return;
2246
2247	int startIndex = (int)rays.size();
2248	// Debug<<"B"<<flush;
2249
2250	AxisAlignedBox3 box = GetBBox(leaf);
2251	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2252
2253	float boxSize = Magnitude(box.Diagonal());
2254
2255	const float smoothRange = 0.1f;
2256	if (smoothRange != 0.0f) {
2257	box.Scale(1.0f + smoothRange);
2258	dirBox.Scale(1.0f + smoothRange);
2259	}
2260
2261	// Debug<<"R"<<flush;
2262
2263	for (int i=0; i < numberOfRays; i++) {
2264	// Debug<<i<<flush;
2265	float r[5];
2266
2267	r[0] = RandomValue(0.0f, 1.0f);
2268	r[1] = RandomValue(0.0f, 1.0f);
2269	r[2] = RandomValue(0.0f, 1.0f);
2270	r[3] = RandomValue(0.0f, 1.0f);
2271	r[4] = RandomValue(0.0f, 1.0f);
2272
2273	GenerateLeafRay(r, rays, box, dirBox);
2274	}
2275
2276	// Debug<<"D"<<flush;
2277
2278	float density = 1.0f;
2279	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2280	for (int i=startIndex; i < rays.size(); i++) {
2281	rays[i].mPdf = p;
2282	}
2283	}
2284
2285	void
2286	RssTree::GenerateLeafRay(
2287	float *params,
2288	SimpleRayContainer &rays,
2289	const AxisAlignedBox3 &sbox,
2290	const AxisAlignedBox3 &sDirBBox
2291	)
2292	{
2293	Vector3 origin = sbox.GetPoint(Vector3(params[0], params[1], params[2]));
2294
2295	Vector3 dirVector = sDirBBox.GetPoint(Vector3(params[3], params[4], 0.0f));
2296
2297	Vector3 direction = VssRay::GetInvDirParam(dirVector.x, dirVector.y);
2298
2299	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2300	float dist = 0.0f;
2301	float minT, maxT;
2302
2303	float boxSize = Magnitude(sbox.Size());
2304	// compute interection of the ray with the box
2305	#if 1
2306	if (sbox.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2307	dist = (maxT + 0.0f*boxSize);
2308	#else
2309	dist = 0.5f*boxSize;
2310	#endif
2311
2312	origin += direction*dist;
2313
2314	//Debug<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2315	rays.push_back(SimpleRay(origin,
2316	direction,
2317	SamplingStrategy::RSS_BASED_DISTRIBUTION,
2318	1.0f
2319	));
2320	}
2321
2322	void
2323	RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
2324	const int numberOfRays,
2325	SimpleRayContainer &rays)
2326	{
2327
2328	if (numberOfRays == 0)
2329	return;
2330
2331
2332	int nrays = (int)leaf->rays.size();
2333	int startIndex = (int)rays.size();
2334
2335	AxisAlignedBox3 box = GetBBox(leaf);
2336	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
2337
2338
2339	const float smoothRange = 0.1f;
2340	if (smoothRange != 0.0f) {
2341	box.Scale(1.0f + smoothRange);
2342	dirBox.Scale(1.0f + smoothRange);
2343	}
2344
2345	Exporter *exporter = NULL;
2346	VssRayContainer selectedRays;
2347	static int counter = 0;
2348	if (counter < 0) {
2349	char filename[128];
2350	sprintf(filename, "selected-rays%04d.x3d", counter);
2351	exporter = Exporter::GetExporter(filename);
2352	// exporter->SetWireframe();
2353	// exporter->ExportKdTree(*mKdTree);
2354	exporter->SetWireframe();
2355	// exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
2356	exporter->SetWireframe();
2357	exporter->ExportBox(box);
2358	exporter->SetFilled();
2359	counter++;
2360	}
2361
2362
2363	float extendedConvexCombinationProb = 0.0f; //0.7f
2364	// float silhouetteCheckPercentage = 1.0f; //0.5f
2365	float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
2366	float silhouetteObjectCheckPercentage = 0.8f;
2367
2368	// int numberOfTries = numberOfRays*4;
2369	int numberOfTries = numberOfRays;
2370	int generated = 0;
2371
2372	vector<Intersectable *> pvs(0);
2373	if (silhouetteCheckPercentage != 0.0f) {
2374	Intersectable::NewMail();
2375	pvs.reserve(leaf->GetPvsSize());
2376	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
2377	ri != leaf->rays.end();
2378	ri++) {
2379	Intersectable object = (ri).GetObject();
2380	if (object) {
2381	if (!object->Mailed()) {
2382	object->Mail();
2383	object->mCounter = 1;
2384	pvs.push_back(object);
2385	if (exporter)
2386	exporter->ExportIntersectable(object);
2387	} else {
2388	object->mCounter++;
2389	}
2390	}
2391	}
2392	// sort objects based on their mCounter
2393	sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
2394	}
2395
2396	for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
2397	bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
2398	extendedConvexCombinationProb));
2399
2400	Vector3 origin, direction;
2401	// generate half of convex combination and half of random rays
2402	if (useExtendedConvexCombination) {
2403	// pickup 2 random rays
2404	int r1, r2;
2405	PickEdgeRays(leaf, r1, r2);
2406
2407
2408	Vector3 o1 = leaf->rays[r1].GetOrigin();
2409	Vector3 o2 = leaf->rays[r2].GetOrigin();
2410
2411	const float overlap = 0.0f;
2412	float w1, w2;
2413	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2414	origin = w1o1 + w2o2;
2415	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
2416	direction =
2417	w1*leaf->rays[r1].GetDir() +
2418	w2*leaf->rays[r2].GetDir();
2419
2420	// shift the origin a little bit
2421	direction.Normalize();
2422
2423	} else {
2424	Vector3 dirVector;
2425
2426	Vector3 pVector;
2427	Vector3 dVector;
2428
2429	pVector = Vector3(RandomValue(0.0f, 1.0f),
2430	RandomValue(0.0f, 1.0f),
2431	RandomValue(0.0f, 1.0f));
2432
2433	dVector = Vector3(RandomValue(0.0f, 1.0f),
2434	RandomValue(0.0f, 1.0f),
2435	0.0f);
2436
2437
2438	origin = box.GetPoint(pVector);
2439	dirVector = dirBox.GetPoint(dVector);
2440
2441	direction = VssRay::GetInvDirParam(dirVector.x,dirVector.y);
2442	}
2443
2444
2445	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
2446	float dist = 0.0f;
2447	float minT, maxT;
2448	// compute interection of the ray with the box
2449
2450	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
2451	dist = maxT;
2452
2453
2454	origin += direction*dist;
2455
2456	bool intersects = false;
2457
2458	if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
2459	if (exporter) {
2460	VssRay ray = new VssRay(origin, origin + 100direction, NULL, NULL);
2461	selectedRays.push_back(ray);
2462	}
2463
2464	// check whether the ray does not intersect already visible objects
2465	Ray traversalRay;
2466	traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
2467	for(vector<Intersectable *>::const_iterator oi = pvs.begin();
2468	oi != pvs.end();
2469	oi++) {
2470	Intersectable object = oi;
2471	// do not test every object
2472	if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
2473	object->CastRay(traversalRay) ) {
2474	intersects = true;
2475	break;
2476	}
2477	}
2478	}
2479
2480
2481	if (!intersects) {
2482	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
2483	rays.push_back(SimpleRay(origin,
2484	direction,
2485	SamplingStrategy::RSS_BASED_DISTRIBUTION,
2486	1.0f
2487	));
2488	generated++;
2489	}
2490	}
2491
2492	// cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
2493	// assign pdfs to the generated rays
2494	// float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
2495	float density = 1.0f;
2496
2497	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
2498	for (int i=startIndex; i < rays.size(); i++) {
2499	rays[i].mPdf = p;
2500	}
2501
2502
2503	if (exporter) {
2504	exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
2505	delete exporter;
2506	CLEAR_CONTAINER(selectedRays);
2507	}
2508
2509	}
2510
2511	int
2512	RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
2513	const float ratio)
2514	{
2515	int i;
2516	int j;
2517
2518	for (j=0, i=0; i < leaf->rays.size(); i++) {
2519	if (Random(1.0f) < ratio) {
2520	// copy a valid sample
2521	if (i!=j)
2522	leaf->rays[j] = leaf->rays[i];
2523	j++;
2524	} else {
2525	// delete the ray
2526	leaf->rays[i].mRay->Unref();
2527	if (leaf->rays[i].mRay->RefCount() != 0) {
2528	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
2529	exit(1);
2530	}
2531	delete leaf->rays[i].mRay;
2532	}
2533	}
2534
2535
2536	leaf->rays.resize(j);
2537	int removed = (i-j);
2538	stat.rayRefs -= removed;
2539	return removed;
2540	}
2541
2542	int
2543	RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
2544	const float ratio)
2545	{
2546	int i;
2547
2548	if (leaf->rays.size() == 0)
2549	return 0;
2550
2551	sort(leaf->rays.begin(),
2552	leaf->rays.end(),
2553	RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
2554
2555	int desired = int(ratio*leaf->rays.size());
2556	int removed = (int)leaf->rays.size() - desired;
2557
2558	for (i=desired; i < (int)leaf->rays.size(); i++) {
2559	// delete the ray
2560	leaf->rays[i].mRay->Unref();
2561	if (leaf->rays[i].mRay->RefCount() != 0) {
2562	cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
2563	cerr<<"desired="<<desired<<"i="<<i<<" size="<<(int)leaf->rays.size()<<endl;
2564	exit(1);
2565	}
2566	delete leaf->rays[i].mRay;
2567	}
2568
2569	leaf->rays.resize(desired);
2570	stat.rayRefs -= removed;
2571	return removed;
2572	}
2573
2574
2575	int
2576	RssTree::PruneRays(
2577	const int maxRays
2578	)
2579	{
2580
2581	stack<RssTreeNode *> tstack;
2582	int prunned = 0;
2583
2584	// prune more rays to amortize the procedure
2585	int desired = (int)(maxRays * 0.8f);
2586
2587	Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;
2588
2589	// prune random rays from each leaf so that the ratio's remain the same
2590	PushRoots(tstack);
2591
2592	if ( maxRays >= stat.rayRefs )
2593	return 0;
2594
2595	float ratio = desired/(float)stat.rayRefs;
2596
2597	while (!tstack.empty()) {
2598	RssTreeNode *node = tstack.top();
2599	tstack.pop();
2600
2601	if (node->IsLeaf()) {
2602	RssTreeLeaf leaf = (RssTreeLeaf )node;
2603	//prunned += PruneRaysRandom(leaf, ratio);
2604	prunned += PruneRaysContribution(leaf, ratio);
2605	} else {
2606	RssTreeInterior in = (RssTreeInterior )node;
2607	// both nodes for directional splits
2608	tstack.push(in->front);
2609	tstack.push(in->back);
2610	}
2611	}
2612
2613	Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
2614
2615	return prunned;
2616	}
2617
2618	int
2619	RssTree::GenerateRays(const int numberOfRays,
2620	const int numberOfLeaves,
2621	SimpleRayContainer &rays)
2622	{
2623
2624	if (1) {
2625	float param[5];
2626
2627	for (int i=0; i < numberOfRays; i++) {
2628	mHalton.GetNext(5, param);
2629	GenerateRay(param, rays);
2630	}
2631
2632	return rays.size();
2633	}
2634
2635	vector<RssTreeLeaf *> leaves;
2636
2637	Debug<<"Collecting leaves..."<<endl<<flush;
2638	CollectLeaves(leaves);
2639	Debug<<"done."<<endl<<flush;
2640
2641
2642	if (numberOfLeaves != leaves.size()) {
2643	Debug<<"sorting leaves..."<<endl<<flush;
2644	sort(leaves.begin(),
2645	leaves.end(),
2646	GreaterContribution);
2647	Debug<<"done."<<endl<<flush;
2648	}
2649
2650
2651	Debug<<"Evaluating leaf importance..."<<endl<<flush;
2652
2653	float sumContrib = 0.0;
2654	int i;
2655	int k = 0;
2656	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2657	if (ValidLeaf(leaves[i])) {
2658	float c = leaves[i]->GetImportance();
2659	sumContrib += c;
2660	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
2661	k++;
2662	}
2663
2664	Debug<<"done."<<endl<<flush;
2665
2666	float avgContrib = sumContrib/numberOfLeaves;
2667
2668	// always generate at leat n ray per leaf
2669	int fixedPerLeaf = 1;
2670	// int fixedPerLeaf = 0;
2671	int fixed = fixedPerLeaf*(int)leaves.size();
2672
2673	// check if the number of fixed is not too big
2674
2675	if (fixed >= numberOfRays/2) {
2676	fixedPerLeaf = 0;
2677	fixed = 0;
2678	}
2679
2680	int iGenerated = numberOfRays - fixed;
2681	float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
2682
2683	k = 0;
2684
2685	Debug<<"generating leaf rays..."<<endl<<flush;
2686
2687	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
2688	if (ValidLeaf(leaves[i])) {
2689	RssTreeLeaf *leaf = leaves[i];
2690	float c = leaf->GetImportance();
2691
2692	int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
2693	GenerateLeafRaysSimple(leaf, num, rays);
2694
2695	k++;
2696	}
2697
2698	Debug<<"done."<<endl<<flush;
2699
2700	return (int)rays.size();
2701	}
2702
2703
2704	float
2705	RssTree::GetAvgPvsSize()
2706	{
2707	stack<RssTreeNode *> tstack;
2708	PushRoots(tstack);
2709
2710	int sumPvs = 0;
2711	int leaves = 0;
2712	while (!tstack.empty()) {
2713	RssTreeNode *node = tstack.top();
2714	tstack.pop();
2715
2716	if (node->IsLeaf()) {
2717	RssTreeLeaf leaf = (RssTreeLeaf )node;
2718	// update pvs size
2719	UpdatePvsSize(leaf);
2720	sumPvs += leaf->GetPvsSize();
2721	leaves++;
2722	} else {
2723	RssTreeInterior in = (RssTreeInterior )node;
2724	// both nodes for directional splits
2725	tstack.push(in->front);
2726	tstack.push(in->back);
2727	}
2728	}
2729
2730
2731	return sumPvs/(float)leaves;
2732	}
2733
2734	float weightAbsContributions = ABS_CONTRIBUTION_WEIGHT;
2735	// if small very high importance of the last sample
2736	// if 1.0f then weighs = 1 1/2 1/3 1/4
2737	float passSampleWeightDecay = 1.0f;
2738	//float passSampleWeightDecay = 1000.0f;
2739
2740	float
2741	RssTree::GetSampleWeight(const int pass)
2742	{
2743	int passDiff = mCurrentPass - pass;
2744	float weight;
2745	if (0)
2746	weight = 1.0f/(passDiff + passSampleWeightDecay);
2747	else
2748	switch (passDiff) {
2749	case 0:
2750	// weight = 1.0f;
2751	// break;
2752	default:
2753	weight = 1.0f;
2754	break;
2755	// case 1:
2756	// weight = 0.5f;
2757	// break;
2758	// case 2:
2759	// weight = 0.25f;
2760	// break;
2761	// case 3:
2762	// weight = 0.12f;
2763	// break;
2764	// case 4:
2765	// weight = 0.06f;
2766	// break;
2767	// default:
2768	// weight = 0.03f;
2769	// break;
2770	}
2771	return weight;
2772	}
2773
2774	void
2775	RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
2776	float &weight,
2777	float &contribution)
2778	{
2779	VssRay *ray = info.mRay;
2780	weight = GetSampleWeight(mCurrentPass);
2781	contribution =
2782	weightAbsContributions*ray->mPvsContribution +
2783	(1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
2784	// store the computed value
2785	info.mRay->mWeightedPvsContribution = weight*contribution;
2786	}
2787
2788
2789	void
2790	RssTree::ComputeImportance(RssTreeLeaf *leaf)
2791	{
2792	#if 1
2793	if (leaf->mTotalRays)
2794	leaf->mImportance = leaf->mTotalContribution / leaf->mTotalRays;
2795	else
2796	leaf->mImportance = Limits::Small;
2797	return;
2798	#endif
2799
2800	if (0)
2801	leaf->mImportance = leaf->GetAvgRayContribution();
2802	else {
2803	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
2804	it_end = leaf->rays.end();
2805
2806	float sumContributions = 0.0f;
2807	float sumRelContributions = 0.0f;
2808	float sumWeights = 0.0f;
2809	float maxContribution = 0.0f;
2810	float maxRelContribution = 0.0f;
2811	float sumLength = 0;
2812
2813	for (; it != it_end; ++it) {
2814	VssRay ray = (it).mRay;
2815	float sumLength = ray->Length();
2816
2817	float weight = GetSampleWeight(ray->mPass);
2818
2819	float c1 = weight*ray->mPvsContribution;
2820	float c2 = weight*ray->mRelativePvsContribution;
2821	sumContributions += c1;
2822	if (c1 > maxContribution)
2823	maxContribution = c1;
2824
2825	//$$ 20.7. changed to sqr to pronouce higher contribution so that
2826	// they do not get smoothed!!
2827	//sumRelContributions += weight*ray->mRelativePvsContribution;
2828
2829	sumRelContributions += c2;
2830	if (c2 > maxRelContribution)
2831	maxRelContribution = c2;
2832
2833	//sumWeights += weight;
2834	sumWeights += 1.0f;
2835	}
2836
2837	float distImportance = 0.0f;
2838
2839	float spatialArea = GetBBox(leaf).SurfaceArea();
2840	float dirArea = GetDirBBox(leaf).SurfaceArea();
2841	if (leaf->rays.size()) {
2842	// compute average length of the ray
2843	float avgLength = sumLength/leaf->rays.size();
2844	// compute estimated area of the visible surface corresponding to these rays
2845	float ss = spatialArea/2.0f;
2846	float sd = dirArea;
2847	float s = ss + avgLength(2sqrt(sssd) + avgLengthsd);
2848
2849	distImportance = s/leaf->rays.size();
2850	}
2851
2852	// $$
2853	// sumWeights = leaf->mTotalRays;
2854
2855	//$$ test 30.11. 2006
2856	// normalize per volume of the cell not the number of rays!
2857	// sumWeights = spatialArea*dirArea;
2858
2859	if (sumWeights != 0.0f) {
2860	leaf->mImportance =
2861	sqr(((weightAbsContributions*sumContributions +
2862	(1.0f - weightAbsContributions)*sumRelContributions)/sumWeights));
2863
2864	// leaf->mImportance =
2865	// (weightAbsContributions*maxContribution +
2866	// (1.0f - weightAbsContributions)*maxRelContribution)/sumWeights;
2867
2868	} else
2869	leaf->mImportance = 0.0f;
2870
2871	}
2872
2873	// return GetAvgRayContribution()*mImportance;
2874	//return GetAvgRayContribution();
2875	}
2876
2877
2878	float
2879	RssTreeNode::GetImportance() const
2880	{
2881	// return sqr(mImportance);
2882	return mImportance;
2883	}
2884
2885
2886	float
2887	RssTreeLeaf::ComputePvsEntropy() const
2888	{
2889	int samples = 0;
2890	Intersectable::NewMail();
2891	// set all object as belonging to the fron pvs
2892	for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
2893	ri != rays.end();
2894	ri++)
2895	if ((*ri).mRay->IsActive()) {
2896	Intersectable object = (ri).GetObject();
2897	if (object) {
2898	if (!object->Mailed()) {
2899	object->Mail();
2900	object->mCounter = 1;
2901	} else
2902	object->mCounter++;
2903	samples++;
2904	}
2905	}
2906
2907	float entropy = 0.0f;
2908
2909	if (samples > 1) {
2910	Intersectable::NewMail();
2911	for(RayInfoContainer::const_iterator ri = rays.begin();
2912	ri != rays.end();
2913	ri++)
2914	if ((*ri).mRay->IsActive()) {
2915	Intersectable object = (ri).GetObject();
2916	if (object) {
2917	if (!object->Mailed()) {
2918	object->Mail();
2919	float p = object->mCounter/(float)samples;
2920	entropy -= p*log(p);
2921	}
2922	}
2923	}
2924	entropy = entropy/log((float)samples);
2925	}
2926	else
2927	entropy = 1.0f;
2928
2929	return entropy;
2930	}
2931
2932	float
2933	RssTreeLeaf::ComputeRayLengthEntropy() const
2934	{
2935	// get sum of all ray lengths
2936	// consider only passing rays or originating rays
2937	float sum = 0.0f;
2938	int samples = 0;
2939	int i=0;
2940	for(RayInfoContainer::const_iterator ri = rays.begin();
2941	ri != rays.end();
2942	ri++)
2943	if ((*ri).mRay->IsActive()) {
2944	// float s;
2945	// if (i == 0)
2946	// s = 200;
2947	// else
2948	// s = 100;
2949	// i++;
2950
2951	sum += (*ri).mRay->GetSize();
2952	samples++;
2953	}
2954
2955
2956	float entropy = 0.0f;
2957
2958	if (samples > 1) {
2959	i = 0;
2960	for(RayInfoContainer::const_iterator ri = rays.begin();
2961	ri != rays.end();
2962	ri++)
2963	if ((*ri).mRay->IsActive()) {
2964	// float s;
2965	// if (i==0)
2966	// s = 200;
2967	// else
2968	// s = 100;
2969	// i++;
2970	// float p = s/sum;
2971	float p = (*ri).mRay->GetSize()/sum;
2972	entropy -= p*log(p);
2973	}
2974	entropy = entropy/log((float)samples);
2975	} else
2976	entropy = 1.0f;
2977
2978	return entropy;
2979	}
2980
2981
2982
2983	float
2984	RssTreeLeaf::ComputeEntropyImportance() const
2985	{
2986
2987	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2988	return 1.0f - ComputePvsEntropy();
2989	}
2990
2991	float
2992	RssTreeLeaf::ComputeRayContributionImportance() const
2993	{
2994	// mPvsEntropy = ComputePvsEntropy();
2995	// mRayLengthEntropy = ComputeRayLengthEntropy();
2996
2997	// mEntropy = 1.0f - ComputeRayLengthEntropy();
2998	return 1.0f - ComputePvsEntropy();
2999	}
3000
3001
3002	AxisAlignedBox3
3003	RssTree::GetShrankedBBox(const RssTreeNode *node) const
3004	{
3005	if (node->parent == NULL)
3006	return bbox;
3007
3008	if (!node->IsLeaf())
3009	return ((RssTreeInterior *)node)->bbox;
3010
3011	// evaluate bounding box from the ray origins
3012	AxisAlignedBox3 box;
3013	box.Initialize();
3014	RssTreeLeaf leaf = (RssTreeLeaf )node;
3015	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
3016	ri != leaf->rays.end();
3017	ri++)
3018	if ((*ri).mRay->IsActive()) {
3019	box.Include((*ri).GetOrigin());
3020	}
3021	return box;
3022	}
3023
3024
3025	int
3026	RssTree::CollectRays(VssRayContainer &rays,
3027	const int number)
3028	{
3029	VssRayContainer allRays;
3030	CollectRays(allRays);
3031
3032
3033	int desired = min(number, (int)allRays.size());
3034	float prob = desired/(float)allRays.size();
3035	//int skip = allRays.size()/desired;
3036	while (rays.size() < desired) {
3037	VssRayContainer::const_iterator it = allRays.begin();
3038	for (; it != allRays.end() && rays.size() < desired; it++) {
3039	if (Random(1.0f) < prob)
3040	rays.push_back(*it);
3041	}
3042	}
3043	return (int)rays.size();
3044	}
3045
3046
3047	int
3048	RssTree::CollectRays(VssRayContainer &rays
3049	)
3050	{
3051	VssRay::NewMail();
3052
3053	stack<RssTreeNode *> tstack;
3054	PushRoots(tstack);
3055
3056	while (!tstack.empty()) {
3057	RssTreeNode *node = tstack.top();
3058	tstack.pop();
3059
3060	if (node->IsLeaf()) {
3061	RssTreeLeaf leaf = (RssTreeLeaf )node;
3062	// update pvs size
3063	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
3064	for (;it != leaf->rays.end(); ++it)
3065	if (!(*it).mRay->Mailed()) {
3066	(*it).mRay->Mail();
3067	rays.push_back((*it).mRay);
3068	}
3069	} else {
3070	RssTreeInterior in = (RssTreeInterior )node;
3071	// both nodes for directional splits
3072	tstack.push(in->front);
3073	tstack.push(in->back);
3074	}
3075	}
3076
3077	return (int)rays.size();
3078	}
3079
3080
3081	RssTreeNode *
3082	RssTree::GetRoot(Intersectable *object) const
3083	{
3084	if (mPerObjectTree && object) {
3085	int id = object->GetId()-1;
3086	if (id < 0 \|\| id >= mRoots.size()) {
3087	Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<(int)mRoots.size()<<
3088	endl<<flush;
3089	id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
3090	}
3091	return mRoots[id];
3092	} else
3093	return mRoots[0];
3094	}
3095
3096	void
3097	RssTree::PushRoots(priority_queue<TraversalData> &st) const
3098	{
3099	for (int i=0; i < mRoots.size(); i++)
3100	st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
3101	}
3102
3103	void
3104	RssTree::PushRoots(stack<RssTreeNode *> &st) const
3105	{
3106	for (int i=0; i < mRoots.size(); i++)
3107	st.push(mRoots[i]);
3108	}
3109
3110	void
3111	RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info) const
3112	{
3113	for (int i=0; i < mRoots.size(); i++)
3114	st.push(RayTraversalData(mRoots[i], info));
3115	}
3116
3117
3118
3119
3120
3121	}

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

Download in other formats: