Context Navigation

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

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

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

Download in other formats: