Context Navigation

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

Revision 1027, 81.1 KB checked in by mattausch, 18 years ago (diff)

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

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

Download in other formats: