Context Navigation

source: GTP/trunk/Lib/Vis/Preprocessing/src/RssTree.cpp @ 1548

Revision 1548, 69.9 KB checked in by mattausch, 18 years ago (diff)

Rev	Line
[447]	1	// ================================================================
	2	// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
	3	// ****************************************************************
	4	/**
	5	The KD tree based LSDS
	6	*/
	7	// Initial coding by
	8	/**
	9	@author Jiri Bittner
	10	*/
	11
	12	// Standard headers
	13	#include <stack>
	14	#include <queue>
	15	#include <algorithm>
	16	#include <fstream>
	17	#include <string>
	18
	19	#include "RssTree.h"
	20
	21	#include "Environment.h"
	22	#include "VssRay.h"
	23	#include "Intersectable.h"
	24	#include "Ray.h"
	25	#include "Containers.h"
[492]	26	#include "ViewCell.h"
	27	#include "Exporter.h"
	28	#include "Preprocessor.h"
	29	#include "SceneGraph.h"
[447]	30
[860]	31
[863]	32	namespace GtpVisibilityPreprocessor {
[860]	33
	34
[447]	35	#define DEBUG_SPLIT_COST 0
	36	#define DEBUG_SPLITS 0
	37
[1199]	38	#define EVAL_VIEWCELLS 0
	39
[447]	40	// Static variables
	41	int
	42	RssTreeLeaf::mailID = 0;
	43
	44	inline void
	45	AddObject2Pvs(Intersectable *object,
	46	const int side,
	47	int &pvsBack,
	48	int &pvsFront)
	49	{
	50
	51	if (!object)
	52	return;
	53
	54	if (side <= 0) {
	55	if (!object->Mailed() && !object->Mailed(2)) {
	56	pvsBack++;
	57	if (object->Mailed(1))
	58	object->Mail(2);
	59	else
	60	object->Mail();
	61	}
	62	}
	63
	64	if (side >= 0) {
	65	if (!object->Mailed(1) && !object->Mailed(2)) {
	66	pvsFront++;
	67	if (object->Mailed())
	68	object->Mail(2);
	69	else
	70	object->Mail(1);
	71	}
	72	}
	73	}
	74
[492]	75	inline void
	76	AddViewcells2Pvs(const ViewCellContainer &viewcells,
	77	const int side,
	78	int &viewcellsBack,
	79	int &viewcellsFront)
	80	{
	81	ViewCellContainer::const_iterator it = viewcells.begin();
	82
	83	for (; it != viewcells.end(); ++it) {
	84	ViewCell viewcell = it;
	85	if (side <= 0) {
	86	if (!viewcell->Mailed() && !viewcell->Mailed(2)) {
	87	viewcellsBack++;
	88	if (viewcell->Mailed(1))
	89	viewcell->Mail(2);
	90	else
	91	viewcell->Mail();
	92	}
	93	}
	94
	95	if (side >= 0) {
	96	if (!viewcell->Mailed(1) && !viewcell->Mailed(2)) {
	97	viewcellsFront++;
	98	if (viewcell->Mailed())
	99	viewcell->Mail(2);
	100	else
	101	viewcell->Mail(1);
	102	}
	103	}
	104	}
	105	}
	106
	107
[447]	108	// Constructor
	109	RssTree::RssTree()
	110	{
[1004]	111	Environment::GetSingleton()->GetIntValue("RssTree.maxDepth", termMaxDepth);
	112	Environment::GetSingleton()->GetIntValue("RssTree.minPvs", termMinPvs);
	113	Environment::GetSingleton()->GetIntValue("RssTree.minRays", termMinRays);
	114	Environment::GetSingleton()->GetFloatValue("RssTree.maxRayContribution", termMaxRayContribution);
	115	Environment::GetSingleton()->GetFloatValue("RssTree.maxCostRatio", termMaxCostRatio);
[447]	116
[1004]	117	Environment::GetSingleton()->GetFloatValue("RssTree.minSize", termMinSize);
[447]	118	termMinSize = sqr(termMinSize);
	119
[1004]	120	Environment::GetSingleton()->GetFloatValue("RssTree.refDirBoxMaxSize", refDirBoxMaxSize);
[447]	121	refDirBoxMaxSize = sqr(refDirBoxMaxSize);
	122
[1004]	123	Environment::GetSingleton()->GetFloatValue("RssTree.epsilon", epsilon);
	124	Environment::GetSingleton()->GetFloatValue("RssTree.ct_div_ci", ct_div_ci);
[447]	125
[1004]	126	Environment::GetSingleton()->GetFloatValue("RssTree.maxTotalMemory", maxTotalMemory);
	127	Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
[447]	128
[1004]	129	Environment::GetSingleton()->GetFloatValue("RssTree.maxStaticMemory", maxStaticMemory);
[447]	130
	131
	132
[1004]	133	Environment::GetSingleton()->GetIntValue("RssTree.accessTimeThreshold", accessTimeThreshold);
[447]	134	//= 1000;
[1004]	135	Environment::GetSingleton()->GetIntValue("RssTree.minCollapseDepth", minCollapseDepth);
[447]	136	// int minCollapseDepth = 4;
	137
	138	// pRefDirThresh = cos(0.5M_PI - M_PIrefDirAngle/180.0);
	139	// cosRefDir = cos(M_PI*refDirAngle/180.0);
	140	// sinRefDir = sin(M_PI*refDirAngle/180.0);
	141
	142
	143	// split type
	144	char sname[128];
[1004]	145	Environment::GetSingleton()->GetStringValue("RssTree.splitType", sname);
[447]	146	string name(sname);
	147
	148	if (name.compare("regular") == 0)
	149	splitType = ESplitRegular;
	150	else
	151	if (name.compare("heuristic") == 0)
	152	splitType = ESplitHeuristic;
	153	else
	154	if (name.compare("hybrid") == 0)
	155	splitType = ESplitHybrid;
	156	else {
	157	cerr<<"Invalid RssTree split type "<<name<<endl;
	158	exit(1);
	159	}
[1199]	160
	161	Environment::GetSingleton()->GetIntValue("RssTree.hybridDepth", mHybridDepth);
	162
[1004]	163	Environment::GetSingleton()->GetBoolValue("RssTree.randomize", randomize);
	164	Environment::GetSingleton()->GetBoolValue("RssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
[447]	165
[1004]	166	Environment::GetSingleton()->GetBoolValue("RssTree.interleaveDirSplits", mInterleaveDirSplits);
	167	Environment::GetSingleton()->GetIntValue("RssTree.dirSplitDepth", mDirSplitDepth);
[447]	168
[1004]	169	Environment::GetSingleton()->GetBoolValue("RssTree.importanceBasedCost", mImportanceBasedCost);
[492]	170
[1004]	171	Environment::GetSingleton()->GetIntValue("RssTree.maxRays", mMaxRays);
[492]	172
[1004]	173	Environment::GetSingleton()->GetBoolValue("RssTree.perObjectTree", mPerObjectTree);
[447]	174
[516]	175	// mRoots;
	176
[447]	177	splitCandidates = new vector<SortableEntry>;
	178	}
	179
	180
	181	RssTree::~RssTree()
	182	{
[516]	183	for (int i=0; i < mRoots.size(); i++)
	184	if (mRoots[i])
	185	delete mRoots[i];
[447]	186	}
	187
	188
	189
	190
	191	void
	192	RssStatistics::Print(ostream &app) const
	193	{
[492]	194	app << "###### RssTree statistics ######\n";
[447]	195
	196	app << "#N_RAYS ( Number of rays )\n"
	197	<< rays <<endl;
	198
	199	app << "#N_INITPVS ( Initial PVS size )\n"
	200	<< initialPvsSize <<endl;
	201
	202	app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
	203
	204	app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
	205
[667]	206	app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
[447]	207	for (int i=0; i<7; i++)
	208	app << splits[i] <<" ";
	209	app <<endl;
	210
	211	app << "#N_RAYREFS ( Number of rayRefs )\n" <<
	212	rayRefs << "\n";
	213
	214	app << "#N_RAYRAYREFS ( Number of rayRefs / ray )\n" <<
	215	rayRefs/(double)rays << "\n";
	216
	217	app << "#N_LEAFRAYREFS ( Number of rayRefs / leaf )\n" <<
	218	rayRefs/(double)Leaves() << "\n";
	219
	220	app << "#N_MAXRAYREFS ( Max number of rayRefs / leaf )\n" <<
	221	maxRayRefs << "\n";
	222
	223
	224	// app << setprecision(4);
	225
	226	app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
	227	maxDepthNodes*100/(double)Leaves()<<endl;
	228
[551]	229	app << "#N_PMINPVSLEAVES ( Percentage of leaves with minPvs )\n"<<
[447]	230	minPvsNodes*100/(double)Leaves()<<endl;
	231
[551]	232	app << "#N_PMINRAYSLEAVES ( Percentage of leaves with minRays )\n"<<
[447]	233	minRaysNodes*100/(double)Leaves()<<endl;
	234
	235	app << "#N_PMINSIZELEAVES ( Percentage of leaves with minSize )\n"<<
	236	minSizeNodes*100/(double)Leaves()<<endl;
	237
	238	app << "#N_PMAXRAYCONTRIBLEAVES ( Percentage of leaves with maximal ray contribution )\n"<<
	239	maxRayContribNodes*100/(double)Leaves()<<endl;
	240
	241	app << "#N_PMAXCOSTRATIOLEAVES ( Percentage of leaves with max cost ratio )\n"<<
	242	maxCostRatioNodes*100/(double)Leaves()<<endl;
	243
	244	app << "#N_ADDED_RAYREFS (Number of dynamically added ray references )\n"<<
	245	addedRayRefs<<endl;
	246
	247	app << "#N_REMOVED_RAYREFS (Number of dynamically removed ray references )\n"<<
	248	removedRayRefs<<endl;
	249
	250	// app << setprecision(4);
	251
	252	app << "#N_CTIME ( Construction time [s] )\n"
	253	<< Time() << " \n";
	254
[492]	255	app << "###### END OF RssTree statistics ######\n";
[447]	256
	257	}
	258
	259
	260	void
[492]	261	RssTree::UpdatePvsSize(RssTreeLeaf *leaf)
[447]	262	{
[492]	263	if (!leaf->mValidPvs) {
[447]	264	Intersectable::NewMail();
	265	int pvsSize = 0;
[492]	266	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	267	ri != leaf->rays.end();
[447]	268	ri++)
	269	if ((*ri).mRay->IsActive()) {
	270	Intersectable *object;
[492]	271	object = (*ri).GetObject();
[447]	272	if (object && !object->Mailed()) {
	273	pvsSize++;
	274	object->Mail();
	275	}
	276	}
[492]	277	leaf->SetPvsSize(pvsSize);
	278
	279	ComputeImportance(leaf);
[447]	280	}
	281	}
	282
	283	bool
	284	RssTree::ClipRay(
	285	RssTreeNode::RayInfo &rayInfo,
	286	const AxisAlignedBox3 &box
	287	)
	288	{
	289	float tmin, tmax;
	290	static Ray ray;
[492]	291	cerr<<"Clip not reimplented yet...Quiting\n";
	292	exit(1);
	293	ray.Init(rayInfo.GetOrigin(), rayInfo.GetDir(), Ray::LINE_SEGMENT);
	294
	295	if (!box.ComputeMinMaxT(ray, &tmin, &tmax) \|\| tmin>=tmax)
[447]	296	return false;
	297
	298	// now check if the ray origin lies inside the box
	299	if ( tmax < rayInfo.mRay->GetSize() ) {
	300	// this ray does not leave the box
	301	rayInfo.mRay->SetupEndPoints(
	302	ray.Extrap(tmax),
	303	rayInfo.mRay->mTermination
	304	);
	305	return true;
	306	}
	307
	308	return false;
	309	}
	310
	311
	312	void
	313	RssTree::Construct(
[516]	314	ObjectContainer &objects,
	315	VssRayContainer &rays
[447]	316	// forced bounding box is only used when computing from-box
	317	// visibility
[516]	318	// AxisAlignedBox3 *forcedBoundingBox
[447]	319	)
	320	{
[516]	321	cout<<"Constructing rss tree"<<endl<<flush;
[447]	322	stat.Start();
	323
	324	maxMemory = maxStaticMemory;
	325
[516]	326	// if (root)
	327	// delete root;
[492]	328
[516]	329	if (mPerObjectTree) {
	330	// get max id from the rays
	331	int i;
	332	mRoots.resize(objects.size());
	333	for (i = 0; i < objects.size(); i++) {
	334	RssTreeLeaf *leaf = new RssTreeLeaf(NULL, 0);
	335	// leaf->bbox.Initialize();
	336	leaf->dirBBox.Initialize();
	337	leaf->dirBBox.SetMin(2, 0.0f);
	338	leaf->dirBBox.SetMax(2, 1.0f);
	339	mRoots[i] = leaf;
	340	}
[1112]	341
[516]	342	// init spatial bounding boxes
	343	for (i = 0; i < objects.size(); i++) {
	344	GetRoot(objects[i])->bbox = objects[i]->GetBox();
	345	}
	346	stat.nodes = i;
[569]	347	stat.leaves = i;
[516]	348	} else {
	349	mRoots.resize(1);
	350	RssTreeLeaf *leaf = new RssTreeLeaf(NULL, rays.size());
	351	leaf->bbox.Initialize();
	352	leaf->dirBBox.Initialize();
	353	leaf->dirBBox.SetMin(2, 0.0f);
	354	leaf->dirBBox.SetMax(2, 1.0f);
	355	mRoots[0] = leaf;
	356	stat.nodes = 1;
[569]	357	stat.leaves = 1;
[516]	358	}
[447]	359
	360	for(VssRayContainer::const_iterator ri = rays.begin();
[516]	361	ri != rays.end();
	362	ri++) {
	363
	364	RssTreeNode::RayInfo info(*ri);
	365
	366	// first construct a leaf that will get subdivide
[569]	367	RssTreeLeaf leaf = (RssTreeLeaf ) GetRoot(info.GetSourceObject());
[447]	368
[516]	369	leaf->AddRay(info);
	370
	371	// leaf bbox contains bbox of origins only
[1112]	372
	373	// include both origin and terminatin in the global bbox
	374	#if USE_ORIGIN
[516]	375	leaf->bbox.Include((*ri)->GetOrigin());
	376	bbox.Include((*ri)->GetOrigin());
[1112]	377	#else
[516]	378	bbox.Include((*ri)->GetTermination());
[1112]	379	leaf->bbox.Include((*ri)->GetTermination());
	380	#endif
[516]	381	Vector3 dVec = Vector3(
	382	(*ri)->GetDirParametrization(0),
	383	(*ri)->GetDirParametrization(1),
	384	0
	385	);
	386
	387	leaf->dirBBox.Include(dVec);
	388	dirBBox.Include(dVec);
[447]	389	}
[516]	390
[492]	391	// make the z axis (unused) a unit size
	392	// important for volume computation
[516]	393
[492]	394
[516]	395	// if ( forcedBoundingBox )
	396	// bbox = *forcedBoundingBox;
	397
[447]	398	cout<<"Bbox = "<<bbox<<endl;
	399	cout<<"Dirr Bbox = "<<dirBBox<<endl;
[516]	400
	401	stat.rays = rays.size();
	402	stat.initialPvsSize = 0;
	403	for (int i=0; i < mRoots.size(); i++) {
	404	RssTreeLeaf leaf = (RssTreeLeaf )mRoots[i];
[1112]	405	UpdatePvsSize(leaf);
[516]	406	stat.initialPvsSize += leaf->GetPvsSize();
	407	mRoots[i] = Subdivide(TraversalData(leaf, GetBBox(leaf), 0));
	408	}
	409
[447]	410	if (splitCandidates) {
	411	// force realease of this vector
	412	delete splitCandidates;
	413	splitCandidates = new vector<SortableEntry>;
	414	}
	415
	416	stat.Stop();
	417	stat.Print(cout);
	418	cout<<"#Total memory="<<GetMemUsage()<<endl;
	419
[492]	420	// this rotine also updates importances etc...
[447]	421	}
	422
	423	int
	424	RssTree::UpdateSubdivision()
	425	{
	426	priority_queue<TraversalData> tStack;
	427	// stack<TraversalData> tStack;
	428
[516]	429	// tStack.push(TraversalData(root, bbox, 0));
	430	PushRoots(tStack);
	431
[447]	432	AxisAlignedBox3 backBox;
	433	AxisAlignedBox3 frontBox;
	434
	435	maxMemory = maxTotalMemory;
	436	int subdivided = 0;
	437	int lastMem = 0;
	438	while (!tStack.empty()) {
	439
	440	float mem = GetMemUsage();
	441
	442	if ( lastMem/10 != ((int)mem)/10) {
	443	cout<<mem<<" MB"<<endl;
	444	}
	445	lastMem = (int)mem;
	446
	447	if ( mem > maxMemory ) {
	448	// count statistics on unprocessed leafs
	449	while (!tStack.empty()) {
	450	// EvaluateLeafStats(tStack.top());
	451	tStack.pop();
	452	}
	453	break;
	454	}
	455
	456	TraversalData data = tStack.top();
	457	tStack.pop();
	458
	459	if (data.node->IsLeaf()) {
	460	RssTreeNode node = SubdivideNode((RssTreeLeaf ) data.node,
	461	data.bbox,
	462	backBox,
	463	frontBox
	464	);
	465	if (!node->IsLeaf()) {
	466	subdivided++;
[492]	467
[1112]	468
[447]	469	RssTreeInterior interior = (RssTreeInterior ) node;
	470	// push the children on the stack
	471	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
	472	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
	473	} else {
	474	// EvaluateLeafStats(data);
	475	}
	476	} else {
	477	RssTreeInterior interior = (RssTreeInterior ) data.node;
	478	tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
	479	tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
	480	}
	481	}
	482	return subdivided;
	483	}
	484
	485
	486	RssTreeNode *
	487	RssTree::Subdivide(const TraversalData &tdata)
	488	{
	489	RssTreeNode *result = NULL;
	490
	491	priority_queue<TraversalData> tStack;
	492	// stack<TraversalData> tStack;
	493
	494	tStack.push(tdata);
	495
	496	AxisAlignedBox3 backBox;
	497	AxisAlignedBox3 frontBox;
	498
	499
	500	int lastMem = 0;
	501	while (!tStack.empty()) {
	502
	503	float mem = GetMemUsage();
	504
	505	if ( lastMem/10 != ((int)mem)/10) {
	506	cout<<mem<<" MB"<<endl;
	507	}
	508	lastMem = (int)mem;
	509
	510	if ( mem > maxMemory ) {
	511	// count statistics on unprocessed leafs
	512	while (!tStack.empty()) {
	513	EvaluateLeafStats(tStack.top());
	514	tStack.pop();
	515	}
	516	break;
	517	}
	518
	519	TraversalData data = tStack.top();
	520	tStack.pop();
	521
	522	#if DEBUG_SPLITS
	523	Debug<<"#Splitting node"<<endl;
	524	data.node->Print(Debug);
	525	#endif
	526	RssTreeNode node = SubdivideNode((RssTreeLeaf ) data.node,
	527	data.bbox,
	528	backBox,
	529	frontBox
	530	);
	531
	532
	533	if (result == NULL)
	534	result = node;
	535
	536	if (!node->IsLeaf()) {
	537
	538	RssTreeInterior interior = (RssTreeInterior ) node;
	539	// push the children on the stack
	540	tStack.push(TraversalData(interior->back, backBox, data.depth+1));
	541	tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
	542
	543	#if DEBUG_SPLITS
	544	Debug<<"#New nodes"<<endl;
	545	interior->back->Print(Debug);
	546	interior->front->Print(Debug);
	547	Debug<<"#####################################"<<endl;
	548	#endif
	549
	550	} else {
	551	EvaluateLeafStats(data);
	552	}
	553	}
	554
	555	return result;
	556	}
	557
	558
	559	// returns selected plane for subdivision
	560	int
	561	RssTree::SelectPlane(
	562	RssTreeLeaf *leaf,
	563	const AxisAlignedBox3 &box,
[492]	564	SplitInfo &info
[447]	565	)
	566	{
	567
[492]	568	BestCostRatio(leaf,
	569	info);
	570
[447]	571	#if DEBUG_SPLIT_COST
[492]	572	Debug<<"Split Info:"<<endl;
	573	Debug<<"axis="<<info.axis<<" ratio="<<info.costRatio<<endl;
[1199]	574	Debug<<"rays="<<info.rays<<
	575	" rays back="<<info.raysBack<<
	576	" rays front="<<info.raysFront<<endl;
	577	Debug<<"c="<<info.contribution<<
	578	" c back="<<info.contributionBack<<
	579	" c front="<<info.contributionFront<<endl;
	580
	581	// Debug<<"viewcells="<<info.viewCells<<
	582	// " viewcells back="<<info.viewCellsBack<<
	583	// " viewcells back="<<info.viewCellsFront<<endl;
[447]	584	#endif
[492]	585
	586	if (info.costRatio > termMaxCostRatio) {
[447]	587	// cout<<"Too big cost ratio "<<costRatio<<endl;
	588	stat.maxCostRatioNodes++;
	589	return -1;
	590	}
	591
[1199]	592	#if 0
[447]	593	cout<<
	594	"pvs="<<leaf->mPvsSize<<
	595	" rays="<<leaf->rays.size()<<
	596	" rc="<<leaf->GetAvgRayContribution()<<
[492]	597	" axis="<<info.axis<<endl;
[447]	598	#endif
[492]	599
	600	return info.axis;
[447]	601	}
	602
	603
[492]	604	void
[447]	605	RssTree::GetCostRatio(
	606	RssTreeLeaf *leaf,
[492]	607	SplitInfo &info
[447]	608	)
	609	{
	610
	611	AxisAlignedBox3 box;
	612	float minBox, maxBox;
	613
[492]	614	if (info.axis < 3) {
[447]	615	box = GetBBox(leaf);
[492]	616	minBox = box.Min(info.axis);
	617	maxBox = box.Max(info.axis);
[447]	618	} else {
	619	box = GetDirBBox(leaf);
[492]	620	minBox = box.Min(info.axis-3);
	621	maxBox = box.Max(info.axis-3);
[447]	622	}
	623
	624	float sizeBox = maxBox - minBox;
	625
	626	int pvsSize = leaf->GetPvsSize();
	627
[492]	628	if (!mImportanceBasedCost) {
	629	const int costMethod = 0;
	630
	631	switch (costMethod) {
	632	case 0: {
	633	float sum = info.pvsBack(info.position - minBox) + info.pvsFront(maxBox - info.position);
	634	float newCost = ct_div_ci + sum/sizeBox;
[1500]	635	float oldCost = (float)pvsSize;
[492]	636	info.costRatio = newCost/oldCost;
	637	break;
	638	}
[1199]	639	case 6: {
	640	float sum = info.raysBack(info.position - minBox) + info.raysFront(maxBox - info.position);
	641	float newCost = ct_div_ci + sum/sizeBox;
[1548]	642	float oldCost = (float)info.rays;
[1199]	643	info.costRatio = newCost/oldCost;
	644	break;
	645	}
	646	case 7: {
	647	float benefit =
	648	(sqr(info.contributionBack) + sqr(info.contributionFront))/
	649	(2.0f*sqr(info.contribution));
	650	info.costRatio = 1.0f/(benefit + Limits::Small);
	651	break;
	652	}
	653
	654
[492]	655	case 1: {
	656	float newContrib =
[1112]	657	(info.contributionBack*(info.position - minBox) +
	658	+
	659	info.contributionFront*(maxBox - info.position))/sizeBox;
	660	float oldContrib = info.contribution;
	661	info.costRatio = newContrib/oldContrib;
	662	// cout<<info.contribution<<" "<<info.contributionBack<<" "<<info.contributionFront<<endl;
[492]	663	break;
	664	}
	665	case 2: {
	666	float sum =
	667	info.viewCellsBack*(info.position - minBox) +
	668	info.viewCellsFront*(maxBox - info.position);
	669	float newCost = ct_div_ci + sum/sizeBox;
[1548]	670	float oldCost = (float)info.viewCells;
[492]	671	info.costRatio = newCost/oldCost;
	672	break;
	673	}
	674	case 3: {
[1548]	675	float newCost = (float)(info.raysBackinfo.pvsBack + info.raysFrontinfo.pvsFront);
[1189]	676	float oldCost = (float)leaf->rays.size()*pvsSize;
[492]	677	info.costRatio = newCost/oldCost;
	678	}
	679	}
[447]	680	} else {
[492]	681	const int costMethod = 1;
[447]	682	// importance based cost
[492]	683	switch (costMethod) {
	684	case 0: {
	685	break;
	686	}
	687	case 1: {
	688	float newContrib =
	689	sqr(info.pvsBack/(info.raysBack + Limits::Small)) +
	690	sqr(info.pvsFront/(info.raysFront + Limits::Small));
	691	float oldContrib = sqr(leaf->GetAvgRayContribution());
	692	info.costRatio = oldContrib/newContrib;
	693	break;
	694	}
	695	case 2: {
	696	float newCost = (info.pvsBack + info.pvsFront)*0.5f;
[1548]	697	float oldCost = (float)pvsSize;
[492]	698	info.costRatio = newCost/oldCost;
	699	break;
	700	}
	701	case 3: {
[1548]	702	float newCost = (float)abs(info.raysBack - info.raysFront);
[1189]	703	float oldCost = (float)leaf->rays.size();
[492]	704	info.costRatio = newCost/oldCost;
	705	break;
	706	}
[447]	707	}
[492]	708	}
[447]	709	}
	710
	711
[492]	712	void
[447]	713	RssTree::EvalCostRatio(
	714	RssTreeLeaf *leaf,
[492]	715	SplitInfo &info
[447]	716	)
	717	{
[1199]	718	info.rays = 0;
[492]	719	info.raysBack = 0;
	720	info.raysFront = 0;
	721	info.pvsFront = 0;
	722	info.pvsBack = 0;
	723	info.viewCells = 0;
	724	info.viewCellsFront = 0;
	725	info.viewCellsBack = 0;
[447]	726
	727
[1199]	728
[492]	729	float sumWeights = Limits::Small;
	730	float sumWeightsBack = Limits::Small;
	731	float sumWeightsFront = Limits::Small;
	732
	733	float sumContribution = 0.0f;
	734	float sumContributionBack = 0.0f;
	735	float sumContributionFront = 0.0f;
[447]	736
[492]	737	Intersectable::NewMail();
[447]	738
[1199]	739	#if EVAL_VIEWCELLS
[492]	740	// count the numebr of viewcells first
	741	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	742	ri != leaf->rays.end();
	743	ri++)
	744	if ((*ri).mRay->IsActive()) {
[1199]	745
[492]	746	ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
	747	for (; it != (*ri).mRay->mViewCells.end(); ++it) {
	748	if (!(*it)->Mailed()) {
	749	(*it)->Mail();
	750	info.viewCells++;
	751	}
	752	}
	753	}
[1199]	754	#endif
[492]	755
	756	Intersectable::NewMail(3);
	757
	758	// this is the main ray classification loop!
	759	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	760	ri != leaf->rays.end();
	761	ri++)
	762	if ((*ri).mRay->IsActive()) {
	763	int side;
	764
	765
	766	// determine the side of this ray with respect to the plane
	767	side = (*ri).ComputeRaySide(info.axis, info.position);
	768
	769	float weight, contribution;
	770
	771	GetRayContribution(*ri, weight, contribution);
	772	sumWeights += weight;
	773	sumContribution += contribution;
[1199]	774	info.rays++;
[492]	775	if (side <= 0) {
	776	info.raysBack++;
	777	sumWeightsBack += weight;
	778	sumContributionBack += contribution;
	779	}
	780
	781	if (side >= 0) {
	782	info.raysFront++;
	783	sumWeightsFront += weight;
	784	sumContributionFront += contribution;
	785	}
	786
	787	AddObject2Pvs((*ri).GetObject(), side, info.pvsBack, info.pvsFront);
[1199]	788
	789	#if EVAL_VIEWCELLS
[492]	790	AddViewcells2Pvs((*ri).mRay->mViewCells,
	791	side,
	792	info.viewCellsBack,
	793	info.viewCellsFront);
[1199]	794	#endif
[492]	795	}
	796
[1112]	797	if (sumWeights!=0.0f)
	798	info.contribution = sumContribution/sumWeights;
	799	else
	800	info.contribution = 0.0f;
	801
	802	if (sumWeightsBack!=0.0f)
	803	info.contributionBack = sumContributionBack/sumWeightsBack;
	804	else
	805	info.contributionBack = 0.0f;
	806
	807	if (sumWeightsFront!=0.0f)
	808	info.contributionFront = sumContributionFront/sumWeightsFront;
	809	else
	810	info.contributionFront = 0.0f;
[1199]	811
	812	// info.costRatio = 0.1 + Random(0.5f);
	813	// return;
	814
[492]	815	GetCostRatio(
	816	leaf,
	817	info);
	818
[1112]	819	//cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
[447]	820	// float oldCost = leaf->rays.size();
[492]	821
[1112]	822	// cout<<"ratio="<<ratio<<endl;
[447]	823	}
	824
[492]	825	void
[447]	826	RssTree::BestCostRatio(
	827	RssTreeLeaf *leaf,
[492]	828	SplitInfo &info
[447]	829	)
	830	{
[492]	831	SplitInfo nInfo[5];
[447]	832	int bestAxis = -1;
[492]	833
[447]	834	AxisAlignedBox3 sBox = GetBBox(leaf);
	835	AxisAlignedBox3 dBox = GetDirBBox(leaf);
	836	// int sAxis = box.Size().DrivingAxis();
	837	int sAxis = sBox.Size().DrivingAxis();
	838	int dAxis = dBox.Size().DrivingAxis() + 3;
[492]	839
[447]	840	float dirSplitBoxSize = 0.01f;
	841	bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
	842
[1199]	843	bool useCyclingAxis = true;
	844
	845	if (splitType == ESplitHeuristic \|\|
	846	(splitType == ESplitHybrid && leaf->depth > mHybridDepth))
	847	useCyclingAxis = false;
[447]	848
[1199]	849	int cyclingAxis = 0;
	850
	851	if (leaf->parent)
	852	cyclingAxis = (leaf->parent->axis + 1)%5;
	853
	854	if (0 && useCyclingAxis) {
	855	int axis = cyclingAxis;
	856	info.position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
	857	info.axis = axis;
	858	info.costRatio = 0.5f; // not true but good for speeding up the subdivision at the top levels!
	859	info.raysBack = info.raysFront = leaf->rays.size()/2;
	860	return;
	861	}
	862
	863
[492]	864	for (int axis = 0; axis < 5; axis++)
[447]	865	if (
	866	(axis < 3 && (leaf->depth < mDirSplitDepth \|\| mInterleaveDirSplits)) \|\|
	867	(axis >= 3 && (leaf->depth >= mDirSplitDepth))
	868	) {
[492]	869	nInfo[axis].axis = axis;
[1199]	870	if ((!useCyclingAxis && (!mSplitUseOnlyDrivingAxis \|\| axis == sAxis \|\| axis == dAxis))
	871	\|\| (useCyclingAxis && (axis == cyclingAxis))
	872	) {
[447]	873
	874	if (splitType == ESplitRegular) {
	875	if (axis < 3)
[492]	876	nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
[447]	877	else
[492]	878	nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
	879	EvalCostRatio(leaf,
	880	nInfo[axis]);
[447]	881	} else
	882	if (splitType == ESplitHeuristic) {
[492]	883	EvalCostRatioHeuristic(
	884	leaf,
	885	nInfo[axis]
	886	);
[447]	887	} else
	888	if (splitType == ESplitHybrid) {
[1199]	889	if (leaf->depth > mHybridDepth)
[492]	890	EvalCostRatioHeuristic(
	891	leaf,
	892	nInfo[axis]
	893	);
[447]	894	else {
	895	if (axis < 3)
[492]	896	nInfo[axis].position = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
[447]	897	else
[492]	898	nInfo[axis].position = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
[447]	899
[492]	900	EvalCostRatio(leaf,
	901	nInfo[axis]
	902	);
[447]	903	}
	904	} else {
	905	cerr<<"RssTree: Unknown split heuristics\n";
	906	exit(1);
	907	}
	908
	909	if ( bestAxis == -1)
	910	bestAxis = axis;
	911	else
[492]	912	if ( nInfo[axis].costRatio < nInfo[bestAxis].costRatio )
[447]	913	bestAxis = axis;
	914	}
	915	}
	916
[492]	917	info = nInfo[bestAxis];
[447]	918	}
	919
	920
[492]	921	void
[447]	922	RssTree::EvalCostRatioHeuristic(
	923	RssTreeLeaf *leaf,
[492]	924	SplitInfo &info
[447]	925	)
	926	{
	927	AxisAlignedBox3 box;
	928	float minBox, maxBox;
	929
[492]	930	if (info.axis < 3) {
[447]	931	box = GetBBox(leaf);
[492]	932	minBox = box.Min(info.axis);
	933	maxBox = box.Max(info.axis);
[447]	934	} else {
	935	box = GetDirBBox(leaf);
[492]	936	minBox = box.Min(info.axis-3);
	937	maxBox = box.Max(info.axis-3);
[447]	938	}
	939
[1233]	940	SortSubdivisionCandidates(leaf, info.axis);
[447]	941
	942	// go through the lists, count the number of objects left and right
	943	// and evaluate the following cost funcion:
	944	// C = ct_div_ci + (qlrl + qrrr)/queries
[492]	945
	946	SplitInfo currInfo;
	947	currInfo.axis = info.axis;
	948
	949	currInfo.raysBack = 0;
	950	currInfo.raysFront = leaf->rays.size();
	951
	952	currInfo.pvsBack = 0;
	953	currInfo.pvsFront = leaf->GetPvsSize();
	954
[1199]	955
	956	float sumWeights = Limits::Small;
	957	float sumWeightsBack = Limits::Small;
	958	float sumWeightsFront = Limits::Small;
	959
	960	float sumContribution = 0.0f;
	961	float sumContributionBack = 0.0f;
	962	float sumContributionFront = 0.0f;
	963
	964
[447]	965	float sizeBox = maxBox - minBox;
[492]	966
[485]	967	float minBand = minBox + 0.1f*(maxBox - minBox);
	968	float maxBand = minBox + 0.9f*(maxBox - minBox);
[447]	969
[492]	970	// best cost ratio
	971	info.costRatio = 1e20f;
	972
	973	currInfo.viewCells = 0;
[1199]	974	currInfo.rays = 0;
[492]	975
[447]	976	Intersectable::NewMail();
	977	// set all object as belonging to the fron pvs
	978	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	979	ri != leaf->rays.end();
	980	ri++)
	981	if ((*ri).mRay->IsActive()) {
[1199]	982
	983	float weight, contribution;
	984	GetRayContribution(*ri, weight, contribution);
	985
	986	sumWeights += weight;
	987	sumContribution += contribution;
	988	currInfo.rays++;
	989
[492]	990	Intersectable object = (ri).GetObject();
[447]	991	if (object)
	992	if (!object->Mailed()) {
	993	object->Mail();
	994	object->mCounter = 1;
	995	} else
	996	object->mCounter++;
[492]	997
[1199]	998	#if EVAL_VIEWCELLS
[492]	999	// and do the same for all viewcells
	1000	ViewCellContainer::const_iterator it = (*ri).mRay->mViewCells.begin();
	1001
	1002	for (; it != (*ri).mRay->mViewCells.end(); ++it) {
	1003	ViewCell viewcell = it;
	1004	if (!viewcell->Mailed()) {
	1005	currInfo.viewCells++;
	1006	viewcell->Mail();
	1007	viewcell->mCounter = 1;
	1008	} else
	1009	viewcell->mCounter++;
	1010	}
[1199]	1011	#endif
[447]	1012	}
[492]	1013
[1199]	1014	if (sumWeights!=0.0f)
	1015	currInfo.contribution = sumContribution/sumWeights;
	1016	else
	1017	currInfo.contribution = 0.0f;
	1018
	1019	sumWeightsFront = sumWeights;
	1020	sumContributionFront = sumContribution;
	1021
	1022	sumWeightsBack = 0;
	1023	sumContributionBack = 0;
	1024
[492]	1025	currInfo.viewCellsBack = 0;
	1026	currInfo.viewCellsFront = currInfo.viewCells;
[1199]	1027
[447]	1028	Intersectable::NewMail();
[1199]	1029
[447]	1030	for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
	1031	ci < splitCandidates->end();
	1032	ci++) {
	1033	switch ((*ci).type) {
	1034	case SortableEntry::ERayMin: {
[492]	1035	currInfo.raysFront--;
	1036	currInfo.raysBack++;
[1199]	1037
[492]	1038	RssTreeNode::RayInfo rayInfo = (RssTreeNode::RayInfo ) (*ci).data;
[1199]	1039
	1040	float weight, contribution;
	1041	GetRayContribution(*rayInfo, weight, contribution);
	1042
	1043	sumWeightsBack += weight;
	1044	sumContributionBack += contribution;
	1045
	1046	sumWeightsFront -= weight;
	1047	sumContributionFront -= contribution;
	1048
[492]	1049	Intersectable *object = rayInfo->GetObject();
[447]	1050	if (object) {
	1051	if (!object->Mailed()) {
	1052	object->Mail();
[492]	1053	currInfo.pvsBack++;
[447]	1054	}
	1055	if (--object->mCounter == 0)
[492]	1056	currInfo.pvsFront--;
[447]	1057	}
[1199]	1058
	1059	#if EVAL_VIEWCELLS
[492]	1060	ViewCellContainer::const_iterator it = rayInfo->mRay->mViewCells.begin();
	1061	for (; it != rayInfo->mRay->mViewCells.end(); ++it) {
	1062	ViewCell viewcell = it;
	1063	if (!viewcell->Mailed()) {
	1064	viewcell->Mail();
	1065	currInfo.viewCellsBack++;
	1066	}
	1067	if (--viewcell->mCounter == 0)
	1068	currInfo.viewCellsFront--;
	1069	}
[1199]	1070	#endif
	1071
[447]	1072	break;
	1073	}
	1074	}
	1075
	1076	float position = (*ci).value;
	1077
	1078	if (position > minBand && position < maxBand) {
[492]	1079	currInfo.position = position;
[1199]	1080
[492]	1081
[1199]	1082	if (sumWeightsBack!=0.0f)
	1083	info.contributionBack = sumContributionBack/sumWeightsBack;
	1084	else
	1085	info.contributionBack = 0.0f;
	1086
	1087	if (sumWeightsFront!=0.0f)
	1088	info.contributionFront = sumContributionFront/sumWeightsFront;
	1089	else
	1090	info.contributionFront = 0.0f;
	1091
[492]	1092	GetCostRatio(
	1093	leaf,
	1094	currInfo);
	1095
[447]	1096
	1097	// cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
	1098	// cout<<"cost= "<<sum<<endl;
[492]	1099
	1100	if (currInfo.costRatio < info.costRatio) {
	1101	info = currInfo;
[447]	1102	}
	1103	}
	1104	}
	1105
[492]	1106
[447]	1107	// cout<<"===================="<<endl;
	1108	// cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
	1109	// <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
	1110	}
	1111
	1112	void
[1233]	1113	RssTree::SortSubdivisionCandidates(
[447]	1114	RssTreeLeaf *node,
	1115	const int axis
	1116	)
	1117	{
	1118
	1119	splitCandidates->clear();
	1120
	1121	int requestedSize = 2*(node->rays.size());
	1122	// creates a sorted split candidates array
	1123	if (splitCandidates->capacity() > 500000 &&
	1124	requestedSize < (int)(splitCandidates->capacity()/10) ) {
	1125
	1126	delete splitCandidates;
	1127	splitCandidates = new vector<SortableEntry>;
	1128	}
	1129
	1130	splitCandidates->reserve(requestedSize);
	1131
	1132	// insert all queries
	1133	for(RssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
	1134	ri < node->rays.end();
	1135	ri++) {
	1136	if ((*ri).mRay->IsActive()) {
	1137	if (axis < 3) {
	1138	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
[492]	1139	(*ri).GetOrigin(axis),
	1140	(void )&(ri))
[447]	1141	);
	1142	} else {
[492]	1143	float pos = (*ri).GetDirParametrization(axis-3);
[447]	1144	splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
	1145	pos,
[492]	1146	(void )&(ri))
[447]	1147	);
	1148	}
	1149	}
	1150	}
	1151
	1152	stable_sort(splitCandidates->begin(), splitCandidates->end());
	1153	}
	1154
	1155
	1156	void
	1157	RssTree::EvaluateLeafStats(const TraversalData &data)
	1158	{
	1159
	1160	// the node became a leaf -> evaluate stats for leafs
	1161	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
	1162
	1163	if (data.depth >= termMaxDepth)
	1164	stat.maxDepthNodes++;
	1165
	1166	// if ( (int)(leaf->rays.size()) < termMinCost)
	1167	// stat.minCostNodes++;
[556]	1168	if ( leaf->GetPvsSize() <= termMinPvs)
[447]	1169	stat.minPvsNodes++;
	1170
[556]	1171	if ( leaf->GetPvsSize() <= termMinRays)
[447]	1172	stat.minRaysNodes++;
	1173
	1174	if (leaf->GetAvgRayContribution() > termMaxRayContribution )
	1175	stat.maxRayContribNodes++;
	1176
	1177	if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
	1178	stat.minSizeNodes++;
	1179	}
	1180
	1181	if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
	1182	stat.maxRayRefs = leaf->rays.size();
	1183
	1184	}
	1185
	1186	bool
	1187	RssTree::TerminationCriteriaSatisfied(RssTreeLeaf *leaf)
	1188	{
[556]	1189	return ( (leaf->GetPvsSize() <= termMinPvs) \|\|
	1190	(leaf->rays.size() <= termMinRays) \|\|
[447]	1191	// (leaf->GetAvgRayContribution() > termMaxRayContribution ) \|\|
	1192	(leaf->depth >= termMaxDepth) \|\|
[464]	1193	(SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
[447]	1194	);
	1195	}
	1196
	1197
	1198	RssTreeNode *
	1199	RssTree::SubdivideNode(
	1200	RssTreeLeaf *leaf,
	1201	const AxisAlignedBox3 &box,
	1202	AxisAlignedBox3 &backBBox,
	1203	AxisAlignedBox3 &frontBBox
	1204	)
	1205	{
	1206
	1207	if (TerminationCriteriaSatisfied(leaf)) {
	1208	#if 0
	1209	if (leaf->depth >= termMaxDepth) {
	1210	cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
	1211	cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
	1212	}
	1213	#endif
	1214
	1215	return leaf;
	1216	}
	1217
[492]	1218	SplitInfo info;
[447]	1219
	1220	// select subdivision axis
[492]	1221	int axis = SelectPlane( leaf,
	1222	box,
	1223	info
	1224	);
	1225	// Debug<<"axis="<<axis<<" depth="<<(int)leaf->depth<<" rb="<<raysBack<<" rf="<<raysFront<<" pvsb="<<pvsBack<<" pvsf="<<pvsFront<<endl;
[447]	1226
	1227	if (axis == -1) {
	1228	return leaf;
	1229	}
	1230
	1231	stat.nodes+=2;
[569]	1232	stat.leaves += 1;
[447]	1233	stat.splits[axis]++;
	1234
	1235	// add the new nodes to the tree
	1236	RssTreeInterior *node = new RssTreeInterior(leaf->parent);
	1237
	1238	node->axis = axis;
[492]	1239	node->position = info.position;
[516]	1240	// node->bbox = box;
	1241	// node->dirBBox = GetDirBBox(leaf);
[447]	1242
[492]	1243
	1244	RssTreeLeaf *back = new RssTreeLeaf(node, info.raysBack);
	1245	RssTreeLeaf *front = new RssTreeLeaf(node, info.raysFront);
[447]	1246
[516]	1247
[492]	1248	// update halton generator
	1249	back->halton.index = leaf->halton.index;
	1250	front->halton.index = leaf->halton.index;
	1251
[447]	1252	// replace a link from node's parent
	1253	if ( leaf->parent )
	1254	leaf->parent->ReplaceChildLink(leaf, node);
	1255	// and setup child links
	1256	node->SetupChildLinks(back, front);
[516]	1257
	1258	back->bbox = leaf->bbox;
	1259	front->bbox = leaf->bbox;
	1260	back->dirBBox = leaf->dirBBox;
	1261	front->dirBBox = leaf->dirBBox;
	1262
[447]	1263	if (axis <= RssTreeNode::SPLIT_Z) {
[516]	1264	back->bbox.SetMax(axis, info.position);
	1265	front->bbox.SetMin(axis, info.position);
	1266	} else {
	1267	back->dirBBox.SetMax(axis-3, info.position);
	1268	front->dirBBox.SetMin(axis-3, info.position);
[492]	1269	}
[516]	1270
[492]	1271	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1272	ri != leaf->rays.end();
	1273	ri++) {
	1274	if ((*ri).mRay->IsActive()) {
	1275
	1276	// first unref ray from the former leaf
	1277	(*ri).mRay->Unref();
	1278
	1279	// Debug << "computed t: " << (*ri).mRay->mT << endl;
	1280	// determine the side of this ray with respect to the plane
	1281	int side = node->ComputeRaySide(*ri);
	1282
	1283	if (side == 1)
	1284	front->AddRay(*ri);
	1285	else
	1286	back->AddRay(*ri);
	1287	} else
	1288	(*ri).mRay->Unref();
[447]	1289	}
	1290
[492]	1291	// distribute the total number of rays according to the distribution
	1292	// of rays which remained
	1293	// front->mTotalRays = front->rays.size()*leaf->mTotalRays/leaf->rays.size();
	1294	// back->mTotalRays = back->rays.size()*leaf->mTotalRays/leaf->rays.size();
	1295
[447]	1296	#if 0
	1297	front->SetPvsSize(pvsFront);
	1298	back->SetPvsSize(pvsBack);
	1299	// compute entropy as well
	1300	front->ComputeEntropyImportance();
	1301	back->ComputeEntropyImportance();
	1302	#else
[492]	1303	UpdatePvsSize(front);
	1304	UpdatePvsSize(back);
[447]	1305	#endif
	1306
	1307
	1308	// update stats
[485]	1309	stat.rayRefs -= (int)leaf->rays.size();
[1199]	1310	stat.rayRefs += back->rays.size() + front->rays.size();
[447]	1311
	1312
	1313	delete leaf;
	1314	return node;
	1315	}
	1316
	1317
	1318
	1319
	1320
	1321
	1322	int
	1323	RssTree::ReleaseMemory(const int time)
	1324	{
	1325	stack<RssTreeNode *> tstack;
	1326
	1327	// find a node in the tree which subtree will be collapsed
	1328	int maxAccessTime = time - accessTimeThreshold;
	1329	int released;
	1330
[516]	1331	PushRoots(tstack);
[447]	1332
	1333	while (!tstack.empty()) {
	1334	RssTreeNode *node = tstack.top();
	1335	tstack.pop();
	1336
	1337
	1338	if (!node->IsLeaf()) {
	1339	RssTreeInterior in = (RssTreeInterior )node;
	1340	// cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
	1341	if (in->depth >= minCollapseDepth &&
	1342	in->lastAccessTime <= maxAccessTime) {
	1343	released = CollapseSubtree(node, time);
	1344	break;
	1345	}
	1346
	1347	if (in->back->GetAccessTime() <
	1348	in->front->GetAccessTime()) {
	1349	tstack.push(in->front);
	1350	tstack.push(in->back);
	1351	} else {
	1352	tstack.push(in->back);
	1353	tstack.push(in->front);
	1354	}
	1355	}
	1356	}
	1357
	1358	while (tstack.empty()) {
	1359	// could find node to collaps...
	1360	// cout<<"Could not find a node to release "<<endl;
	1361	break;
	1362	}
	1363
	1364	return released;
	1365	}
	1366
	1367
	1368
	1369
	1370	RssTreeNode *
	1371	RssTree::SubdivideLeaf(
	1372	RssTreeLeaf *leaf
	1373	)
	1374	{
	1375	RssTreeNode *node = leaf;
	1376
	1377	AxisAlignedBox3 leafBBox = GetBBox(leaf);
	1378
	1379	static int pass = 0;
	1380	pass ++;
	1381
	1382	// check if we should perform a dynamic subdivision of the leaf
	1383	if (!TerminationCriteriaSatisfied(leaf)) {
	1384
	1385	// memory check and realese...
	1386	if (GetMemUsage() > maxTotalMemory) {
	1387	ReleaseMemory( pass );
	1388	}
	1389
	1390	AxisAlignedBox3 backBBox, frontBBox;
	1391
	1392	// subdivide the node
	1393	node =
	1394	SubdivideNode(leaf,
	1395	leafBBox,
	1396	backBBox,
	1397	frontBBox
	1398	);
	1399	}
	1400
	1401	return node;
	1402	}
	1403
	1404
	1405
	1406	void
[492]	1407	RssTree::UpdateRays(
	1408	VssRayContainer &remove,
[447]	1409	VssRayContainer &add
	1410	)
	1411	{
	1412	RssTreeLeaf::NewMail();
	1413
	1414	// schedule rays for removal
	1415	for(VssRayContainer::const_iterator ri = remove.begin();
	1416	ri != remove.end();
	1417	ri++) {
	1418	(*ri)->ScheduleForRemoval();
	1419	}
	1420
	1421	int inactive=0;
	1422
	1423	for(VssRayContainer::const_iterator ri = remove.begin();
	1424	ri != remove.end();
	1425	ri++) {
	1426	if ((*ri)->ScheduledForRemoval())
	1427	// RemoveRay(*ri, NULL, false);
	1428	// !!! BUG - with true it does not work correctly - aggreated delete
	1429	RemoveRay(*ri, NULL, true);
	1430	else
	1431	inactive++;
	1432	}
	1433
	1434
	1435	// cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
	1436
	1437	for(VssRayContainer::const_iterator ri = add.begin();
	1438	ri != add.end();
	1439	ri++) {
	1440	RssTreeNode::RayInfo info(*ri);
[516]	1441	// if (mForcedBoundingBox==NULL \|\| ClipRay(info, bbox))
	1442	AddRay(info);
[447]	1443	}
[492]	1444
	1445	stat.rayRefs += add.size() - remove.size();
	1446
	1447	UpdateTreeStatistics();
	1448	// check whether the tree should be prunned
	1449	if (stat.rayRefs > mMaxRays) {
	1450	PruneRays(mMaxRays);
	1451	// UpdateTreeStatistics();
	1452	}
	1453
[447]	1454	}
	1455
[492]	1456
[447]	1457
[492]	1458
[447]	1459	void
	1460	RssTree::RemoveRay(VssRay *ray,
	1461	vector<RssTreeLeaf > affectedLeaves,
	1462	const bool removeAllScheduledRays
	1463	)
	1464	{
	1465
	1466	stack<RayTraversalData> tstack;
	1467
[516]	1468	PushRoots(tstack, RssTreeLeaf::RayInfo(ray));
[447]	1469
	1470	RayTraversalData data;
	1471
	1472	// cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
	1473
	1474	while (!tstack.empty()) {
	1475	data = tstack.top();
	1476	tstack.pop();
	1477
	1478	if (!data.node->IsLeaf()) {
	1479	// split the set of rays in two groups intersecting the
	1480	// two subtrees
	1481
	1482	TraverseInternalNode(data, tstack);
	1483
	1484	} else {
	1485	// remove the ray from the leaf
	1486	// find the ray in the leaf and swap it with the last ray...
	1487	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
	1488
	1489	if (!leaf->Mailed()) {
	1490	leaf->Mail();
	1491	if (affectedLeaves)
	1492	affectedLeaves->push_back(leaf);
	1493
	1494	if (removeAllScheduledRays) {
[485]	1495	int tail = (int)leaf->rays.size()-1;
[447]	1496
	1497	for (int i=0; i < (int)(leaf->rays.size()); i++) {
	1498	if (leaf->rays[i].mRay->ScheduledForRemoval()) {
	1499	// find a ray to replace it with
	1500	while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
	1501	stat.removedRayRefs++;
	1502	leaf->rays[tail].mRay->Unref();
	1503	leaf->rays.pop_back();
	1504	tail--;
	1505	}
	1506
	1507	if (tail < i)
	1508	break;
	1509
	1510	stat.removedRayRefs++;
	1511	leaf->rays[i].mRay->Unref();
	1512	leaf->rays[i] = leaf->rays[tail];
	1513	leaf->rays.pop_back();
	1514	tail--;
	1515	}
	1516	}
	1517	}
	1518	}
	1519
	1520	if (!removeAllScheduledRays)
	1521	for (int i=0; i < (int)leaf->rays.size(); i++) {
	1522	if (leaf->rays[i].mRay == ray) {
	1523	stat.removedRayRefs++;
	1524	ray->Unref();
	1525	leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
	1526	leaf->rays.pop_back();
	1527	// check this ray again
	1528	break;
	1529	}
	1530	}
	1531
	1532	}
	1533	}
	1534
	1535	if (ray->RefCount() != 0) {
	1536	cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
	1537	exit(1);
	1538	}
	1539
	1540	}
	1541
	1542
	1543	void
	1544	RssTree::AddRay(RssTreeNode::RayInfo &info)
	1545	{
	1546
	1547	stack<RayTraversalData> tstack;
[516]	1548
[569]	1549	RssTreeNode *root = GetRoot(info.GetSourceObject());
[447]	1550	tstack.push(RayTraversalData(root, info));
	1551
	1552	RayTraversalData data;
	1553
	1554	while (!tstack.empty()) {
	1555	data = tstack.top();
	1556	tstack.pop();
	1557
	1558	if (!data.node->IsLeaf()) {
	1559	TraverseInternalNode(data, tstack);
	1560	} else {
	1561	// remove the ray from the leaf
	1562	// find the ray in the leaf and swap it with the last ray...
	1563	RssTreeLeaf leaf = (RssTreeLeaf )data.node;
	1564	leaf->AddRay(data.rayData);
	1565	stat.addedRayRefs++;
	1566	}
	1567	}
	1568	}
	1569
	1570	void
	1571	RssTree::TraverseInternalNode(
	1572	RayTraversalData &data,
	1573	stack<RayTraversalData> &tstack)
	1574	{
	1575	RssTreeInterior in = (RssTreeInterior ) data.node;
	1576
[492]	1577	// determine the side of this ray with respect to the plane
	1578	int side = in->ComputeRaySide(data.rayData
	1579	);
[447]	1580
[492]	1581	if (side == 1)
	1582	tstack.push(RayTraversalData(in->front, data.rayData));
	1583	else
	1584	tstack.push(RayTraversalData(in->back, data.rayData));
	1585
[447]	1586	}
	1587
	1588
	1589	int
	1590	RssTree::CollapseSubtree(RssTreeNode *sroot, const int time)
	1591	{
	1592	// first count all rays in the subtree
	1593	// use mail 1 for this purpose
	1594	stack<RssTreeNode *> tstack;
	1595	int rayCount = 0;
	1596	int totalRayCount = 0;
	1597	int collapsedNodes = 0;
	1598
	1599	#if DEBUG_COLLAPSE
	1600	cout<<"Collapsing subtree"<<endl;
	1601	cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
	1602	cout<<"depth="<<(int)sroot->depth<<endl;
	1603	#endif
	1604
	1605	// tstat.collapsedSubtrees++;
	1606	// tstat.collapseDepths += (int)sroot->depth;
	1607	// tstat.collapseAccessTimes += time - sroot->GetAccessTime();
	1608
	1609	tstack.push(sroot);
	1610	VssRay::NewMail();
	1611
	1612	while (!tstack.empty()) {
	1613	collapsedNodes++;
	1614	RssTreeNode *node = tstack.top();
	1615	tstack.pop();
	1616
	1617	if (node->IsLeaf()) {
	1618	RssTreeLeaf leaf = (RssTreeLeaf ) node;
	1619	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1620	ri != leaf->rays.end();
	1621	ri++) {
	1622
	1623	totalRayCount++;
	1624	if ((ri).mRay->IsActive() && !(ri).mRay->Mailed()) {
	1625	(*ri).mRay->Mail();
	1626	rayCount++;
	1627	}
	1628	}
	1629	} else {
	1630	tstack.push(((RssTreeInterior *)node)->back);
	1631	tstack.push(((RssTreeInterior *)node)->front);
	1632	}
	1633	}
	1634
	1635	VssRay::NewMail();
	1636
	1637	// create a new node that will hold the rays
	1638	RssTreeLeaf *newLeaf = new RssTreeLeaf( sroot->parent, rayCount );
	1639	if ( newLeaf->parent )
	1640	newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
	1641
	1642
	1643	tstack.push( sroot );
	1644
	1645	while (!tstack.empty()) {
	1646
	1647	RssTreeNode *node = tstack.top();
	1648	tstack.pop();
	1649
	1650	if (node->IsLeaf()) {
	1651	RssTreeLeaf leaf = (RssTreeLeaf ) node;
	1652
	1653	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1654	ri != leaf->rays.end();
	1655	ri++) {
	1656
	1657	// unref this ray from the old node
	1658
	1659	if ((*ri).mRay->IsActive()) {
	1660	(*ri).mRay->Unref();
	1661	if (!(*ri).mRay->Mailed()) {
	1662	(*ri).mRay->Mail();
	1663	newLeaf->AddRay(*ri);
	1664	}
	1665	} else
	1666	(*ri).mRay->Unref();
	1667
	1668	}
	1669	} else {
	1670	tstack.push(((RssTreeInterior *)node)->back);
	1671	tstack.push(((RssTreeInterior *)node)->front);
	1672	}
	1673	}
	1674
	1675	// delete the node and all its children
	1676	delete sroot;
	1677
	1678	// for(RssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
	1679	// ri != newLeaf->rays.end();
	1680	// ri++)
	1681	// (*ri).ray->UnMail(2);
	1682
	1683
	1684	#if DEBUG_COLLAPSE
	1685	cout<<"Total memory before="<<GetMemUsage()<<endl;
	1686	#endif
	1687
	1688	stat.nodes -= collapsedNodes - 1;
[569]	1689	stat.leaves -= collapsedNodes/2 - 1;
[447]	1690	stat.rayRefs -= totalRayCount - rayCount;
	1691
	1692	#if DEBUG_COLLAPSE
	1693	cout<<"collapsed nodes"<<collapsedNodes<<endl;
	1694	cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
	1695	cout<<"Total memory after="<<GetMemUsage()<<endl;
	1696	cout<<"================================"<<endl;
	1697	#endif
	1698
	1699	// tstat.collapsedNodes += collapsedNodes;
	1700	// tstat.collapsedRays += totalRayCount - rayCount;
	1701
	1702	return totalRayCount - rayCount;
	1703	}
	1704
	1705
	1706	int
	1707	RssTree::GetPvsSize(const AxisAlignedBox3 &box) const
	1708	{
	1709	stack<RssTreeNode *> tstack;
[516]	1710	PushRoots(tstack);
[447]	1711
	1712	Intersectable::NewMail();
	1713	int pvsSize = 0;
	1714
	1715	while (!tstack.empty()) {
	1716	RssTreeNode *node = tstack.top();
	1717	tstack.pop();
	1718
	1719
	1720	if (node->IsLeaf()) {
	1721	RssTreeLeaf leaf = (RssTreeLeaf )node;
	1722	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1723	ri != leaf->rays.end();
	1724	ri++)
	1725	if ((*ri).mRay->IsActive()) {
	1726	Intersectable *object;
[492]	1727	object = (*ri).GetObject();
[447]	1728	if (object && !object->Mailed()) {
	1729	pvsSize++;
	1730	object->Mail();
	1731	}
	1732	}
	1733	} else {
	1734	RssTreeInterior in = (RssTreeInterior )node;
	1735	if (in->axis < 3) {
	1736	if (box.Max(in->axis) >= in->position )
	1737	tstack.push(in->front);
	1738
	1739	if (box.Min(in->axis) <= in->position )
	1740	tstack.push(in->back);
	1741	} else {
	1742	// both nodes for directional splits
	1743	tstack.push(in->front);
	1744	tstack.push(in->back);
	1745	}
	1746	}
	1747	}
	1748	return pvsSize;
	1749	}
	1750
	1751	int
	1752	RssTree::CollectPvs(const AxisAlignedBox3 &box,
	1753	ObjectContainer &pvs
	1754	) const
	1755	{
	1756	stack<RssTreeNode *> tstack;
[516]	1757	PushRoots(tstack);
[447]	1758
	1759	Intersectable::NewMail();
	1760
	1761	while (!tstack.empty()) {
	1762	RssTreeNode *node = tstack.top();
	1763	tstack.pop();
	1764
	1765
	1766	if (node->IsLeaf()) {
	1767	RssTreeLeaf leaf = (RssTreeLeaf )node;
	1768	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1769	ri != leaf->rays.end();
	1770	ri++)
	1771	if ((*ri).mRay->IsActive()) {
	1772	Intersectable *object;
[492]	1773	object = (*ri).GetObject();
[447]	1774	if (object && !object->Mailed()) {
	1775	pvs.push_back(object);
	1776	object->Mail();
	1777	}
	1778	}
	1779	} else {
	1780	RssTreeInterior in = (RssTreeInterior )node;
	1781	if (in->axis < 3) {
	1782	if (box.Max(in->axis) >= in->position )
	1783	tstack.push(in->front);
	1784
	1785	if (box.Min(in->axis) <= in->position )
	1786	tstack.push(in->back);
	1787	} else {
	1788	// both nodes for directional splits
	1789	tstack.push(in->front);
	1790	tstack.push(in->back);
	1791	}
	1792	}
	1793	}
	1794	return pvs.size();
	1795	}
	1796
	1797	void
[492]	1798	RssTree::UpdateTreeStatistics()
[447]	1799	{
	1800	stack<RssTreeNode *> tstack;
[516]	1801	PushRoots(tstack);
[447]	1802
	1803	float sumPvsSize = 0.0f;
	1804	float sumRayContribution = 0.0f;
[492]	1805	float sumWeightedRayContribution = 0.0f;
[447]	1806	float sumImportance = 0.0f;
	1807	float sumPvsEntropy = 0.0f;
	1808	float sumRayLengthEntropy = 0.0f;
[463]	1809	float sumRays = 0.0f;
[447]	1810
	1811	int leaves = 0;
	1812
	1813	while (!tstack.empty()) {
	1814	RssTreeNode *node = tstack.top();
	1815	tstack.pop();
	1816
	1817	if (node->IsLeaf()) {
	1818	leaves++;
	1819	RssTreeLeaf leaf = (RssTreeLeaf )node;
[492]	1820	UpdatePvsSize(leaf);
[463]	1821
[447]	1822	sumPvsSize += leaf->GetPvsSize();
	1823	sumRayContribution += leaf->GetAvgRayContribution();
[492]	1824
	1825
	1826	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
	1827	for (;it != leaf->rays.end(); ++it) {
	1828	float weight, contribution;
	1829	GetRayContribution(*it, weight, contribution);
	1830	sumWeightedRayContribution += weight*contribution;
	1831	}
	1832
	1833	// sumPvsEntropy += leaf->mPvsEntropy;
	1834	// sumRayLengthEntropy += leaf->mRayLengthEntropy;
[463]	1835	sumRays += leaf->rays.size();
	1836
[447]	1837	float imp = leaf->GetImportance();
	1838
[492]	1839	// if (imp > 1.0f)
	1840	// cout<<"warning imp > 1.0f:"<<imp<<endl;
[447]	1841
	1842	sumImportance += imp;
	1843
	1844
	1845	} else {
	1846	RssTreeInterior in = (RssTreeInterior )node;
	1847	// both nodes for directional splits
	1848	tstack.push(in->front);
	1849	tstack.push(in->back);
	1850	}
	1851	}
	1852
[492]	1853	stat.avgPvsSize = sumPvsSize/(float)leaves;
	1854	stat.avgRays = sumRays/(float)leaves;
	1855	stat.avgRayContribution = sumRayContribution/(float)leaves;
	1856	// avgPvsEntropy = sumPvsEntropy/(float)leaves;
	1857	// avgRayLengthEntropy = sumRayLengthEntropy/(float)leaves;
	1858	stat.avgImportance = sumImportance/(float)leaves;
	1859	stat.avgWeightedRayContribution = sumWeightedRayContribution/(float)sumRays;
	1860	stat.rayRefs = (int)sumRays;
[447]	1861	}
	1862
	1863
[492]	1864
[447]	1865	int
[492]	1866	RssTree::GenerateRays(
	1867	const float ratioPerLeaf,
[447]	1868	SimpleRayContainer &rays)
	1869	{
	1870	stack<RssTreeNode *> tstack;
[516]	1871	PushRoots(tstack);
[447]	1872
	1873	while (!tstack.empty()) {
	1874	RssTreeNode *node = tstack.top();
	1875	tstack.pop();
	1876
	1877	if (node->IsLeaf()) {
	1878	RssTreeLeaf leaf = (RssTreeLeaf )node;
	1879	float c = leaf->GetImportance();
[1189]	1880	int num = int(c*ratioPerLeaf + 0.5f);
[447]	1881	// cout<<num<<" ";
	1882
	1883	for (int i=0; i < num; i++) {
	1884	Vector3 origin = GetBBox(leaf).GetRandomPoint();
	1885	Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
	1886	Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
	1887	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	1888	rays.push_back(SimpleRay(origin, direction));
	1889	}
	1890
	1891	} else {
	1892	RssTreeInterior in = (RssTreeInterior )node;
	1893	// both nodes for directional splits
	1894	tstack.push(in->front);
	1895	tstack.push(in->back);
	1896	}
	1897	}
	1898
	1899	return rays.size();
	1900	}
	1901
	1902	void
	1903	RssTree::CollectLeaves(vector<RssTreeLeaf *> &leaves)
	1904	{
	1905	stack<RssTreeNode *> tstack;
[516]	1906	PushRoots(tstack);
[447]	1907
	1908	while (!tstack.empty()) {
	1909	RssTreeNode *node = tstack.top();
	1910	tstack.pop();
	1911
	1912	if (node->IsLeaf()) {
	1913	RssTreeLeaf leaf = (RssTreeLeaf )node;
	1914	leaves.push_back(leaf);
	1915	} else {
	1916	RssTreeInterior in = (RssTreeInterior )node;
	1917	// both nodes for directional splits
	1918	tstack.push(in->front);
	1919	tstack.push(in->back);
	1920	}
	1921	}
	1922	}
	1923
	1924	bool
	1925	RssTree::ValidLeaf(RssTreeLeaf *leaf) const
	1926	{
[463]	1927	return true;
	1928	//return leaf->rays.size() > termMinRays/4;
[447]	1929	}
	1930
	1931
	1932	void
[464]	1933	RssTree::PickEdgeRays(RssTreeLeaf *leaf,
	1934	int &r1,
	1935	int &r2
	1936	)
	1937	{
	1938	int nrays = leaf->rays.size();
	1939
	1940	if (nrays == 2) {
	1941	r1 = 0;
	1942	r2 = 1;
	1943	return;
	1944	}
[492]	1945
	1946	#if 1
[464]	1947	int tries = min(20, nrays);
	1948
	1949	while (--tries) {
	1950	r1 = Random(nrays);
	1951	r2 = Random(nrays);
[492]	1952	if (leaf->rays[r1].GetObject() != leaf->rays[r2].GetObject())
[464]	1953	break;
	1954	}
	1955
	1956	if (r1 == r2)
	1957	r2 = (r1+1)%leaf->rays.size();
[492]	1958
	1959	#else
	1960	// pick a random ray
	1961	int base = Random(nrays);
	1962	RssTreeNode::RayInfo *baseRay = &leaf->rays[base];
	1963	Intersectable *baseObject = baseRay->GetObject();
	1964
	1965	// and a random 5D derivative which will be used to find the minimal projected distances
	1966	Vector3 spatialDerivative;
	1967	Vector3 directionalDerivative;
	1968
	1969	while (1) {
	1970	spatialDerivative = Vector3(RandomValue(-1.0f, 1.0f),
	1971	RandomValue(-1.0f, 1.0f),
	1972	RandomValue(-1.0f, 1.0f));
	1973	float m = Magnitude(spatialDerivative);
	1974	if (m != 0) {
	1975	spatialDerivative /= m*Magnitude(GetBBox(leaf).Size());
	1976	break;
	1977	}
	1978	}
	1979
	1980	while (1) {
	1981	directionalDerivative = Vector3(RandomValue(-1.0f, 1.0f),
	1982	RandomValue(-1.0f, 1.0f),
	1983	0.0f);
	1984	float m = Magnitude(directionalDerivative);
	1985	if (m != 0) {
	1986	directionalDerivative /= m*Magnitude(GetDirBBox(leaf).Size());
	1987	break;
	1988	}
	1989	}
	1990
	1991	// now find the furthest sample from the same object and the closest from a different object
	1992	int i = 0;
	1993	float minDist = MAX_FLOAT;
	1994	float maxDist = -MAX_FLOAT;
	1995	r1 = base;
	1996	r2 = base;
	1997	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	1998	ri != leaf->rays.end();
	1999	ri++, i++) {
	2000	float dist = RssTreeNode::RayInfo::SqrDistance5D(*baseRay,
	2001	*ri,
	2002	spatialDerivative,
	2003	directionalDerivative
	2004	);
	2005
	2006	if ((*ri).GetObject() == baseObject) {
	2007	if (dist > maxDist) {
	2008	maxDist = dist;
	2009	r1 = i;
	2010	}
	2011	} else {
	2012	if (dist > 0.0f && dist < minDist) {
	2013	minDist = dist;
	2014	r2 = i;
	2015	}
	2016	}
	2017	}
	2018
	2019	#endif
[464]	2020	}
	2021
	2022
[492]	2023	struct RayNeighbor {
	2024	int rayInfo;
	2025	float distance;
	2026	RayNeighbor():rayInfo(0), distance(MAX_FLOAT) {}
	2027	};
	2028
[464]	2029	void
[492]	2030	RssTree::FindSilhouetteRays(RssTreeLeaf *leaf,
	2031	vector<RssTreeLeaf::SilhouetteRays> &rays
	2032	)
	2033	{
	2034	// for every leaf find its neares neighbor from a different object
	2035	vector<RayNeighbor> neighbors(leaf->rays.size());
	2036	int i, j;
	2037	for (i=0; i < leaf->rays.size(); i++)
	2038	for (j=0; j < leaf->rays.size(); j++)
	2039	if (i!=j) {
	2040	float d = RssTreeLeaf::RayInfo::Distance(leaf->rays[i], leaf->rays[j]);
	2041	if (d < neighbors[i].distance) {
	2042	neighbors[i].distance = d;
	2043	neighbors[i].rayInfo = j;
	2044	}
	2045	}
	2046
	2047	// now check which are the pairs of nearest neighbors
	2048	for (i=0; i < leaf->rays.size(); i++) {
	2049	int j = neighbors[i].rayInfo;
	2050	if (neighbors[j].rayInfo == i) {
	2051	// this is a silhouette edge pair
	2052	if (i < j) {
	2053	// generate an silhoutte ray pair
	2054	rays.push_back(RssTreeLeaf::SilhouetteRays(&leaf->rays[i],
	2055	&leaf->rays[j]));
	2056	}
	2057	} else {
	2058	// this is not a silhouette ray - delete???
	2059	}
	2060	}
	2061
	2062	}
	2063
	2064
	2065	bool
	2066	GetRandomTripple(vector<RssTreeLeaf::RayInfo *> &rays,
	2067	const int index,
	2068	int &i1,
	2069	int &i2,
	2070	int &i3)
	2071	{
	2072	int found = 0;
	2073	int indices[3];
	2074
	2075	int size = rays.size();
	2076	// use russian roulete selection for the tripple
	2077	// number of free positions for the bullet
	2078	int positions = size - 1;
	2079	int i;
	2080	for (i=0; i < size; i++) {
	2081	if (rays[i]->mRay->Mailed())
	2082	positions--;
	2083	}
	2084
	2085	if (positions < 3)
	2086	return false;
	2087
	2088	for (i=0; i < size; i++) {
	2089	if (i != index && !rays[i]->mRay->Mailed()) {
	2090	float p = (3 - found)/(float)positions;
	2091	if (Random(1.0f) < p) {
	2092	indices[found] = i;
	2093	found++;
	2094	}
	2095	}
	2096	positions--;
	2097	}
[500]	2098	return true; // corr. matt
[492]	2099	}
	2100
	2101	bool
	2102	RssTree::IsRayConvexCombination(const RssTreeNode::RayInfo &ray,
	2103	const RssTreeNode::RayInfo &r1,
	2104	const RssTreeNode::RayInfo &r2,
	2105	const RssTreeNode::RayInfo &r3)
	2106	{
	2107
	2108
	2109	return false;
	2110	}
	2111
	2112	int
	2113	RssTree::RemoveInteriorRays(
	2114	RssTreeLeaf *leaf
	2115	)
	2116	{
	2117	#if 1
	2118	// first collect all objects refered in this leaf
	2119	map<Intersectable , vector<RssTreeLeaf::RayInfo > > rayMap;
	2120
	2121	RssTreeLeaf::RayInfoContainer::iterator it = leaf->rays.begin();
	2122	for (; it != leaf->rays.end(); ++it) {
	2123	Intersectable object = (it).GetObject();
	2124
	2125	rayMap[object].push_back(&(*it));
	2126	// vector<RayInfo > data = rayMap.Find(object);
	2127	// if (data) {
	2128	// data->push_back(&(*it));
	2129	// } else {
	2130	// // rayMap[object] = vector<RayInfo *>;
	2131	// rayMap[object].push_back(&(*it));
	2132	// }
	2133	}
	2134
	2135	// now go through all objects
	2136	map<Intersectable , vector<RssTreeLeaf::RayInfo > >::iterator mi;
	2137
	2138	// infos of mailed rays are scheduled for removal
	2139	VssRay::NewMail();
	2140	for (mi = rayMap.begin(); mi != rayMap.end(); ++ mi) {
	2141	vector<RssTreeLeaf::RayInfo > &rays = (mi).second;
	2142	vector<RssTreeLeaf::RayInfo *>::iterator ri = rays.begin();
	2143	int rayIndex = 0;
	2144
	2145	for (; ri != rays.end(); ++ri, ++rayIndex) {
	2146	RssTreeNode::RayInfo ray = ri;
	2147	int tries = rays.size();
	2148	for (int i = 0; i < tries; i++) {
	2149	int r1, r2, r3;
	2150	GetRandomTripple(rays,
	2151	rayIndex,
	2152	r1,
	2153	r2,
	2154	r3);
	2155	if (IsRayConvexCombination(*ray,
	2156	*rays[r1],
	2157	*rays[r2],
	2158	*rays[r3])) {
	2159	ray->mRay->Mail();
	2160	}
	2161	}
	2162	}
	2163	}
	2164
	2165
	2166	#endif
	2167	return 0;
	2168	}
	2169
	2170	void
[447]	2171	RssTree::GenerateLeafRays(RssTreeLeaf *leaf,
	2172	const int numberOfRays,
	2173	SimpleRayContainer &rays)
	2174	{
[463]	2175
[492]	2176	if (numberOfRays == 0)
	2177	return;
	2178
[507]	2179
[447]	2180	int nrays = leaf->rays.size();
[534]	2181	int startIndex = rays.size();
	2182
[492]	2183	AxisAlignedBox3 box = GetBBox(leaf);
	2184	AxisAlignedBox3 dirBox = GetDirBBox(leaf);
	2185
[507]	2186
[1112]	2187	const float smoothRange = 0.1f;
[534]	2188	if (smoothRange != 0.0f) {
	2189	box.Scale(1.0f + smoothRange);
	2190	dirBox.Scale(1.0f + smoothRange);
[492]	2191	}
	2192
[534]	2193
[492]	2194	Exporter *exporter = NULL;
	2195	VssRayContainer selectedRays;
	2196	static int counter = 0;
	2197	if (counter < 0) {
	2198	char filename[128];
	2199	sprintf(filename, "selected-rays%04d.x3d", counter);
	2200	exporter = Exporter::GetExporter(filename);
	2201	// exporter->SetWireframe();
	2202	// exporter->ExportKdTree(*mKdTree);
	2203	exporter->SetWireframe();
[1328]	2204	// exporter->ExportScene(preprocessor->mSceneGraph->GetRoot());
[492]	2205	exporter->SetWireframe();
	2206	exporter->ExportBox(box);
	2207	exporter->SetFilled();
	2208	counter++;
	2209	}
	2210
[1199]	2211
	2212	float extendedConvexCombinationProb = 0.0f; //0.7f
	2213	// float silhouetteCheckPercentage = 1.0f; //0.5f
	2214	float silhouetteCheckPercentage = 0.0f; //0.9f; // 0.5f;
	2215	float silhouetteObjectCheckPercentage = 0.8f;
	2216
[1112]	2217	// int numberOfTries = numberOfRays*4;
	2218	int numberOfTries = numberOfRays;
[492]	2219	int generated = 0;
	2220
	2221	vector<Intersectable *> pvs(0);
[1199]	2222	if (silhouetteCheckPercentage != 0.0f) {
	2223	Intersectable::NewMail();
	2224	pvs.reserve(leaf->GetPvsSize());
	2225	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	2226	ri != leaf->rays.end();
	2227	ri++) {
	2228	Intersectable object = (ri).GetObject();
	2229	if (object) {
	2230	if (!object->Mailed()) {
	2231	object->Mail();
	2232	object->mCounter = 1;
	2233	pvs.push_back(object);
	2234	if (exporter)
	2235	exporter->ExportIntersectable(object);
	2236	} else {
	2237	object->mCounter++;
	2238	}
[556]	2239	}
[492]	2240	}
[1199]	2241	// sort objects based on their mCounter
	2242	sort(pvs.begin(), pvs.end(), Intersectable::GreaterCounter);
[492]	2243	}
[556]	2244
[492]	2245	for (int i=0; generated < numberOfRays && i < numberOfTries; i++) {
[549]	2246	bool useExtendedConvexCombination = ((nrays >= 2) && (Random(1.0f) <
	2247	extendedConvexCombinationProb));
[492]	2248
[447]	2249	Vector3 origin, direction;
	2250	// generate half of convex combination and half of random rays
	2251	if (useExtendedConvexCombination) {
[492]	2252	// pickup 2 random rays
[464]	2253	int r1, r2;
	2254	PickEdgeRays(leaf, r1, r2);
[492]	2255
[463]	2256
	2257	Vector3 o1 = leaf->rays[r1].GetOrigin();
	2258	Vector3 o2 = leaf->rays[r2].GetOrigin();
	2259
[492]	2260	const float overlap = 0.0f;
[464]	2261	float w1, w2;
	2262	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
	2263	origin = w1o1 + w2o2;
[492]	2264	GenerateExtendedConvexCombinationWeights2(w1, w2, overlap);
[447]	2265	direction =
[492]	2266	w1*leaf->rays[r1].GetDir() +
	2267	w2*leaf->rays[r2].GetDir();
[447]	2268	} else {
[492]	2269	Vector3 dirVector;
	2270
	2271	Vector3 pVector;
	2272	Vector3 dVector;
	2273
	2274	bool useHalton = true;
	2275
	2276	if (useHalton) {
	2277	// generate a random 5D vector
	2278	pVector = Vector3(leaf->halton.GetNumber(1),
	2279	leaf->halton.GetNumber(2),
	2280	leaf->halton.GetNumber(3));
	2281
	2282	dVector = Vector3(leaf->halton.GetNumber(4),
	2283	leaf->halton.GetNumber(5),
	2284	0.0f);
[534]	2285
[492]	2286	leaf->halton.GenerateNext();
	2287	} else {
	2288	pVector = Vector3(RandomValue(0.0f, 1.0f),
	2289	RandomValue(0.0f, 1.0f),
	2290	RandomValue(0.0f, 1.0f));
	2291
	2292	dVector = Vector3(RandomValue(0.0f, 1.0f),
	2293	RandomValue(0.0f, 1.0f),
	2294	0.0f);
	2295	}
	2296
[534]	2297	origin = box.GetPoint(pVector);
	2298	dirVector = dirBox.GetPoint(dVector);
	2299
[447]	2300	direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
	2301	}
[492]	2302
	2303	// shift the origin a little bit
	2304	direction.Normalize();
	2305
	2306	// float dist = Min(avgLength*0.5f, Magnitude(GetBBox(leaf).Size()));
	2307	float dist = 0.0f;
	2308	float minT, maxT;
	2309	// compute interection of the ray with the box
	2310
	2311	if (box.ComputeMinMaxT(origin, direction, &minT, &maxT) && minT < maxT)
	2312	dist = maxT;
	2313
	2314
	2315	origin += direction*dist;
	2316
	2317	bool intersects = false;
	2318
[1112]	2319	if (i >= numberOfRays*(1.0f - silhouetteCheckPercentage)) {
[492]	2320	if (exporter) {
	2321	VssRay ray = new VssRay(origin, origin + 100direction, NULL, NULL);
	2322	selectedRays.push_back(ray);
	2323	}
[567]	2324
[492]	2325	// check whether the ray does not intersect already visible objects
	2326	Ray traversalRay;
	2327	traversalRay.Init(origin, direction, Ray::LOCAL_RAY);
	2328	for(vector<Intersectable *>::const_iterator oi = pvs.begin();
	2329	oi != pvs.end();
	2330	oi++) {
	2331	Intersectable object = oi;
[567]	2332	// do not test every object
[1112]	2333	if ( Random(1.0f) <= silhouetteObjectCheckPercentage &&
	2334	object->CastRay(traversalRay) ) {
[492]	2335	intersects = true;
	2336	break;
	2337	}
	2338	}
	2339	}
	2340
	2341
	2342	if (!intersects) {
	2343	//cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
	2344	rays.push_back(SimpleRay(origin, direction));
	2345	generated++;
	2346	}
[447]	2347	}
[492]	2348
[534]	2349	// cout<<"desired="<<numberOfRays<<" tries="<<i<<endl;
[556]	2350	// assign pdfs to the generated rays
[1112]	2351	// float density = generated/(box.SurfaceArea()*dirBox.GetVolume());
	2352	float density = 1.0f;
	2353
	2354	float p = 1.0f/density; //(box.GetVolume()*dirBox.GetVolume())/i;
[537]	2355	for (i=startIndex; i < rays.size(); i++) {
[556]	2356	rays[i].mPdf = p;
[537]	2357	}
	2358
	2359
[492]	2360	if (exporter) {
	2361	exporter->ExportRays(selectedRays, RgbColor(1, 0, 0));
	2362	delete exporter;
	2363	CLEAR_CONTAINER(selectedRays);
	2364	}
[534]	2365
[447]	2366	}
	2367
	2368	int
[492]	2369	RssTree::PruneRays(RssTreeLeaf *leaf,
	2370	const float contributionThreshold)
	2371	{
	2372	int i;
	2373	int j;
	2374
	2375	for (j=0, i=0; i < leaf->rays.size(); i++) {
	2376
	2377	if (leaf->rays[i].mRay->mWeightedPvsContribution > contributionThreshold) {
	2378	// copy a valid sample
	2379	if (i!=j)
	2380	leaf->rays[j] = leaf->rays[i];
	2381	j++;
	2382	} else {
	2383	// delete the ray
	2384	leaf->rays[i].mRay->Unref();
	2385	if (leaf->rays[i].mRay->RefCount() != 0) {
	2386	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
	2387	exit(1);
	2388	}
	2389	delete leaf->rays[i].mRay;
	2390	}
	2391	}
	2392
	2393
	2394	leaf->rays.resize(j);
	2395	int removed = (i-j);
	2396	stat.rayRefs -= removed;
	2397	return removed;
	2398	}
	2399
	2400	int
	2401	RssTree::PruneRaysRandom(RssTreeLeaf *leaf,
	2402	const float ratio)
	2403	{
	2404	int i;
	2405	int j;
	2406
	2407	for (j=0, i=0; i < leaf->rays.size(); i++) {
	2408
	2409	if (Random(1.0f) < ratio) {
	2410	// copy a valid sample
	2411	if (i!=j)
	2412	leaf->rays[j] = leaf->rays[i];
	2413	j++;
	2414	} else {
	2415	// delete the ray
	2416	leaf->rays[i].mRay->Unref();
	2417	if (leaf->rays[i].mRay->RefCount() != 0) {
	2418	cerr<<"Error: refcount!=0, but"<<leaf->rays[j].mRay->RefCount()<<endl;
	2419	exit(1);
	2420	}
	2421	delete leaf->rays[i].mRay;
	2422	}
	2423	}
	2424
	2425
	2426	leaf->rays.resize(j);
	2427	int removed = (i-j);
	2428	stat.rayRefs -= removed;
	2429	return removed;
	2430	}
	2431
	2432	int
	2433	RssTree::PruneRaysContribution(RssTreeLeaf *leaf,
	2434	const float ratio)
	2435	{
	2436	int i;
	2437
	2438	if (leaf->rays.size() == 0)
	2439	return 0;
	2440
	2441	sort(leaf->rays.begin(),
	2442	leaf->rays.end(),
	2443	RssTreeLeaf::RayInfo::GreaterWeightedPvsContribution);
	2444
[1189]	2445	int desired = int(ratio*leaf->rays.size());
	2446	int removed = (int)leaf->rays.size() - desired;
[492]	2447
[1189]	2448	for (i=desired; i < (int)leaf->rays.size(); i++) {
[492]	2449	// delete the ray
	2450	leaf->rays[i].mRay->Unref();
	2451	if (leaf->rays[i].mRay->RefCount() != 0) {
	2452	cerr<<"Error: refcount!=0, but"<<leaf->rays[i].mRay->RefCount()<<endl;
[1112]	2453	cerr<<"desired="<<desired<<"i="<<i<<" size="<<leaf->rays.size()<<endl;
[492]	2454	exit(1);
	2455	}
	2456	delete leaf->rays[i].mRay;
	2457	}
	2458
	2459	leaf->rays.resize(desired);
	2460	stat.rayRefs -= removed;
	2461	return removed;
	2462	}
	2463
	2464
	2465	int
	2466	RssTree::PruneRays(
[1199]	2467	const int maxRays
[492]	2468	)
	2469	{
[1199]	2470	bool globalPrunning = false;
[492]	2471
	2472	stack<RssTreeNode *> tstack;
	2473	int prunned = 0;
	2474
[1199]	2475	// prune more rays to amortize the procedure
[1501]	2476	int desired = (int)(maxRays * 0.8f);
[492]	2477
[1199]	2478	Debug<<"Prunning rays...\nOriginal size "<<stat.rayRefs<<endl<<flush;
	2479
[492]	2480	if (globalPrunning) {
	2481	VssRayContainer allRays;
	2482	allRays.reserve(stat.rayRefs);
	2483	CollectRays(allRays);
	2484	sort(allRays.begin(),
	2485	allRays.end(),
	2486	GreaterWeightedPvsContribution);
	2487
[1199]	2488	if ( maxRays >= allRays.size() )
[492]	2489	return 0;
[1199]	2490
[492]	2491	float contributionThreshold = allRays[desired]->mWeightedPvsContribution;
	2492
[516]	2493	PushRoots(tstack);
[492]	2494
	2495	while (!tstack.empty()) {
	2496	RssTreeNode *node = tstack.top();
	2497	tstack.pop();
	2498
	2499	if (node->IsLeaf()) {
	2500	RssTreeLeaf leaf = (RssTreeLeaf )node;
	2501	prunned += PruneRays(leaf, contributionThreshold);
	2502	} else {
	2503	RssTreeInterior in = (RssTreeInterior )node;
	2504	// both nodes for directional splits
	2505	tstack.push(in->front);
	2506	tstack.push(in->back);
	2507	}
	2508	}
	2509	} else {
	2510	// prune random rays from each leaf so that the ratio's remain the same
[516]	2511	PushRoots(tstack);
[1199]	2512
	2513	if ( maxRays >= stat.rayRefs )
	2514	return 0;
	2515
[492]	2516	float ratio = desired/(float)stat.rayRefs;
	2517
	2518	while (!tstack.empty()) {
	2519	RssTreeNode *node = tstack.top();
	2520	tstack.pop();
	2521
	2522	if (node->IsLeaf()) {
	2523	RssTreeLeaf leaf = (RssTreeLeaf )node;
[567]	2524	// prunned += PruneRaysRandom(leaf, ratio);
[492]	2525	prunned += PruneRaysContribution(leaf, ratio);
	2526	} else {
	2527	RssTreeInterior in = (RssTreeInterior )node;
	2528	// both nodes for directional splits
	2529	tstack.push(in->front);
	2530	tstack.push(in->back);
	2531	}
	2532	}
	2533
	2534
	2535
	2536
	2537	}
	2538
	2539
	2540
[1199]	2541	Debug<<"Remained "<<stat.rayRefs<<" rays"<<endl<<flush;
[492]	2542
	2543	return prunned;
	2544	}
	2545
	2546	int
[447]	2547	RssTree::GenerateRays(const int numberOfRays,
	2548	const int numberOfLeaves,
	2549	SimpleRayContainer &rays)
	2550	{
	2551
[492]	2552
[447]	2553	vector<RssTreeLeaf *> leaves;
	2554
	2555	CollectLeaves(leaves);
	2556
[1112]	2557	if (numberOfLeaves != leaves.size())
	2558	sort(leaves.begin(),
	2559	leaves.end(),
	2560	GreaterContribution);
[447]	2561
	2562
	2563	float sumContrib = 0.0;
	2564	int i;
	2565	int k = 0;
	2566	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	2567	if (ValidLeaf(leaves[i])) {
	2568	float c = leaves[i]->GetImportance();
	2569	sumContrib += c;
	2570	// cout<<"ray contrib "<<i<<" : "<<c<<endl;
	2571	k++;
	2572	}
	2573
	2574	float avgContrib = sumContrib/numberOfLeaves;
[492]	2575
	2576	// always generate at leat n ray per leaf
[1112]	2577	int fixedPerLeaf = 0;
	2578	// int fixedPerLeaf = 0;
	2579	int fixed = fixedPerLeaf*leaves.size();
[492]	2580	int iGenerated = numberOfRays;
	2581	float ratioPerLeaf = iGenerated /(avgContrib*numberOfLeaves);
	2582
[447]	2583	k = 0;
	2584	for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
	2585	if (ValidLeaf(leaves[i])) {
	2586	k++;
	2587	RssTreeLeaf *leaf = leaves[i];
	2588	float c = leaf->GetImportance();
[492]	2589
	2590	int num = fixedPerLeaf + (int)(c*ratioPerLeaf + 0.5f);
[447]	2591	GenerateLeafRays(leaf, num, rays);
	2592	}
	2593
	2594	return rays.size();
	2595	}
	2596
	2597
	2598	float
	2599	RssTree::GetAvgPvsSize()
	2600	{
	2601	stack<RssTreeNode *> tstack;
[516]	2602	PushRoots(tstack);
[447]	2603
	2604	int sumPvs = 0;
	2605	int leaves = 0;
	2606	while (!tstack.empty()) {
	2607	RssTreeNode *node = tstack.top();
	2608	tstack.pop();
	2609
	2610	if (node->IsLeaf()) {
	2611	RssTreeLeaf leaf = (RssTreeLeaf )node;
	2612	// update pvs size
[492]	2613	UpdatePvsSize(leaf);
[447]	2614	sumPvs += leaf->GetPvsSize();
	2615	leaves++;
	2616	} else {
	2617	RssTreeInterior in = (RssTreeInterior )node;
	2618	// both nodes for directional splits
	2619	tstack.push(in->front);
	2620	tstack.push(in->back);
	2621	}
	2622	}
	2623
	2624
	2625	return sumPvs/(float)leaves;
	2626	}
	2627
[492]	2628	float weightAbsContributions = 0.0f;
	2629	// if small very high importance of the last sample
	2630	// if 1.0f then weighs = 1 1/2 1/3 1/4
[567]	2631	float passSampleWeightDecay = 1.0f;
[1112]	2632	//float passSampleWeightDecay = 0.0001f;
[492]	2633
	2634
[447]	2635	float
[492]	2636	RssTree::GetSampleWeight(const int pass)
[447]	2637	{
[492]	2638	int passDiff = mCurrentPass - pass;
	2639	float weight;
	2640	if (1)
[1112]	2641	weight = 1.0f/sqr(passDiff + passSampleWeightDecay);
[492]	2642	else
	2643	switch (passDiff) {
[1233]	2644	case 0:
[1112]	2645	// weight = 1.0f;
	2646	// break;
	2647	default:
[492]	2648	weight = 1.0f;
	2649	break;
	2650	// case 1:
	2651	// weight = 0.5f;
	2652	// break;
	2653	// case 2:
	2654	// weight = 0.25f;
	2655	// break;
	2656	// case 3:
	2657	// weight = 0.12f;
	2658	// break;
	2659	// case 4:
	2660	// weight = 0.06f;
	2661	// break;
	2662	// default:
	2663	// weight = 0.03f;
	2664	// break;
	2665	}
	2666	return weight;
	2667	}
[447]	2668
[492]	2669	void
[1112]	2670	RssTree::GetRayContribution(const RssTreeNode::RayInfo &info,
	2671	float &weight,
	2672	float &contribution)
	2673	{
	2674	VssRay *ray = info.mRay;
	2675	weight = GetSampleWeight(mCurrentPass);
	2676	contribution =
	2677	weightAbsContributions*ray->mPvsContribution +
	2678	(1.0f - weightAbsContributions)*ray->mRelativePvsContribution;
	2679	// store the computed value
	2680	info.mRay->mWeightedPvsContribution = weight*contribution;
	2681	}
	2682
	2683
	2684	void
[492]	2685	RssTree::ComputeImportance(RssTreeLeaf *leaf)
	2686	{
	2687	if (0)
	2688	leaf->mImportance = leaf->GetAvgRayContribution();
	2689	else {
	2690	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin(),
	2691	it_end = leaf->rays.end();
	2692
	2693	float sumContributions = 0.0f;
	2694	float sumRelContributions = 0.0f;
	2695	float sumWeights = 0.0f;
	2696	for (; it != it_end; ++it) {
	2697	VssRay ray = (it).mRay;
	2698	float weight = GetSampleWeight(ray->mPass);
	2699
	2700	sumContributions += weight*ray->mPvsContribution;
[1199]	2701	//$$ 20.7. changed to sqr to pronouce higher contribution so that they do not get smoothed!!
	2702	//sumRelContributions += weight*ray->mRelativePvsContribution;
	2703
	2704	sumRelContributions += sqr(weight*ray->mRelativePvsContribution);
[492]	2705
[1112]	2706	//sumWeights += weight;
[492]	2707	sumWeights += 1.0f;
	2708	}
	2709	// $$
	2710	// sumWeights = leaf->mTotalRays;
	2711
	2712	if (sumWeights != 0.0f)
	2713	leaf->mImportance =
	2714	(weightAbsContributions*sumContributions +
	2715	(1.0f - weightAbsContributions)*sumRelContributions)/sumWeights;
	2716	else
	2717	leaf->mImportance = 0.0f;
	2718
[447]	2719	}
[492]	2720
	2721	// return GetAvgRayContribution()*mImportance;
	2722	//return GetAvgRayContribution();
[447]	2723	}
	2724
	2725
	2726	float
[492]	2727	RssTreeLeaf::GetImportance() const
[447]	2728	{
[1112]	2729	//return sqr(mImportance);
[492]	2730	return mImportance;
	2731	}
	2732
	2733
	2734	float
	2735	RssTreeLeaf::ComputePvsEntropy() const
	2736	{
[447]	2737	int samples = 0;
	2738	Intersectable::NewMail();
	2739	// set all object as belonging to the fron pvs
[492]	2740	for(RssTreeNode::RayInfoContainer::const_iterator ri = rays.begin();
[447]	2741	ri != rays.end();
	2742	ri++)
	2743	if ((*ri).mRay->IsActive()) {
[492]	2744	Intersectable object = (ri).GetObject();
[447]	2745	if (object) {
	2746	if (!object->Mailed()) {
	2747	object->Mail();
	2748	object->mCounter = 1;
	2749	} else
	2750	object->mCounter++;
	2751	samples++;
	2752	}
	2753	}
	2754
	2755	float entropy = 0.0f;
	2756
	2757	if (samples > 1) {
	2758	Intersectable::NewMail();
	2759	for(RayInfoContainer::const_iterator ri = rays.begin();
	2760	ri != rays.end();
	2761	ri++)
	2762	if ((*ri).mRay->IsActive()) {
[492]	2763	Intersectable object = (ri).GetObject();
[447]	2764	if (object) {
	2765	if (!object->Mailed()) {
	2766	object->Mail();
	2767	float p = object->mCounter/(float)samples;
	2768	entropy -= p*log(p);
	2769	}
	2770	}
	2771	}
	2772	entropy = entropy/log((float)samples);
	2773	}
	2774	else
	2775	entropy = 1.0f;
	2776
	2777	return entropy;
	2778	}
	2779
	2780	float
[492]	2781	RssTreeLeaf::ComputeRayLengthEntropy() const
[447]	2782	{
	2783	// get sum of all ray lengths
	2784	// consider only passing rays or originating rays
	2785	float sum = 0.0f;
	2786	int samples = 0;
	2787	int i=0;
	2788	for(RayInfoContainer::const_iterator ri = rays.begin();
	2789	ri != rays.end();
	2790	ri++)
	2791	if ((*ri).mRay->IsActive()) {
	2792	// float s;
	2793	// if (i == 0)
	2794	// s = 200;
	2795	// else
	2796	// s = 100;
	2797	// i++;
	2798
	2799	sum += (*ri).mRay->GetSize();
	2800	samples++;
	2801	}
	2802
	2803
	2804	float entropy = 0.0f;
	2805
	2806	if (samples > 1) {
	2807	i = 0;
	2808	for(RayInfoContainer::const_iterator ri = rays.begin();
	2809	ri != rays.end();
	2810	ri++)
	2811	if ((*ri).mRay->IsActive()) {
	2812	// float s;
	2813	// if (i==0)
	2814	// s = 200;
	2815	// else
	2816	// s = 100;
	2817	// i++;
	2818	// float p = s/sum;
	2819	float p = (*ri).mRay->GetSize()/sum;
	2820	entropy -= p*log(p);
	2821	}
	2822	entropy = entropy/log((float)samples);
	2823	} else
	2824	entropy = 1.0f;
	2825
	2826	return entropy;
	2827	}
	2828
	2829
	2830
[492]	2831	float
	2832	RssTreeLeaf::ComputeEntropyImportance() const
[447]	2833	{
	2834
	2835	// mEntropy = 1.0f - ComputeRayLengthEntropy();
[492]	2836	return 1.0f - ComputePvsEntropy();
	2837	}
	2838
	2839	float
	2840	RssTreeLeaf::ComputeRayContributionImportance() const
	2841	{
	2842	// mPvsEntropy = ComputePvsEntropy();
	2843	// mRayLengthEntropy = ComputeRayLengthEntropy();
[447]	2844
[492]	2845	// mEntropy = 1.0f - ComputeRayLengthEntropy();
	2846	return 1.0f - ComputePvsEntropy();
[447]	2847	}
	2848
	2849
	2850	AxisAlignedBox3
	2851	RssTree::GetShrankedBBox(const RssTreeNode *node) const
	2852	{
	2853	if (node->parent == NULL)
	2854	return bbox;
	2855
	2856	if (!node->IsLeaf())
	2857	return ((RssTreeInterior *)node)->bbox;
	2858
	2859	// evaluate bounding box from the ray origins
	2860	AxisAlignedBox3 box;
	2861	box.Initialize();
	2862	RssTreeLeaf leaf = (RssTreeLeaf )node;
	2863	for(RssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
	2864	ri != leaf->rays.end();
	2865	ri++)
	2866	if ((*ri).mRay->IsActive()) {
	2867	box.Include((*ri).GetOrigin());
	2868	}
	2869	return box;
	2870	}
[466]	2871
	2872
	2873	int
	2874	RssTree::CollectRays(VssRayContainer &rays,
	2875	const int number)
	2876	{
	2877	VssRayContainer allRays;
	2878	CollectRays(allRays);
[492]	2879
[466]	2880
	2881	int desired = min(number, (int)allRays.size());
	2882	float prob = desired/(float)allRays.size();
	2883	while (rays.size() < desired) {
	2884	VssRayContainer::const_iterator it = allRays.begin();
	2885	for (; it != allRays.end() && rays.size() < desired; it++) {
	2886	if (Random(1.0f) < prob)
	2887	rays.push_back(*it);
	2888	}
	2889	}
[485]	2890	return (int)rays.size();
[466]	2891	}
	2892
	2893
	2894	int
	2895	RssTree::CollectRays(VssRayContainer &rays
	2896	)
	2897	{
	2898	VssRay::NewMail();
	2899
	2900	stack<RssTreeNode *> tstack;
[516]	2901	PushRoots(tstack);
[466]	2902
	2903	while (!tstack.empty()) {
	2904	RssTreeNode *node = tstack.top();
	2905	tstack.pop();
	2906
	2907	if (node->IsLeaf()) {
	2908	RssTreeLeaf leaf = (RssTreeLeaf )node;
	2909	// update pvs size
	2910	RssTreeNode::RayInfoContainer::const_iterator it = leaf->rays.begin();
	2911	for (;it != leaf->rays.end(); ++it)
	2912	if (!(*it).mRay->Mailed()) {
	2913	(*it).mRay->Mail();
	2914	rays.push_back((*it).mRay);
	2915	}
	2916	} else {
	2917	RssTreeInterior in = (RssTreeInterior )node;
	2918	// both nodes for directional splits
	2919	tstack.push(in->front);
	2920	tstack.push(in->back);
	2921	}
	2922	}
	2923
[485]	2924	return (int)rays.size();
[466]	2925	}
[492]	2926
[516]	2927
	2928	RssTreeNode *
	2929	RssTree::GetRoot(Intersectable *object) const
	2930	{
	2931	if (mPerObjectTree && object) {
[569]	2932	int id = object->GetId()-1;
	2933	if (id < 0 \|\| id >= mRoots.size()) {
	2934	Debug<<"Error: object Id out of range, Id="<<id<<" roots.size()="<<mRoots.size()<<
	2935	endl<<flush;
[579]	2936	id = (int)mRoots.size()-1; // $$ last tree is used by all unsigned objects
[569]	2937	}
[516]	2938	return mRoots[id];
	2939	} else
	2940	return mRoots[0];
	2941	}
	2942
	2943	void
	2944	RssTree::PushRoots(priority_queue<TraversalData> &st) const
	2945	{
	2946	for (int i=0; i < mRoots.size(); i++)
	2947	st.push(TraversalData(mRoots[i], GetBBox(mRoots[i]), 0));
	2948	}
	2949
	2950	void
	2951	RssTree::PushRoots(stack<RssTreeNode *> &st) const
	2952	{
	2953	for (int i=0; i < mRoots.size(); i++)
	2954	st.push(mRoots[i]);
	2955	}
	2956
	2957	void
	2958	RssTree::PushRoots(stack<RayTraversalData> &st, RssTreeNode::RayInfo &info) const
	2959	{
	2960	for (int i=0; i < mRoots.size(); i++)
	2961	st.push(RayTraversalData(mRoots[i], info));
	2962	}
[860]	2963
	2964
[1112]	2965	}

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

Download in other formats: