Context Navigation

source: trunk/VUT/GtpVisibilityPreprocessor/src/VspBspTree.cpp @ 568

Revision 568, 86.9 KB checked in by mattausch, 18 years ago (diff)
debug version for testing shuffling and visualizing neighbors

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

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

Download in other formats: