Context Navigation

VssTree.cpp @ 2575

Revision 2575, 47.3 KB checked in by bittner, 17 years ago (diff)
big merge: preparation for havran ray caster, check if everything works

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

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

Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/VssTree.cpp @ 2575

Download in other formats: