Context Navigation

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

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

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

Download in other formats: