Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/VspBspTree.cpp @ 685

Revision 685, 87.5 KB checked in by bittner, 18 years ago (diff)

Line
1	#include <stack>
2	#include <time.h>
3	#include <iomanip>
4
5	#include "Plane3.h"
6	#include "VspBspTree.h"
7	#include "Mesh.h"
8	#include "common.h"
9	#include "ViewCell.h"
10	#include "Environment.h"
11	#include "Polygon3.h"
12	#include "Ray.h"
13	#include "AxisAlignedBox3.h"
14	#include "Exporter.h"
15	#include "Plane3.h"
16	#include "ViewCellBsp.h"
17	#include "ViewCellsManager.h"
18	#include "Beam.h"
19
20	#define USE_FIXEDPOINT_T 0
21
22
23	//-- static members
24
25
26	int VspBspTree::sFrontId = 0;
27	int VspBspTree::sBackId = 0;
28	int VspBspTree::sFrontAndBackId = 0;
29
30
31
32	// pvs penalty can be different from pvs size
33	inline float EvalPvsPenalty(const int pvs,
34	const int lower,
35	const int upper)
36	{
37	// clamp to minmax values
38	if (pvs < lower)
39	return (float)lower;
40	if (pvs > upper)
41	return (float)upper;
42
43	return (float)pvs;
44	}
45
46
47
48
49	/******************************************************************************/
50	/* class VspBspTree implementation */
51	/******************************************************************************/
52
53
54	VspBspTree::VspBspTree():
55	mRoot(NULL),
56	mUseAreaForPvs(false),
57	mCostNormalizer(Limits::Small),
58	mViewCellsManager(NULL),
59	mOutOfBoundsCell(NULL),
60	mStoreRays(false),
61	mRenderCostWeight(0.5),
62	mUseRandomAxis(false),
63	mTimeStamp(1)
64	{
65	bool randomize = false;
66	environment->GetBoolValue("VspBspTree.Construction.randomize", randomize);
67	if (randomize)
68	Randomize(); // initialise random generator for heuristics
69
70	//-- termination criteria for autopartition
71	environment->GetIntValue("VspBspTree.Termination.maxDepth", mTermMaxDepth);
72	environment->GetIntValue("VspBspTree.Termination.minPvs", mTermMinPvs);
73	environment->GetIntValue("VspBspTree.Termination.minRays", mTermMinRays);
74	environment->GetFloatValue("VspBspTree.Termination.minProbability", mTermMinProbability);
75	environment->GetFloatValue("VspBspTree.Termination.maxRayContribution", mTermMaxRayContribution);
76	environment->GetFloatValue("VspBspTree.Termination.minAccRayLenght", mTermMinAccRayLength);
77	environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
78	environment->GetIntValue("VspBspTree.Termination.missTolerance", mTermMissTolerance);
79	environment->GetIntValue("VspBspTree.Termination.maxViewCells", mMaxViewCells);
80
81	//-- max cost ratio for early tree termination
82	environment->GetFloatValue("VspBspTree.Termination.maxCostRatio", mTermMaxCostRatio);
83
84	environment->GetFloatValue("VspBspTree.Termination.minGlobalCostRatio", mTermMinGlobalCostRatio);
85	environment->GetIntValue("VspBspTree.Termination.globalCostMissTolerance", mTermGlobalCostMissTolerance);
86
87	// HACK//mTermMinPolygons = 25;
88
89	//-- factors for bsp tree split plane heuristics
90	environment->GetFloatValue("VspBspTree.Factor.pvs", mPvsFactor);
91	environment->GetFloatValue("VspBspTree.Termination.ct_div_ci", mCtDivCi);
92
93
94	//-- partition criteria
95	environment->GetIntValue("VspBspTree.maxPolyCandidates", mMaxPolyCandidates);
96	environment->GetIntValue("VspBspTree.maxRayCandidates", mMaxRayCandidates);
97	environment->GetIntValue("VspBspTree.splitPlaneStrategy", mSplitPlaneStrategy);
98
99	environment->GetFloatValue("VspBspTree.Construction.epsilon", mEpsilon);
100	environment->GetIntValue("VspBspTree.maxTests", mMaxTests);
101
102	// if only the driving axis is used for axis aligned split
103	environment->GetBoolValue("VspBspTree.splitUseOnlyDrivingAxis", mOnlyDrivingAxis);
104
105	//-- termination criteria for axis aligned split
106	environment->GetFloatValue("VspBspTree.Termination.AxisAligned.maxRayContribution",
107	mTermMaxRayContriForAxisAligned);
108	environment->GetIntValue("VspBspTree.Termination.AxisAligned.minRays",
109	mTermMinRaysForAxisAligned);
110
111	//environment->GetFloatValue("VspBspTree.maxTotalMemory", mMaxTotalMemory);
112	environment->GetFloatValue("VspBspTree.maxStaticMemory", mMaxMemory);
113
114	environment->GetFloatValue("VspBspTree.Construction.renderCostWeight", mRenderCostWeight);
115	environment->GetBoolValue("VspBspTree.usePolygonSplitIfAvailable", mUsePolygonSplitIfAvailable);
116
117	environment->GetBoolValue("VspBspTree.useCostHeuristics", mUseCostHeuristics);
118	environment->GetBoolValue("VspBspTree.useSplitCostQueue", mUseSplitCostQueue);
119	environment->GetBoolValue("VspBspTree.simulateOctree", mSimulateOctree);
120	environment->GetBoolValue("VspBspTree.useRandomAxis", mUseRandomAxis);
121	environment->GetBoolValue("VspBspTree.useBreathFirstSplits", mBreathFirstSplits);
122
123	environment->GetBoolValue("ViewCells.PostProcess.emptyViewCellsMerge", mEmptyViewCellsMergeAllowed);
124
125	char subdivisionStatsLog[100];
126	environment->GetStringValue("VspBspTree.subdivisionStats", subdivisionStatsLog);
127	mSubdivisionStats.open(subdivisionStatsLog);
128
129	//-- debug output
130
131	Debug << "***** VSP BSP options ****** " << endl;
132	Debug << "max depth: " << mTermMaxDepth << endl;
133	Debug << "min PVS: " << mTermMinPvs << endl;
134	Debug << "min probabiliy: " << mTermMinProbability << endl;
135	Debug << "min rays: " << mTermMinRays << endl;
136	Debug << "max ray contri: " << mTermMaxRayContribution << endl;
137	Debug << "max cost ratio: " << mTermMaxCostRatio << endl;
138	Debug << "miss tolerance: " << mTermMissTolerance << endl;
139	Debug << "max view cells: " << mMaxViewCells << endl;
140	Debug << "max polygon candidates: " << mMaxPolyCandidates << endl;
141	//Debug << "max plane candidates: " << mMaxRayCandidates << endl;
142	Debug << "randomize: " << randomize << endl;
143
144	Debug << "using area for pvs: " << mUseAreaForPvs << endl;
145	Debug << "render cost weight: " << mRenderCostWeight << endl;
146	Debug << "min global cost ratio: " << mTermMinGlobalCostRatio << endl;
147	Debug << "global cost miss tolerance: " << mTermGlobalCostMissTolerance << endl;
148	Debug << "only driving axis: " << mOnlyDrivingAxis << endl;
149	Debug << "max memory: " << mMaxMemory << endl;
150	Debug << "use poly split if available: " << mUsePolygonSplitIfAvailable << endl;
151	Debug << "use cost heuristics: " << mUseCostHeuristics << endl;
152	Debug << "use split cost queue: " << mUseSplitCostQueue << endl;
153	Debug << "subdivision stats log: " << subdivisionStatsLog << endl;
154	Debug << "use random axis: " << mUseRandomAxis << endl;
155	Debug << "breath first splits: " << mBreathFirstSplits << endl;
156	Debug << "empty view cells merge: " << mEmptyViewCellsMergeAllowed << endl;
157
158	Debug << "octree: " << mSimulateOctree << endl;
159
160	Debug << "Split plane strategy: ";
161
162	if (mSplitPlaneStrategy & RANDOM_POLYGON)
163	{
164	Debug << "random polygon ";
165	}
166	if (mSplitPlaneStrategy & AXIS_ALIGNED)
167	{
168	Debug << "axis aligned ";
169	}
170	if (mSplitPlaneStrategy & LEAST_RAY_SPLITS)
171	{
172	mCostNormalizer += mLeastRaySplitsFactor;
173	Debug << "least ray splits ";
174	}
175	if (mSplitPlaneStrategy & BALANCED_RAYS)
176	{
177	mCostNormalizer += mBalancedRaysFactor;
178	Debug << "balanced rays ";
179	}
180	if (mSplitPlaneStrategy & PVS)
181	{
182	mCostNormalizer += mPvsFactor;
183	Debug << "pvs";
184	}
185
186
187	mSplitCandidates = new vector<SortableEntry>;
188
189	Debug << endl;
190	}
191
192
193	BspViewCell *VspBspTree::GetOutOfBoundsCell()
194	{
195	return mOutOfBoundsCell;
196	}
197
198
199	BspViewCell *VspBspTree::GetOrCreateOutOfBoundsCell()
200	{
201	if (!mOutOfBoundsCell)
202	{
203	mOutOfBoundsCell = new BspViewCell();
204	mOutOfBoundsCell->SetId(-1);
205	mOutOfBoundsCell->SetValid(false);
206	}
207
208	return mOutOfBoundsCell;
209	}
210
211
212	const BspTreeStatistics &VspBspTree::GetStatistics() const
213	{
214	return mBspStats;
215	}
216
217
218	VspBspTree::~VspBspTree()
219	{
220	DEL_PTR(mRoot);
221	DEL_PTR(mSplitCandidates);
222	}
223
224
225	int VspBspTree::AddMeshToPolygons(Mesh *mesh,
226	PolygonContainer &polys,
227	MeshInstance *parent)
228	{
229	FaceContainer::const_iterator fi;
230
231	// copy the face data to polygons
232	for (fi = mesh->mFaces.begin(); fi != mesh->mFaces.end(); ++ fi)
233	{
234	Polygon3 poly = new Polygon3((fi), mesh);
235
236	if (poly->Valid(mEpsilon))
237	{
238	poly->mParent = parent; // set parent intersectable
239	polys.push_back(poly);
240	}
241	else
242	DEL_PTR(poly);
243	}
244	return (int)mesh->mFaces.size();
245	}
246
247
248	int VspBspTree::AddToPolygonSoup(const ViewCellContainer &viewCells,
249	PolygonContainer &polys,
250	int maxObjects)
251	{
252	int limit = (maxObjects > 0) ?
253	Min((int)viewCells.size(), maxObjects) : (int)viewCells.size();
254
255	int polysSize = 0;
256
257	for (int i = 0; i < limit; ++ i)
258	{
259	if (viewCells[i]->GetMesh()) // copy the mesh data to polygons
260	{
261	mBox.Include(viewCells[i]->GetBox()); // add to BSP tree aabb
262	polysSize +=
263	AddMeshToPolygons(viewCells[i]->GetMesh(),
264	polys,
265	viewCells[i]);
266	}
267	}
268
269	return polysSize;
270	}
271
272
273	int VspBspTree::AddToPolygonSoup(const ObjectContainer &objects,
274	PolygonContainer &polys,
275	int maxObjects)
276	{
277	int limit = (maxObjects > 0) ?
278	Min((int)objects.size(), maxObjects) : (int)objects.size();
279
280	for (int i = 0; i < limit; ++i)
281	{
282	Intersectable object = objects[i];//it;
283	Mesh *mesh = NULL;
284
285	switch (object->Type()) // extract the meshes
286	{
287	case Intersectable::MESH_INSTANCE:
288	mesh = dynamic_cast<MeshInstance *>(object)->GetMesh();
289	break;
290	case Intersectable::VIEW_CELL:
291	mesh = dynamic_cast<ViewCell *>(object)->GetMesh();
292	break;
293	// TODO: handle transformed mesh instances
294	default:
295	Debug << "intersectable type not supported" << endl;
296	break;
297	}
298
299	if (mesh) // copy the mesh data to polygons
300	{
301	mBox.Include(object->GetBox()); // add to BSP tree aabb
302	AddMeshToPolygons(mesh, polys, NULL);
303	}
304	}
305
306	return (int)polys.size();
307	}
308
309
310	void VspBspTree::Construct(const VssRayContainer &sampleRays,
311	AxisAlignedBox3 *forcedBoundingBox)
312	{
313	mBspStats.nodes = 1;
314	mBox.Initialize(); // initialise BSP tree bounding box
315
316	if (forcedBoundingBox)
317	mBox = *forcedBoundingBox;
318
319	PolygonContainer polys;
320	RayInfoContainer *rays = new RayInfoContainer();
321
322	VssRayContainer::const_iterator rit, rit_end = sampleRays.end();
323
324	long startTime = GetTime();
325
326	cout << "Extracting polygons from rays ... ";
327
328	Intersectable::NewMail();
329
330	int numObj = 0;
331
332	//-- extract polygons intersected by the rays
333	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
334	{
335	VssRay ray = rit;
336
337	if ((mBox.IsInside(ray->mTermination) \|\| !forcedBoundingBox) &&
338	ray->mTerminationObject &&
339	!ray->mTerminationObject->Mailed())
340	{
341	ray->mTerminationObject->Mail();
342	MeshInstance obj = dynamic_cast<MeshInstance >(ray->mTerminationObject);
343	AddMeshToPolygons(obj->GetMesh(), polys, obj);
344	++ numObj;
345
346	//-- compute bounding box
347	if (!forcedBoundingBox)
348	mBox.Include(ray->mTermination);
349	}
350
351	if ((mBox.IsInside(ray->mOrigin) \|\| !forcedBoundingBox) &&
352	ray->mOriginObject &&
353	!ray->mOriginObject->Mailed())
354	{
355	ray->mOriginObject->Mail();
356	MeshInstance obj = dynamic_cast<MeshInstance >(ray->mOriginObject);
357	AddMeshToPolygons(obj->GetMesh(), polys, obj);
358	++ numObj;
359
360	//-- compute bounding box
361	if (!forcedBoundingBox)
362	mBox.Include(ray->mOrigin);
363	}
364	}
365
366	Debug << "maximal pvs (i.e., pvs still considered as valid) : "
367	<< mViewCellsManager->GetMaxPvsSize() << endl;
368
369	//-- store rays
370	for (rit = sampleRays.begin(); rit != rit_end; ++ rit)
371	{
372	VssRay ray = rit;
373
374	float minT, maxT;
375
376	static Ray hray;
377	hray.Init(*ray);
378
379	// TODO: not very efficient to implictly cast between rays types
380	if (mBox.GetRaySegment(hray, minT, maxT))
381	{
382	float len = ray->Length();
383
384	if (!len)
385	len = Limits::Small;
386
387	rays->push_back(RayInfo(ray, minT / len, maxT / len));
388	}
389	}
390
391	// normalize
392	if (mUseAreaForPvs)
393	mTermMinProbability *= mBox.SurfaceArea();
394	else
395	mTermMinProbability *= mBox.GetVolume();
396
397	// throw out unnecessary polygons
398	PreprocessPolygons(polys);
399
400	mBspStats.polys = (int)polys.size();
401	mGlobalCostMisses = 0;
402
403	cout << "finished" << endl;
404
405	Debug << "\nPolygon extraction: " << (int)polys.size() << " polys extracted from "
406	<< (int)sampleRays.size() << " rays in "
407	<< TimeDiff(startTime, GetTime())*1e-3 << " secs" << endl << endl;
408
409	// use split cost priority queue
410	if (mUseSplitCostQueue)
411	{
412	ConstructWithSplitQueue(polys, rays);
413	}
414	else
415	{
416	Construct(polys, rays);
417	}
418
419	// clean up polygons
420	CLEAR_CONTAINER(polys);
421	}
422
423
424	// TODO: return memory usage in MB
425	float VspBspTree::GetMemUsage() const
426	{
427	return (float)
428	(sizeof(VspBspTree) +
429	mBspStats.Leaves() * sizeof(BspLeaf) +
430	mCreatedViewCells * sizeof(BspViewCell) +
431	mBspStats.pvs * sizeof(ObjectPvsData) +
432	mBspStats.Interior() * sizeof(BspInterior) +
433	mBspStats.accumRays * sizeof(RayInfo)) / (1024.0f * 1024.0f);
434	}
435
436
437
438	void VspBspTree::Construct(const PolygonContainer &polys, RayInfoContainer *rays)
439	{
440	VspBspTraversalQueue tQueue;
441
442	mRoot = new BspLeaf();
443
444	// constrruct root node geometry
445	BspNodeGeometry *geom = new BspNodeGeometry();
446	ConstructGeometry(mRoot, *geom);
447
448	const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();
449
450	VspBspTraversalData tData(mRoot,
451	new PolygonContainer(polys),
452	0,
453	rays,
454	ComputePvsSize(*rays),
455	prop,
456	geom);
457
458	EvalPriority(tData);
459
460
461	if (mSimulateOctree)
462	tData.mAxis = 0;
463
464	// first node is kd node, i.e. an axis aligned box
465	if (0)
466	tData.mIsKdNode = true;
467	else
468	tData.mIsKdNode = false;
469
470	tQueue.push(tData);
471
472
473	mTotalCost = tData.mPvs * tData.mProbability / mBox.GetVolume();
474	mTotalPvsSize = tData.mPvs;
475
476	mSubdivisionStats
477	<< "#ViewCells\n1\n" << endl
478	<< "#RenderCostDecrease\n0\n" << endl
479	<< "#TotalRenderCost\n" << mTotalCost << endl
480	<< "#AvgRenderCost\n" << mTotalPvsSize << endl;
481
482	Debug << "total cost: " << mTotalCost << endl;
483
484
485	mBspStats.Start();
486	cout << "Constructing vsp bsp tree ... \n";
487
488	long startTime = GetTime();
489	int nLeaves = 500;
490	int nViewCells = 500;
491
492	// used for intermediate time measurements and progress
493	long interTime = GetTime();
494
495	mOutOfMemory = false;
496
497	mCreatedViewCells = 0;
498
499	while (!tQueue.empty())
500	{
501	tData = tQueue.top();
502	tQueue.pop();
503
504	if (0 && !mOutOfMemory)
505	{
506	float mem = GetMemUsage();
507
508	if (mem > mMaxMemory)
509	{
510	mOutOfMemory = true;
511	Debug << "memory limit reached: " << mem << endl;
512	}
513	}
514
515	// subdivide leaf node
516	BspNode *r = Subdivide(tQueue, tData);
517
518	if (r == mRoot)
519	Debug << "VSP BSP tree construction time spent at root: "
520	<< TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
521
522	if (mBspStats.Leaves() >= nLeaves)
523	{
524	nLeaves += 500;
525
526	cout << "leaves=" << mBspStats.Leaves() << endl;
527	Debug << "needed "
528	<< TimeDiff(interTime, GetTime())*1e-3
529	<< " secs to create 500 view cells" << endl;
530	interTime = GetTime();
531	}
532
533	if (mCreatedViewCells >= nViewCells)
534	{
535	nViewCells += 500;
536
537	cout << "generated " << mCreatedViewCells << " viewcells" << endl;
538	}
539	}
540
541	Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
542	cout << "finished\n";
543
544	mBspStats.Stop();
545	}
546
547
548
549	void VspBspTree::ConstructWithSplitQueue(const PolygonContainer &polys,
550	RayInfoContainer *rays)
551	{
552	VspBspSplitQueue tQueue;
553
554	mRoot = new BspLeaf();
555
556	// constrruct root node geometry
557	BspNodeGeometry *geom = new BspNodeGeometry();
558	ConstructGeometry(mRoot, *geom);
559
560	const float prop = mUseAreaForPvs ? geom->GetArea() : geom->GetVolume();
561
562	VspBspTraversalData tData(mRoot,
563	new PolygonContainer(polys),
564	0,
565	rays,
566	ComputePvsSize(*rays),
567	prop,
568	geom);
569
570
571
572	// compute first split candidate
573	VspBspSplitCandidate splitCandidate;
574	EvalSplitCandidate(tData, splitCandidate);
575
576	tQueue.push(splitCandidate);
577
578	mTotalCost = tData.mPvs * tData.mProbability / mBox.GetVolume();
579	mTotalPvsSize = tData.mPvs;
580
581	mSubdivisionStats
582	<< "#ViewCells\n1\n" << endl
583	<< "#RenderCostDecrease\n0\n" << endl
584	<< "#dummy\n0\n" << endl
585	<< "#TotalRenderCost\n" << mTotalCost << endl
586	<< "#AvgRenderCost\n" << mTotalPvsSize << endl;
587
588	Debug << "total cost: " << mTotalCost << endl;
589
590
591	mBspStats.Start();
592	cout << "Constructing vsp bsp tree ... \n";
593
594	long startTime = GetTime();
595	int nLeaves = 500;
596	int nViewCells = 500;
597
598	// used for intermediate time measurements and progress
599	long interTime = GetTime();
600
601	mOutOfMemory = false;
602
603	mCreatedViewCells = 0;
604
605	while (!tQueue.empty())
606	{
607	splitCandidate = tQueue.top();
608	tQueue.pop();
609
610	// cost ratio of cost decrease / totalCost
611	float costRatio = splitCandidate.GetCost() / mTotalCost;
612
613	//Debug << "cost ratio: " << costRatio << endl;
614
615	if (costRatio < mTermMinGlobalCostRatio)
616	++ mGlobalCostMisses;
617
618	if (0 && !mOutOfMemory)
619	{
620	float mem = GetMemUsage();
621
622	if (mem > mMaxMemory)
623	{
624	mOutOfMemory = true;
625	Debug << "memory limit reached: " << mem << endl;
626	}
627	}
628
629	// subdivide leaf node
630	BspNode *r = Subdivide(tQueue, splitCandidate);
631
632	if (r == mRoot)
633	Debug << "VSP BSP tree construction time spent at root: "
634	<< TimeDiff(startTime, GetTime())*1e-3 << "s" << endl;
635
636	if (mBspStats.Leaves() >= nLeaves)
637	{
638	nLeaves += 500;
639
640	cout << "leaves=" << mBspStats.Leaves() << endl;
641	Debug << "needed "
642	<< TimeDiff(interTime, GetTime())*1e-3
643	<< " secs to create 500 view cells" << endl;
644	interTime = GetTime();
645	}
646
647	if (mCreatedViewCells == nViewCells)
648	{
649	nViewCells += 500;
650
651	cout << "generated " << mCreatedViewCells << " viewcells" << endl;
652	}
653	}
654
655	Debug << "Used Memory: " << GetMemUsage() << " MB" << endl << endl;
656	cout << "finished\n";
657
658	mBspStats.Stop();
659	}
660
661
662	bool VspBspTree::LocalTerminationCriteriaMet(const VspBspTraversalData &data) const
663	{
664	return
665	(((int)data.mRays->size() <= mTermMinRays) \|\|
666	(data.mPvs <= mTermMinPvs) \|\|
667	(data.mProbability <= mTermMinProbability) \|\|
668	(data.GetAvgRayContribution() > mTermMaxRayContribution) \|\|
669	(data.mDepth >= mTermMaxDepth));
670	}
671
672
673	bool VspBspTree::GlobalTerminationCriteriaMet(const VspBspTraversalData &data) const
674	{
675	return
676	(mOutOfMemory
677	\|\| (mBspStats.Leaves() >= mMaxViewCells)
678	\|\| (mGlobalCostMisses >= mTermGlobalCostMissTolerance)
679	);
680	}
681
682
683
684	BspNode *VspBspTree::Subdivide(VspBspTraversalQueue &tQueue,
685	VspBspTraversalData &tData)
686	{
687	BspNode *newNode = tData.mNode;
688
689	if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
690	{
691	PolygonContainer coincident;
692
693	VspBspTraversalData tFrontData;
694	VspBspTraversalData tBackData;
695
696	if (mSimulateOctree)
697	{ // choose axes in circular motion
698	tFrontData.mAxis = (tData.mAxis + 1) % 3;
699	tBackData.mAxis = (tData.mAxis + 1) % 3;
700	}
701
702	// create new interior node and two leaf nodes
703	// or return leaf as it is (if maxCostRatio missed)
704	int splitAxis;
705	bool splitFurther = true;
706	int maxCostMisses = tData.mMaxCostMisses;
707
708	Plane3 splitPlane;
709	BspLeaf leaf = dynamic_cast<BspLeaf >(tData.mNode);
710
711	if (!SelectPlane(splitPlane, leaf, tData, tFrontData, tBackData, splitAxis))
712	{
713	++ maxCostMisses;
714
715	if (maxCostMisses > mTermMissTolerance)
716	{
717	// terminate branch because of max cost
718	++ mBspStats.maxCostNodes;
719	splitFurther = false;
720	}
721	}
722
723	if (splitFurther) //-- continue subdivision
724	{
725	newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);
726
727	if (splitAxis < 3)
728	++ mBspStats.splits[splitAxis];
729	else
730	++ mBspStats.polySplits;
731
732	tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && splitAxis < 3);
733
734	// how often was max cost ratio missed in this branch?
735	tFrontData.mMaxCostMisses = maxCostMisses;
736	tBackData.mMaxCostMisses = maxCostMisses;
737
738	EvalPriority(tFrontData);
739	EvalPriority(tBackData);
740
741	if (1)
742	{
743
744	float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
745	float cBack = (float)tBackData.mPvs * tBackData.mProbability;
746	float cData = (float)tData.mPvs * tData.mProbability;;
747
748	float costDecr = (cFront + cBack - cData) / mBox.GetVolume();
749
750	mTotalCost += costDecr;
751	mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
752
753	mSubdivisionStats
754	<< "#ViewCells\n" << mBspStats.Leaves() << endl
755	<< "#RenderCostDecrease\n" << -costDecr << endl
756	<< "#TotalRenderCost\n" << mTotalCost << endl
757	<< "#AvgRenderCost\n" << (float)mTotalPvsSize / (float)mBspStats.Leaves() << endl;
758	}
759
760	// push the children on the stack
761	tQueue.push(tFrontData);
762	tQueue.push(tBackData);
763
764	// delete old leaf node
765	DEL_PTR(tData.mNode);
766	}
767	}
768
769	//-- terminate traversal and create new view cell
770	if (newNode->IsLeaf())
771	{
772	BspLeaf leaf = dynamic_cast<BspLeaf >(newNode);
773	BspViewCell *viewCell = new BspViewCell();
774
775	leaf->SetViewCell(viewCell);
776
777	//-- update pvs
778	int conSamp = 0;
779	float sampCon = 0.0f;
780	AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
781
782	mBspStats.contributingSamples += conSamp;
783	mBspStats.sampleContributions +=(int) sampCon;
784
785	//-- store additional info
786	if (mStoreRays)
787	{
788	RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
789	for (it = tData.mRays->begin(); it != it_end; ++ it)
790	{
791	(*it).mRay->Ref();
792	leaf->mVssRays.push_back((*it).mRay);
793	}
794	}
795
796	// should I check here?
797	if (0 && !mViewCellsManager->CheckValidity(viewCell, 0, mViewCellsManager->GetMaxPvsSize()))
798	{
799	viewCell->SetValid(false);
800	leaf->SetTreeValid(false);
801	PropagateUpValidity(leaf);
802
803	++ mBspStats.invalidLeaves;
804	}
805
806	viewCell->mLeaf = leaf;
807
808	if (mUseAreaForPvs)
809	viewCell->SetArea(tData.mProbability);
810	else
811	viewCell->SetVolume(tData.mProbability);
812
813	leaf->mProbability = tData.mProbability;
814
815	EvaluateLeafStats(tData);
816	}
817
818	//-- cleanup
819	tData.Clear();
820
821	return newNode;
822	}
823
824
825	BspNode *VspBspTree::Subdivide(VspBspSplitQueue &tQueue,
826	VspBspSplitCandidate &splitCandidate)
827	{
828	VspBspTraversalData &tData = splitCandidate.mParentData;
829
830	BspNode *newNode = tData.mNode;
831
832	if (!LocalTerminationCriteriaMet(tData) && !GlobalTerminationCriteriaMet(tData))
833	{
834	PolygonContainer coincident;
835
836	VspBspTraversalData tFrontData;
837	VspBspTraversalData tBackData;
838
839	//-- continue subdivision
840	// create new interior node and two leaf node
841	const Plane3 splitPlane = splitCandidate.mSplitPlane;
842
843	newNode = SubdivideNode(splitPlane, tData, tFrontData, tBackData, coincident);
844
845	const int splitAxis = splitCandidate.mSplitAxis;
846	const int maxCostMisses = splitCandidate.mMaxCostMisses;
847
848	if (splitAxis < 3)
849	++ mBspStats.splits[splitAxis];
850	else
851	++ mBspStats.polySplits;
852
853	tFrontData.mIsKdNode = tBackData.mIsKdNode = (tData.mIsKdNode && splitAxis < 3);
854
855	// how often was max cost ratio missed in this branch?
856	tFrontData.mMaxCostMisses = maxCostMisses;
857	tBackData.mMaxCostMisses = maxCostMisses;
858
859
860	if (1)
861	{
862	float cFront = (float)tFrontData.mPvs * tFrontData.mProbability;
863	float cBack = (float)tBackData.mPvs * tBackData.mProbability;
864	float cData = (float)tData.mPvs * tData.mProbability;
865
866
867	float costDecr =
868	(cFront + cBack - cData) / mBox.GetVolume();
869
870	mTotalCost += costDecr;
871	mTotalPvsSize += tFrontData.mPvs + tBackData.mPvs - tData.mPvs;
872
873	mSubdivisionStats
874	<< "#ViewCells\n" << mBspStats.Leaves() << endl
875	<< "#RenderCostDecrease\n" << -costDecr << endl
876	<< "#dummy\n" << splitCandidate.GetCost() << endl
877	<< "#TotalRenderCost\n" << mTotalCost << endl
878	<< "#AvgRenderCost\n" << (float)mTotalPvsSize / (float)mBspStats.Leaves() << endl;
879	}
880
881
882	//-- push the new split candidates on the stack
883	VspBspSplitCandidate frontCandidate;
884	VspBspSplitCandidate backCandidate;
885
886	EvalSplitCandidate(tFrontData, frontCandidate);
887	EvalSplitCandidate(tBackData, backCandidate);
888
889	tQueue.push(frontCandidate);
890	tQueue.push(backCandidate);
891
892	// delete old leaf node
893	DEL_PTR(tData.mNode);
894	}
895
896
897	//-- terminate traversal and create new view cell
898	if (newNode->IsLeaf())
899	{
900	BspLeaf leaf = dynamic_cast<BspLeaf >(newNode);
901	BspViewCell *viewCell = new BspViewCell();
902
903	leaf->SetViewCell(viewCell);
904
905	//-- update pvs
906	int conSamp = 0;
907	float sampCon = 0.0f;
908	AddToPvs(leaf, *tData.mRays, sampCon, conSamp);
909
910	mBspStats.contributingSamples += conSamp;
911	mBspStats.sampleContributions +=(int) sampCon;
912
913	//-- store additional info
914	if (mStoreRays)
915	{
916	RayInfoContainer::const_iterator it, it_end = tData.mRays->end();
917	for (it = tData.mRays->begin(); it != it_end; ++ it)
918	{
919	(*it).mRay->Ref();
920	leaf->mVssRays.push_back((*it).mRay);
921	}
922	}
923
924	// should I check here?
925	viewCell->mLeaf = leaf;
926
927	if (mUseAreaForPvs)
928	viewCell->SetArea(tData.mProbability);
929	else
930	viewCell->SetVolume(tData.mProbability);
931
932	leaf->mProbability = tData.mProbability;
933
934	EvaluateLeafStats(tData);
935	}
936
937	//-- cleanup
938	tData.Clear();
939
940	return newNode;
941	}
942
943
944	void VspBspTree::EvalPriority(VspBspTraversalData &tData) const
945	{
946	tData.mPriority = mBreathFirstSplits ?
947	(float)-tData.mDepth : tData.mPvs * tData.mProbability;
948	//cout << "priority: " << tData.mPriority << endl;
949	}
950
951
952	void VspBspTree::EvalSplitCandidate(VspBspTraversalData &tData,
953	VspBspSplitCandidate &splitData)
954	{
955	VspBspTraversalData frontData;
956	VspBspTraversalData backData;
957
958	BspLeaf leaf = dynamic_cast<BspLeaf >(tData.mNode);
959
960	// compute locally best split plane
961	bool success = SelectPlane(splitData.mSplitPlane, leaf, tData,
962	frontData, backData, splitData.mSplitAxis);
963
964	// TODO: reuse
965	DEL_PTR(frontData.mGeometry);
966	DEL_PTR(backData.mGeometry);
967
968	// compute global decrease in render cost
969	splitData.mRenderCost = EvalRenderCostDecrease(splitData.mSplitPlane, tData);
970	splitData.mParentData = tData;
971	splitData.mMaxCostMisses = success ? tData.mMaxCostMisses : tData.mMaxCostMisses + 1;
972	}
973
974
975	BspInterior *VspBspTree::SubdivideNode(const Plane3 &splitPlane,
976	VspBspTraversalData &tData,
977	VspBspTraversalData &frontData,
978	VspBspTraversalData &backData,
979	PolygonContainer &coincident)
980	{
981	BspLeaf leaf = dynamic_cast<BspLeaf >(tData.mNode);
982
983	//-- the front and back traversal data is filled with the new values
984	frontData.mDepth = tData.mDepth + 1;
985	frontData.mPolygons = new PolygonContainer();
986	frontData.mRays = new RayInfoContainer();
987
988	backData.mDepth = tData.mDepth + 1;
989	backData.mPolygons = new PolygonContainer();
990	backData.mRays = new RayInfoContainer();
991
992
993	//-- subdivide rays
994	SplitRays(splitPlane,
995	*tData.mRays,
996	*frontData.mRays,
997	*backData.mRays);
998
999
1000	// compute pvs
1001	frontData.mPvs = ComputePvsSize(*frontData.mRays);
1002	backData.mPvs = ComputePvsSize(*backData.mRays);
1003
1004	// split front and back node geometry and compute area
1005
1006	// if geometry was not already computed
1007	if (!frontData.mGeometry && !backData.mGeometry)
1008	{
1009	frontData.mGeometry = new BspNodeGeometry();
1010	backData.mGeometry = new BspNodeGeometry();
1011
1012	tData.mGeometry->SplitGeometry(*frontData.mGeometry,
1013	*backData.mGeometry,
1014	splitPlane,
1015	mBox,
1016	//0.0f);
1017	mEpsilon);
1018
1019	if (mUseAreaForPvs)
1020	{
1021	frontData.mProbability = frontData.mGeometry->GetArea();
1022	backData.mProbability = backData.mGeometry->GetArea();
1023	}
1024	else
1025	{
1026	frontData.mProbability = frontData.mGeometry->GetVolume();
1027	backData.mProbability = tData.mProbability - frontData.mProbability;
1028
1029	if (frontData.mProbability < -0.00001)
1030	Debug << "here2 error f: " << frontData.mProbability << endl;
1031	if (backData.mProbability < -0.00001)
1032	Debug << "here2 error b: " << backData.mProbability << endl;
1033
1034	// clamp because of precision issues
1035	if (0)
1036	{
1037	if (frontData.mProbability < 0) frontData.mProbability = 0;
1038	if (backData.mProbability < 0) backData.mProbability = 0;
1039	}
1040	}
1041	}
1042
1043
1044	// subdivide polygons
1045	SplitPolygons(splitPlane,
1046	*tData.mPolygons,
1047	*frontData.mPolygons,
1048	*backData.mPolygons,
1049	coincident);
1050
1051
1052
1053	///////////////////////////////////////
1054	// subdivide further
1055
1056	mBspStats.nodes += 2;
1057
1058
1059	BspInterior *interior = new BspInterior(splitPlane);
1060
1061	#ifdef _DEBUG
1062	Debug << interior << endl;
1063	#endif
1064
1065	//-- create front and back leaf
1066
1067	BspInterior *parent = leaf->GetParent();
1068
1069	// replace a link from node's parent
1070	if (parent)
1071	{
1072	parent->ReplaceChildLink(leaf, interior);
1073	interior->SetParent(parent);
1074	}
1075	else // new root
1076	{
1077	mRoot = interior;
1078	}
1079
1080	// and setup child links
1081	interior->SetupChildLinks(new BspLeaf(interior), new BspLeaf(interior));
1082
1083	frontData.mNode = interior->GetFront();
1084	backData.mNode = interior->GetBack();
1085
1086	interior->mTimeStamp = mTimeStamp ++;
1087
1088	//DEL_PTR(leaf);
1089	return interior;
1090	}
1091
1092
1093	void VspBspTree::AddToPvs(BspLeaf *leaf,
1094	const RayInfoContainer &rays,
1095	float &sampleContributions,
1096	int &contributingSamples)
1097	{
1098	sampleContributions = 0;
1099	contributingSamples = 0;
1100
1101	RayInfoContainer::const_iterator it, it_end = rays.end();
1102
1103	ViewCell *vc = leaf->GetViewCell();
1104
1105	// add contributions from samples to the PVS
1106	for (it = rays.begin(); it != it_end; ++ it)
1107	{
1108	float sc = 0.0f;
1109	VssRay ray = (it).mRay;
1110	bool madeContrib = false;
1111	float contribution;
1112	if (ray->mTerminationObject) {
1113	if (vc->GetPvs().AddSample(ray->mTerminationObject, ray->mPdf, contribution))
1114	madeContrib = true;
1115	sc += contribution;
1116	}
1117
1118	if (ray->mOriginObject) {
1119	if (vc->GetPvs().AddSample(ray->mOriginObject, ray->mPdf, contribution))
1120	madeContrib = true;
1121	sc += contribution;
1122	}
1123
1124	sampleContributions += sc;
1125	if (madeContrib)
1126	++ contributingSamples;
1127
1128	//leaf->mVssRays.push_back(ray);
1129	}
1130	}
1131
1132
1133	void VspBspTree::SortSplitCandidates(const RayInfoContainer &rays, const int axis)
1134	{
1135	mSplitCandidates->clear();
1136
1137	int requestedSize = 2 * (int)(rays.size());
1138	// creates a sorted split candidates array
1139	if (mSplitCandidates->capacity() > 500000 &&
1140	requestedSize < (int)(mSplitCandidates->capacity() / 10) )
1141	{
1142	delete mSplitCandidates;
1143	mSplitCandidates = new vector<SortableEntry>;
1144	}
1145
1146	mSplitCandidates->reserve(requestedSize);
1147
1148	// insert all queries
1149	for (RayInfoContainer::const_iterator ri = rays.begin(); ri < rays.end(); ++ ri)
1150	{
1151	const bool positive = (*ri).mRay->HasPosDir(axis);
1152	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
1153	(ri).ExtrapOrigin(axis), (ri).mRay));
1154
1155	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
1156	(ri).ExtrapTermination(axis), (ri).mRay));
1157	}
1158
1159	stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
1160	}
1161
1162
1163	float VspBspTree::BestCostRatioHeuristics(const RayInfoContainer &rays,
1164	const AxisAlignedBox3 &box,
1165	const int pvsSize,
1166	const int &axis,
1167	float &position)
1168	{
1169	SortSplitCandidates(rays, axis);
1170
1171	// go through the lists, count the number of objects left and right
1172	// and evaluate the following cost funcion:
1173	// C = ct_div_ci + (qlrl + qrrr)/queries
1174
1175	int pvsl = 0, pvsr = pvsSize;
1176
1177	int pvsBack = pvsl;
1178	int pvsFront = pvsr;
1179
1180	float minBox = box.Min(axis);
1181	float maxBox = box.Max(axis);
1182	float sizeBox = maxBox - minBox;
1183
1184	float minBand = minBox + 0.1f * (maxBox - minBox);
1185	float maxBand = minBox + 0.9f * (maxBox - minBox);
1186
1187	float sum = (float)pvsSize * sizeBox;
1188	float minSum = 1e20f;
1189
1190	Intersectable::NewMail();
1191
1192	RayInfoContainer::const_iterator ri, ri_end = rays.end();
1193
1194	// set all object as belonging to the front pvs
1195	for(ri = rays.begin(); ri != ri_end; ++ ri)
1196	{
1197	Intersectable oObject = (ri).mRay->mOriginObject;
1198	Intersectable tObject = (ri).mRay->mTerminationObject;
1199
1200	if (oObject)
1201	{
1202	if (!oObject->Mailed())
1203	{
1204	oObject->Mail();
1205	oObject->mCounter = 1;
1206	}
1207	else
1208	{
1209	++ oObject->mCounter;
1210	}
1211	}
1212	if (tObject)
1213	{
1214	if (!tObject->Mailed())
1215	{
1216	tObject->Mail();
1217	tObject->mCounter = 1;
1218	}
1219	else
1220	{
1221	++ tObject->mCounter;
1222	}
1223	}
1224	}
1225
1226	Intersectable::NewMail();
1227
1228	vector<SortableEntry>::const_iterator ci, ci_end = mSplitCandidates->end();
1229
1230	for (ci = mSplitCandidates->begin(); ci < ci_end; ++ ci)
1231	{
1232	VssRay *ray;
1233	ray = (*ci).ray;
1234
1235	Intersectable *oObject = ray->mOriginObject;
1236	Intersectable *tObject = ray->mTerminationObject;
1237
1238
1239	switch ((*ci).type)
1240	{
1241	case SortableEntry::ERayMin:
1242	{
1243	if (oObject && !oObject->Mailed())
1244	{
1245	oObject->Mail();
1246	++ pvsl;
1247	}
1248	if (tObject && !tObject->Mailed())
1249	{
1250	tObject->Mail();
1251	++ pvsl;
1252	}
1253	break;
1254	}
1255	case SortableEntry::ERayMax:
1256	{
1257	if (oObject)
1258	{
1259	if (-- oObject->mCounter == 0)
1260	-- pvsr;
1261	}
1262
1263	if (tObject)
1264	{
1265	if (-- tObject->mCounter == 0)
1266	-- pvsr;
1267	}
1268
1269	break;
1270	}
1271	}
1272
1273	// Note: sufficient to compare size of bounding boxes of front and back side?
1274	if ((ci).value > minBand && (ci).value < maxBand)
1275	{
1276	sum = pvsl * ((ci).value - minBox) + pvsr (maxBox - (*ci).value);
1277
1278	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
1279	// cout<<"cost= "<<sum<<endl;
1280
1281	if (sum < minSum)
1282	{
1283	minSum = sum;
1284	position = (*ci).value;
1285
1286	pvsBack = pvsl;
1287	pvsFront = pvsr;
1288	}
1289	}
1290	}
1291
1292	// -- compute cost
1293
1294	const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1295	const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1296
1297	const float pOverall = sizeBox;
1298
1299	const float pBack = position - minBox;
1300	const float pFront = maxBox - position;
1301
1302	const float penaltyOld = EvalPvsPenalty(pvsSize, lowerPvsLimit, upperPvsLimit);
1303	const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1304	const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1305
1306	const float oldRenderCost = penaltyOld * pOverall;
1307	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1308
1309	float ratio = mPvsFactor * newRenderCost / (oldRenderCost + Limits::Small);
1310
1311	//Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
1312	// <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
1313
1314	return ratio;
1315	}
1316
1317
1318	float VspBspTree::SelectAxisAlignedPlane(Plane3 &plane,
1319	const VspBspTraversalData &tData,
1320	int &axis,
1321	BspNodeGeometry **frontGeom,
1322	BspNodeGeometry **backGeom,
1323	float &pFront,
1324	float &pBack,
1325	const bool useKdSplit)
1326	{
1327	float nPosition[3];
1328	float nCostRatio[3];
1329	float nProbFront[3];
1330	float nProbBack[3];
1331
1332	BspNodeGeometry *nFrontGeom[3];
1333	BspNodeGeometry *nBackGeom[3];
1334
1335	for (int i = 0; i < 3; ++ i)
1336	{
1337	nFrontGeom[i] = NULL;
1338	nBackGeom[i] = NULL;
1339	}
1340
1341	int bestAxis = -1;
1342
1343	// create bounding box of node geometry
1344	AxisAlignedBox3 box;
1345	box.Initialize();
1346
1347	//TODO: for kd split geometry already is box => only take minmax vertices
1348	if (1)
1349	{
1350	PolygonContainer::const_iterator it, it_end = tData.mGeometry->GetPolys().end();
1351
1352	for(it = tData.mGeometry->GetPolys().begin(); it < it_end; ++ it)
1353	box.Include((it));
1354	}
1355	else
1356	{
1357	RayInfoContainer::const_iterator ri, ri_end = tData.mRays->end();
1358
1359	for(ri = tData.mRays->begin(); ri < ri_end; ++ ri)
1360	box.Include((*ri).ExtrapTermination());
1361	}
1362
1363	int sAxis = 0;
1364
1365	bool useSpecialAxis = mOnlyDrivingAxis \|\| mUseRandomAxis \|\| mSimulateOctree;
1366
1367	// use some kind of specialised fixed axis
1368	if (mOnlyDrivingAxis)
1369	sAxis = box.Size().DrivingAxis();
1370	else if (mUseRandomAxis)
1371	sAxis = Random(3);
1372	else if (mSimulateOctree)
1373	sAxis = tData.mAxis;
1374
1375	//Debug << "use special axis: " << useSpecialAxis << endl;
1376	//Debug << "axis: " << sAxis << " drivingaxis: " << box.Size().DrivingAxis();
1377
1378	for (axis = 0; axis < 3; ++ axis)
1379	{
1380	if (!useSpecialAxis \|\| (axis == sAxis))
1381	{
1382	if (!mUseCostHeuristics)
1383	{
1384	nFrontGeom[axis] = new BspNodeGeometry();
1385	nBackGeom[axis] = new BspNodeGeometry();
1386
1387	nPosition[axis] = (box.Min()[axis] + box.Max()[axis]) * 0.5f;
1388	Vector3 normal(0,0,0); normal[axis] = 1.0f;
1389
1390	// allows faster split because we have axis aligned kd tree boxes
1391	if (0 && useKdSplit)
1392	{
1393	nCostRatio[axis] = EvalAxisAlignedSplitCost(tData,
1394	box,
1395	axis,
1396	nPosition[axis],
1397	nProbFront[axis],
1398	nProbBack[axis]);
1399
1400	Vector3 pos;
1401
1402	pos = box.Max(); pos[axis] = nPosition[axis];
1403	AxisAlignedBox3 bBox(box.Min(), pos);
1404	// TODO
1405	#if 0
1406	bBox.ExtractPolys(nBackGeom[axis]->GetPolys());
1407	#endif
1408	pos = box.Min(); pos[axis] = nPosition[axis];
1409	AxisAlignedBox3 fBox(pos, box.Max());
1410	// TODO
1411	#if 0
1412	fBox.ExtractPolys(nFrontGeom[axis]->GetPolys());
1413	#endif
1414	}
1415	else
1416	{
1417	nCostRatio[axis] =
1418	EvalSplitPlaneCost(Plane3(normal, nPosition[axis]),
1419	tData, nFrontGeom[axis], nBackGeom[axis],
1420	nProbFront[axis], nProbBack[axis]);
1421	}
1422	}
1423	else
1424	{
1425	nCostRatio[axis] =
1426	BestCostRatioHeuristics(*tData.mRays,
1427	box,
1428	tData.mPvs,
1429	axis,
1430	nPosition[axis]);
1431	}
1432
1433	if (bestAxis == -1)
1434	{
1435	bestAxis = axis;
1436	}
1437	else if (nCostRatio[axis] < nCostRatio[bestAxis])
1438	{
1439	bestAxis = axis;
1440	}
1441
1442	}
1443	}
1444
1445	//-- assign values
1446	axis = bestAxis;
1447	pFront = nProbFront[bestAxis];
1448	pBack = nProbBack[bestAxis];
1449
1450	// assign best split nodes geometry
1451	*frontGeom = nFrontGeom[bestAxis];
1452	*backGeom = nBackGeom[bestAxis];
1453
1454	// and delete other geometry
1455	DEL_PTR(nFrontGeom[(bestAxis + 1) % 3]);
1456	DEL_PTR(nBackGeom[(bestAxis + 2) % 3]);
1457
1458	//-- split plane
1459	Vector3 normal(0,0,0); normal[bestAxis] = 1;
1460	plane = Plane3(normal, nPosition[bestAxis]);
1461
1462	return nCostRatio[bestAxis];
1463	}
1464
1465
1466	bool VspBspTree::SelectPlane(Plane3 &bestPlane,
1467	BspLeaf *leaf,
1468	VspBspTraversalData &data,
1469	VspBspTraversalData &frontData,
1470	VspBspTraversalData &backData,
1471	int &splitAxis)
1472	{
1473	// simplest strategy: just take next polygon
1474	if (mSplitPlaneStrategy & RANDOM_POLYGON)
1475	{
1476	if (!data.mPolygons->empty())
1477	{
1478	const int randIdx =
1479	(int)RandomValue(0, (Real)((int)data.mPolygons->size() - 1));
1480	Polygon3 nextPoly = (data.mPolygons)[randIdx];
1481
1482	bestPlane = nextPoly->GetSupportingPlane();
1483	return true;
1484	}
1485	}
1486
1487	//-- use heuristics to find appropriate plane
1488
1489	// intermediate plane
1490	Plane3 plane;
1491	float lowestCost = MAX_FLOAT;
1492
1493	// decides if the first few splits should be only axisAligned
1494	const bool onlyAxisAligned =
1495	(mSplitPlaneStrategy & AXIS_ALIGNED) &&
1496	((int)data.mRays->size() > mTermMinRaysForAxisAligned) &&
1497	((int)data.GetAvgRayContribution() < mTermMaxRayContriForAxisAligned);
1498
1499	const int limit = onlyAxisAligned ? 0 :
1500	Min((int)data.mPolygons->size(), mMaxPolyCandidates);
1501
1502	float candidateCost;
1503
1504	int maxIdx = (int)data.mPolygons->size();
1505
1506	for (int i = 0; i < limit; ++ i)
1507	{
1508	// the already taken candidates are stored behind maxIdx
1509	// => assure that no index is taken twice
1510	const int candidateIdx = (int)RandomValue(0, (Real)(-- maxIdx));
1511	Polygon3 poly = (data.mPolygons)[candidateIdx];
1512
1513	// swap candidate to the end to avoid testing same plane
1514	std::swap((data.mPolygons)[maxIdx], (data.mPolygons)[candidateIdx]);
1515	//Polygon3 poly = (data.mPolygons)[(int)RandomValue(0, (int)polys.size() - 1)];
1516
1517	// evaluate current candidate
1518	BspNodeGeometry fGeom, bGeom;
1519	float fArea, bArea;
1520	plane = poly->GetSupportingPlane();
1521	candidateCost = EvalSplitPlaneCost(plane, data, fGeom, bGeom, fArea, bArea);
1522
1523	if (candidateCost < lowestCost)
1524	{
1525	bestPlane = plane;
1526	lowestCost = candidateCost;
1527	}
1528	}
1529
1530	//-- evaluate axis aligned splits
1531	int axis;
1532	BspNodeGeometry fGeom, bGeom;
1533	float pFront, pBack;
1534
1535	candidateCost = 99999999.0f;
1536
1537	// option: axis aligned split only if no polygon available
1538	if (!mUsePolygonSplitIfAvailable \|\| data.mPolygons->empty())
1539	{
1540	candidateCost = SelectAxisAlignedPlane(plane,
1541	data,
1542	axis,
1543	&fGeom,
1544	&bGeom,
1545	pFront,
1546	pBack,
1547	data.mIsKdNode);
1548	}
1549
1550	splitAxis = 3;
1551
1552	if (candidateCost < lowestCost)
1553	{
1554	bestPlane = plane;
1555	lowestCost = candidateCost;
1556	splitAxis = axis;
1557
1558	// assign already computed values
1559	// we can do this because we always save the
1560	// computed values from the axis aligned splits
1561
1562	if (fGeom && bGeom)
1563	{
1564	frontData.mGeometry = fGeom;
1565	backData.mGeometry = bGeom;
1566
1567	frontData.mProbability = pFront;
1568	backData.mProbability = pBack;
1569	}
1570	}
1571	else
1572	{
1573	DEL_PTR(fGeom);
1574	DEL_PTR(bGeom);
1575	}
1576
1577
1578	// if (lowestCost > 10)
1579	// Debug << "warning!! lowest cost: " << lowestCost << endl;
1580
1581	#ifdef _DEBUG
1582	Debug << "plane lowest cost: " << lowestCost << endl;
1583	#endif
1584
1585	if (lowestCost > mTermMaxCostRatio)
1586	{
1587	return false;
1588	}
1589
1590	return true;
1591	}
1592
1593
1594	Plane3 VspBspTree::ChooseCandidatePlane(const RayInfoContainer &rays) const
1595	{
1596	const int candidateIdx = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1597
1598	const Vector3 minPt = rays[candidateIdx].ExtrapOrigin();
1599	const Vector3 maxPt = rays[candidateIdx].ExtrapTermination();
1600
1601	const Vector3 pt = (maxPt + minPt) * 0.5;
1602	const Vector3 normal = Normalize(rays[candidateIdx].mRay->GetDir());
1603
1604	return Plane3(normal, pt);
1605	}
1606
1607
1608	Plane3 VspBspTree::ChooseCandidatePlane2(const RayInfoContainer &rays) const
1609	{
1610	Vector3 pt[3];
1611
1612	int idx[3];
1613	int cmaxT = 0;
1614	int cminT = 0;
1615	bool chooseMin = false;
1616
1617	for (int j = 0; j < 3; ++ j)
1618	{
1619	idx[j] = (int)RandomValue(0, (Real)((int)rays.size() * 2 - 1));
1620
1621	if (idx[j] >= (int)rays.size())
1622	{
1623	idx[j] -= (int)rays.size();
1624
1625	chooseMin = (cminT < 2);
1626	}
1627	else
1628	chooseMin = (cmaxT < 2);
1629
1630	RayInfo rayInf = rays[idx[j]];
1631	pt[j] = chooseMin ? rayInf.ExtrapOrigin() : rayInf.ExtrapTermination();
1632	}
1633
1634	return Plane3(pt[0], pt[1], pt[2]);
1635	}
1636
1637
1638	Plane3 VspBspTree::ChooseCandidatePlane3(const RayInfoContainer &rays) const
1639	{
1640	Vector3 pt[3];
1641
1642	int idx1 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1643	int idx2 = (int)RandomValue(0, (Real)((int)rays.size() - 1));
1644
1645	// check if rays different
1646	if (idx1 == idx2)
1647	idx2 = (idx2 + 1) % (int)rays.size();
1648
1649	const RayInfo ray1 = rays[idx1];
1650	const RayInfo ray2 = rays[idx2];
1651
1652	// normal vector of the plane parallel to both lines
1653	const Vector3 norm = Normalize(CrossProd(ray1.mRay->GetDir(), ray2.mRay->GetDir()));
1654
1655	// vector from line 1 to line 2
1656	const Vector3 vd = ray2.ExtrapOrigin() - ray1.ExtrapOrigin();
1657
1658	// project vector on normal to get distance
1659	const float dist = DotProd(vd, norm);
1660
1661	// point on plane lies halfway between the two planes
1662	const Vector3 planePt = ray1.ExtrapOrigin() + norm * dist * 0.5;
1663
1664	return Plane3(norm, planePt);
1665	}
1666
1667
1668	inline void VspBspTree::GenerateUniqueIdsForPvs()
1669	{
1670	Intersectable::NewMail(); sBackId = Intersectable::sMailId;
1671	Intersectable::NewMail(); sFrontId = Intersectable::sMailId;
1672	Intersectable::NewMail(); sFrontAndBackId = Intersectable::sMailId;
1673	}
1674
1675
1676	float VspBspTree::EvalRenderCostDecrease(const Plane3 &candidatePlane,
1677	const VspBspTraversalData &data) const
1678	{
1679	int pvsFront = 0;
1680	int pvsBack = 0;
1681	int totalPvs = 0;
1682
1683	// probability that view point lies in back / front node
1684	float pOverall = data.mProbability;
1685	float pFront = 0;
1686	float pBack = 0;
1687
1688
1689	// create unique ids for pvs heuristics
1690	GenerateUniqueIdsForPvs();
1691
1692	for (int i = 0; i < data.mRays->size(); ++ i)
1693	{
1694	RayInfo rayInf = (*data.mRays)[i];
1695
1696	float t;
1697	VssRay *ray = rayInf.mRay;
1698	const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
1699
1700	// find front and back pvs for origing and termination object
1701	AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1702	AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1703	}
1704
1705
1706	BspNodeGeometry geomFront;
1707	BspNodeGeometry geomBack;
1708
1709	// construct child geometry with regard to the candidate split plane
1710	data.mGeometry->SplitGeometry(geomFront,
1711	geomBack,
1712	candidatePlane,
1713	mBox,
1714	//0.0f);
1715	mEpsilon);
1716
1717	if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
1718	{
1719	pFront = geomFront.GetVolume();
1720	pBack = pOverall - pFront;
1721
1722	if ((pFront < 0.0) \|\| (pBack < 0.0))
1723	{
1724	Debug << "vol f " << pFront << " b " << pBack << " p " << pOverall << endl;
1725	Debug << "real vol f " << pFront << " b " << geomBack.GetVolume() << " p " << pOverall << endl;
1726	Debug << "polys f " << geomFront.Size() << " b " << geomBack.Size() << " data " << data.mGeometry->Size() << endl;
1727	}
1728
1729	// clamp because of possible precision issues
1730	if (0)
1731	{
1732	if (pFront < 0) pFront = 0;
1733	if (pBack < 0) pBack = 0;
1734	}
1735	}
1736	else
1737	{
1738	pFront = geomFront.GetArea();
1739	pBack = geomBack.GetArea();
1740	}
1741
1742
1743	// -- pvs rendering heuristics
1744	const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1745	const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1746
1747	// only render cost heuristics or combined with standard deviation
1748	const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
1749	const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1750	const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1751
1752	const float oldRenderCost = pOverall * penaltyOld;
1753	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1754
1755	//Debug << "decrease: " << oldRenderCost - newRenderCost << endl;
1756	return (oldRenderCost - newRenderCost) / mBox.GetVolume();
1757	}
1758
1759
1760	float VspBspTree::EvalSplitPlaneCost(const Plane3 &candidatePlane,
1761	const VspBspTraversalData &data,
1762	BspNodeGeometry &geomFront,
1763	BspNodeGeometry &geomBack,
1764	float &pFront,
1765	float &pBack) const
1766	{
1767	float cost = 0;
1768
1769	int totalPvs = 0;
1770	int pvsFront = 0;
1771	int pvsBack = 0;
1772
1773	// probability that view point lies in back / front node
1774	float pOverall = 0;
1775	pFront = 0;
1776	pBack = 0;
1777
1778
1779	int limit;
1780	bool useRand;
1781
1782	// create unique ids for pvs heuristics
1783	GenerateUniqueIdsForPvs();
1784
1785	// choose test rays randomly if too much
1786	if ((int)data.mRays->size() > mMaxTests)
1787	{
1788	useRand = true;
1789	limit = mMaxTests;
1790	}
1791	else
1792	{
1793	useRand = false;
1794	limit = (int)data.mRays->size();
1795	}
1796
1797	for (int i = 0; i < limit; ++ i)
1798	{
1799	const int testIdx = useRand ?
1800	(int)RandomValue(0, (Real)((int)data.mRays->size() - 1)) : i;
1801	RayInfo rayInf = (*data.mRays)[testIdx];
1802
1803	float t;
1804	VssRay *ray = rayInf.mRay;
1805	const int cf = rayInf.ComputeRayIntersection(candidatePlane, t);
1806
1807	// find front and back pvs for origing and termination object
1808	AddObjToPvs(ray->mTerminationObject, cf, pvsFront, pvsBack, totalPvs);
1809	AddObjToPvs(ray->mOriginObject, cf, pvsFront, pvsBack, totalPvs);
1810	}
1811
1812	// construct child geometry with regard to the candidate split plane
1813	bool splitSuccessFull = data.mGeometry->SplitGeometry(geomFront,
1814	geomBack,
1815	candidatePlane,
1816	mBox,
1817	//0.0f);
1818	mEpsilon);
1819
1820	pOverall = data.mProbability;
1821
1822	if (!mUseAreaForPvs) // use front and back cell areas to approximate volume
1823	{
1824	pFront = geomFront.GetVolume();
1825	pBack = pOverall - pFront;
1826
1827	// clamp because of possible precision issues
1828	if (1 &&
1829	(!splitSuccessFull \|\| (pFront <= 0) \|\| (pBack <= 0) \|\|
1830	!geomFront.Valid() \|\| !geomBack.Valid()))
1831	{
1832	Debug << "error f: " << pFront << " b: " << pBack << endl;
1833	return 999;
1834	}
1835	}
1836	else
1837	{
1838	pFront = geomFront.GetArea();
1839	pBack = geomBack.GetArea();
1840	}
1841
1842
1843	// -- pvs rendering heuristics
1844	const int lowerPvsLimit = mViewCellsManager->GetMinPvsSize();
1845	const int upperPvsLimit = mViewCellsManager->GetMaxPvsSize();
1846
1847	// only render cost heuristics or combined with standard deviation
1848	const float penaltyOld = EvalPvsPenalty(totalPvs, lowerPvsLimit, upperPvsLimit);
1849	const float penaltyFront = EvalPvsPenalty(pvsFront, lowerPvsLimit, upperPvsLimit);
1850	const float penaltyBack = EvalPvsPenalty(pvsBack, lowerPvsLimit, upperPvsLimit);
1851
1852	const float oldRenderCost = pOverall * penaltyOld;
1853	const float newRenderCost = penaltyFront * pFront + penaltyBack * pBack;
1854
1855	float oldCost, newCost;
1856
1857	// only render cost
1858	if (1)
1859	{
1860	oldCost = oldRenderCost;
1861	newCost = newRenderCost;
1862	}
1863	else // also considering standard deviation
1864	{
1865	// standard deviation is difference of back and front pvs
1866	const float expectedCost = 0.5f * (penaltyFront + penaltyBack);
1867
1868	const float newDeviation = 0.5f *
1869	fabs(penaltyFront - expectedCost) + fabs(penaltyBack - expectedCost);
1870
1871	const float oldDeviation = penaltyOld;
1872
1873	newCost = mRenderCostWeight * newRenderCost + (1.0f - mRenderCostWeight) * newDeviation;
1874	oldCost = mRenderCostWeight * oldRenderCost + (1.0f - mRenderCostWeight) * oldDeviation;
1875	}
1876
1877	cost += mPvsFactor * newCost / (oldCost + Limits::Small);
1878
1879
1880	#ifdef _DEBUG
1881	Debug << "totalpvs: " << data.mPvs << " ptotal: " << pOverall
1882	<< " frontpvs: " << pvsFront << " pFront: " << pFront
1883	<< " backpvs: " << pvsBack << " pBack: " << pBack << endl << endl;
1884	Debug << "cost: " << cost << endl;
1885	#endif
1886
1887	return cost;
1888	}
1889
1890
1891	float VspBspTree::EvalAxisAlignedSplitCost(const VspBspTraversalData &data,
1892	const AxisAlignedBox3 &box,
1893	const int axis,
1894	const float &position,
1895	float &pFront,
1896	float &pBack) const
1897	{
1898	int pvsTotal = 0;
1899	int pvsFront = 0;
1900	int pvsBack = 0;
1901
1902	// create unique ids for pvs heuristics
1903	GenerateUniqueIdsForPvs();
1904
1905	const int pvsSize = data.mPvs;
1906
1907	RayInfoContainer::const_iterator rit, rit_end = data.mRays->end();
1908
1909	// this is the main ray classification loop!
1910	for(rit = data.mRays->begin(); rit != rit_end; ++ rit)
1911	{
1912	//if (!(*rit).mRay->IsActive()) continue;
1913
1914	// determine the side of this ray with respect to the plane
1915	float t;
1916	const int side = (*rit).ComputeRayIntersection(axis, position, t);
1917
1918	AddObjToPvs((*rit).mRay->mTerminationObject, side, pvsFront, pvsBack, pvsTotal);
1919	AddObjToPvs((*rit).mRay->mOriginObject, side, pvsFront, pvsBack, pvsTotal);
1920	}
1921
1922	//-- pvs heuristics
1923	float pOverall;
1924
1925	//-- compute heurstics
1926	// we take simplified computation for mid split
1927
1928	pOverall = data.mProbability;
1929
1930	if (!mUseAreaForPvs)
1931	{ // volume
1932	pBack = pFront = pOverall * 0.5f;
1933
1934	#if 0
1935	// box length substitute for probability
1936	const float minBox = box.Min(axis);
1937	const float maxBox = box.Max(axis);
1938
1939	pBack = position - minBox;
1940	pFront = maxBox - position;
1941	pOverall = maxBox - minBox;
1942	#endif
1943	}
1944	else //-- area substitute for probability
1945	{
1946	const int axis2 = (axis + 1) % 3;
1947	const int axis3 = (axis + 2) % 3;
1948
1949	const float faceArea =
1950	(box.Max(axis2) - box.Min(axis2)) *
1951	(box.Max(axis3) - box.Min(axis3));
1952
1953	pBack = pFront = pOverall * 0.5f + faceArea;
1954	}
1955
1956	#ifdef _DEBUG
1957	Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
1958	Debug << pFront << " " << pBack << " " << pOverall << endl;
1959	#endif
1960
1961
1962	const float newCost = pvsBack * pBack + pvsFront * pFront;
1963	const float oldCost = (float)pvsSize * pOverall + Limits::Small;
1964
1965	return (mCtDivCi + newCost) / oldCost;
1966	}
1967
1968
1969	void VspBspTree::AddObjToPvs(Intersectable *obj,
1970	const int cf,
1971	int &frontPvs,
1972	int &backPvs,
1973	int &totalPvs) const
1974	{
1975	if (!obj)
1976	return;
1977	// new object
1978	if ((obj->mMailbox != sFrontId) &&
1979	(obj->mMailbox != sBackId) &&
1980	(obj->mMailbox != sFrontAndBackId))
1981	{
1982	++ totalPvs;
1983	}
1984
1985	// TODO: does this really belong to no pvs?
1986	//if (cf == Ray::COINCIDENT) return;
1987
1988	// object belongs to both PVS
1989	if (cf >= 0)
1990	{
1991	if ((obj->mMailbox != sFrontId) &&
1992	(obj->mMailbox != sFrontAndBackId))
1993	{
1994	++ frontPvs;
1995
1996	if (obj->mMailbox == sBackId)
1997	obj->mMailbox = sFrontAndBackId;
1998	else
1999	obj->mMailbox = sFrontId;
2000	}
2001	}
2002
2003	if (cf <= 0)
2004	{
2005	if ((obj->mMailbox != sBackId) &&
2006	(obj->mMailbox != sFrontAndBackId))
2007	{
2008	++ backPvs;
2009
2010	if (obj->mMailbox == sFrontId)
2011	obj->mMailbox = sFrontAndBackId;
2012	else
2013	obj->mMailbox = sBackId;
2014	}
2015	}
2016	}
2017
2018
2019	void VspBspTree::CollectLeaves(vector<BspLeaf *> &leaves,
2020	const bool onlyUnmailed,
2021	const int maxPvsSize) const
2022	{
2023	stack<BspNode *> nodeStack;
2024	nodeStack.push(mRoot);
2025
2026	while (!nodeStack.empty())
2027	{
2028	BspNode *node = nodeStack.top();
2029	nodeStack.pop();
2030
2031	if (node->IsLeaf())
2032	{
2033	// test if this leaf is in valid view space
2034	BspLeaf leaf = dynamic_cast<BspLeaf >(node);
2035	if (leaf->TreeValid() &&
2036	(!onlyUnmailed \|\| !leaf->Mailed()) &&
2037	((maxPvsSize < 0) \|\| (leaf->GetViewCell()->GetPvs().GetSize() <= maxPvsSize)))
2038	{
2039	leaves.push_back(leaf);
2040	}
2041	}
2042	else
2043	{
2044	BspInterior interior = dynamic_cast<BspInterior >(node);
2045
2046	nodeStack.push(interior->GetBack());
2047	nodeStack.push(interior->GetFront());
2048	}
2049	}
2050	}
2051
2052
2053	AxisAlignedBox3 VspBspTree::GetBoundingBox() const
2054	{
2055	return mBox;
2056	}
2057
2058
2059	BspNode *VspBspTree::GetRoot() const
2060	{
2061	return mRoot;
2062	}
2063
2064
2065	void VspBspTree::EvaluateLeafStats(const VspBspTraversalData &data)
2066	{
2067	// the node became a leaf -> evaluate stats for leafs
2068	BspLeaf leaf = dynamic_cast<BspLeaf >(data.mNode);
2069
2070	// store maximal and minimal depth
2071	if (data.mDepth > mBspStats.maxDepth)
2072	mBspStats.maxDepth = data.mDepth;
2073
2074	if (data.mPvs > mBspStats.maxPvs)
2075	mBspStats.maxPvs = data.mPvs;
2076
2077	mBspStats.pvs += data.mPvs;
2078
2079	if (data.mDepth < mBspStats.minDepth)
2080	mBspStats.minDepth = data.mDepth;
2081
2082	if (data.mDepth >= mTermMaxDepth)
2083	++ mBspStats.maxDepthNodes;
2084
2085	// accumulate rays to compute rays / leaf
2086	mBspStats.accumRays += (int)data.mRays->size();
2087
2088	if (data.mPvs < mTermMinPvs)
2089	++ mBspStats.minPvsNodes;
2090
2091	if ((int)data.mRays->size() < mTermMinRays)
2092	++ mBspStats.minRaysNodes;
2093
2094	if (data.GetAvgRayContribution() > mTermMaxRayContribution)
2095	++ mBspStats.maxRayContribNodes;
2096
2097	if (data.mProbability <= mTermMinProbability)
2098	++ mBspStats.minProbabilityNodes;
2099
2100	// accumulate depth to compute average depth
2101	mBspStats.accumDepth += data.mDepth;
2102
2103	++ mCreatedViewCells;
2104
2105	#ifdef _DEBUG
2106	Debug << "BSP stats: "
2107	<< "Depth: " << data.mDepth << " (max: " << mTermMaxDepth << "), "
2108	<< "PVS: " << data.mPvs << " (min: " << mTermMinPvs << "), "
2109	// << "Area: " << data.mProbability << " (min: " << mTermMinProbability << "), "
2110	<< "#rays: " << (int)data.mRays->size() << " (max: " << mTermMinRays << "), "
2111	<< "#pvs: " << leaf->GetViewCell()->GetPvs().GetSize() << "=, "
2112	<< "#avg ray contrib (pvs): " << (float)data.mPvs / (float)data.mRays->size() << endl;
2113	#endif
2114	}
2115
2116
2117	int VspBspTree::CastRay(Ray &ray)
2118	{
2119	int hits = 0;
2120
2121	stack<BspRayTraversalData> tQueue;
2122
2123	float maxt, mint;
2124
2125	if (!mBox.GetRaySegment(ray, mint, maxt))
2126	return 0;
2127
2128	Intersectable::NewMail();
2129	ViewCell::NewMail();
2130	Vector3 entp = ray.Extrap(mint);
2131	Vector3 extp = ray.Extrap(maxt);
2132
2133	BspNode *node = mRoot;
2134	BspNode *farChild = NULL;
2135
2136	while (1)
2137	{
2138	if (!node->IsLeaf())
2139	{
2140	BspInterior in = dynamic_cast<BspInterior >(node);
2141
2142	Plane3 splitPlane = in->GetPlane();
2143	const int entSide = splitPlane.Side(entp);
2144	const int extSide = splitPlane.Side(extp);
2145
2146	if (entSide < 0)
2147	{
2148	node = in->GetBack();
2149
2150	if(extSide <= 0) // plane does not split ray => no far child
2151	continue;
2152
2153	farChild = in->GetFront(); // plane splits ray
2154
2155	} else if (entSide > 0)
2156	{
2157	node = in->GetFront();
2158
2159	if (extSide >= 0) // plane does not split ray => no far child
2160	continue;
2161
2162	farChild = in->GetBack(); // plane splits ray
2163	}
2164	else // ray and plane are coincident
2165	{
2166	// WHAT TO DO IN THIS CASE ?
2167	//break;
2168	node = in->GetFront();
2169	continue;
2170	}
2171
2172	// push data for far child
2173	tQueue.push(BspRayTraversalData(farChild, extp, maxt));
2174
2175	// find intersection of ray segment with plane
2176	float t;
2177	extp = splitPlane.FindIntersection(ray.GetLoc(), extp, &t);
2178	maxt *= t;
2179
2180	} else // reached leaf => intersection with view cell
2181	{
2182	BspLeaf leaf = dynamic_cast<BspLeaf >(node);
2183
2184	if (!leaf->GetViewCell()->Mailed())
2185	{
2186	//ray.bspIntersections.push_back(Ray::VspBspIntersection(maxt, leaf));
2187	leaf->GetViewCell()->Mail();
2188	++ hits;
2189	}
2190
2191	//-- fetch the next far child from the stack
2192	if (tQueue.empty())
2193	break;
2194
2195	entp = extp;
2196	mint = maxt; // NOTE: need this?
2197
2198	if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
2199	break;
2200
2201	BspRayTraversalData &s = tQueue.top();
2202
2203	node = s.mNode;
2204	extp = s.mExitPoint;
2205	maxt = s.mMaxT;
2206
2207	tQueue.pop();
2208	}
2209	}
2210
2211	return hits;
2212	}
2213
2214
2215	void VspBspTree::CollectViewCells(ViewCellContainer &viewCells, bool onlyValid) const
2216	{
2217	ViewCell::NewMail();
2218
2219	CollectViewCells(mRoot, onlyValid, viewCells, true);
2220	}
2221
2222
2223	void VspBspTree::CollapseViewCells()
2224	{
2225	// TODO
2226	#if VC_HISTORY
2227	stack<BspNode *> nodeStack;
2228
2229	if (!mRoot)
2230	return;
2231
2232	nodeStack.push(mRoot);
2233
2234	while (!nodeStack.empty())
2235	{
2236	BspNode *node = nodeStack.top();
2237	nodeStack.pop();
2238
2239	if (node->IsLeaf())
2240	{
2241	BspViewCell viewCell = dynamic_cast<BspLeaf >(node)->GetViewCell();
2242
2243	if (!viewCell->GetValid())
2244	{
2245	BspViewCell viewCell = dynamic_cast<BspLeaf >(node)->GetViewCell();
2246
2247	ViewCellContainer leaves;
2248	mViewCellsTree->CollectLeaves(viewCell, leaves);
2249
2250	ViewCellContainer::const_iterator it, it_end = leaves.end();
2251
2252	for (it = leaves.begin(); it != it_end; ++ it)
2253	{
2254	BspLeaf l = dynamic_cast<BspViewCell >(*it)->mLeaf;
2255	l->SetViewCell(GetOrCreateOutOfBoundsCell());
2256	++ mBspStats.invalidLeaves;
2257	}
2258
2259	// add to unbounded view cell
2260	GetOrCreateOutOfBoundsCell()->GetPvs().AddPvs(viewCell->GetPvs());
2261	DEL_PTR(viewCell);
2262	}
2263	}
2264	else
2265	{
2266	BspInterior interior = dynamic_cast<BspInterior >(node);
2267
2268	nodeStack.push(interior->GetFront());
2269	nodeStack.push(interior->GetBack());
2270	}
2271	}
2272
2273	Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
2274	#endif
2275	}
2276
2277
2278	void VspBspTree::CollectRays(VssRayContainer &rays)
2279	{
2280	vector<BspLeaf *> leaves;
2281
2282	vector<BspLeaf *>::const_iterator lit, lit_end = leaves.end();
2283
2284	for (lit = leaves.begin(); lit != lit_end; ++ lit)
2285	{
2286	BspLeaf leaf = lit;
2287	VssRayContainer::const_iterator rit, rit_end = leaf->mVssRays.end();
2288
2289	for (rit = leaf->mVssRays.begin(); rit != rit_end; ++ rit)
2290	rays.push_back(*rit);
2291	}
2292	}
2293
2294
2295	void VspBspTree::ValidateTree()
2296	{
2297	stack<BspNode *> nodeStack;
2298
2299	if (!mRoot)
2300	return;
2301
2302	nodeStack.push(mRoot);
2303
2304	mBspStats.invalidLeaves = 0;
2305	while (!nodeStack.empty())
2306	{
2307	BspNode *node = nodeStack.top();
2308	nodeStack.pop();
2309
2310	if (node->IsLeaf())
2311	{
2312	BspLeaf leaf = dynamic_cast<BspLeaf >(node);
2313
2314	if (!leaf->GetViewCell()->GetValid())
2315	++ mBspStats.invalidLeaves;
2316
2317	// validity flags don't match => repair
2318	if (leaf->GetViewCell()->GetValid() != leaf->TreeValid())
2319	{
2320	leaf->SetTreeValid(leaf->GetViewCell()->GetValid());
2321	PropagateUpValidity(leaf);
2322	}
2323	}
2324	else
2325	{
2326	BspInterior interior = dynamic_cast<BspInterior >(node);
2327
2328	nodeStack.push(interior->GetFront());
2329	nodeStack.push(interior->GetBack());
2330	}
2331	}
2332
2333	Debug << "invalid leaves: " << mBspStats.invalidLeaves << endl;
2334	}
2335
2336
2337
2338	void VspBspTree::CollectViewCells(BspNode *root,
2339	bool onlyValid,
2340	ViewCellContainer &viewCells,
2341	bool onlyUnmailed) const
2342	{
2343	stack<BspNode *> nodeStack;
2344
2345	if (!root)
2346	return;
2347
2348	nodeStack.push(root);
2349
2350	while (!nodeStack.empty())
2351	{
2352	BspNode *node = nodeStack.top();
2353	nodeStack.pop();
2354
2355	if (node->IsLeaf())
2356	{
2357	if (!onlyValid \|\| node->TreeValid())
2358	{
2359	ViewCell *viewCell =
2360	mViewCellsTree->GetActiveViewCell(dynamic_cast<BspLeaf *>(node)->GetViewCell());
2361
2362	if (!onlyUnmailed \|\| !viewCell->Mailed())
2363	{
2364	viewCell->Mail();
2365	viewCells.push_back(viewCell);
2366	}
2367	}
2368	}
2369	else
2370	{
2371	BspInterior interior = dynamic_cast<BspInterior >(node);
2372
2373	nodeStack.push(interior->GetFront());
2374	nodeStack.push(interior->GetBack());
2375	}
2376	}
2377
2378	}
2379
2380
2381	void VspBspTree::PreprocessPolygons(PolygonContainer &polys)
2382	{
2383	// preprocess: throw out polygons coincident to the view space box (not needed)
2384	PolygonContainer boxPolys;
2385	mBox.ExtractPolys(boxPolys);
2386	vector<Plane3> boxPlanes;
2387
2388	PolygonContainer::iterator pit, pit_end = boxPolys.end();
2389
2390	// extract planes of box
2391	// TODO: can be done more elegantly than first extracting polygons
2392	// and take their planes
2393	for (pit = boxPolys.begin(); pit != pit_end; ++ pit)
2394	{
2395	boxPlanes.push_back((*pit)->GetSupportingPlane());
2396	}
2397
2398	pit_end = polys.end();
2399
2400	for (pit = polys.begin(); pit != pit_end; ++ pit)
2401	{
2402	vector<Plane3>::const_iterator bit, bit_end = boxPlanes.end();
2403
2404	for (bit = boxPlanes.begin(); (bit != bit_end) && (*pit); ++ bit)
2405	{
2406	const int cf = (pit)->ClassifyPlane(bit, mEpsilon);
2407
2408	if (cf == Polygon3::COINCIDENT)
2409	{
2410	DEL_PTR(*pit);
2411	//Debug << "coincident!!" << endl;
2412	}
2413	}
2414	}
2415
2416	// remove deleted entries
2417	for (int i = 0; i < (int)polys.size(); ++ i)
2418	{
2419	while (!polys[i] && (i < (int)polys.size()))
2420	{
2421	swap(polys[i], polys.back());
2422	polys.pop_back();
2423	}
2424	}
2425	}
2426
2427
2428	float VspBspTree::AccumulatedRayLength(const RayInfoContainer &rays) const
2429	{
2430	float len = 0;
2431
2432	RayInfoContainer::const_iterator it, it_end = rays.end();
2433
2434	for (it = rays.begin(); it != it_end; ++ it)
2435	len += (*it).SegmentLength();
2436
2437	return len;
2438	}
2439
2440
2441	int VspBspTree::SplitRays(const Plane3 &plane,
2442	RayInfoContainer &rays,
2443	RayInfoContainer &frontRays,
2444	RayInfoContainer &backRays) const
2445	{
2446	int splits = 0;
2447
2448	RayInfoContainer::const_iterator it, it_end = rays.end();
2449
2450	for (it = rays.begin(); it != it_end; ++ it)
2451	{
2452	RayInfo bRay = *it;
2453
2454	VssRay *ray = bRay.mRay;
2455	float t;
2456
2457	// get classification and receive new t
2458	const int cf = bRay.ComputeRayIntersection(plane, t);
2459
2460	switch (cf)
2461	{
2462	case -1:
2463	backRays.push_back(bRay);
2464	break;
2465	case 1:
2466	frontRays.push_back(bRay);
2467	break;
2468	case 0:
2469	{
2470	//-- split ray
2471	// test if start point behind or in front of plane
2472	const int side = plane.Side(bRay.ExtrapOrigin());
2473
2474	if (side <= 0)
2475	{
2476	backRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2477	frontRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2478	}
2479	else
2480	{
2481	frontRays.push_back(RayInfo(ray, bRay.GetMinT(), t));
2482	backRays.push_back(RayInfo(ray, t, bRay.GetMaxT()));
2483	}
2484	}
2485	break;
2486	default:
2487	Debug << "Should not come here" << endl;
2488	break;
2489	}
2490	}
2491
2492	return splits;
2493	}
2494
2495
2496	void VspBspTree::ExtractHalfSpaces(BspNode *n, vector<Plane3> &halfSpaces) const
2497	{
2498	BspNode *lastNode;
2499
2500	do
2501	{
2502	lastNode = n;
2503
2504	// want to get planes defining geometry of this node => don't take
2505	// split plane of node itself
2506	n = n->GetParent();
2507
2508	if (n)
2509	{
2510	BspInterior interior = dynamic_cast<BspInterior >(n);
2511	Plane3 halfSpace = dynamic_cast<BspInterior *>(interior)->GetPlane();
2512
2513	if (interior->GetBack() != lastNode)
2514	halfSpace.ReverseOrientation();
2515
2516	halfSpaces.push_back(halfSpace);
2517	}
2518	}
2519	while (n);
2520	}
2521
2522
2523	void VspBspTree::ConstructGeometry(BspNode *n,
2524	BspNodeGeometry &geom) const
2525	{
2526	vector<Plane3> halfSpaces;
2527	ExtractHalfSpaces(n, halfSpaces);
2528
2529	PolygonContainer candidatePolys;
2530	vector<Plane3> candidatePlanes;
2531
2532	// bounded planes are added to the polygons
2533	for (int i = 0; i < (int)halfSpaces.size(); ++ i)
2534	{
2535	Polygon3 *p = GetBoundingBox().CrossSection(halfSpaces[i]);
2536
2537	if (p->Valid(mEpsilon))
2538	{
2539	candidatePolys.push_back(p);
2540	candidatePlanes.push_back(halfSpaces[i]);
2541	}
2542	}
2543
2544	// add faces of bounding box (also could be faces of the cell)
2545	for (int i = 0; i < 6; ++ i)
2546	{
2547	VertexContainer vertices;
2548
2549	for (int j = 0; j < 4; ++ j)
2550	vertices.push_back(mBox.GetFace(i).mVertices[j]);
2551
2552	Polygon3 *poly = new Polygon3(vertices);
2553
2554	candidatePolys.push_back(poly);
2555	candidatePlanes.push_back(poly->GetSupportingPlane());
2556	}
2557
2558	for (int i = 0; i < (int)candidatePolys.size(); ++ i)
2559	{
2560	// polygon is split by all other planes
2561	for (int j = 0; (j < (int)halfSpaces.size()) && candidatePolys[i]; ++ j)
2562	{
2563	if (i == j) // polygon and plane are coincident
2564	continue;
2565
2566	VertexContainer splitPts;
2567	Polygon3 frontPoly, backPoly;
2568
2569	const int cf =
2570	candidatePolys[i]->ClassifyPlane(halfSpaces[j],
2571	mEpsilon);
2572
2573	switch (cf)
2574	{
2575	case Polygon3::SPLIT:
2576	frontPoly = new Polygon3();
2577	backPoly = new Polygon3();
2578
2579	candidatePolys[i]->Split(halfSpaces[j],
2580	*frontPoly,
2581	*backPoly,
2582	mEpsilon);
2583
2584	DEL_PTR(candidatePolys[i]);
2585
2586	if (backPoly->Valid(mEpsilon))
2587	candidatePolys[i] = backPoly;
2588	else
2589	DEL_PTR(backPoly);
2590
2591	// outside, don't need this
2592	DEL_PTR(frontPoly);
2593	break;
2594	// polygon outside of halfspace
2595	case Polygon3::FRONT_SIDE:
2596	DEL_PTR(candidatePolys[i]);
2597	break;
2598	// just take polygon as it is
2599	case Polygon3::BACK_SIDE:
2600	case Polygon3::COINCIDENT:
2601	default:
2602	break;
2603	}
2604	}
2605
2606	if (candidatePolys[i])
2607	{
2608	geom.Add(candidatePolys[i], candidatePlanes[i]);
2609	// geom.mPolys.push_back(candidates[i]);
2610	}
2611	}
2612	}
2613
2614
2615	void VspBspTree::ConstructGeometry(ViewCell *vc,
2616	BspNodeGeometry &vcGeom) const
2617	{
2618	ViewCellContainer leaves;
2619
2620	mViewCellsTree->CollectLeaves(vc, leaves);
2621
2622	ViewCellContainer::const_iterator it, it_end = leaves.end();
2623
2624	for (it = leaves.begin(); it != it_end; ++ it)
2625	{
2626	BspLeaf l = dynamic_cast<BspViewCell >(*it)->mLeaf;
2627
2628	ConstructGeometry(l, vcGeom);
2629	}
2630	}
2631
2632
2633	typedef pair<BspNode , BspNodeGeometry > bspNodePair;
2634
2635
2636	int VspBspTree::FindNeighbors(BspNode n, vector<BspLeaf > &neighbors,
2637	const bool onlyUnmailed) const
2638	{
2639	stack<bspNodePair> nodeStack;
2640
2641	BspNodeGeometry nodeGeom;
2642	ConstructGeometry(n, nodeGeom);
2643
2644	// split planes from the root to this node
2645	// needed to verify that we found neighbor leaf
2646	// TODO: really needed?
2647	vector<Plane3> halfSpaces;
2648	ExtractHalfSpaces(n, halfSpaces);
2649
2650
2651	BspNodeGeometry *rgeom = new BspNodeGeometry();
2652	ConstructGeometry(mRoot, *rgeom);
2653
2654	nodeStack.push(bspNodePair(mRoot, rgeom));
2655
2656	while (!nodeStack.empty())
2657	{
2658	BspNode *node = nodeStack.top().first;
2659	BspNodeGeometry *geom = nodeStack.top().second;
2660
2661	nodeStack.pop();
2662
2663	if (node->IsLeaf())
2664	{
2665	// test if this leaf is in valid view space
2666	if (node->TreeValid() &&
2667	(node != n) &&
2668	(!onlyUnmailed \|\| !node->Mailed()))
2669	{
2670	bool isAdjacent = true;
2671
2672	if (1)
2673	{
2674	// test all planes of current node if still adjacent
2675	for (int i = 0; (i < halfSpaces.size()) && isAdjacent; ++ i)
2676	{
2677	const int cf =
2678	Polygon3::ClassifyPlane(geom->GetPolys(),
2679	halfSpaces[i],
2680	mEpsilon);
2681
2682	if (cf == Polygon3::FRONT_SIDE)
2683	{
2684	isAdjacent = false;
2685	}
2686	}
2687	}
2688	else
2689	{
2690	// TODO: why is this wrong??
2691	// test all planes of current node if still adjacent
2692	for (int i = 0; (i < nodeGeom.Size()) && isAdjacent; ++ i)
2693	{
2694	Polygon3 *poly = nodeGeom.GetPolys()[i];
2695
2696	const int cf =
2697	Polygon3::ClassifyPlane(geom->GetPolys(),
2698	poly->GetSupportingPlane(),
2699	mEpsilon);
2700
2701	if (cf == Polygon3::FRONT_SIDE)
2702	{
2703	isAdjacent = false;
2704	}
2705	}
2706	}
2707	// neighbor was found
2708	if (isAdjacent)
2709	{
2710	neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2711	}
2712	}
2713	}
2714	else
2715	{
2716	BspInterior interior = dynamic_cast<BspInterior >(node);
2717
2718	const int cf = Polygon3::ClassifyPlane(nodeGeom.GetPolys(),
2719	interior->GetPlane(),
2720	mEpsilon);
2721
2722	BspNode *front = interior->GetFront();
2723	BspNode *back = interior->GetBack();
2724
2725	BspNodeGeometry *fGeom = new BspNodeGeometry();
2726	BspNodeGeometry *bGeom = new BspNodeGeometry();
2727
2728	geom->SplitGeometry(*fGeom,
2729	*bGeom,
2730	interior->GetPlane(),
2731	mBox,
2732	//0.0000001f);
2733	mEpsilon);
2734
2735	if (cf == Polygon3::BACK_SIDE)
2736	{
2737	nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
2738	DEL_PTR(fGeom);
2739	}
2740	else
2741	{
2742	if (cf == Polygon3::FRONT_SIDE)
2743	{
2744	nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
2745	DEL_PTR(bGeom);
2746	}
2747	else
2748	{ // random decision
2749	nodeStack.push(bspNodePair(front, fGeom));
2750	nodeStack.push(bspNodePair(back, bGeom));
2751	}
2752	}
2753	}
2754
2755	DEL_PTR(geom);
2756	}
2757
2758	return (int)neighbors.size();
2759	}
2760
2761
2762
2763	int VspBspTree::FindApproximateNeighbors(BspNode n, vector<BspLeaf > &neighbors,
2764	const bool onlyUnmailed) const
2765	{
2766	stack<bspNodePair> nodeStack;
2767
2768	BspNodeGeometry nodeGeom;
2769	ConstructGeometry(n, nodeGeom);
2770
2771	// split planes from the root to this node
2772	// needed to verify that we found neighbor leaf
2773	// TODO: really needed?
2774	vector<Plane3> halfSpaces;
2775	ExtractHalfSpaces(n, halfSpaces);
2776
2777
2778	BspNodeGeometry *rgeom = new BspNodeGeometry();
2779	ConstructGeometry(mRoot, *rgeom);
2780
2781	nodeStack.push(bspNodePair(mRoot, rgeom));
2782
2783	while (!nodeStack.empty())
2784	{
2785	BspNode *node = nodeStack.top().first;
2786	BspNodeGeometry *geom = nodeStack.top().second;
2787
2788	nodeStack.pop();
2789
2790	if (node->IsLeaf())
2791	{
2792	// test if this leaf is in valid view space
2793	if (node->TreeValid() &&
2794	(node != n) &&
2795	(!onlyUnmailed \|\| !node->Mailed()))
2796	{
2797	bool isAdjacent = true;
2798
2799	// neighbor was found
2800	if (isAdjacent)
2801	{
2802	neighbors.push_back(dynamic_cast<BspLeaf *>(node));
2803	}
2804	}
2805	}
2806	else
2807	{
2808	BspInterior interior = dynamic_cast<BspInterior >(node);
2809
2810	const int cf = Polygon3::ClassifyPlane(nodeGeom.GetPolys(),
2811	interior->GetPlane(),
2812	mEpsilon);
2813
2814	BspNode *front = interior->GetFront();
2815	BspNode *back = interior->GetBack();
2816
2817	BspNodeGeometry *fGeom = new BspNodeGeometry();
2818	BspNodeGeometry *bGeom = new BspNodeGeometry();
2819
2820	geom->SplitGeometry(*fGeom,
2821	*bGeom,
2822	interior->GetPlane(),
2823	mBox,
2824	//0.0000001f);
2825	mEpsilon);
2826
2827	if (cf == Polygon3::BACK_SIDE)
2828	{
2829	nodeStack.push(bspNodePair(interior->GetBack(), bGeom));
2830	DEL_PTR(fGeom);
2831	}
2832	else
2833	{
2834	if (cf == Polygon3::FRONT_SIDE)
2835	{
2836	nodeStack.push(bspNodePair(interior->GetFront(), fGeom));
2837	DEL_PTR(bGeom);
2838	}
2839	else
2840	{ // random decision
2841	nodeStack.push(bspNodePair(front, fGeom));
2842	nodeStack.push(bspNodePair(back, bGeom));
2843	}
2844	}
2845	}
2846
2847	DEL_PTR(geom);
2848	}
2849
2850	return (int)neighbors.size();
2851	}
2852
2853
2854
2855	BspLeaf *VspBspTree::GetRandomLeaf(const Plane3 &halfspace)
2856	{
2857	stack<BspNode *> nodeStack;
2858	nodeStack.push(mRoot);
2859
2860	int mask = rand();
2861
2862	while (!nodeStack.empty())
2863	{
2864	BspNode *node = nodeStack.top();
2865	nodeStack.pop();
2866
2867	if (node->IsLeaf())
2868	{
2869	return dynamic_cast<BspLeaf *>(node);
2870	}
2871	else
2872	{
2873	BspInterior interior = dynamic_cast<BspInterior >(node);
2874	BspNode *next;
2875	BspNodeGeometry geom;
2876
2877	// todo: not very efficient: constructs full cell everytime
2878	ConstructGeometry(interior, geom);
2879
2880	const int cf =
2881	Polygon3::ClassifyPlane(geom.GetPolys(), halfspace, mEpsilon);
2882
2883	if (cf == Polygon3::BACK_SIDE)
2884	next = interior->GetFront();
2885	else
2886	if (cf == Polygon3::FRONT_SIDE)
2887	next = interior->GetFront();
2888	else
2889	{
2890	// random decision
2891	if (mask & 1)
2892	next = interior->GetBack();
2893	else
2894	next = interior->GetFront();
2895	mask = mask >> 1;
2896	}
2897
2898	nodeStack.push(next);
2899	}
2900	}
2901
2902	return NULL;
2903	}
2904
2905	BspLeaf *VspBspTree::GetRandomLeaf(const bool onlyUnmailed)
2906	{
2907	stack<BspNode *> nodeStack;
2908
2909	nodeStack.push(mRoot);
2910
2911	int mask = rand();
2912
2913	while (!nodeStack.empty())
2914	{
2915	BspNode *node = nodeStack.top();
2916	nodeStack.pop();
2917
2918	if (node->IsLeaf())
2919	{
2920	if ( (!onlyUnmailed \|\| !node->Mailed()) )
2921	return dynamic_cast<BspLeaf *>(node);
2922	}
2923	else
2924	{
2925	BspInterior interior = dynamic_cast<BspInterior >(node);
2926
2927	// random decision
2928	if (mask & 1)
2929	nodeStack.push(interior->GetBack());
2930	else
2931	nodeStack.push(interior->GetFront());
2932
2933	mask = mask >> 1;
2934	}
2935	}
2936
2937	return NULL;
2938	}
2939
2940	int VspBspTree::ComputePvsSize(const RayInfoContainer &rays) const
2941	{
2942	int pvsSize = 0;
2943
2944	RayInfoContainer::const_iterator rit, rit_end = rays.end();
2945
2946	Intersectable::NewMail();
2947
2948	for (rit = rays.begin(); rit != rays.end(); ++ rit)
2949	{
2950	VssRay ray = (rit).mRay;
2951
2952	if (ray->mOriginObject)
2953	{
2954	if (!ray->mOriginObject->Mailed())
2955	{
2956	ray->mOriginObject->Mail();
2957	++ pvsSize;
2958	}
2959	}
2960	if (ray->mTerminationObject)
2961	{
2962	if (!ray->mTerminationObject->Mailed())
2963	{
2964	ray->mTerminationObject->Mail();
2965	++ pvsSize;
2966	}
2967	}
2968	}
2969
2970	return pvsSize;
2971	}
2972
2973	float VspBspTree::GetEpsilon() const
2974	{
2975	return mEpsilon;
2976	}
2977
2978
2979	int VspBspTree::SplitPolygons(const Plane3 &plane,
2980	PolygonContainer &polys,
2981	PolygonContainer &frontPolys,
2982	PolygonContainer &backPolys,
2983	PolygonContainer &coincident) const
2984	{
2985	int splits = 0;
2986
2987	PolygonContainer::const_iterator it, it_end = polys.end();
2988
2989	for (it = polys.begin(); it != polys.end(); ++ it)
2990	{
2991	Polygon3 poly = it;
2992
2993	// classify polygon
2994	const int cf = poly->ClassifyPlane(plane, mEpsilon);
2995
2996	switch (cf)
2997	{
2998	case Polygon3::COINCIDENT:
2999	coincident.push_back(poly);
3000	break;
3001	case Polygon3::FRONT_SIDE:
3002	frontPolys.push_back(poly);
3003	break;
3004	case Polygon3::BACK_SIDE:
3005	backPolys.push_back(poly);
3006	break;
3007	case Polygon3::SPLIT:
3008	backPolys.push_back(poly);
3009	frontPolys.push_back(poly);
3010	++ splits;
3011	break;
3012	default:
3013	Debug << "SHOULD NEVER COME HERE\n";
3014	break;
3015	}
3016	}
3017
3018	return splits;
3019	}
3020
3021
3022	int VspBspTree::CastLineSegment(const Vector3 &origin,
3023	const Vector3 &termination,
3024	vector<ViewCell *> &viewcells)
3025	{
3026	int hits = 0;
3027	stack<BspRayTraversalData> tQueue;
3028
3029	float mint = 0.0f, maxt = 1.0f;
3030
3031	Intersectable::NewMail();
3032	ViewCell::NewMail();
3033
3034	Vector3 entp = origin;
3035	Vector3 extp = termination;
3036
3037	BspNode *node = mRoot;
3038	BspNode *farChild = NULL;
3039
3040	float t;
3041	//const float thresh = 1 ? 1e-6f : 0.0f;
3042	const float thresh = 0.01f;
3043
3044	while (1)
3045	{
3046	if (!node->IsLeaf())
3047	{
3048	BspInterior in = dynamic_cast<BspInterior >(node);
3049
3050	Plane3 splitPlane = in->GetPlane();
3051
3052	const int entSide = splitPlane.Side(entp, thresh);
3053	const int extSide = splitPlane.Side(extp, thresh);
3054
3055	if (entSide < 0)
3056	{
3057	node = in->GetBack();
3058
3059	// plane does not split ray => no far child
3060	if (extSide <= 0)
3061	continue;
3062
3063	farChild = in->GetFront(); // plane splits ray
3064	}
3065	else if (entSide > 0)
3066	{
3067	node = in->GetFront();
3068
3069	if (extSide >= 0) // plane does not split ray => no far child
3070	continue;
3071
3072	farChild = in->GetBack(); // plane splits ray
3073	}
3074	else // one of the ray end points on the plane
3075	{ // NOTE: what to do if ray is coincident with plane?
3076	if (extSide < 0)
3077	node = in->GetBack();
3078	else
3079	node = in->GetFront();
3080
3081	continue; // no far child
3082	}
3083
3084	// push data for far child
3085	tQueue.push(BspRayTraversalData(farChild, extp));
3086
3087	// find intersection of ray segment with plane
3088	extp = splitPlane.FindIntersection(origin, extp, &t);
3089	}
3090	else
3091	{
3092	// reached leaf => intersection with view cell
3093	BspLeaf leaf = dynamic_cast<BspLeaf >(node);
3094	ViewCell *viewCell;
3095
3096	viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
3097	//viewCell = leaf->GetViewCell();
3098
3099	if (!viewCell->Mailed())
3100	{
3101	viewcells.push_back(viewCell);
3102	viewCell->Mail();
3103	++ hits;
3104	}
3105
3106	//-- fetch the next far child from the stack
3107	if (tQueue.empty())
3108	break;
3109
3110	entp = extp;
3111
3112	BspRayTraversalData &s = tQueue.top();
3113
3114	node = s.mNode;
3115	extp = s.mExitPoint;
3116
3117	tQueue.pop();
3118	}
3119	}
3120
3121	return hits;
3122	}
3123
3124
3125
3126
3127	int VspBspTree::TreeDistance(BspNode n1, BspNode n2) const
3128	{
3129	std::deque<BspNode *> path1;
3130	BspNode *p1 = n1;
3131
3132	// create path from node 1 to root
3133	while (p1)
3134	{
3135	if (p1 == n2) // second node on path
3136	return (int)path1.size();
3137
3138	path1.push_front(p1);
3139	p1 = p1->GetParent();
3140	}
3141
3142	int depth = n2->GetDepth();
3143	int d = depth;
3144
3145	BspNode *p2 = n2;
3146
3147	// compare with same depth
3148	while (1)
3149	{
3150	if ((d < (int)path1.size()) && (p2 == path1[d]))
3151	return (depth - d) + ((int)path1.size() - 1 - d);
3152
3153	-- d;
3154	p2 = p2->GetParent();
3155	}
3156
3157	return 0; // never come here
3158	}
3159
3160
3161	BspNode VspBspTree::CollapseTree(BspNode node, int &collapsed)
3162	{
3163	// TODO
3164	#if VC_HISTORY
3165	if (node->IsLeaf())
3166	return node;
3167
3168	BspInterior interior = dynamic_cast<BspInterior >(node);
3169
3170	BspNode *front = CollapseTree(interior->GetFront(), collapsed);
3171	BspNode *back = CollapseTree(interior->GetBack(), collapsed);
3172
3173	if (front->IsLeaf() && back->IsLeaf())
3174	{
3175	BspLeaf frontLeaf = dynamic_cast<BspLeaf >(front);
3176	BspLeaf backLeaf = dynamic_cast<BspLeaf >(back);
3177
3178	//-- collapse tree
3179	if (frontLeaf->GetViewCell() == backLeaf->GetViewCell())
3180	{
3181	BspViewCell *vc = frontLeaf->GetViewCell();
3182
3183	BspLeaf *leaf = new BspLeaf(interior->GetParent(), vc);
3184	leaf->SetTreeValid(frontLeaf->TreeValid());
3185
3186	// replace a link from node's parent
3187	if (leaf->GetParent())
3188	leaf->GetParent()->ReplaceChildLink(node, leaf);
3189	else
3190	mRoot = leaf;
3191
3192	++ collapsed;
3193	delete interior;
3194
3195	return leaf;
3196	}
3197	}
3198	#endif
3199	return node;
3200	}
3201
3202
3203	int VspBspTree::CollapseTree()
3204	{
3205	int collapsed = 0;
3206	//TODO
3207	#if VC_HISTORY
3208	(void)CollapseTree(mRoot, collapsed);
3209
3210	// revalidate leaves
3211	RepairViewCellsLeafLists();
3212	#endif
3213	return collapsed;
3214	}
3215
3216
3217	void VspBspTree::RepairViewCellsLeafLists()
3218	{
3219	// TODO
3220	#if VC_HISTORY
3221	// list not valid anymore => clear
3222	stack<BspNode *> nodeStack;
3223	nodeStack.push(mRoot);
3224
3225	ViewCell::NewMail();
3226
3227	while (!nodeStack.empty())
3228	{
3229	BspNode *node = nodeStack.top();
3230	nodeStack.pop();
3231
3232	if (node->IsLeaf())
3233	{
3234	BspLeaf leaf = dynamic_cast<BspLeaf >(node);
3235
3236	BspViewCell *viewCell = leaf->GetViewCell();
3237
3238	if (!viewCell->Mailed())
3239	{
3240	viewCell->mLeaves.clear();
3241	viewCell->Mail();
3242	}
3243
3244	viewCell->mLeaves.push_back(leaf);
3245
3246	}
3247	else
3248	{
3249	BspInterior interior = dynamic_cast<BspInterior >(node);
3250
3251	nodeStack.push(interior->GetFront());
3252	nodeStack.push(interior->GetBack());
3253	}
3254	}
3255	// TODO
3256	#endif
3257	}
3258
3259
3260	typedef pair<BspNode , BspNodeGeometry > bspNodePair;
3261
3262
3263	int VspBspTree::CastBeam(Beam &beam)
3264	{
3265	stack<bspNodePair> nodeStack;
3266	BspNodeGeometry *rgeom = new BspNodeGeometry();
3267	ConstructGeometry(mRoot, *rgeom);
3268
3269	nodeStack.push(bspNodePair(mRoot, rgeom));
3270
3271	ViewCell::NewMail();
3272
3273	while (!nodeStack.empty())
3274	{
3275	BspNode *node = nodeStack.top().first;
3276	BspNodeGeometry *geom = nodeStack.top().second;
3277	nodeStack.pop();
3278
3279	AxisAlignedBox3 box;
3280	box.Initialize();
3281	geom->IncludeInBox(box);
3282
3283	const int side = beam.ComputeIntersection(box);
3284
3285	switch (side)
3286	{
3287	case -1:
3288	CollectViewCells(node, true, beam.mViewCells, true);
3289	break;
3290	case 0:
3291
3292	if (node->IsLeaf())
3293	{
3294	BspLeaf leaf = dynamic_cast<BspLeaf >(node);
3295
3296	if (!leaf->GetViewCell()->Mailed() && leaf->TreeValid())
3297	{
3298	leaf->GetViewCell()->Mail();
3299	beam.mViewCells.push_back(leaf->GetViewCell());
3300	}
3301	}
3302	else
3303	{
3304	BspInterior interior = dynamic_cast<BspInterior >(node);
3305
3306	BspNode *first = interior->GetFront();
3307	BspNode *second = interior->GetBack();
3308
3309	BspNodeGeometry *firstGeom = new BspNodeGeometry();
3310	BspNodeGeometry *secondGeom = new BspNodeGeometry();
3311
3312	geom->SplitGeometry(*firstGeom,
3313	*secondGeom,
3314	interior->GetPlane(),
3315	mBox,
3316	//0.0000001f);
3317	mEpsilon);
3318
3319	// decide on the order of the nodes
3320	if (DotProd(beam.mPlanes[0].mNormal,
3321	interior->GetPlane().mNormal) > 0)
3322	{
3323	swap(first, second);
3324	swap(firstGeom, secondGeom);
3325	}
3326
3327	nodeStack.push(bspNodePair(first, firstGeom));
3328	nodeStack.push(bspNodePair(second, secondGeom));
3329	}
3330
3331	break;
3332	default:
3333	// default: cull
3334	break;
3335	}
3336
3337	DEL_PTR(geom);
3338
3339	}
3340
3341	return (int)beam.mViewCells.size();
3342	}
3343
3344
3345	void VspBspTree::SetViewCellsManager(ViewCellsManager *vcm)
3346	{
3347	mViewCellsManager = vcm;
3348	}
3349
3350
3351	int VspBspTree::CollectMergeCandidates(const vector<BspLeaf *> leaves,
3352	vector<MergeCandidate> &candidates)
3353	{
3354	BspLeaf::NewMail();
3355
3356	vector<BspLeaf *>::const_iterator it, it_end = leaves.end();
3357
3358	int numCandidates = 0;
3359
3360	// find merge candidates and push them into queue
3361	for (it = leaves.begin(); it != it_end; ++ it)
3362	{
3363	BspLeaf leaf = it;
3364
3365	// the same leaves must not be part of two merge candidates
3366	leaf->Mail();
3367
3368	vector<BspLeaf *> neighbors;
3369	if (0)
3370	FindNeighbors(leaf, neighbors, true);
3371	else
3372	FindApproximateNeighbors(leaf, neighbors, true);
3373	vector<BspLeaf *>::const_iterator nit, nit_end = neighbors.end();
3374
3375	// TODO: test if at least one ray goes from one leaf to the other
3376	for (nit = neighbors.begin(); nit != nit_end; ++ nit)
3377	{
3378	if ((*nit)->GetViewCell() != leaf->GetViewCell())
3379	{
3380	MergeCandidate mc(leaf->GetViewCell(), (*nit)->GetViewCell());
3381
3382	// dont't merge view cells if they are empty and not front and back leaf of the same parent
3383	// in the tree?
3384	if (mEmptyViewCellsMergeAllowed \|\|
3385	!leaf->GetViewCell()->GetPvs().Empty() \|\| !(*nit)->GetViewCell()->GetPvs().Empty() \|\|
3386	leaf->IsSibling(*nit))
3387	{
3388	candidates.push_back(mc);
3389	}
3390
3391	++ numCandidates;
3392	if ((numCandidates % 1000) == 0)
3393	{
3394	cout << "collected " << numCandidates << " merge candidates" << endl;
3395	}
3396	}
3397	}
3398	}
3399
3400	Debug << "merge queue: " << (int)candidates.size() << endl;
3401	Debug << "leaves in queue: " << numCandidates << endl;
3402
3403
3404	return (int)leaves.size();
3405	}
3406
3407
3408	int VspBspTree::CollectMergeCandidates(const VssRayContainer &rays,
3409	vector<MergeCandidate> &candidates)
3410	{
3411	ViewCell::NewMail();
3412	long startTime = GetTime();
3413
3414	map<BspLeaf , vector<BspLeaf> > neighborMap;
3415	ViewCellContainer::const_iterator iit;
3416
3417	int numLeaves = 0;
3418
3419	BspLeaf::NewMail();
3420
3421	for (int i = 0; i < (int)rays.size(); ++ i)
3422	{
3423	VssRay *ray = rays[i];
3424
3425	// traverse leaves stored in the rays and compare and
3426	// merge consecutive leaves (i.e., the neighbors in the tree)
3427	if (ray->mViewCells.size() < 2)
3428	continue;
3429	//TODO viewcellhierarchy
3430	iit = ray->mViewCells.begin();
3431	BspViewCell bspVc = dynamic_cast<BspViewCell >(*(iit ++));
3432	BspLeaf *leaf = bspVc->mLeaf;
3433
3434	// traverse intersections
3435	// consecutive leaves are neighbors => add them to queue
3436	for (; iit != ray->mViewCells.end(); ++ iit)
3437	{
3438	// next pair
3439	BspLeaf *prevLeaf = leaf;
3440	bspVc = dynamic_cast<BspViewCell >(iit);
3441	leaf = bspVc->mLeaf;
3442
3443	// view space not valid or same view cell
3444	if (!leaf->TreeValid() \|\| !prevLeaf->TreeValid() \|\|
3445	(leaf->GetViewCell() == prevLeaf->GetViewCell()))
3446	continue;
3447
3448	vector<BspLeaf *> &neighbors = neighborMap[leaf];
3449
3450	bool found = false;
3451
3452	// both leaves inserted in queue already =>
3453	// look if double pair already exists
3454	if (leaf->Mailed() && prevLeaf->Mailed())
3455	{
3456	vector<BspLeaf *>::const_iterator it, it_end = neighbors.end();
3457
3458	for (it = neighbors.begin(); !found && (it != it_end); ++ it)
3459	if (*it == prevLeaf)
3460	found = true; // already in queue
3461	}
3462
3463	if (!found)
3464	{
3465	// this pair is not in map yet
3466	// => insert into the neighbor map and the queue
3467	neighbors.push_back(prevLeaf);
3468	neighborMap[prevLeaf].push_back(leaf);
3469
3470	leaf->Mail();
3471	prevLeaf->Mail();
3472
3473	MergeCandidate mc(leaf->GetViewCell(), prevLeaf->GetViewCell());
3474
3475	candidates.push_back(mc);
3476
3477	if (((int)candidates.size() % 1000) == 0)
3478	{
3479	cout << "collected " << (int)candidates.size() << " merge candidates" << endl;
3480	}
3481	}
3482	}
3483	}
3484
3485	Debug << "neighbormap size: " << (int)neighborMap.size() << endl;
3486	Debug << "merge queue: " << (int)candidates.size() << endl;
3487	Debug << "leaves in queue: " << numLeaves << endl;
3488
3489
3490	//-- collect the leaves which haven't been found by ray casting
3491	if (0)
3492	{
3493	cout << "finding additional merge candidates using geometry" << endl;
3494	vector<BspLeaf *> leaves;
3495	CollectLeaves(leaves, true);
3496	Debug << "found " << (int)leaves.size() << " new leaves" << endl << endl;
3497	CollectMergeCandidates(leaves, candidates);
3498	}
3499
3500	return numLeaves;
3501	}
3502
3503
3504
3505
3506	ViewCell *VspBspTree::GetViewCell(const Vector3 &point)
3507	{
3508	if (mRoot == NULL)
3509	return NULL;
3510
3511	stack<BspNode *> nodeStack;
3512	nodeStack.push(mRoot);
3513
3514	ViewCell *viewcell = NULL;
3515
3516	while (!nodeStack.empty()) {
3517	BspNode *node = nodeStack.top();
3518	nodeStack.pop();
3519
3520	if (node->IsLeaf()) {
3521	viewcell = dynamic_cast<BspLeaf *>(node)->GetViewCell();
3522	break;
3523	} else {
3524
3525	BspInterior interior = dynamic_cast<BspInterior >(node);
3526
3527	// random decision
3528	if (interior->GetPlane().Side(point) < 0)
3529	nodeStack.push(interior->GetBack());
3530	else
3531	nodeStack.push(interior->GetFront());
3532	}
3533	}
3534
3535	return viewcell;
3536	}
3537
3538
3539	bool VspBspTree::ViewPointValid(const Vector3 &viewPoint) const
3540	{
3541	BspNode *node = mRoot;
3542
3543	while (1)
3544	{
3545	// early exit
3546	if (node->TreeValid())
3547	return true;
3548
3549	if (node->IsLeaf())
3550	return false;
3551
3552	BspInterior in = dynamic_cast<BspInterior >(node);
3553
3554	if (in->GetPlane().Side(viewPoint) <= 0)
3555	{
3556	node = in->GetBack();
3557	}
3558	else
3559	{
3560	node = in->GetFront();
3561	}
3562	}
3563
3564	// should never come here
3565	return false;
3566	}
3567
3568
3569	void VspBspTree::PropagateUpValidity(BspNode *node)
3570	{
3571	const bool isValid = node->TreeValid();
3572
3573	// propagative up invalid flag until only invalid nodes exist over this node
3574	if (!isValid)
3575	{
3576	while (!node->IsRoot() && node->GetParent()->TreeValid())
3577	{
3578	node = node->GetParent();
3579	node->SetTreeValid(false);
3580	}
3581	}
3582	else
3583	{
3584	// propagative up valid flag until one of the subtrees is invalid
3585	while (!node->IsRoot() && !node->TreeValid())
3586	{
3587	node = node->GetParent();
3588	BspInterior interior = dynamic_cast<BspInterior >(node);
3589
3590	// the parent is valid iff both leaves are valid
3591	node->SetTreeValid(interior->GetBack()->TreeValid() &&
3592	interior->GetFront()->TreeValid());
3593	}
3594	}
3595	}
3596
3597
3598	bool VspBspTree::Export(ofstream &stream)
3599	{
3600	ExportNode(mRoot, stream);
3601
3602	return true;
3603	}
3604
3605
3606	void VspBspTree::ExportNode(BspNode *node, ofstream &stream)
3607	{
3608	if (node->IsLeaf())
3609	{
3610	BspLeaf leaf = dynamic_cast<BspLeaf >(node);
3611	ViewCell *viewCell = mViewCellsTree->GetActiveViewCell(leaf->GetViewCell());
3612
3613	int id = -1;
3614	if (viewCell != mOutOfBoundsCell)
3615	id = viewCell->GetId();
3616
3617	stream << "<Leaf viewCellId=\"" << id << "\" />" << endl;
3618	}
3619	else
3620	{
3621	BspInterior interior = dynamic_cast<BspInterior >(node);
3622
3623	Plane3 plane = interior->GetPlane();
3624	stream << "<Interior plane=\"" << plane.mNormal.x << " "
3625	<< plane.mNormal.y << " " << plane.mNormal.z << " "
3626	<< plane.mD << "\">" << endl;
3627
3628	ExportNode(interior->GetBack(), stream);
3629	ExportNode(interior->GetFront(), stream);
3630
3631	stream << "</Interior>" << endl;
3632	}
3633	}

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

Download in other formats: