Context Navigation

source: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp @ 453

Revision 453, 44.0 KB checked in by mattausch, 19 years ago (diff)

Rev	Line
[408]	1
	2	// ================================================================
	3	// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
	4	// ****************************************************************
	5	/**
	6	The KD tree based LSDS
	7	*/
	8	// Initial coding by
	9	/**
	10	@author Jiri Bittner
	11	*/
	12
	13	// Standard headers
[406]	14	#include <stack>
[408]	15	#include <queue>
[406]	16	#include <algorithm>
[408]	17	#include <fstream>
	18	#include <string>
	19
	20	#include "VspKdTree.h"
	21
[406]	22	#include "Environment.h"
[408]	23	#include "VssRay.h"
	24	#include "Intersectable.h"
	25	#include "Ray.h"
[420]	26	#include "RayInfo.h"
[453]	27	#include "ViewCellsManager.h"
	28	#include "ViewCell.h"
[406]	29
[418]	30	// Static variables
	31	int VspKdTreeLeaf::mailID = 0;
[406]	32
[428]	33	#define USE_FIXEDPOINT_T 0
	34
[420]	35	/// Adds object to the pvs of the front and back node
[409]	36	inline void AddObject2Pvs(Intersectable *object,
	37	const int side,
	38	int &pvsBack,
	39	int &pvsFront)
[406]	40	{
[408]	41	if (!object)
	42	return;
	43
[409]	44	if (side <= 0)
	45	{
	46	if (!object->Mailed() && !object->Mailed(2))
	47	{
	48	++ pvsBack;
	49
[408]	50	if (object->Mailed(1))
	51	object->Mail(2);
	52	else
	53	object->Mail();
	54	}
	55	}
	56
[428]	57	if (side >= 0)
[409]	58	{
	59	if (!object->Mailed(1) && !object->Mailed(2))
	60	{
	61	++ pvsFront;
[428]	62
[408]	63	if (object->Mailed())
	64	object->Mail(2);
	65	else
	66	object->Mail(1);
	67	}
	68	}
[406]	69	}
	70
[418]	71	/**************************************************************/
[420]	72	/* class VspKdTreeNode implementation */
	73	/**************************************************************/
	74
	75	// Inline constructor
	76	VspKdTreeNode::VspKdTreeNode(VspKdTreeInterior *p):
	77	mParent(p), mAxis(-1), mDepth(p ? p->mDepth + 1 : 0)
	78	{}
	79
	80	VspKdTreeNode::~VspKdTreeNode()
	81	{};
	82
	83	inline VspKdTreeInterior *VspKdTreeNode::GetParent() const
	84	{
	85	return mParent;
	86	}
	87
	88	inline void VspKdTreeNode::SetParent(VspKdTreeInterior *p)
	89	{
	90	mParent = p;
	91	}
	92
	93	bool VspKdTreeNode::IsLeaf() const
	94	{
	95	return mAxis == -1;
	96	}
	97
	98	int VspKdTreeNode::GetAccessTime()
	99	{
	100	return 0x7FFFFFF;
	101	}
	102
	103	/**************************************************************/
[419]	104	/* VspKdTreeInterior implementation */
[418]	105	/**************************************************************/
[408]	106
[418]	107	VspKdTreeInterior::VspKdTreeInterior(VspKdTreeInterior *p):
	108	VspKdTreeNode(p), mBack(NULL), mFront(NULL), mAccesses(0), mLastAccessTime(-1)
	109	{
	110	}
	111
	112	int VspKdTreeInterior::GetAccessTime()
	113	{
[419]	114	return mLastAccessTime;
[418]	115	}
	116
	117	void VspKdTreeInterior::SetupChildLinks(VspKdTreeNode b, VspKdTreeNode f)
	118	{
	119	mBack = b;
	120	mFront = f;
[420]	121	b->SetParent(this);
	122	f->SetParent(this);
[418]	123	}
	124
[423]	125	void VspKdTreeInterior::ReplaceChildLink(VspKdTreeNode *oldChild,
	126	VspKdTreeNode *newChild)
[418]	127	{
	128	if (mBack == oldChild)
	129	mBack = newChild;
	130	else
	131	mFront = newChild;
	132	}
	133
	134	int VspKdTreeInterior::Type() const
	135	{
[419]	136	return EInterior;
[418]	137	}
	138
	139	VspKdTreeInterior::~VspKdTreeInterior()
	140	{
	141	DEL_PTR(mBack);
	142	DEL_PTR(mFront);
	143	}
	144
	145	void VspKdTreeInterior::Print(ostream &s) const
	146	{
[423]	147	switch (mAxis)
	148	{
	149	case 0:
	150	s << "x "; break;
	151	case 1:
	152	s << "y "; break;
	153	case 2:
	154	s << "z "; break;
	155	}
	156
[419]	157	s << mPosition << " ";
[418]	158
	159	mBack->Print(s);
	160	mFront->Print(s);
	161	}
	162
	163	int VspKdTreeInterior::ComputeRayIntersection(const RayInfo &rayData, float &t)
	164	{
[419]	165	return rayData.ComputeRayIntersection(mAxis, mPosition, t);
[418]	166	}
	167
[419]	168	VspKdTreeNode *VspKdTreeInterior::GetBack() const
	169	{
	170	return mBack;
	171	}
[418]	172
[419]	173	VspKdTreeNode *VspKdTreeInterior::GetFront() const
[406]	174	{
[419]	175	return mFront;
	176	}
[406]	177
[420]	178
[425]	179	/**************************************************************/
	180	/* class VspKdTreeLeaf implementation */
	181	/**************************************************************/
[420]	182	VspKdTreeLeaf::VspKdTreeLeaf(VspKdTreeInterior *p, const int nRays):
[452]	183	VspKdTreeNode(p), mRays(), mPvsSize(0), mValidPvs(false), mViewCell(NULL)
[420]	184	{
	185	mRays.reserve(nRays);
	186	}
	187
	188	VspKdTreeLeaf::~VspKdTreeLeaf()
	189	{}
	190
	191	int VspKdTreeLeaf::Type() const
	192	{
	193	return ELeaf;
	194	}
	195
	196	void VspKdTreeLeaf::Print(ostream &s) const
	197	{
	198	s << endl << "L: r = " << (int)mRays.size() << endl;
	199	};
	200
	201	void VspKdTreeLeaf::AddRay(const RayInfo &data)
	202	{
	203	mValidPvs = false;
	204	mRays.push_back(data);
	205	data.mRay->Ref();
	206	}
	207
	208	int VspKdTreeLeaf::GetPvsSize() const
	209	{
	210	return mPvsSize;
	211	}
	212
	213	void VspKdTreeLeaf::SetPvsSize(const int s)
	214	{
	215	mPvsSize = s;
	216	}
	217
	218	void VspKdTreeLeaf::Mail()
	219	{
	220	mMailbox = mailID;
	221	}
	222
	223	void VspKdTreeLeaf::NewMail()
	224	{
	225	++ mailID;
	226	}
	227
	228	bool VspKdTreeLeaf::Mailed() const
	229	{
	230	return mMailbox == mailID;
	231	}
	232
	233	bool VspKdTreeLeaf::Mailed(const int mail)
	234	{
	235	return mMailbox >= mailID + mail;
	236	}
	237
	238	float VspKdTreeLeaf::GetAvgRayContribution() const
	239	{
	240	return GetPvsSize() / ((float)mRays.size() + Limits::Small);
	241	}
	242
	243	float VspKdTreeLeaf::GetSqrRayContribution() const
	244	{
	245	return sqr(GetAvgRayContribution());
	246	}
	247
[426]	248	RayInfoContainer &VspKdTreeLeaf::GetRays()
	249	{
	250	return mRays;
	251	}
[428]	252
	253	void VspKdTreeLeaf::ExtractPvs(ObjectContainer &objects) const
	254	{
	255	RayInfoContainer::const_iterator it, it_end = mRays.end();
	256
	257	for (it = mRays.begin(); it != it_end; ++ it)
	258	{
	259	if ((*it).mRay->mTerminationObject)
	260	objects.push_back((*it).mRay->mTerminationObject);
	261	if ((*it).mRay->mOriginObject)
	262	objects.push_back((*it).mRay->mOriginObject);
	263	}
	264	}
	265
	266	void VspKdTreeLeaf::GetRays(VssRayContainer &rays)
	267	{
	268	RayInfoContainer::const_iterator it, it_end = mRays.end();
	269
	270	for (it = mRays.begin(); it != mRays.end(); ++ it)
	271	rays.push_back((*it).mRay);
	272	}
	273
[453]	274	ViewCell *VspKdTreeLeaf::GetViewCell()
	275	{
	276	return mViewCell;
	277	}
	278
[420]	279	/*********************************************************/
	280	/* class VspKdTree implementation */
	281	/*********************************************************/
	282
	283
[428]	284	VspKdTree::VspKdTree(): mOnlyDrivingAxis(true)
[419]	285	{
[421]	286	environment->GetIntValue("VspKdTree.Termination.maxDepth", mTermMaxDepth);
	287	environment->GetIntValue("VspKdTree.Termination.minPvs", mTermMinPvs);
	288	environment->GetIntValue("VspKdTree.Termination.minRays", mTermMinRays);
	289	environment->GetFloatValue("VspKdTree.Termination.maxRayContribution", mTermMaxRayContribution);
	290	environment->GetFloatValue("VspKdTree.Termination.maxCostRatio", mTermMaxCostRatio);
[422]	291	environment->GetFloatValue("VspKdTree.Termination.minSize", mTermMinSize);
[421]	292
	293	mTermMinSize = sqr(mTermMinSize);
	294
[422]	295	environment->GetFloatValue("VspKdTree.epsilon", mEpsilon);
	296	environment->GetFloatValue("VspKdTree.ct_div_ci", mCtDivCi);
[408]	297
[422]	298	environment->GetFloatValue("VspKdTree.maxTotalMemory", mMaxTotalMemory);
	299	environment->GetFloatValue("VspKdTree.maxStaticMemory", mMaxStaticMemory);
	300
	301	environment->GetIntValue("VspKdTree.accessTimeThreshold", mAccessTimeThreshold);
	302	environment->GetIntValue("VspKdTree.minCollapseDepth", mMinCollapseDepth);
[408]	303
[409]	304	// split type
[408]	305	char sname[128];
	306	environment->GetStringValue("VspKdTree.splitType", sname);
[409]	307	string name(sname);
[423]	308
	309	Debug << "======= vsp kd tree options ========" << endl;
	310	Debug << "max depth: "<< mTermMaxDepth << endl;
	311	Debug << "min pvs: "<< mTermMinPvs << endl;
	312	Debug << "min rays: "<< mTermMinRays << endl;
	313	Debug << "max ray contribution: "<< mTermMaxRayContribution << endl;
	314	Debug << "max cost ratio: "<< mTermMaxCostRatio << endl;
	315	Debug << "min size: "<<mTermMinSize << endl;
	316
[409]	317	if (name.compare("regular") == 0)
[423]	318	{
	319	Debug << "using regular split" << endl;
[409]	320	splitType = ESplitRegular;
[423]	321	}
[409]	322	else
	323	{
	324	if (name.compare("heuristic") == 0)
[423]	325	{
	326	Debug << "using heuristic split" << endl;
[409]	327	splitType = ESplitHeuristic;
[423]	328	}
[409]	329	else
	330	{
	331	cerr << "Invalid VspKdTree split type " << name << endl;
[408]	332	exit(1);
	333	}
[409]	334	}
[408]	335
[418]	336	mRoot = NULL;
[409]	337
[422]	338	mSplitCandidates = new vector<SortableEntry>;
[406]	339	}
	340
[408]	341
	342	VspKdTree::~VspKdTree()
[406]	343	{
[419]	344	DEL_PTR(mRoot);
[449]	345	DEL_PTR(mSplitCandidates);
[408]	346	}
[406]	347
[418]	348	void VspKdStatistics::Print(ostream &app) const
[408]	349	{
[426]	350	app << "===== VspKdTree statistics ===============\n";
[408]	351
[426]	352	app << "#N_RAYS ( Number of rays )\n"
	353	<< rays << endl;
[410]	354
[426]	355	app << "#N_INITPVS ( Initial PVS size )\n"
	356	<< initialPvsSize << endl;
[406]	357
[426]	358	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
[406]	359
[426]	360	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
[406]	361
[426]	362	app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
	363
	364	for (int i=0; i<7; i++)
	365	app << splits[i] <<" ";
	366	app << endl;
[406]	367
[426]	368	app << "#N_RAYREFS ( Number of rayRefs )\n" <<
	369	rayRefs << "\n";
[408]	370
[426]	371	app << "#N_RAYRAYREFS ( Number of rayRefs / ray )\n" <<
	372	rayRefs / (double)rays << "\n";
[408]	373
[426]	374	app << "#N_LEAFRAYREFS ( Number of rayRefs / leaf )\n" <<
	375	rayRefs / (double)Leaves() << "\n";
[408]	376
[426]	377	app << "#N_MAXRAYREFS ( Max number of rayRefs / leaf )\n" <<
	378	maxRayRefs << "\n";
[408]	379
	380	// app << setprecision(4);
	381
[426]	382	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n" <<
	383	maxDepthNodes * 100 / (double)Leaves() << endl;
[408]	384
[426]	385	app << "#N_PMINPVSLEAVES ( Percentage of leaves with mininimal PVS )\n" <<
	386	minPvsNodes * 100 / (double)Leaves() << endl;
[408]	387
[426]	388	app << "#N_PMINRAYSLEAVES ( Percentage of leaves with minimal number of rays)\n" <<
	389	minRaysNodes * 100 / (double)Leaves() << endl;
	390
[408]	391	app << "#N_PMINSIZELEAVES ( Percentage of leaves with minSize )\n"<<
[426]	392	minSizeNodes * 100 / (double)Leaves() << endl;
[408]	393
[426]	394	app << "#N_PMAXRAYCONTRIBLEAVES ( Percentage of leaves with maximal ray contribution )\n" <<
	395	maxRayContribNodes * 100 / (double)Leaves() << endl;
[408]	396
[426]	397	app << "#N_ADDED_RAYREFS ( Number of dynamically added ray references )\n"<<
	398	addedRayRefs << endl;
[408]	399
[426]	400	app << "#N_REMOVED_RAYREFS ( Number of dynamically removed ray references )\n"<<
	401	removedRayRefs << endl;
[408]	402
[426]	403	// app << setprecision(4);
[408]	404
[426]	405	app << "#N_MAXPVS ( Maximal PVS size / leaf)\n"
	406	<< maxPvsSize << endl;
[408]	407
[426]	408	app << "#N_CTIME ( Construction time [s] )\n"
	409	<< Time() << " \n";
[408]	410
[426]	411	app << "===== END OF VspKdTree statistics ==========\n";
[406]	412	}
	413
[408]	414
	415	void
	416	VspKdTreeLeaf::UpdatePvsSize()
[406]	417	{
[409]	418	if (!mValidPvs)
	419	{
[408]	420	Intersectable::NewMail();
	421	int pvsSize = 0;
[420]	422	for(RayInfoContainer::iterator ri = mRays.begin();
	423	ri != mRays.end(); ++ ri)
[409]	424	{
	425	if ((*ri).mRay->IsActive())
	426	{
[408]	427	Intersectable *object;
	428	#if BIDIRECTIONAL_RAY
	429	object = (*ri).mRay->mOriginObject;
[409]	430
	431	if (object && !object->Mailed())
	432	{
	433	++ pvsSize;
[408]	434	object->Mail();
	435	}
	436	#endif
	437	object = (*ri).mRay->mTerminationObject;
[409]	438	if (object && !object->Mailed())
	439	{
	440	++ pvsSize;
[408]	441	object->Mail();
	442	}
	443	}
[409]	444	}
	445
[408]	446	mPvsSize = pvsSize;
	447	mValidPvs = true;
	448	}
	449	}
[406]	450
	451
[408]	452	void
[440]	453	VspKdTree::Construct(const VssRayContainer &rays,
[409]	454	AxisAlignedBox3 *forcedBoundingBox)
[452]	455	{
[414]	456	mStat.Start();
[452]	457
[422]	458	mMaxMemory = mMaxStaticMemory;
[419]	459	DEL_PTR(mRoot);
[408]	460
[419]	461	mRoot = new VspKdTreeLeaf(NULL, (int)rays.size());
[408]	462
[423]	463	// first construct a leaf that will get subdivided
[419]	464	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(mRoot);
[409]	465
[414]	466	mStat.nodes = 1;
[448]	467	mBox.Initialize();
[452]	468
[448]	469	//-- compute bounding box
[409]	470	if (forcedBoundingBox)
[419]	471	mBox = *forcedBoundingBox;
[448]	472	else
	473	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
	474	{
[452]	475	if ((*ri)->mOriginObject)
	476	mBox.Include((*ri)->GetOrigin());
	477	if ((*ri)->mTerminationObject)
	478	mBox.Include((*ri)->GetTermination());
[448]	479	}
[408]	480
[422]	481	Debug << "bbox = " << mBox << endl;
[448]	482
	483	for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
	484	{
	485	//leaf->AddRay(RayInfo(*ri));
	486	VssRay ray = ri;
	487
	488	float minT, maxT;
	489
	490	// TODO: not very efficient to implictly cast between rays types ...
	491	if (mBox.GetRaySegment(*ray, minT, maxT))
	492	{
	493	float len = ray->Length();
	494
	495	if (!len)
	496	len = Limits::Small;
	497
	498	leaf->AddRay(RayInfo(ray, minT / len, maxT / len));
	499	}
	500	}
[419]	501
[420]	502	mStat.rays = (int)leaf->mRays.size();
[408]	503	leaf->UpdatePvsSize();
[409]	504
[414]	505	mStat.initialPvsSize = leaf->GetPvsSize();
[409]	506
	507	// Subdivide();
[419]	508	mRoot = Subdivide(TraversalData(leaf, mBox, 0));
[408]	509
[422]	510	if (mSplitCandidates)
[409]	511	{
[422]	512	// force release of this vector
	513	delete mSplitCandidates;
	514	mSplitCandidates = new vector<SortableEntry>;
[409]	515	}
[408]	516
[414]	517	mStat.Stop();
[408]	518
[414]	519	mStat.Print(cout);
[423]	520	Debug << "#Total memory=" << GetMemUsage() << endl;
[408]	521	}
	522
	523
	524
[409]	525	VspKdTreeNode *VspKdTree::Subdivide(const TraversalData &tdata)
[408]	526	{
[409]	527	VspKdTreeNode *result = NULL;
[408]	528
[409]	529	priority_queue<TraversalData> tStack;
[423]	530	//stack<TraversalData> tStack;
[406]	531
[409]	532	tStack.push(tdata);
[406]	533
[409]	534	AxisAlignedBox3 backBox;
	535	AxisAlignedBox3 frontBox;
[406]	536
[408]	537	int lastMem = 0;
[423]	538
[409]	539	while (!tStack.empty())
	540	{
[408]	541	float mem = GetMemUsage();
	542
[423]	543	if (lastMem / 10 != ((int)mem) / 10)
[409]	544	{
[425]	545	Debug << mem << " MB" << endl;
[408]	546	}
	547	lastMem = (int)mem;
	548
[423]	549	if (1 && mem > mMaxMemory)
[409]	550	{
[423]	551	Debug << "memory limit reached: " << mem << endl;
[409]	552	// count statistics on unprocessed leafs
	553	while (!tStack.empty())
	554	{
[408]	555	EvaluateLeafStats(tStack.top());
	556	tStack.pop();
[409]	557	}
	558	break;
	559	}
[408]	560
[409]	561	TraversalData data = tStack.top();
[423]	562	tStack.pop();
[409]	563
[419]	564	VspKdTreeNode node = SubdivideNode((VspKdTreeLeaf ) data.mNode,
	565	data.mBox, backBox, frontBox);
[409]	566	if (result == NULL)
	567	result = node;
[406]	568
[409]	569	if (!node->IsLeaf())
	570	{
[423]	571	VspKdTreeInterior interior = dynamic_cast<VspKdTreeInterior >(node);
[408]	572
[409]	573	// push the children on the stack
[419]	574	tStack.push(TraversalData(interior->GetBack(), backBox, data.mDepth + 1));
	575	tStack.push(TraversalData(interior->GetFront(), frontBox, data.mDepth + 1));
[409]	576	}
	577	else
	578	{
	579	EvaluateLeafStats(data);
	580	}
	581	}
	582
	583	return result;
[408]	584	}
[406]	585
[408]	586
	587	// returns selected plane for subdivision
[424]	588	int VspKdTree::SelectPlane(VspKdTreeLeaf *leaf,
	589	const AxisAlignedBox3 &box,
	590	float &position,
	591	int &raysBack,
	592	int &raysFront,
	593	int &pvsBack,
	594	int &pvsFront)
[408]	595	{
[410]	596	int minDirDepth = 6;
[428]	597	int axis = 0;
	598	float costRatio = 0;
[408]	599
[422]	600	if (splitType == ESplitRegular)
	601	{
[408]	602	costRatio = BestCostRatioRegular(leaf,
[410]	603	axis,
	604	position,
	605	raysBack,
	606	raysFront,
	607	pvsBack,
	608	pvsFront);
	609	}
[423]	610	else if (splitType == ESplitHeuristic)
[410]	611	{
	612	costRatio = BestCostRatioHeuristic(leaf,
	613	axis,
	614	position,
	615	raysBack,
	616	raysFront,
	617	pvsBack,
	618	pvsFront);
[423]	619	}
	620	else
	621	{
	622	cerr << "VspKdTree: Unknown split heuristics\n";
	623	exit(1);
	624	}
[408]	625
[428]	626	if (costRatio > mTermMaxCostRatio)
[410]	627	{
[425]	628	Debug << "Too big cost ratio " << costRatio << endl;
[408]	629	return -1;
	630	}
	631
[428]	632	if (0)
	633	Debug << "pvs=" << leaf->mPvsSize
	634	<< " rays=" << (int)leaf->mRays.size()
	635	<< " rc=" << leaf->GetAvgRayContribution()
	636	<< " axis=" << axis << endl;
[408]	637
	638	return axis;
[406]	639	}
	640
	641
[408]	642
[425]	643	float VspKdTree::EvalCostRatio(VspKdTreeLeaf *leaf,
	644	const int axis,
	645	const float position,
	646	int &raysBack,
	647	int &raysFront,
	648	int &pvsBack,
	649	int &pvsFront)
[406]	650	{
[408]	651	raysBack = 0;
	652	raysFront = 0;
	653	pvsFront = 0;
	654	pvsBack = 0;
	655
	656	float newCost;
	657
	658	Intersectable::NewMail(3);
	659
	660	// eval pvs size
	661	int pvsSize = leaf->GetPvsSize();
	662
	663	Intersectable::NewMail(3);
	664
[419]	665	// this is the main ray classification loop!
[420]	666	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[424]	667	ri != leaf->mRays.end(); ++ ri)
[419]	668	{
	669	if (!(*ri).mRay->IsActive())
	670	continue;
[410]	671
[419]	672	// determine the side of this ray with respect to the plane
	673	int side = (ri).ComputeRayIntersection(axis, position, (ri).mRay->mT);
	674	// (*ri).mRay->mSide = side;
	675
	676	if (side <= 0)
	677	++ raysBack;
[408]	678
[419]	679	if (side >= 0)
	680	++ raysFront;
[408]	681
[419]	682	AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
	683	}
[408]	684
[428]	685	if (0)
[425]	686	{
[428]	687	const float sum = float(pvsBack + pvsFront);
[425]	688	const float oldCost = (float)pvsSize * 2;
	689
	690	return sum / oldCost;
	691	}
[428]	692	else
	693	{
[425]	694	AxisAlignedBox3 box = GetBBox(leaf);
	695
	696	float minBox = box.Min(axis);
	697	float maxBox = box.Max(axis);
	698	float sizeBox = maxBox - minBox;
[408]	699
[425]	700	// float sum = raysBack(position - minBox) + raysFront(maxBox - position);
[428]	701	const float sum = pvsBack * (position - minBox) + pvsFront * (maxBox - position);
[408]	702
[425]	703	newCost = mCtDivCi + sum / sizeBox;
[419]	704
[425]	705	//Debug << axis << " " << pvsSize << " " << pvsBack << " " << pvsFront << endl;
	706	// float oldCost = leaf->mRays.size();
[428]	707	const float oldCost = (float)pvsSize;
[425]	708
[452]	709	return newCost / oldCost;
[425]	710	}
[406]	711	}
	712
[422]	713	float VspKdTree::BestCostRatioRegular(VspKdTreeLeaf *leaf,
	714	int &axis,
	715	float &position,
	716	int &raysBack,
	717	int &raysFront,
	718	int &pvsBack,
	719	int &pvsFront)
[408]	720	{
[422]	721	int nRaysBack[3], nRaysFront[3];
	722	int nPvsBack[3], nPvsFront[3];
	723
	724	float nPosition[3];
	725	float nCostRatio[3];
[408]	726	int bestAxis = -1;
	727
	728	AxisAlignedBox3 sBox = GetBBox(leaf);
[410]	729
[428]	730	const int sAxis = sBox.Size().DrivingAxis();
[408]	731
[422]	732	for (axis = 0; axis < 3; ++ axis)
[410]	733	{
[422]	734	if (!mOnlyDrivingAxis \|\| axis == sAxis)
[410]	735	{
[452]	736
[424]	737	nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis]) * 0.5f;
[419]	738
[408]	739	nCostRatio[axis] = EvalCostRatio(leaf,
[410]	740	axis,
	741	nPosition[axis],
	742	nRaysBack[axis],
	743	nRaysFront[axis],
	744	nPvsBack[axis],
	745	nPvsFront[axis]);
	746	if (bestAxis == -1)
[452]	747	{
[408]	748	bestAxis = axis;
[452]	749	}
	750	else if (nCostRatio[axis] < nCostRatio[bestAxis])
	751	{
	752	/*Debug << "pvs front " << nPvsBack[axis]
	753	<< " pvs back " << nPvsFront[axis]
	754	<< " overall pvs " << leaf->GetPvsSize() << endl;*/
	755	bestAxis = axis;
	756	}
	757
[408]	758	}
	759	}
	760
[428]	761	//-- assign best axis
[408]	762	axis = bestAxis;
	763	position = nPosition[bestAxis];
	764
	765	raysBack = nRaysBack[bestAxis];
	766	raysFront = nRaysFront[bestAxis];
	767
	768	pvsBack = nPvsBack[bestAxis];
	769	pvsFront = nPvsFront[bestAxis];
	770
	771	return nCostRatio[bestAxis];
	772	}
	773
[410]	774	float VspKdTree::BestCostRatioHeuristic(VspKdTreeLeaf *leaf,
	775	int &axis,
	776	float &position,
	777	int &raysBack,
	778	int &raysFront,
	779	int &pvsBack,
	780	int &pvsFront)
[406]	781	{
[408]	782	AxisAlignedBox3 box = GetBBox(leaf);
	783	// AxisAlignedBox3 dirBox = GetDirBBox(node);
	784
	785	axis = box.Size().DrivingAxis();
	786
	787	SortSplitCandidates(leaf, axis);
[406]	788
[408]	789	// go through the lists, count the number of objects left and right
	790	// and evaluate the following cost funcion:
	791	// C = ct_div_ci + (qlrl + qrrr)/queries
	792
[420]	793	int rl = 0, rr = (int)leaf->mRays.size();
	794	int pl = 0, pr = leaf->GetPvsSize();
[424]	795
[408]	796	float minBox = box.Min(axis);
	797	float maxBox = box.Max(axis);
	798	float sizeBox = maxBox - minBox;
	799
[419]	800	float minBand = minBox + 0.1f*(maxBox - minBox);
	801	float maxBand = minBox + 0.9f*(maxBox - minBox);
[408]	802
	803	float sum = rr*sizeBox;
[411]	804	float minSum = 1e20f;
[408]	805
	806	Intersectable::NewMail();
[424]	807
[408]	808	// set all object as belonging to the fron pvs
[420]	809	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
[424]	810	ri != leaf->mRays.end(); ++ ri)
[410]	811	{
	812	if ((*ri).mRay->IsActive())
	813	{
[408]	814	Intersectable object = (ri).mRay->mTerminationObject;
[410]	815
[408]	816	if (object)
[410]	817	if (!object->Mailed())
	818	{
[408]	819	object->Mail();
	820	object->mCounter = 1;
[410]	821	}
	822	else
	823	++ object->mCounter;
[408]	824	}
[410]	825	}
[408]	826
	827	Intersectable::NewMail();
	828
[422]	829	for (vector<SortableEntry>::const_iterator ci = mSplitCandidates->begin();
	830	ci < mSplitCandidates->end(); ++ ci)
[410]	831	{
[408]	832	VssRay *ray;
[410]	833
	834	switch ((*ci).type)
	835	{
	836	case SortableEntry::ERayMin:
	837	{
	838	++ rl;
	839	ray = (VssRay ) (ci).data;
	840	Intersectable *object = ray->mTerminationObject;
	841	if (object && !object->Mailed())
	842	{
	843	object->Mail();
	844	++ pl;
	845	}
	846	break;
	847	}
	848	case SortableEntry::ERayMax:
	849	{
	850	-- rr;
	851
	852	ray = (VssRay ) (ci).data;
	853	Intersectable *object = ray->mTerminationObject;
	854
	855	if (object)
	856	{
	857	if (-- object->mCounter == 0)
	858	-- pr;
	859	}
	860
	861	break;
	862	}
[408]	863	}
	864
[410]	865	if ((ci).value > minBand && (ci).value < maxBand)
	866	{
[408]	867	sum = pl((ci).value - minBox) + pr(maxBox - (ci).value);
[410]	868
	869	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
	870	// cout<<"cost= "<<sum<<endl;
[408]	871
[410]	872	if (sum < minSum)
	873	{
[408]	874	minSum = sum;
	875	position = (*ci).value;
[410]	876
[408]	877	raysBack = rl;
	878	raysFront = rr;
[410]	879
[408]	880	pvsBack = pl;
	881	pvsFront = pr;
	882
[410]	883	}
	884	}
	885	}
	886
[419]	887	float oldCost = (float)leaf->GetPvsSize();
[422]	888	float newCost = mCtDivCi + minSum / sizeBox;
[410]	889	float ratio = newCost / oldCost;
[408]	890
[422]	891	//Debug << "costRatio=" << ratio << " pos=" << position << " t=" << (position - minBox) / (maxBox - minBox)
	892	// <<"\t q=(" << queriesBack << "," << queriesFront << ")\t r=(" << raysBack << "," << raysFront << ")" << endl;
[410]	893
[408]	894	return ratio;
	895	}
	896
[410]	897	void VspKdTree::SortSplitCandidates(VspKdTreeLeaf *node,
	898	const int axis)
[408]	899	{
[422]	900	mSplitCandidates->clear();
[408]	901
[420]	902	int requestedSize = 2 * (int)(node->mRays.size());
[410]	903	// creates a sorted split candidates array
[422]	904	if (mSplitCandidates->capacity() > 500000 &&
	905	requestedSize < (int)(mSplitCandidates->capacity()/10) )
[410]	906	{
[449]	907	DEL_PTR(mSplitCandidates);
[422]	908	mSplitCandidates = new vector<SortableEntry>;
[410]	909	}
[408]	910
[422]	911	mSplitCandidates->reserve(requestedSize);
[408]	912
[410]	913	// insert all queries
[420]	914	for(RayInfoContainer::const_iterator ri = node->mRays.begin();
	915	ri < node->mRays.end(); ++ ri)
[410]	916	{
	917	bool positive = (*ri).mRay->HasPosDir(axis);
[422]	918	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
	919	(ri).ExtrapOrigin(axis), (void )&*ri));
[408]	920
[422]	921	mSplitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
	922	(ri).ExtrapTermination(axis), (void )&*ri));
[410]	923	}
[408]	924
[422]	925	stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
[406]	926	}
	927
[408]	928
[410]	929	void VspKdTree::EvaluateLeafStats(const TraversalData &data)
[406]	930	{
[410]	931	// the node became a leaf -> evaluate stats for leafs
[419]	932	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(data.mNode);
[408]	933
[426]	934	if (leaf->GetPvsSize() > mStat.maxPvsSize)
	935	mStat.maxPvsSize = leaf->GetPvsSize();
	936
	937	if ((int)(leaf->mRays.size()) > mStat.maxRayRefs)
	938	mStat.maxRayRefs = (int)leaf->mRays.size();
	939
	940
[421]	941	if (data.mDepth >= mTermMaxDepth)
[413]	942	++ mStat.maxDepthNodes;
[406]	943
[421]	944	if (leaf->GetPvsSize() < mTermMinPvs)
[413]	945	++ mStat.minPvsNodes;
[408]	946
[426]	947	if ((int)leaf->GetRays().size() < mTermMinRays)
[413]	948	++ mStat.minRaysNodes;
[408]	949
[426]	950	if (leaf->GetAvgRayContribution() > mTermMaxRayContribution)
[413]	951	++ mStat.maxRayContribNodes;
[408]	952
[421]	953	if (SqrMagnitude(data.mBox.Size()) <= mTermMinSize)
[413]	954	++ mStat.minSizeNodes;
[408]	955	}
	956
	957
[421]	958	inline bool VspKdTree::TerminationCriteriaMet(const VspKdTreeLeaf *leaf,
	959	const AxisAlignedBox3 &box) const
	960	{
[423]	961	return ((leaf->GetPvsSize() < mTermMinPvs) \|\|
[426]	962	(leaf->mRays.size() < mTermMinRays) \|\|
	963	(leaf->GetAvgRayContribution() > mTermMaxRayContribution ) \|\|
[423]	964	(leaf->mDepth >= mTermMaxDepth) \|\|
	965	(SqrMagnitude(box.Size()) <= mTermMinSize));
[421]	966	}
[408]	967
[421]	968
[410]	969	VspKdTreeNode VspKdTree::SubdivideNode(VspKdTreeLeaf leaf,
	970	const AxisAlignedBox3 &box,
	971	AxisAlignedBox3 &backBBox,
	972	AxisAlignedBox3 &frontBBox)
[408]	973	{
[421]	974	if (TerminationCriteriaMet(leaf, box))
[410]	975	{
[428]	976	if (1 && leaf->mDepth >= mTermMaxDepth)
[422]	977	{
[428]	978	Debug << "Warning: max depth reached depth=" << (int)leaf->mDepth<<" rays=" << (int)leaf->mRays.size() << endl;
	979	Debug << "Bbox: " << GetBBox(leaf) << endl;
[408]	980	}
[428]	981	//Debug << "depth: " << (int)leaf->mDepth << " pvs: " << leaf->GetPvsSize() << " rays: " << leaf->mRays.size() << endl;
	982
[408]	983	return leaf;
	984	}
	985
	986	float position;
	987	// first count ray sides
[410]	988	int raysBack;
	989	int raysFront;
[408]	990	int pvsBack;
	991	int pvsFront;
	992
[410]	993	// select subdivision axis
[425]	994	const int axis = SelectPlane(leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
[406]	995
[426]	996	//Debug << "rays back=" << raysBack << " rays front=" << raysFront << " pvs back=" << pvsBack << " pvs front=" << pvsFront << endl;
[408]	997
[410]	998	if (axis == -1)
	999	{
	1000	return leaf;
	1001	}
[406]	1002
[413]	1003	mStat.nodes += 2;
	1004	//++ mStat.splits[axis];
[408]	1005
[410]	1006	// add the new nodes to the tree
[419]	1007	VspKdTreeInterior *node = new VspKdTreeInterior(leaf->mParent);
[406]	1008
[419]	1009	node->mAxis = axis;
	1010	node->mPosition = position;
	1011	node->mBox = box;
	1012
[410]	1013	backBBox = box;
	1014	frontBBox = box;
[406]	1015
[410]	1016	VspKdTreeLeaf *back = new VspKdTreeLeaf(node, raysBack);
[408]	1017	back->SetPvsSize(pvsBack);
[410]	1018	VspKdTreeLeaf *front = new VspKdTreeLeaf(node, raysFront);
[408]	1019	front->SetPvsSize(pvsFront);
[406]	1020
[410]	1021	// replace a link from node's parent
[419]	1022	if (leaf->mParent)
	1023	leaf->mParent->ReplaceChildLink(leaf, node);
[410]	1024	// and setup child links
	1025	node->SetupChildLinks(back, front);
[408]	1026
[410]	1027	if (axis <= VspKdTreeNode::SPLIT_Z)
	1028	{
[408]	1029	backBBox.SetMax(axis, position);
[410]	1030	frontBBox.SetMin(axis, position);
[408]	1031
[420]	1032	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
	1033	ri != leaf->mRays.end(); ++ ri)
[410]	1034	{
	1035	if ((*ri).mRay->IsActive())
	1036	{
[408]	1037	// first unref ray from the former leaf
	1038	(*ri).mRay->Unref();
[406]	1039
[408]	1040	// determine the side of this ray with respect to the plane
	1041	int side = node->ComputeRayIntersection(ri, (ri).mRay->mT);
	1042
[410]	1043	if (side == 0)
	1044	{
	1045	if ((*ri).mRay->HasPosDir(axis))
	1046	{
[420]	1047	back->AddRay(RayInfo((*ri).mRay,
	1048	(*ri).mMinT,
	1049	(*ri).mRay->mT));
	1050	front->AddRay(RayInfo((*ri).mRay,
	1051	(*ri).mRay->mT,
	1052	(*ri).mMaxT));
[410]	1053	}
	1054	else
	1055	{
[420]	1056	back->AddRay(RayInfo((*ri).mRay,
	1057	(*ri).mRay->mT,
	1058	(*ri).mMaxT));
	1059	front->AddRay(RayInfo((*ri).mRay,
	1060	(*ri).mMinT,
	1061	(*ri).mRay->mT));
[408]	1062	}
[410]	1063	}
	1064	else
[408]	1065	if (side == 1)
	1066	front->AddRay(*ri);
	1067	else
	1068	back->AddRay(*ri);
[410]	1069	} else
[408]	1070	(*ri).mRay->Unref();
[410]	1071	}
	1072	}
	1073	else
	1074	{
	1075	// rays front/back
[406]	1076
[420]	1077	for (RayInfoContainer::iterator ri = leaf->mRays.begin();
	1078	ri != leaf->mRays.end(); ++ ri)
[410]	1079	{
	1080	if ((*ri).mRay->IsActive())
	1081	{
[408]	1082	// first unref ray from the former leaf
	1083	(*ri).mRay->Unref();
	1084
	1085	int side;
	1086	if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
	1087	side = 1;
	1088	else
	1089	side = -1;
	1090
	1091	if (side == 1)
	1092	front->AddRay(*ri);
	1093	else
[410]	1094	back->AddRay(*ri);
	1095	}
	1096	else
[408]	1097	(*ri).mRay->Unref();
[410]	1098	}
	1099	}
[408]	1100
[410]	1101	// update stats
[420]	1102	mStat.rayRefs -= (int)leaf->mRays.size();
[414]	1103	mStat.rayRefs += raysBack + raysFront;
[408]	1104
[420]	1105	DEL_PTR(leaf);
	1106
[410]	1107	return node;
[406]	1108	}
	1109
[410]	1110	int VspKdTree::ReleaseMemory(const int time)
	1111	{
	1112	stack<VspKdTreeNode *> tstack;
[406]	1113
[410]	1114	// find a node in the tree which subtree will be collapsed
[422]	1115	int maxAccessTime = time - mAccessTimeThreshold;
[420]	1116	int released = 0;
[406]	1117
[419]	1118	tstack.push(mRoot);
[406]	1119
[423]	1120	while (!tstack.empty())
[410]	1121	{
	1122	VspKdTreeNode *node = tstack.top();
	1123	tstack.pop();
[408]	1124
[410]	1125	if (!node->IsLeaf())
	1126	{
[419]	1127	VspKdTreeInterior in = dynamic_cast<VspKdTreeInterior >(node);
[410]	1128	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
[422]	1129	if (in->mDepth >= mMinCollapseDepth && in->mLastAccessTime <= maxAccessTime)
[410]	1130	{
	1131	released = CollapseSubtree(node, time);
	1132	break;
	1133	}
[420]	1134	if (in->GetBack()->GetAccessTime() < in->GetFront()->GetAccessTime())
[410]	1135	{
[419]	1136	tstack.push(in->GetFront());
	1137	tstack.push(in->GetBack());
[410]	1138	}
	1139	else
	1140	{
[419]	1141	tstack.push(in->GetBack());
	1142	tstack.push(in->GetFront());
[410]	1143	}
	1144	}
	1145	}
[406]	1146
[410]	1147	while (tstack.empty())
	1148	{
	1149	// could find node to collaps...
	1150	// cout<<"Could not find a node to release "<<endl;
	1151	break;
	1152	}
[408]	1153
[410]	1154	return released;
[408]	1155	}
[406]	1156
	1157
[410]	1158	VspKdTreeNode VspKdTree::SubdivideLeaf(VspKdTreeLeaf leaf,
	1159	const float sizeThreshold)
[408]	1160	{
[410]	1161	VspKdTreeNode *node = leaf;
[406]	1162
[410]	1163	AxisAlignedBox3 leafBBox = GetBBox(leaf);
[406]	1164
[408]	1165	static int pass = 0;
[410]	1166	++ pass;
[408]	1167
[410]	1168	// check if we should perform a dynamic subdivision of the leaf
[420]	1169	if (// leaf->mRays.size() > (unsigned)termMinCost &&
[421]	1170	(leaf->GetPvsSize() >= mTermMinPvs) &&
[410]	1171	(SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
	1172	{
[423]	1173	// memory check and release
[422]	1174	if (GetMemUsage() > mMaxTotalMemory)
[410]	1175	ReleaseMemory(pass);
	1176
	1177	AxisAlignedBox3 backBBox, frontBBox;
[406]	1178
[410]	1179	// subdivide the node
	1180	node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
	1181	}
	1182	return node;
[406]	1183	}
	1184
	1185
	1186
[423]	1187	void VspKdTree::UpdateRays(VssRayContainer &remove,
	1188	VssRayContainer &add)
[406]	1189	{
[410]	1190	VspKdTreeLeaf::NewMail();
[406]	1191
[410]	1192	// schedule rays for removal
	1193	for(VssRayContainer::const_iterator ri = remove.begin();
	1194	ri != remove.end(); ++ ri)
	1195	{
	1196	(*ri)->ScheduleForRemoval();
	1197	}
[406]	1198
[410]	1199	int inactive = 0;
[406]	1200
[410]	1201	for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
	1202	{
	1203	if ((*ri)->ScheduledForRemoval())
	1204	// RemoveRay(*ri, NULL, false);
	1205	// !!! BUG - with true it does not work correctly - aggreated delete
	1206	RemoveRay(*ri, NULL, true);
	1207	else
	1208	++ inactive;
	1209	}
[406]	1210
[410]	1211	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
[423]	1212	for (VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
[410]	1213	{
	1214	AddRay(*ri);
	1215	}
[406]	1216	}
	1217
	1218
[410]	1219	void VspKdTree::RemoveRay(VssRay *ray,
	1220	vector<VspKdTreeLeaf > affectedLeaves,
	1221	const bool removeAllScheduledRays)
[406]	1222	{
[410]	1223	stack<RayTraversalData> tstack;
[406]	1224
[420]	1225	tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[406]	1226
[410]	1227	RayTraversalData data;
[406]	1228
[410]	1229	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
	1230	while (!tstack.empty())
	1231	{
	1232	data = tstack.top();
	1233	tstack.pop();
[406]	1234
[419]	1235	if (!data.mNode->IsLeaf())
[410]	1236	{
	1237	// split the set of rays in two groups intersecting the
	1238	// two subtrees
	1239	TraverseInternalNode(data, tstack);
	1240	}
	1241	else
	1242	{
	1243	// remove the ray from the leaf
	1244	// find the ray in the leaf and swap it with the last ray...
[419]	1245	VspKdTreeLeaf leaf = (VspKdTreeLeaf )data.mNode;
[408]	1246
[410]	1247	if (!leaf->Mailed())
	1248	{
	1249	leaf->Mail();
	1250
	1251	if (affectedLeaves)
	1252	affectedLeaves->push_back(leaf);
[406]	1253
[410]	1254	if (removeAllScheduledRays)
	1255	{
[420]	1256	int tail = (int)leaf->mRays.size() - 1;
[406]	1257
[420]	1258	for (int i=0; i < (int)(leaf->mRays.size()); ++ i)
[410]	1259	{
[420]	1260	if (leaf->mRays[i].mRay->ScheduledForRemoval())
[410]	1261	{
	1262	// find a ray to replace it with
[420]	1263	while (tail >= i && leaf->mRays[tail].mRay->ScheduledForRemoval())
[410]	1264	{
[414]	1265	++ mStat.removedRayRefs;
[420]	1266	leaf->mRays[tail].mRay->Unref();
	1267	leaf->mRays.pop_back();
[410]	1268
	1269	-- tail;
	1270	}
	1271
	1272	if (tail < i)
	1273	break;
[408]	1274
[414]	1275	++ mStat.removedRayRefs;
[410]	1276
[420]	1277	leaf->mRays[i].mRay->Unref();
	1278	leaf->mRays[i] = leaf->mRays[tail];
	1279	leaf->mRays.pop_back();
[410]	1280	-- tail;
	1281	}
	1282	}
	1283	}
	1284	}
[406]	1285
[410]	1286	if (!removeAllScheduledRays)
[420]	1287	for (int i=0; i < (int)leaf->mRays.size(); i++)
[410]	1288	{
[420]	1289	if (leaf->mRays[i].mRay == ray)
[410]	1290	{
[414]	1291	++ mStat.removedRayRefs;
[410]	1292	ray->Unref();
[420]	1293	leaf->mRays[i] = leaf->mRays[leaf->mRays.size() - 1];
	1294	leaf->mRays.pop_back();
[410]	1295	// check this ray again
	1296	break;
	1297	}
	1298	}
	1299	}
[408]	1300	}
	1301
[410]	1302	if (ray->RefCount() != 0)
	1303	{
	1304	cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
	1305	exit(1);
	1306	}
	1307
[406]	1308	}
	1309
[410]	1310	void VspKdTree::AddRay(VssRay *ray)
[406]	1311	{
[410]	1312	stack<RayTraversalData> tstack;
[406]	1313
[420]	1314	tstack.push(RayTraversalData(mRoot, RayInfo(ray)));
[406]	1315
[410]	1316	RayTraversalData data;
[406]	1317
[410]	1318	while (!tstack.empty())
	1319	{
	1320	data = tstack.top();
	1321	tstack.pop();
[406]	1322
[419]	1323	if (!data.mNode->IsLeaf())
[410]	1324	{
	1325	TraverseInternalNode(data, tstack);
	1326	}
	1327	else
	1328	{
	1329	// remove the ray from the leaf
[420]	1330	// find the ray in the leaf and swap it with the last ray
	1331	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(data.mNode);
[423]	1332
[419]	1333	leaf->AddRay(data.mRayData);
[413]	1334	++ mStat.addedRayRefs;
[410]	1335	}
	1336	}
[406]	1337	}
	1338
[410]	1339	void VspKdTree::TraverseInternalNode(RayTraversalData &data,
	1340	stack<RayTraversalData> &tstack)
[406]	1341	{
[448]	1342	VspKdTreeInterior in = dynamic_cast<VspKdTreeInterior >(data.mNode);
[408]	1343
[419]	1344	if (in->mAxis <= VspKdTreeNode::SPLIT_Z)
[410]	1345	{
	1346	// determine the side of this ray with respect to the plane
[419]	1347	int side = in->ComputeRayIntersection(data.mRayData, data.mRayData.mRay->mT);
[408]	1348
[433]	1349	if (side == 0)
[410]	1350	{
[433]	1351	if (data.mRayData.mRay->HasPosDir(in->mAxis))
	1352	{
	1353	tstack.push(RayTraversalData(in->GetBack(),
	1354	RayInfo(data.mRayData.mRay, data.mRayData.mMinT, data.mRayData.mRay->mT)));
[408]	1355
[433]	1356	tstack.push(RayTraversalData(in->GetFront(),
	1357	RayInfo(data.mRayData.mRay, data.mRayData.mRay->mT, data.mRayData.mMaxT)));
[408]	1358
[433]	1359	}
	1360	else
	1361	{
	1362	tstack.push(RayTraversalData(in->GetBack(),
	1363	RayInfo(data.mRayData.mRay,
	1364	data.mRayData.mRay->mT,
	1365	data.mRayData.mMaxT)));
	1366	tstack.push(RayTraversalData(in->GetFront(),
	1367	RayInfo(data.mRayData.mRay,
	1368	data.mRayData.mMinT,
	1369	data.mRayData.mRay->mT)));
	1370	}
[410]	1371	}
[433]	1372	else
	1373	if (side == 1)
	1374	tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
	1375	else
	1376	tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[410]	1377	}
	1378	else
	1379	{
	1380	// directional split
[419]	1381	if (data.mRayData.mRay->GetDirParametrization(in->mAxis - 3) > in->mPosition)
	1382	tstack.push(RayTraversalData(in->GetFront(), data.mRayData));
[408]	1383	else
[419]	1384	tstack.push(RayTraversalData(in->GetBack(), data.mRayData));
[410]	1385	}
[406]	1386	}
	1387
[408]	1388
[410]	1389	int VspKdTree::CollapseSubtree(VspKdTreeNode *sroot, const int time)
[406]	1390	{
[410]	1391	// first count all rays in the subtree
	1392	// use mail 1 for this purpose
	1393	stack<VspKdTreeNode *> tstack;
	1394
	1395	int rayCount = 0;
	1396	int totalRayCount = 0;
	1397	int collapsedNodes = 0;
[408]	1398
	1399	#if DEBUG_COLLAPSE
[420]	1400	cout << "Collapsing subtree" << endl;
[428]	1401	cout << "accessTime=" << sroot->GetAccessTime() << endl;
[420]	1402	cout << "depth=" << (int)sroot->depth << endl;
[408]	1403	#endif
	1404
[410]	1405	// tstat.collapsedSubtrees++;
	1406	// tstat.collapseDepths += (int)sroot->depth;
	1407	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
	1408	tstack.push(sroot);
	1409	VssRay::NewMail();
[408]	1410
[410]	1411	while (!tstack.empty())
	1412	{
	1413	++ collapsedNodes;
	1414
	1415	VspKdTreeNode *node = tstack.top();
	1416	tstack.pop();
	1417	if (node->IsLeaf())
	1418	{
	1419	VspKdTreeLeaf leaf = (VspKdTreeLeaf ) node;
	1420
[420]	1421	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
	1422	ri != leaf->mRays.end(); ++ ri)
[410]	1423	{
	1424	++ totalRayCount;
[408]	1425
[410]	1426	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed())
	1427	{
[408]	1428	(*ri).mRay->Mail();
[410]	1429	++ rayCount;
[408]	1430	}
[410]	1431	}
	1432	}
	1433	else
	1434	{
[419]	1435	tstack.push(((VspKdTreeInterior *)node)->GetFront());
	1436	tstack.push(((VspKdTreeInterior *)node)->GetBack());
[410]	1437	}
	1438	}
[406]	1439
[410]	1440	VssRay::NewMail();
[406]	1441
[410]	1442	// create a new node that will hold the rays
[419]	1443	VspKdTreeLeaf *newLeaf = new VspKdTreeLeaf(sroot->mParent, rayCount);
[410]	1444
[419]	1445	if (newLeaf->mParent)
	1446	newLeaf->mParent->ReplaceChildLink(sroot, newLeaf);
[406]	1447
[410]	1448	tstack.push(sroot);
[406]	1449
[410]	1450	while (!tstack.empty())
	1451	{
	1452	VspKdTreeNode *node = tstack.top();
	1453	tstack.pop();
[408]	1454
[410]	1455	if (node->IsLeaf())
	1456	{
[419]	1457	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(node);
[408]	1458
[420]	1459	for(RayInfoContainer::iterator ri = leaf->mRays.begin();
	1460	ri != leaf->mRays.end(); ++ ri)
[410]	1461	{
	1462	// unref this ray from the old node
	1463	if ((*ri).mRay->IsActive())
	1464	{
[408]	1465	(*ri).mRay->Unref();
[410]	1466	if (!(*ri).mRay->Mailed())
	1467	{
[408]	1468	(*ri).mRay->Mail();
	1469	newLeaf->AddRay(*ri);
	1470	}
[410]	1471	}
	1472	else
[408]	1473	(*ri).mRay->Unref();
[410]	1474	}
	1475	}
	1476	else
	1477	{
[419]	1478	VspKdTreeInterior *interior =
	1479	dynamic_cast<VspKdTreeInterior *>(node);
	1480	tstack.push(interior->GetBack());
	1481	tstack.push(interior->GetFront());
[410]	1482	}
	1483	}
[406]	1484
[410]	1485	// delete the node and all its children
[420]	1486	DEL_PTR(sroot);
[408]	1487
[420]	1488	// for(VspKdTreeNode::SRayContainer::iterator ri = newleaf->mRays.begin();
	1489	// ri != newleaf->mRays.end(); ++ ri)
[410]	1490	// (*ri).ray->UnMail(2);
[408]	1491
	1492	#if DEBUG_COLLAPSE
[423]	1493	cout << "Total memory before=" << GetMemUsage() << endl;
[408]	1494	#endif
	1495
[414]	1496	mStat.nodes -= collapsedNodes - 1;
	1497	mStat.rayRefs -= totalRayCount - rayCount;
[408]	1498
	1499	#if DEBUG_COLLAPSE
[410]	1500	cout << "collapsed nodes" << collapsedNodes << endl;
	1501	cout << "collapsed rays" << totalRayCount - rayCount << endl;
	1502	cout << "Total memory after=" << GetMemUsage() << endl;
	1503	cout << "================================" << endl;
[408]	1504	#endif
	1505
	1506	// tstat.collapsedNodes += collapsedNodes;
	1507	// tstat.collapsedRays += totalRayCount - rayCount;
[410]	1508	return totalRayCount - rayCount;
[406]	1509	}
	1510
[408]	1511
[410]	1512	int VspKdTree::GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const
[406]	1513	{
[408]	1514	stack<VspKdTreeNode *> tstack;
[419]	1515	tstack.push(mRoot);
[408]	1516
	1517	Intersectable::NewMail();
	1518	int pvsSize = 0;
	1519
[410]	1520	while (!tstack.empty())
	1521	{
	1522	VspKdTreeNode *node = tstack.top();
	1523	tstack.pop();
[408]	1524
[410]	1525	if (node->IsLeaf())
	1526	{
[408]	1527	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
[420]	1528	for (RayInfoContainer::iterator ri = leaf->mRays.begin();
	1529	ri != leaf->mRays.end(); ++ ri)
[410]	1530	{
	1531	if ((*ri).mRay->IsActive())
	1532	{
[408]	1533	Intersectable *object;
	1534	#if BIDIRECTIONAL_RAY
	1535	object = (*ri).mRay->mOriginObject;
[410]	1536
	1537	if (object && !object->Mailed())
	1538	{
	1539	++ pvsSize;
[408]	1540	object->Mail();
	1541	}
	1542	#endif
	1543	object = (*ri).mRay->mTerminationObject;
[410]	1544	if (object && !object->Mailed())
	1545	{
	1546	++ pvsSize;
[408]	1547	object->Mail();
	1548	}
	1549	}
[410]	1550	}
	1551	}
	1552	else
	1553	{
[419]	1554	VspKdTreeInterior in = dynamic_cast<VspKdTreeInterior >(node);
[410]	1555
[419]	1556	if (in->mAxis < 3)
[410]	1557	{
[419]	1558	if (box.Max(in->mAxis) >= in->mPosition)
	1559	tstack.push(in->GetFront());
[408]	1560
[419]	1561	if (box.Min(in->mAxis) <= in->mPosition)
	1562	tstack.push(in->GetBack());
[410]	1563	}
	1564	else
	1565	{
[408]	1566	// both nodes for directional splits
[419]	1567	tstack.push(in->GetFront());
	1568	tstack.push(in->GetBack());
[408]	1569	}
	1570	}
	1571	}
[410]	1572
[408]	1573	return pvsSize;
[406]	1574	}
	1575
[410]	1576	void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
	1577	float &maxRayContribution,
	1578	float &avgRayContribution)
[406]	1579	{
[408]	1580	stack<VspKdTreeNode *> tstack;
[419]	1581	tstack.push(mRoot);
[408]	1582
	1583	minRayContribution = 1.0f;
	1584	maxRayContribution = 0.0f;
[410]	1585
[408]	1586	float sumRayContribution = 0.0f;
	1587	int leaves = 0;
[406]	1588
[410]	1589	while (!tstack.empty())
	1590	{
	1591	VspKdTreeNode *node = tstack.top();
	1592	tstack.pop();
	1593	if (node->IsLeaf())
	1594	{
[408]	1595	leaves++;
	1596	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
	1597	float c = leaf->GetAvgRayContribution();
[410]	1598
[408]	1599	if (c > maxRayContribution)
	1600	maxRayContribution = c;
	1601	if (c < minRayContribution)
	1602	minRayContribution = c;
	1603	sumRayContribution += c;
	1604
[410]	1605	}
	1606	else {
[408]	1607	VspKdTreeInterior in = (VspKdTreeInterior )node;
	1608	// both nodes for directional splits
[419]	1609	tstack.push(in->GetFront());
	1610	tstack.push(in->GetBack());
[408]	1611	}
	1612	}
	1613
[423]	1614	Debug << "sum=" << sumRayContribution << endl;
	1615	Debug << "leaves=" << leaves << endl;
	1616	avgRayContribution = sumRayContribution / (float)leaves;
[406]	1617	}
	1618
	1619
[410]	1620	int VspKdTree::GenerateRays(const float ratioPerLeaf,
	1621	SimpleRayContainer &rays)
[406]	1622	{
[408]	1623	stack<VspKdTreeNode *> tstack;
[419]	1624	tstack.push(mRoot);
[408]	1625
[410]	1626	while (!tstack.empty())
	1627	{
	1628	VspKdTreeNode *node = tstack.top();
	1629	tstack.pop();
[408]	1630
[410]	1631	if (node->IsLeaf())
	1632	{
[423]	1633	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(node);
	1634
[408]	1635	float c = leaf->GetAvgRayContribution();
[411]	1636	int num = (int)(c*ratioPerLeaf + 0.5);
[408]	1637	// cout<<num<<" ";
	1638
[410]	1639	for (int i=0; i < num; i++)
	1640	{
[408]	1641	Vector3 origin = GetBBox(leaf).GetRandomPoint();
[419]	1642	/*Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
[408]	1643	Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
	1644	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
[419]	1645	rays.push_back(SimpleRay(origin, direction));*/
[408]	1646	}
	1647
[410]	1648	}
	1649	else
	1650	{
[419]	1651	VspKdTreeInterior *in =
	1652	dynamic_cast<VspKdTreeInterior *>(node);
[408]	1653	// both nodes for directional splits
[419]	1654	tstack.push(in->GetFront());
	1655	tstack.push(in->GetBack());
[408]	1656	}
[406]	1657	}
	1658
[411]	1659	return (int)rays.size();
[406]	1660	}
	1661
	1662
[410]	1663	float VspKdTree::GetAvgPvsSize()
[406]	1664	{
[408]	1665	stack<VspKdTreeNode *> tstack;
[419]	1666	tstack.push(mRoot);
[406]	1667
[408]	1668	int sumPvs = 0;
	1669	int leaves = 0;
[410]	1670
	1671	while (!tstack.empty())
	1672	{
	1673	VspKdTreeNode *node = tstack.top();
	1674	tstack.pop();
[408]	1675
[410]	1676	if (node->IsLeaf())
	1677	{
[408]	1678	VspKdTreeLeaf leaf = (VspKdTreeLeaf )node;
	1679	// update pvs size
	1680	leaf->UpdatePvsSize();
	1681	sumPvs += leaf->GetPvsSize();
	1682	leaves++;
[410]	1683	}
	1684	else
	1685	{
[408]	1686	VspKdTreeInterior in = (VspKdTreeInterior )node;
	1687	// both nodes for directional splits
[419]	1688	tstack.push(in->GetFront());
	1689	tstack.push(in->GetBack());
[406]	1690	}
	1691	}
[408]	1692
[410]	1693	return sumPvs / (float)leaves;
[418]	1694	}
	1695
	1696	VspKdTreeNode *VspKdTree::GetRoot() const
	1697	{
	1698	return mRoot;
	1699	}
	1700
[419]	1701	AxisAlignedBox3 VspKdTree::GetBBox(VspKdTreeNode *node) const
[418]	1702	{
[419]	1703	if (node->mParent == NULL)
	1704	return mBox;
[418]	1705
	1706	if (!node->IsLeaf())
[419]	1707	return (dynamic_cast<VspKdTreeInterior *>(node))->mBox;
[418]	1708
[419]	1709	if (node->mParent->mAxis >= 3)
	1710	return node->mParent->mBox;
[418]	1711
[419]	1712	AxisAlignedBox3 box(node->mParent->mBox);
[448]	1713	if (node->mParent->GetFront() == node)
[419]	1714	box.SetMin(node->mParent->mAxis, node->mParent->mPosition);
[418]	1715	else
[419]	1716	box.SetMax(node->mParent->mAxis, node->mParent->mPosition);
	1717
	1718	return box;
[418]	1719	}
	1720
	1721	int VspKdTree::GetRootPvsSize() const
	1722	{
	1723	return GetPvsSize(mRoot, mBox);
	1724	}
[419]	1725
	1726	const VspKdStatistics &VspKdTree::GetStatistics() const
	1727	{
	1728	return mStat;
[420]	1729	}
	1730
	1731	void VspKdTree::AddRays(VssRayContainer &add)
	1732	{
	1733	VssRayContainer remove;
	1734	UpdateRays(remove, add);
	1735	}
	1736
	1737	// return memory usage in MB
	1738	float VspKdTree::GetMemUsage() const
	1739	{
	1740	return
	1741	(sizeof(VspKdTree) +
	1742	mStat.Leaves() * sizeof(VspKdTreeLeaf) +
	1743	mStat.Interior() * sizeof(VspKdTreeInterior) +
	1744	mStat.rayRefs * sizeof(RayInfo)) / (1024.0f * 1024.0f);
	1745	}
	1746
	1747	float VspKdTree::GetRayMemUsage() const
	1748	{
	1749	return mStat.rays * (sizeof(VssRay))/(1024.0f * 1024.0f);
[425]	1750	}
	1751
[426]	1752	int VspKdTree::ComputePvsSize(VspKdTreeNode *node,
	1753	const RayInfoContainer &globalRays) const
[425]	1754	{
	1755	int pvsSize = 0;
	1756
	1757	const AxisAlignedBox3 box = GetBBox(node);
	1758
	1759	RayInfoContainer::const_iterator it, it_end = globalRays.end();
	1760
[426]	1761	// warning: implicit conversion from VssRay to Ray
[425]	1762	for (it = globalRays.begin(); it != globalRays.end(); ++ it)
[426]	1763	pvsSize += box.GetMinMaxT((it).mRay, NULL, NULL);
	1764
[425]	1765	return pvsSize;
[428]	1766	}
	1767
	1768	void VspKdTree::CollectLeaves(vector<VspKdTreeLeaf *> &leaves) const
	1769	{
	1770	stack<VspKdTreeNode *> nodeStack;
	1771	nodeStack.push(mRoot);
	1772
	1773	while (!nodeStack.empty())
	1774	{
	1775	VspKdTreeNode *node = nodeStack.top();
	1776
	1777	nodeStack.pop();
	1778
	1779	if (node->IsLeaf())
	1780	{
	1781	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(node);
	1782	leaves.push_back(leaf);
	1783	}
	1784	else
	1785	{
	1786	VspKdTreeInterior interior = dynamic_cast<VspKdTreeInterior >(node);
	1787
[448]	1788	nodeStack.push(interior->GetBack());
	1789	nodeStack.push(interior->GetFront());
[428]	1790	}
	1791	}
[452]	1792	}
	1793
	1794	int VspKdTree::FindNeighbors(VspKdTreeLeaf *n,
	1795	vector<VspKdTreeLeaf *> &neighbors,
	1796	bool onlyUnmailed)
	1797	{
	1798	stack<VspKdTreeNode *> nodeStack;
	1799
	1800	nodeStack.push(mRoot);
	1801
	1802	AxisAlignedBox3 box = GetBBox(n);
	1803
	1804	while (!nodeStack.empty())
	1805	{
	1806	VspKdTreeNode *node = nodeStack.top();
	1807	nodeStack.pop();
	1808
	1809	if (node->IsLeaf())
	1810	{
	1811	VspKdTreeLeaf leaf = dynamic_cast<VspKdTreeLeaf >(node);
	1812
	1813	if (leaf != n && (!onlyUnmailed \|\| !leaf->Mailed()))
	1814	neighbors.push_back(leaf);
	1815	}
	1816	else
	1817	{
	1818	VspKdTreeInterior interior = dynamic_cast<VspKdTreeInterior >(node);
	1819
	1820	if (interior->mPosition > box.Max(interior->mAxis))
	1821	nodeStack.push(interior->mBack);
	1822	else
	1823	{
	1824	if (interior->mPosition < box.Min(interior->mAxis))
	1825	nodeStack.push(interior->mFront);
	1826	else
	1827	{
	1828	// random decision
	1829	nodeStack.push(interior->mBack);
	1830	nodeStack.push(interior->mFront);
	1831	}
	1832	}
	1833	}
	1834	}
	1835
	1836	return (int)neighbors.size();
	1837	}
	1838
	1839	void VspKdTree::MergeLeaves()
	1840	{
	1841	vector<VspKdTreeLeaf *> leaves;
	1842	priority_queue<LeafPair> mergeCandidates;
	1843	vector<VspKdTreeLeaf *> neighbors;
	1844
	1845	CollectLeaves(leaves);
	1846
	1847	VspKdTreeLeaf::NewMail();
	1848
	1849	vector<VspKdTreeLeaf *>::const_iterator it, it_end = leaves.end();
	1850
	1851
	1852	for (it = leaves.begin(); it != it_end; ++ it)
	1853	{
	1854	VspKdTreeLeaf leaf = it;
	1855
	1856	if (!leaf->Mailed())
	1857	{
	1858	FindNeighbors(leaf, neighbors, true);
	1859
	1860	vector<VspKdTreeLeaf *>::const_iterator nit, nit_end = neighbors.end();
	1861
	1862	for (nit = neighbors.begin(); nit != neighbors.end(); ++ nit)
	1863	mergeCandidates.push(LeafPair(leaf, *nit));
	1864
	1865	leaf->Mail();
	1866	}
	1867	}
	1868	}
	1869
	1870
	1871
	1872	/*********************************************************************/
	1873	/* MergeCandidate implementation */
	1874	/*********************************************************************/
	1875
	1876
	1877	MergeCandidate::MergeCandidate(VspKdTreeLeaf l1, VspKdTreeLeaf l2):
	1878	mLeaf1(l1), mLeaf2(l2)
	1879	{}
	1880
[453]	1881	float MergeCandidate::EvaluateMergeCost() const
[452]	1882	{
[453]	1883	// pvs difference
	1884	VspKdViewCell vc1 = dynamic_cast<VspKdViewCell >(mLeaf1->GetViewCell());
	1885	VspKdViewCell vc2 = dynamic_cast<VspKdViewCell >(mLeaf2->GetViewCell());
[452]	1886
[453]	1887	const int diff1 = vc1->GetPvs().Diff(vc2->GetPvs());
	1888
	1889	//const int diff2 = vc2->GetPvs().Diff(vc1->GetPvs());
	1890	//const float simDiff = (float)max(diff1, diff2);
	1891
	1892	const float oldCost =
	1893	vc1->GetSize() * vc1->GetPvs().GetSize() +
	1894	vc2->GetSize() * vc2->GetPvs().GetSize();
	1895
	1896	const float mergedPvsCost =
	1897	(diff1 + vc1->GetPvs().GetSize()) *
	1898	(vc1->GetSize() + vc2->GetSize());
	1899
	1900	return mergedPvsCost / oldCost;
[419]	1901	}

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

Download in other formats: