Context Navigation

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

Revision 667, 46.7 KB checked in by mattausch, 19 years ago (diff)
test version build script for testing. code is set uo for testing also

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

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

Download in other formats: