Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellBsp.cpp @ 1563

Revision 1563, 83.0 KB checked in by mattausch, 18 years ago (diff)
fixed bug with view space box

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

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

Download in other formats: