Context Navigation

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

Revision 1883, 47.3 KB checked in by bittner, 18 years ago (diff)
mixture distribution initial coding

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

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

Download in other formats: