Context Navigation

source: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp @ 409

Revision 409, 35.7 KB checked in by mattausch, 19 years ago (diff)
worked on kd view space partitioning structure worked on render simulation

Line
1
2	// ================================================================
3	// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
4	// ****************************************************************
5	/**
6	The KD tree based LSDS
7	*/
8	// Initial coding by
9	/**
10	@author Jiri Bittner
11	*/
12
13	// Standard headers
14	#include <stack>
15	#include <queue>
16	#include <algorithm>
17	#include <fstream>
18	#include <string>
19
20	#include "VspKdTree.h"
21
22	#include "Environment.h"
23	#include "VssRay.h"
24	#include "Intersectable.h"
25	#include "Ray.h"
26
27
28	#define DEBUG_DIR_SPLIT 0
29
30
31
32	// Static variables
33	int
34	VspKdTreeLeaf::mailID = 0;
35
36	inline void AddObject2Pvs(Intersectable *object,
37	const int side,
38	int &pvsBack,
39	int &pvsFront)
40	{
41	if (!object)
42	return;
43
44	if (side <= 0)
45	{
46	if (!object->Mailed() && !object->Mailed(2))
47	{
48	++ pvsBack;
49
50	if (object->Mailed(1))
51	object->Mail(2);
52	else
53	object->Mail();
54	}
55	}
56
57	if (side >= 0)
58	{
59	if (!object->Mailed(1) && !object->Mailed(2))
60	{
61	++ pvsFront;
62	if (object->Mailed())
63	object->Mail(2);
64	else
65	object->Mail(1);
66	}
67	}
68	}
69
70
71	// Constructor
72	VspKdTree::VspKdTree()
73	{
74	environment->GetIntValue("VspKdTree.maxDepth", termMaxDepth);
75	environment->GetIntValue("VspKdTree.minPvs", termMinPvs);
76	environment->GetIntValue("VspKdTree.minRays", termMinRays);
77	environment->GetFloatValue("VspKdTree.maxRayContribution", termMaxRayContribution);
78	environment->GetFloatValue("VspKdTree.maxCostRatio", termMaxCostRatio);
79
80	environment->GetFloatValue("VspKdTree.minSize", termMinSize);
81	termMinSize = sqr(termMinSize);
82
83	environment->GetFloatValue("VspKdTree.refDirBoxMaxSize", refDirBoxMaxSize);
84	refDirBoxMaxSize = sqr(refDirBoxMaxSize);
85
86	environment->GetFloatValue("VspKdTree.epsilon", epsilon);
87	environment->GetFloatValue("VspKdTree.ct_div_ci", ct_div_ci);
88
89	environment->GetFloatValue("VspKdTree.maxTotalMemory", maxTotalMemory);
90	environment->GetFloatValue("VspKdTree.maxStaticMemory", maxStaticMemory);
91
92
93	float refDirAngle;
94	environment->GetFloatValue("VspKdTree.refDirAngle", refDirAngle);
95
96	environment->GetIntValue("VspKdTree.accessTimeThreshold", accessTimeThreshold);
97	//= 1000;
98	environment->GetIntValue("VspKdTree.minCollapseDepth", minCollapseDepth);
99
100
101	// split type
102	char sname[128];
103	environment->GetStringValue("VspKdTree.splitType", sname);
104	string name(sname);
105
106	if (name.compare("regular") == 0)
107	splitType = ESplitRegular;
108	else
109	{
110	if (name.compare("heuristic") == 0)
111	splitType = ESplitHeuristic;
112	else
113	{
114	cerr << "Invalid VspKdTree split type " << name << endl;
115	exit(1);
116	}
117	}
118
119	environment->GetBoolValue("VspKdTree.randomize", randomize);
120
121	root = NULL;
122
123	splitCandidates = new vector<SortableEntry>;
124	}
125
126
127	VspKdTree::~VspKdTree()
128	{
129	if (root)
130	delete root;
131	}
132
133
134
135
136	void
137	VspKdStatistics::Print(ostream &app) const
138	{
139	app << "===== VspKdTree statistics ===============\n";
140
141	app << "#N_RAYS ( Number of rays )\n"
142	<< rays <<endl;
143
144	app << "#N_INITPVS ( Initial PVS size )\n"
145	<< initialPvsSize <<endl;
146
147	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
148
149	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
150
151	app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
152	for (int i=0; i<7; i++)
153	app << splits[i] <<" ";
154	app <<endl;
155
156	app << "#N_RAYREFS ( Number of rayRefs )\n" <<
157	rayRefs << "\n";
158
159	app << "#N_RAYRAYREFS ( Number of rayRefs / ray )\n" <<
160	rayRefs/(double)rays << "\n";
161
162	app << "#N_LEAFRAYREFS ( Number of rayRefs / leaf )\n" <<
163	rayRefs/(double)Leaves() << "\n";
164
165	app << "#N_MAXRAYREFS ( Max number of rayRefs / leaf )\n" <<
166	maxRayRefs << "\n";
167
168
169	// app << setprecision(4);
170
171	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
172	maxDepthNodes*100/(double)Leaves()<<endl;
173
174	app << "#N_PMINPVSLEAVES ( Percentage of leaves with minCost )\n"<<
175	minPvsNodes*100/(double)Leaves()<<endl;
176
177	app << "#N_PMINRAYSLEAVES ( Percentage of leaves with minCost )\n"<<
178	minRaysNodes*100/(double)Leaves()<<endl;
179
180	app << "#N_PMINSIZELEAVES ( Percentage of leaves with minSize )\n"<<
181	minSizeNodes*100/(double)Leaves()<<endl;
182
183	app << "#N_PMAXRAYCONTRIBLEAVES ( Percentage of leaves with maximal ray contribution )\n"<<
184	maxRayContribNodes*100/(double)Leaves()<<endl;
185
186	app << "#N_ADDED_RAYREFS (Number of dynamically added ray references )\n"<<
187	addedRayRefs<<endl;
188
189	app << "#N_REMOVED_RAYREFS (Number of dynamically removed ray references )\n"<<
190	removedRayRefs<<endl;
191
192	// app << setprecision(4);
193
194	app << "#N_CTIME ( Construction time [s] )\n"
195	<< Time() << " \n";
196
197	app << "===== END OF VspKdTree statistics ==========\n";
198
199	}
200
201
202	void
203	VspKdTreeLeaf::UpdatePvsSize()
204	{
205	if (!mValidPvs)
206	{
207	Intersectable::NewMail();
208	int pvsSize = 0;
209	for(VspKdTreeNode::RayInfoContainer::iterator ri = rays.begin();
210	ri != rays.end(); ++ ri)
211	{
212	if ((*ri).mRay->IsActive())
213	{
214	Intersectable *object;
215	#if BIDIRECTIONAL_RAY
216	object = (*ri).mRay->mOriginObject;
217
218	if (object && !object->Mailed())
219	{
220	++ pvsSize;
221	object->Mail();
222	}
223	#endif
224	object = (*ri).mRay->mTerminationObject;
225	if (object && !object->Mailed())
226	{
227	++ pvsSize;
228	object->Mail();
229	}
230	}
231	}
232
233	mPvsSize = pvsSize;
234	mValidPvs = true;
235	}
236	}
237
238
239	void
240	VspKdTree::Construct(VssRayContainer &rays,
241	AxisAlignedBox3 *forcedBoundingBox)
242	{
243	stat.Start();
244
245	maxMemory = maxStaticMemory;
246
247	if (root)
248	delete root;
249
250	root = new VspKdTreeLeaf(NULL, rays.size());
251
252	// first construct a leaf that will get subdivide
253	VspKdTreeLeaf leaf = (VspKdTreeLeaf ) root;
254
255
256	stat.nodes = 1;
257
258	bbox.Initialize();
259	dirBBox.Initialize();
260
261
262	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
263	{
264	leaf->AddRay(VspKdTreeNode::RayInfo(*ri));
265
266	bbox.Include((*ri)->GetOrigin());
267	bbox.Include((*ri)->GetTermination());
268
269	dirBBox.Include(Vector3((ri)->GetDirParametrization(0), (ri)->GetDirParametrization(1), 0));
270	}
271
272	if (forcedBoundingBox)
273	bbox = *forcedBoundingBox;
274
275	cout<<"Bbox = "<<bbox<<endl;
276	cout<<"Dirr Bbox = "<<dirBBox<<endl;
277
278	stat.rays = leaf->rays.size();
279	leaf->UpdatePvsSize();
280
281	stat.initialPvsSize = leaf->GetPvsSize();
282
283	// Subdivide();
284	root = Subdivide(TraversalData(leaf, bbox, 0));
285
286	if (splitCandidates)
287	{
288	// force realease of this vector
289	delete splitCandidates;
290	splitCandidates = new vector<SortableEntry>;
291	}
292
293	stat.Stop();
294
295	stat.Print(cout);
296	cout<<"#Total memory=" << GetMemUsage() << endl;
297	}
298
299
300
301	VspKdTreeNode *VspKdTree::Subdivide(const TraversalData &tdata)
302	{
303	VspKdTreeNode *result = NULL;
304
305	priority_queue<TraversalData> tStack;
306	// stack<TraversalData> tStack;
307
308	tStack.push(tdata);
309
310	AxisAlignedBox3 backBox;
311	AxisAlignedBox3 frontBox;
312
313	int lastMem = 0;
314	while (!tStack.empty())
315	{
316	float mem = GetMemUsage();
317
318	if ( lastMem/10 != ((int)mem)/10)
319	{
320	cout << mem << " MB" << endl;
321	}
322	lastMem = (int)mem;
323
324	if ( mem > maxMemory )
325	{
326	// count statistics on unprocessed leafs
327	while (!tStack.empty())
328	{
329	EvaluateLeafStats(tStack.top());
330	tStack.pop();
331	}
332	break;
333	}
334
335	TraversalData data = tStack.top();
336	tStack.pop();
337
338
339	VspKdTreeNode node = SubdivideNode((VspKdTreeLeaf ) data.node,
340	data.bbox, backBox, frontBox);
341	if (result == NULL)
342	result = node;
343
344	if (!node->IsLeaf())
345	{
346	VspKdTreeInterior interior = (VspKdTreeInterior ) node;
347
348	// push the children on the stack
349	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
350	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
351	}
352	else
353	{
354	EvaluateLeafStats(data);
355	}
356	}
357
358	return result;
359	}
360
361
362	// returns selected plane for subdivision
363	int
364	VspKdTree::SelectPlane(VspKdTreeLeaf *leaf,
365	const AxisAlignedBox3 &box,
366	float &position,
367	int &raysBack,
368	int &raysFront,
369	int &pvsBack,
370	int &pvsFront)
371	{
372
373	int minDirDepth = 6;
374	int axis;
375	float costRatio;
376
377	if (splitType == ESplitRegular) {
378	costRatio = BestCostRatioRegular(leaf,
379	axis,
380	position,
381	raysBack,
382	raysFront,
383	pvsBack,
384	pvsFront
385	);
386
387	} else {
388	if (splitType == ESplitHeuristic)
389	costRatio = BestCostRatioHeuristic(leaf,
390	axis,
391	position,
392	raysBack,
393	raysFront,
394	pvsBack,
395	pvsFront
396	);
397	else {
398	cerr<<"VspKdTree: Unknown split heuristics\n";
399	exit(1);
400	}
401	}
402
403	if (costRatio > termMaxCostRatio) {
404	// cout<<"Too big cost ratio "<<costRatio<<endl;
405	return -1;
406	}
407
408	#if 0
409	cout<<
410	"pvs="<<leaf->mPvsSize<<
411	" rays="<<leaf->rays.size()<<
412	" rc="<<leaf->GetAvgRayContribution()<<
413	" axis="<<axis<<endl;
414	#endif
415
416	return axis;
417	}
418
419
420
421
422	float
423	VspKdTree::EvalCostRatio(
424	VspKdTreeLeaf *leaf,
425	const int axis,
426	const float position,
427	int &raysBack,
428	int &raysFront,
429	int &pvsBack,
430	int &pvsFront
431	)
432	{
433	raysBack = 0;
434	raysFront = 0;
435	pvsFront = 0;
436	pvsBack = 0;
437
438	float newCost;
439
440	Intersectable::NewMail(3);
441
442	// eval pvs size
443	int pvsSize = leaf->GetPvsSize();
444
445	Intersectable::NewMail(3);
446
447	if (axis <= VspKdTreeNode::SPLIT_Z) {
448	// this is the main ray classification loop!
449	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
450	ri != leaf->rays.end();
451	ri++)
452	if ((*ri).mRay->IsActive()) {
453
454	// determine the side of this ray with respect to the plane
455	int side = (ri).ComputeRayIntersection(axis, position, (ri).mRay->mT);
456	// (*ri).mRay->mSide = side;
457
458	if (side <= 0)
459	raysBack++;
460
461	if (side >= 0)
462	raysFront++;
463
464	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
465	}
466
467	AxisAlignedBox3 box = GetBBox(leaf);
468
469	float minBox = box.Min(axis);
470	float maxBox = box.Max(axis);
471	float sizeBox = maxBox - minBox;
472
473	// float sum = raysBack(position - minBox) + raysFront(maxBox - position);
474	float sum = pvsBack(position - minBox) + pvsFront(maxBox - position);
475
476	newCost = ct_div_ci + sum/sizeBox;
477
478	} else {
479
480	// directional split
481	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
482	ri != leaf->rays.end();
483	ri++)
484	if ((*ri).mRay->IsActive()) {
485
486	// determine the side of this ray with respect to the plane
487	int side;
488	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
489	side = 1;
490	else
491	side = -1;
492
493	if (side <= 0)
494	raysBack++;
495
496	if (side >= 0)
497	raysFront++;
498
499	// (*ri).mRay->mSide = side;
500	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
501
502	}
503
504	AxisAlignedBox3 box = GetDirBBox(leaf);
505
506	float minBox = box.Min(axis-3);
507	float maxBox = box.Max(axis-3);
508	float sizeBox = maxBox - minBox;
509
510	// float sum = raysBack(position - minBox) + raysFront(maxBox - position);
511	float sum =
512	pvsBack*(position - minBox) +
513	pvsFront*(maxBox - position);
514	// float sum = pvsBack + pvsFront;
515	newCost = ct_div_ci + sum/sizeBox;
516	}
517
518	// cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
519	// float oldCost = leaf->rays.size();
520	float oldCost = pvsSize;
521	float ratio = newCost/oldCost;
522	// cout<<"ratio="<<ratio<<endl;
523	return ratio;
524	}
525
526	float
527	VspKdTree::BestCostRatioRegular(
528	VspKdTreeLeaf *leaf,
529	int &axis,
530	float &position,
531	int &raysBack,
532	int &raysFront,
533	int &pvsBack,
534	int &pvsFront
535	)
536	{
537	int nRaysBack[6], nRaysFront[6];
538	int nPvsBack[6], nPvsFront[6];
539	float nPosition[6];
540	float nCostRatio[6];
541	int bestAxis = -1;
542
543	AxisAlignedBox3 sBox = GetBBox(leaf);
544	AxisAlignedBox3 dBox = GetDirBBox(leaf);
545	// int sAxis = box.Size().DrivingAxis();
546	int sAxis = sBox.Size().DrivingAxis();
547	int dAxis = dBox.Size().DrivingAxis() + 3;
548
549	bool onlyDrivingAxis = true;
550
551	for (axis = 0; axis < 6; axis++) {
552	if (!onlyDrivingAxis \|\| axis == sAxis \|\| axis == dAxis) {
553	if (axis < 3)
554	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
555	else
556	nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
557
558	nCostRatio[axis] = EvalCostRatio(leaf,
559	axis,
560	nPosition[axis],
561	nRaysBack[axis],
562	nRaysFront[axis],
563	nPvsBack[axis],
564	nPvsFront[axis]
565	);
566
567	if ( bestAxis == -1)
568	bestAxis = axis;
569	else
570	if ( nCostRatio[axis] < nCostRatio[bestAxis] )
571	bestAxis = axis;
572	}
573	}
574
575	axis = bestAxis;
576	position = nPosition[bestAxis];
577
578	raysBack = nRaysBack[bestAxis];
579	raysFront = nRaysFront[bestAxis];
580
581	pvsBack = nPvsBack[bestAxis];
582	pvsFront = nPvsFront[bestAxis];
583
584	return nCostRatio[bestAxis];
585	}
586
587	float
588	VspKdTree::BestCostRatioHeuristic(
589	VspKdTreeLeaf *leaf,
590	int &axis,
591	float &position,
592	int &raysBack,
593	int &raysFront,
594	int &pvsBack,
595	int &pvsFront
596	)
597	{
598	AxisAlignedBox3 box = GetBBox(leaf);
599	// AxisAlignedBox3 dirBox = GetDirBBox(node);
600
601	axis = box.Size().DrivingAxis();
602
603	SortSplitCandidates(leaf, axis);
604
605	// go through the lists, count the number of objects left and right
606	// and evaluate the following cost funcion:
607	// C = ct_div_ci + (qlrl + qrrr)/queries
608
609	int rl=0, rr = leaf->rays.size();
610	int pl=0, pr = leaf->GetPvsSize();
611	float minBox = box.Min(axis);
612	float maxBox = box.Max(axis);
613	float sizeBox = maxBox - minBox;
614
615	float minBand = minBox + 0.1*(maxBox - minBox);
616	float maxBand = minBox + 0.9*(maxBox - minBox);
617
618	float sum = rr*sizeBox;
619	float minSum = 1e20;
620
621	Intersectable::NewMail();
622	// set all object as belonging to the fron pvs
623	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
624	ri != leaf->rays.end();
625	ri++)
626	if ((*ri).mRay->IsActive()) {
627	Intersectable object = (ri).mRay->mTerminationObject;
628	if (object)
629	if (!object->Mailed()) {
630	object->Mail();
631	object->mCounter = 1;
632	} else
633	object->mCounter++;
634	}
635
636	Intersectable::NewMail();
637
638	for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
639	ci < splitCandidates->end();
640	ci++) {
641	VssRay *ray;
642	switch ((*ci).type) {
643	case SortableEntry::ERayMin: {
644	rl++;
645	ray = (VssRay ) (ci).data;
646	Intersectable *object = ray->mTerminationObject;
647	if (object && !object->Mailed()) {
648	object->Mail();
649	pl++;
650	}
651	break;
652	}
653	case SortableEntry::ERayMax: {
654	rr--;
655	ray = (VssRay ) (ci).data;
656	Intersectable *object = ray->mTerminationObject;
657	if (object) {
658	if (--object->mCounter == 0)
659	pr--;
660	}
661	break;
662	}
663	}
664	if ((ci).value > minBand && (ci).value < maxBand) {
665
666	sum = pl((ci).value - minBox) + pr(maxBox - (ci).value);
667
668	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
669	// cout<<"cost= "<<sum<<endl;
670
671	if (sum < minSum) {
672	minSum = sum;
673	position = (*ci).value;
674
675	raysBack = rl;
676	raysFront = rr;
677
678	pvsBack = pl;
679	pvsFront = pr;
680
681	}
682	}
683	}
684
685	float oldCost = leaf->GetPvsSize();
686	float newCost = ct_div_ci + minSum/sizeBox;
687	float ratio = newCost/oldCost;
688
689	// cout<<"===================="<<endl;
690	// cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
691	// <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
692	return ratio;
693	}
694
695	void
696	VspKdTree::SortSplitCandidates(
697	VspKdTreeLeaf *node,
698	const int axis
699	)
700	{
701
702	splitCandidates->clear();
703
704	int requestedSize = 2*(node->rays.size());
705	// creates a sorted split candidates array
706	if (splitCandidates->capacity() > 500000 &&
707	requestedSize < (int)(splitCandidates->capacity()/10) ) {
708
709	delete splitCandidates;
710	splitCandidates = new vector<SortableEntry>;
711	}
712
713	splitCandidates->reserve(requestedSize);
714
715	// insert all queries
716	for(VspKdTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
717	ri < node->rays.end();
718	ri++) {
719	bool positive = (*ri).mRay->HasPosDir(axis);
720	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
721	SortableEntry::ERayMax,
722	(*ri).ExtrapOrigin(axis),
723	(void )&ri)
724	);
725
726	splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
727	SortableEntry::ERayMin,
728	(*ri).ExtrapTermination(axis),
729	(void )&ri)
730	);
731	}
732
733	stable_sort(splitCandidates->begin(), splitCandidates->end());
734	}
735
736
737	void
738	VspKdTree::EvaluateLeafStats(const TraversalData &data)
739	{
740
741	// the node became a leaf -> evaluate stats for leafs
742	VspKdTreeLeaf leaf = (VspKdTreeLeaf )data.node;
743
744	if (data.depth >= termMaxDepth)
745	stat.maxDepthNodes++;
746
747	// if ( (int)(leaf->rays.size()) < termMinCost)
748	// stat.minCostNodes++;
749	if ( leaf->GetPvsSize() < termMinPvs)
750	stat.minPvsNodes++;
751
752	if ( leaf->GetPvsSize() < termMinRays)
753	stat.minRaysNodes++;
754
755	if (0 && leaf->GetAvgRayContribution() > termMaxRayContribution )
756	stat.maxRayContribNodes++;
757
758	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
759	stat.minSizeNodes++;
760	}
761
762	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
763	stat.maxRayRefs = leaf->rays.size();
764
765	}
766
767
768
769	VspKdTreeNode *
770	VspKdTree::SubdivideNode(
771	VspKdTreeLeaf *leaf,
772	const AxisAlignedBox3 &box,
773	AxisAlignedBox3 &backBBox,
774	AxisAlignedBox3 &frontBBox
775	)
776	{
777
778	if ( (leaf->GetPvsSize() < termMinPvs) \|\| (leaf->rays.size() < termMinRays) \|\|
779	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) \|\|
780	(leaf->depth >= termMaxDepth) \|\| SqrMagnitude(box.Size()) <= termMinSize )
781	{
782
783	#if 0
784	if (leaf->depth >= termMaxDepth) {
785	cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
786	cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
787	}
788	#endif
789
790	return leaf;
791	}
792
793	float position;
794
795	// first count ray sides
796	int raysBack;
797	int raysFront;
798	int pvsBack;
799	int pvsFront;
800
801	// select subdivision axis
802	int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
803
804	// cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
805	// pvsFront<<endl;
806
807	if (axis == -1) {
808	return leaf;
809	}
810
811	stat.nodes+=2;
812	stat.splits[axis]++;
813
814	// add the new nodes to the tree
815	VspKdTreeInterior *node = new VspKdTreeInterior(leaf->parent);
816
817	node->axis = axis;
818	node->position = position;
819	node->bbox = box;
820	node->dirBBox = GetDirBBox(leaf);
821
822	backBBox = box;
823	frontBBox = box;
824
825	VspKdTreeLeaf *back = new VspKdTreeLeaf(node, raysBack);
826	back->SetPvsSize(pvsBack);
827	VspKdTreeLeaf *front = new VspKdTreeLeaf(node, raysFront);
828	front->SetPvsSize(pvsFront);
829
830	// replace a link from node's parent
831	if ( leaf->parent )
832	leaf->parent->ReplaceChildLink(leaf, node);
833	// and setup child links
834	node->SetupChildLinks(back, front);
835
836	if (axis <= VspKdTreeNode::SPLIT_Z) {
837	backBBox.SetMax(axis, position);
838	frontBBox.SetMin(axis, position);
839
840	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
841	ri != leaf->rays.end();
842	ri++) {
843	if ((*ri).mRay->IsActive()) {
844
845	// first unref ray from the former leaf
846	(*ri).mRay->Unref();
847
848	// determine the side of this ray with respect to the plane
849	int side = node->ComputeRayIntersection(ri, (ri).mRay->mT);
850
851	if (side == 0) {
852	if ((*ri).mRay->HasPosDir(axis)) {
853	back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
854	(*ri).mMinT,
855	(*ri).mRay->mT)
856	);
857	front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
858	(*ri).mRay->mT,
859	(*ri).mMaxT));
860	} else {
861	back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
862	(*ri).mRay->mT,
863	(*ri).mMaxT));
864	front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
865	(*ri).mMinT,
866	(*ri).mRay->mT));
867	}
868	} else
869	if (side == 1)
870	front->AddRay(*ri);
871	else
872	back->AddRay(*ri);
873	} else
874	(*ri).mRay->Unref();
875	}
876	} else {
877	// rays front/back
878
879
880	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
881	ri != leaf->rays.end();
882	ri++) {
883	if ((*ri).mRay->IsActive()) {
884	// first unref ray from the former leaf
885	(*ri).mRay->Unref();
886
887	int side;
888	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
889	side = 1;
890	else
891	side = -1;
892
893	if (side == 1)
894	front->AddRay(*ri);
895	else
896	back->AddRay(*ri);
897
898	} else
899	(*ri).mRay->Unref();
900	}
901	}
902
903	// update stats
904	stat.rayRefs -= leaf->rays.size();
905	stat.rayRefs += raysBack + raysFront;
906
907
908	delete leaf;
909	return node;
910	}
911
912
913
914
915
916
917	int
918	VspKdTree::ReleaseMemory(const int time)
919	{
920	stack<VspKdTreeNode *> tstack;
921
922	// find a node in the tree which subtree will be collapsed
923	int maxAccessTime = time - accessTimeThreshold;
924	int released;
925
926	tstack.push(root);
927
928	while (!tstack.empty()) {
929	VspKdTreeNode *node = tstack.top();
930	tstack.pop();
931
932
933	if (!node->IsLeaf()) {
934	VspKdTreeInterior in = (VspKdTreeInterior )node;
935	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
936	if (in->depth >= minCollapseDepth &&
937	in->lastAccessTime <= maxAccessTime) {
938	released = CollapseSubtree(node, time);
939	break;
940	}
941
942	if (in->back->GetAccessTime() <
943	in->front->GetAccessTime()) {
944	tstack.push(in->front);
945	tstack.push(in->back);
946	} else {
947	tstack.push(in->back);
948	tstack.push(in->front);
949	}
950	}
951	}
952
953	while (tstack.empty()) {
954	// could find node to collaps...
955	// cout<<"Could not find a node to release "<<endl;
956	break;
957	}
958
959	return released;
960	}
961
962
963
964
965	VspKdTreeNode *
966	VspKdTree::SubdivideLeaf(
967	VspKdTreeLeaf *leaf,
968	const float sizeThreshold
969	)
970	{
971	VspKdTreeNode *node = leaf;
972
973	AxisAlignedBox3 leafBBox = GetBBox(leaf);
974
975	static int pass = 0;
976	pass ++;
977
978	// check if we should perform a dynamic subdivision of the leaf
979	if (
980	// leaf->rays.size() > (unsigned)termMinCost &&
981	(leaf->GetPvsSize() >= termMinPvs) &&
982	SqrMagnitude(leafBBox.Size()) > sizeThreshold) {
983
984	// memory check and realese...
985	if (GetMemUsage() > maxTotalMemory) {
986	ReleaseMemory( pass );
987	}
988
989	AxisAlignedBox3 backBBox, frontBBox;
990
991	// subdivide the node
992	node =
993	SubdivideNode(leaf,
994	leafBBox,
995	backBBox,
996	frontBBox
997	);
998	}
999	return node;
1000	}
1001
1002
1003
1004	void
1005	VspKdTree::UpdateRays(VssRayContainer &remove,
1006	VssRayContainer &add
1007	)
1008	{
1009	VspKdTreeLeaf::NewMail();
1010
1011	// schedule rays for removal
1012	for(VssRayContainer::const_iterator ri = remove.begin();
1013	ri != remove.end();
1014	ri++) {
1015	(*ri)->ScheduleForRemoval();
1016	}
1017
1018	int inactive=0;
1019
1020	for(VssRayContainer::const_iterator ri = remove.begin();
1021	ri != remove.end();
1022	ri++) {
1023	if ((*ri)->ScheduledForRemoval())
1024	// RemoveRay(*ri, NULL, false);
1025	// !!! BUG - with true it does not work correctly - aggreated delete
1026	RemoveRay(*ri, NULL, true);
1027	else
1028	inactive++;
1029	}
1030
1031
1032	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1033
1034	for(VssRayContainer::const_iterator ri = add.begin();
1035	ri != add.end();
1036	ri++) {
1037	AddRay(*ri);
1038	}
1039
1040
1041	}
1042
1043
1044	void
1045	VspKdTree::RemoveRay(VssRay *ray,
1046	vector<VspKdTreeLeaf > affectedLeaves,
1047	const bool removeAllScheduledRays
1048	)
1049	{
1050
1051	stack<RayTraversalData> tstack;
1052
1053	tstack.push(RayTraversalData(root, VspKdTreeNode::RayInfo(ray)));
1054
1055	RayTraversalData data;
1056
1057	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1058
1059	while (!tstack.empty()) {
1060	data = tstack.top();
1061	tstack.pop();
1062
1063	if (!data.node->IsLeaf()) {
1064	// split the set of rays in two groups intersecting the
1065	// two subtrees
1066
1067	TraverseInternalNode(data, tstack);
1068
1069	} else {
1070	// remove the ray from the leaf
1071	// find the ray in the leaf and swap it with the last ray...
1072	VspKdTreeLeaf leaf = (VspKdTreeLeaf )data.node;
1073
1074	if (!leaf->Mailed()) {
1075	leaf->Mail();
1076	if (affectedLeaves)
1077	affectedLeaves->push_back(leaf);
1078
1079	if (removeAllScheduledRays) {
1080	int tail = leaf->rays.size()-1;
1081
1082	for (int i=0; i < (int)(leaf->rays.size()); i++) {
1083	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1084	// find a ray to replace it with
1085	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1086	stat.removedRayRefs++;
1087	leaf->rays[tail].mRay->Unref();
1088	leaf->rays.pop_back();
1089	tail--;
1090	}
1091
1092	if (tail < i)
1093	break;
1094
1095	stat.removedRayRefs++;
1096	leaf->rays[i].mRay->Unref();
1097	leaf->rays[i] = leaf->rays[tail];
1098	leaf->rays.pop_back();
1099	tail--;
1100	}
1101	}
1102	}
1103	}
1104
1105	if (!removeAllScheduledRays)
1106	for (int i=0; i < (int)leaf->rays.size(); i++) {
1107	if (leaf->rays[i].mRay == ray) {
1108	stat.removedRayRefs++;
1109	ray->Unref();
1110	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1111	leaf->rays.pop_back();
1112	// check this ray again
1113	break;
1114	}
1115	}
1116
1117	}
1118	}
1119
1120	if (ray->RefCount() != 0) {
1121	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1122	exit(1);
1123	}
1124
1125	}
1126
1127
1128	void
1129	VspKdTree::AddRay(VssRay *ray)
1130	{
1131
1132	stack<RayTraversalData> tstack;
1133
1134	tstack.push(RayTraversalData(root, VspKdTreeNode::RayInfo(ray)));
1135
1136	RayTraversalData data;
1137
1138	while (!tstack.empty()) {
1139	data = tstack.top();
1140	tstack.pop();
1141
1142	if (!data.node->IsLeaf()) {
1143	TraverseInternalNode(data, tstack);
1144	} else {
1145	// remove the ray from the leaf
1146	// find the ray in the leaf and swap it with the last ray...
1147	VspKdTreeLeaf leaf = (VspKdTreeLeaf )data.node;
1148	leaf->AddRay(data.rayData);
1149	stat.addedRayRefs++;
1150	}
1151	}
1152	}
1153
1154	void
1155	VspKdTree::TraverseInternalNode(
1156	RayTraversalData &data,
1157	stack<RayTraversalData> &tstack)
1158	{
1159	VspKdTreeInterior in = (VspKdTreeInterior ) data.node;
1160
1161	if (in->axis <= VspKdTreeNode::SPLIT_Z) {
1162
1163	// determine the side of this ray with respect to the plane
1164	int side = in->ComputeRayIntersection(data.rayData,
1165	data.rayData.mRay->mT);
1166
1167
1168	if (side == 0) {
1169	if (data.rayData.mRay->HasPosDir(in->axis)) {
1170	tstack.push(RayTraversalData(in->back,
1171	VspKdTreeNode::RayInfo(data.rayData.mRay,
1172	data.rayData.mMinT,
1173	data.rayData.mRay->mT))
1174	);
1175
1176	tstack.push(RayTraversalData(in->front,
1177	VspKdTreeNode::RayInfo(data.rayData.mRay,
1178	data.rayData.mRay->mT,
1179	data.rayData.mMaxT
1180	))
1181	);
1182
1183	} else {
1184	tstack.push(RayTraversalData(in->back,
1185	VspKdTreeNode::RayInfo(data.rayData.mRay,
1186	data.rayData.mRay->mT,
1187	data.rayData.mMaxT
1188	))
1189	);
1190
1191	tstack.push(RayTraversalData(in->front,
1192	VspKdTreeNode::RayInfo(data.rayData.mRay,
1193	data.rayData.mMinT,
1194	data.rayData.mRay->mT))
1195	);
1196
1197
1198	}
1199	} else
1200	if (side == 1)
1201	tstack.push(RayTraversalData(in->front, data.rayData));
1202	else
1203	tstack.push(RayTraversalData(in->back, data.rayData));
1204	}
1205	else {
1206	// directional split
1207	if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position)
1208	tstack.push(RayTraversalData(in->front, data.rayData));
1209	else
1210	tstack.push(RayTraversalData(in->back, data.rayData));
1211	}
1212	}
1213
1214
1215	int
1216	VspKdTree::CollapseSubtree(VspKdTreeNode *sroot, const int time)
1217	{
1218	// first count all rays in the subtree
1219	// use mail 1 for this purpose
1220	stack<VspKdTreeNode *> tstack;
1221	int rayCount = 0;
1222	int totalRayCount = 0;
1223	int collapsedNodes = 0;
1224
1225	#if DEBUG_COLLAPSE
1226	cout<<"Collapsing subtree"<<endl;
1227	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1228	cout<<"depth="<<(int)sroot->depth<<endl;
1229	#endif
1230
1231	// tstat.collapsedSubtrees++;
1232	// tstat.collapseDepths += (int)sroot->depth;
1233	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1234
1235	tstack.push(sroot);
1236	VssRay::NewMail();
1237
1238	while (!tstack.empty()) {
1239	collapsedNodes++;
1240	VspKdTreeNode *node = tstack.top();
1241	tstack.pop();
1242
1243	if (node->IsLeaf()) {
1244	VspKdTreeLeaf leaf = (VspKdTreeLeaf ) node;
1245	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1246	ri != leaf->rays.end();
1247	ri++) {
1248
1249	totalRayCount++;
1250	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
1251	(*ri).mRay->Mail();
1252	rayCount++;
1253	}
1254	}
1255	} else {
1256	tstack.push(((VspKdTreeInterior *)node)->back);
1257	tstack.push(((VspKdTreeInterior *)node)->front);
1258	}
1259	}
1260
1261	VssRay::NewMail();
1262
1263	// create a new node that will hold the rays
1264	VspKdTreeLeaf *newLeaf = new VspKdTreeLeaf( sroot->parent, rayCount );
1265	if ( newLeaf->parent )
1266	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1267
1268
1269	tstack.push( sroot );
1270
1271	while (!tstack.empty()) {
1272
1273	VspKdTreeNode *node = tstack.top();
1274	tstack.pop();
1275
1276	if (node->IsLeaf()) {
1277	VspKdTreeLeaf leaf = (VspKdTreeLeaf ) node;
1278
1279	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1280	ri != leaf->rays.end();
1281	ri++) {
1282
1283	// unref this ray from the old node
1284
1285	if ((*ri).mRay->IsActive()) {
1286	(*ri).mRay->Unref();
1287	if (!(*ri).mRay->Mailed()) {
1288	(*ri).mRay->Mail();
1289	newLeaf->AddRay(*ri);
1290	}
1291	} else
1292	(*ri).mRay->Unref();
1293
1294	}
1295	} else {
1296	tstack.push(((VspKdTreeInterior *)node)->back);
1297	tstack.push(((VspKdTreeInterior *)node)->front);
1298	}
1299	}
1300
1301	// delete the node and all its children
1302	delete sroot;
1303
1304	// for(VspKdTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1305	// ri != newLeaf->rays.end();
1306	// ri++)
1307	// (*ri).ray->UnMail(2);
1308
1309
1310	#if DEBUG_COLLAPSE
1311	cout<<"Total memory before="<<GetMemUsage()<<endl;
1312	#endif
1313
1314	stat.nodes -= collapsedNodes - 1;
1315	stat.rayRefs -= totalRayCount - rayCount;
1316
1317	#if DEBUG_COLLAPSE
1318	cout<<"collapsed nodes"<<collapsedNodes<<endl;
1319	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1320	cout<<"Total memory after="<<GetMemUsage()<<endl;
1321	cout<<"================================"<<endl;
1322	#endif
1323
1324	// tstat.collapsedNodes += collapsedNodes;
1325	// tstat.collapsedRays += totalRayCount - rayCount;
1326
1327	return totalRayCount - rayCount;
1328	}
1329
1330
1331	int
1332	VspKdTree::GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const
1333	{
1334	stack<VspKdTreeNode *> tstack;
1335	tstack.push(root);
1336
1337	Intersectable::NewMail();
1338	int pvsSize = 0;
1339
1340	while (!tstack.empty()) {
1341	VspKdTreeNode *node = tstack.top();
1342	tstack.pop();
1343
1344
1345	if (node->IsLeaf()) {
1346	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
1347	for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1348	ri != leaf->rays.end();
1349	ri++)
1350	if ((*ri).mRay->IsActive()) {
1351	Intersectable *object;
1352	#if BIDIRECTIONAL_RAY
1353	object = (*ri).mRay->mOriginObject;
1354	if (object && !object->Mailed()) {
1355	pvsSize++;
1356	object->Mail();
1357	}
1358	#endif
1359	object = (*ri).mRay->mTerminationObject;
1360	if (object && !object->Mailed()) {
1361	pvsSize++;
1362	object->Mail();
1363	}
1364	}
1365	} else {
1366	VspKdTreeInterior in = (VspKdTreeInterior )node;
1367	if (in->axis < 3) {
1368	if (box.Max(in->axis) >= in->position )
1369	tstack.push(in->front);
1370
1371	if (box.Min(in->axis) <= in->position )
1372	tstack.push(in->back);
1373	} else {
1374	// both nodes for directional splits
1375	tstack.push(in->front);
1376	tstack.push(in->back);
1377	}
1378	}
1379	}
1380	return pvsSize;
1381	}
1382
1383	void
1384	VspKdTree::GetRayContributionStatistics(
1385	float &minRayContribution,
1386	float &maxRayContribution,
1387	float &avgRayContribution
1388	)
1389	{
1390	stack<VspKdTreeNode *> tstack;
1391	tstack.push(root);
1392
1393	minRayContribution = 1.0f;
1394	maxRayContribution = 0.0f;
1395	float sumRayContribution = 0.0f;
1396	int leaves = 0;
1397
1398	while (!tstack.empty()) {
1399	VspKdTreeNode *node = tstack.top();
1400	tstack.pop();
1401
1402	if (node->IsLeaf()) {
1403	leaves++;
1404	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
1405	float c = leaf->GetAvgRayContribution();
1406	if (c > maxRayContribution)
1407	maxRayContribution = c;
1408	if (c < minRayContribution)
1409	minRayContribution = c;
1410	sumRayContribution += c;
1411
1412	} else {
1413	VspKdTreeInterior in = (VspKdTreeInterior )node;
1414	// both nodes for directional splits
1415	tstack.push(in->front);
1416	tstack.push(in->back);
1417	}
1418	}
1419
1420	cout<<"sum="<<sumRayContribution<<endl;
1421	cout<<"leaves="<<leaves<<endl;
1422	avgRayContribution = sumRayContribution/(float)leaves;
1423	}
1424
1425
1426	int
1427	VspKdTree::GenerateRays(const float ratioPerLeaf,
1428	SimpleRayContainer &rays)
1429	{
1430	stack<VspKdTreeNode *> tstack;
1431	tstack.push(root);
1432
1433	while (!tstack.empty()) {
1434	VspKdTreeNode *node = tstack.top();
1435	tstack.pop();
1436
1437	if (node->IsLeaf()) {
1438	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
1439	float c = leaf->GetAvgRayContribution();
1440	int num = (c*ratioPerLeaf + 0.5);
1441	// cout<<num<<" ";
1442
1443	for (int i=0; i < num; i++) {
1444	Vector3 origin = GetBBox(leaf).GetRandomPoint();
1445	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1446	Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1447	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1448	rays.push_back(SimpleRay(origin, direction));
1449	}
1450
1451	} else {
1452	VspKdTreeInterior in = (VspKdTreeInterior )node;
1453	// both nodes for directional splits
1454	tstack.push(in->front);
1455	tstack.push(in->back);
1456	}
1457	}
1458
1459	return rays.size();
1460	}
1461
1462
1463	float
1464	VspKdTree::GetAvgPvsSize()
1465	{
1466	stack<VspKdTreeNode *> tstack;
1467	tstack.push(root);
1468
1469	int sumPvs = 0;
1470	int leaves = 0;
1471	while (!tstack.empty()) {
1472	VspKdTreeNode *node = tstack.top();
1473	tstack.pop();
1474
1475	if (node->IsLeaf()) {
1476	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
1477	// update pvs size
1478	leaf->UpdatePvsSize();
1479	sumPvs += leaf->GetPvsSize();
1480	leaves++;
1481	} else {
1482	VspKdTreeInterior in = (VspKdTreeInterior )node;
1483	// both nodes for directional splits
1484	tstack.push(in->front);
1485	tstack.push(in->back);
1486	}
1487	}
1488
1489
1490	return sumPvs/(float)leaves;
1491	}

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

Download in other formats: