Context Navigation

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

Revision 1145, 66.7 KB checked in by mattausch, 18 years ago (diff)
vsposp debug version

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

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

Download in other formats: