Context Navigation

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

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

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

Download in other formats: