Context Navigation

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

Revision 2130, 74.4 KB checked in by mattausch, 18 years ago (diff)
runs also under debug mode now

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

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

Download in other formats: