Context Navigation

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

Revision 643, 74.7 KB checked in by mattausch, 18 years ago (diff)
made demo running

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

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

Download in other formats: