Context Navigation

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

Revision 1328, 69.8 KB checked in by mattausch, 18 years ago (diff)

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

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

Download in other formats: