Index: GTP/trunk/Lib/Vis/Preprocessing/src/run_test2
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/run_test2	(revision 1966)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/run_test2	(revision 1974)
@@ -29,5 +29,5 @@
 VIEWCELLS=../data/vienna/vienna_cropped-gradient-viewcells.xml.gz
 
-PREFIX=../work/plots/osp2-1e5
+PREFIX=../work/plots/osp-1e5
 
 #SCENE=../data/atlanta/atlanta2.x3d
@@ -36,20 +36,20 @@
 
 
-$COMMAND -preprocessor=sampling -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
--preprocessor_visibility_file=$PREFIX-r-reference.xml \
--view_cells_filter_max_size=1 -preprocessor_stats=$PREFIX-r-reference.log \
--preprocessor_histogram_file=$PREFIX-r-reference.hlog
+# $COMMAND -preprocessor=sampling -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
+# -preprocessor_visibility_file=$PREFIX-r-reference.xml \
+# -view_cells_filter_max_size=1 -preprocessor_stats=$PREFIX-r-reference.log \
+# -preprocessor_histogram_file=$PREFIX-r-reference.hlog
 
 # $COMMAND -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
-# -rss_distributions=direction -view_cells_filter_max_size=1 \
-# -preprocessor_visibility_file=$PREFIX-r-reference-global.xml \
-# -preprocessor_stats=$PREFIX-r-reference-global.log \
-# -preprocessor_histogram_file=$PREFIX-r-reference-global.hlog
+#  -rss_distributions=direction -view_cells_filter_max_size=1 \
+#  -preprocessor_visibility_file=$PREFIX-r-reference-global.xml \
+#  -preprocessor_stats=$PREFIX-r-reference-global.log \
+#  -preprocessor_histogram_file=$PREFIX-r-reference-global.hlog
 
-# $COMMAND -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
-# -rss_distributions=mutation+object_direction+spatial -view_cells_filter_max_size=1 \
-# -preprocessor_visibility_file=$PREFIX-i-mixed-b1-n4c.xml \
-# -preprocessor_stats=$PREFIX-i-mixed-b1-n4c.log \
-# -preprocessor_histogram_file=$PREFIX-i-mixed-b1-n4c.hlog
+$COMMAND -preprocessor=combined -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
+ -rss_distributions=mutation+object_direction+spatial -view_cells_filter_max_size=1 \
+ -preprocessor_visibility_file=$PREFIX-i-mixed-b1-n4e.xml \
+ -preprocessor_stats=$PREFIX-i-mixed-b1-n4e.log \
+ -preprocessor_histogram_file=$PREFIX-i-mixed-b1-n4e.hlog
 
 # $COMMAND -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
Index: GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.h	(revision 1942)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.h	(revision 1974)
@@ -92,5 +92,9 @@
 	  float maxPvs;
 	  float avgPvs;
+	  float avgPvsEntries;
+	  
 	  float avgFilteredPvs;
+	  float avgFilteredPvsEntries;
+
 	  float avgFilterContribution;
 	  float avgFilterRadius;
@@ -775,4 +779,5 @@
 
 	bool mUseKdPvs;
+	bool mUseKdPvsAfterFiltering;
 
 	MixtureDistribution *mMixtureDistribution;
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Ray.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Ray.h	(revision 1969)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Ray.h	(revision 1974)
@@ -339,4 +339,11 @@
 	//	mId = sSimpleRayId++;
   }
+
+  float GetParam(const int axis) const {
+	if (axis < 3)
+	  return mOrigin[axis];
+	else
+	  return mDirection[axis-3];
+  }
   
   Vector3 Extrap(const float t) const {
Index: GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h	(revision 1972)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h	(revision 1974)
@@ -60,8 +60,19 @@
 							 const bool pruneInvalidRays = true
 							 ) = 0;
+
+
+  virtual void
+  SortRays(SimpleRayContainer &rays);
+
 	
 
 protected:
-	struct Intersection
+  void
+  _SortRays(SimpleRayContainer &rays,
+			const int l,
+			const int r,
+			const int axis);
+
+  struct Intersection
 	{
 		Intersection(): mObject(NULL), mFaceId(0)
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Environment.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Environment.cpp	(revision 1968)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Environment.cpp	(revision 1974)
@@ -1597,4 +1597,9 @@
 		"false");
 
+	RegisterOption("ViewCells.useKdPvsAfterFiltering",
+				   optBool, 
+				   "view_cells_use_kd_pvs_after_filtering", 
+				   "false");
+
 
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 1972)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 1974)
@@ -1149,6 +1149,13 @@
 					   )
 {
-	const long t1 = GetTime();
-
+
+
+  const long t1 = GetTime();
+
+#if 0
+  mRayCaster->SortRays(rays);
+  cout<<"Rays sorted in "<<TimeDiff(t1, GetTime())<<" s."<<endl;
+#endif
+  
 	SimpleRayContainer::const_iterator rit, rit_end = rays.end();
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/default.env
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/default.env	(revision 1967)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/default.env	(revision 1974)
@@ -36,6 +36,7 @@
 }
 
+
 Preprocessor {
-	totalSamples 10000000
+	totalSamples 250000000
 	samplesPerPass 1000000
 #	initialSamples 2000000
@@ -50,8 +51,8 @@
 	rayCastMethod 1
 	
-	type sampling
+#	type sampling
 #	type vss
 #	type rss
-#	type combined
+	type combined
 #	type render
 	detectEmptyViewSpace true
@@ -196,4 +197,5 @@
 ViewCells {
 	useKdPvs true
+	useKdPvsAfterFiltering true
 	# samples used for view cell construction
 	Construction {
Index: GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.h	(revision 1968)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.h	(revision 1974)
@@ -6,11 +6,13 @@
 
 #include "Halton.h"
-
 namespace GtpVisibilityPreprocessor {
 
-  class VssRay;
+class Vector2;
+class Vector3;
+class VssRay;
 class Preprocessor;
 struct SimpleRay;
 class SimpleRayContainer;
+
 struct VssRayContainer;
 
@@ -238,8 +240,32 @@
 	int mSamples;
 	HaltonSequence mHalton;
-
+	float mImportance;
+	float mCdf;
+
+	float GetSamplingFactor() const { return mSamples/mImportance; }
 	RayEntry() {}
-	RayEntry(VssRay *r):mRay(r), mSamples(0), mHalton() {}
+	RayEntry(VssRay *r):mRay(r), mSamples(0), mHalton(), mImportance(1.0f) {}
   };
+
+
+  Vector3
+  ComputeOriginMutation(const VssRay &ray,
+						const Vector3 &U,
+						const Vector3 &V,
+						const Vector2 vr2,
+						const float radius
+						);
+
+  Vector3
+  ComputeTerminationMutation(const VssRay &ray,
+							 const Vector3 &U,
+							 const Vector3 &V,
+							 const Vector2 vr2,
+							 const float radius
+							 );
+
+  RayEntry &GetEntry(const int index) {
+	return mRays[(mBufferStart+index)%mRays.size()];
+  }
   
   vector<RayEntry> mRays;
Index: GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.cpp	(revision 1972)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.cpp	(revision 1974)
@@ -667,4 +667,148 @@
 
 
+
+
+ViewCellsManager *ViewCellsManager::LoadViewCells(const string &filename, 
+												  ObjectContainer *objects,
+												  bool finalizeViewCells,
+												  BoundingBoxConverter *bconverter)
+												 
+{
+	ViewCellsParser parser;
+	ViewCellsManager *vm = NULL;
+
+	const long startTime = GetTime();
+	bool success = parser.ParseViewCellsFile(filename, &vm, objects, bconverter);
+
+	cout<<"viewcells parsed "<<endl<<flush;
+	
+	if (success)
+	{
+		//vm->ResetViewCells();
+		//hack
+		vm->mViewCells.clear();
+		ViewCellContainer leaves;
+		vm->mViewCellsTree->CollectLeaves(vm->mViewCellsTree->GetRoot(), leaves);
+
+		ViewCellContainer::const_iterator it, it_end = leaves.end();
+
+		for (it = leaves.begin(); it != it_end; ++ it)
+		{
+			vm->mViewCells.push_back(*it);
+		}
+		vm->mViewCellsFinished = true;
+		vm->mMaxPvsSize = (int)objects->size();
+
+		if (finalizeViewCells)
+		{
+			// create the meshes and compute volumes
+			vm->FinalizeViewCells(true);
+			// vm->mViewCellsTree->AssignRandomColors();
+		}
+
+		Debug << (int)vm->mViewCells.size() << " view cells loaded in "
+			  << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
+	}
+	else
+	{
+		Debug << "Error: loading view cells failed!" << endl;
+		DEL_PTR(vm);
+	}
+
+	return vm;
+}
+
+
+bool VspBspViewCellsManager::ExportViewCells(const string filename, 
+											 const bool exportPvs, 
+											 const ObjectContainer &objects)
+{
+	if (!ViewCellsConstructed() || !ViewCellsTreeConstructed())
+	{
+		return false;
+	}
+
+	cout << "exporting view cells to xml ... ";
+
+	OUT_STREAM stream(filename.c_str());
+
+	// for output we need unique ids for each view cell
+	CreateUniqueViewCellIds();
+
+	stream << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>"<<endl;
+	stream << "<VisibilitySolution>" << endl;
+
+	if (exportPvs) 
+	{
+		//-- export bounding boxes
+		stream << "<BoundingBoxes>" << endl;
+
+		if (mUseKdPvs)
+		{
+			vector<KdIntersectable *>::iterator kit, kit_end = GetPreprocessor()->mKdTree->mKdIntersectables.end();
+
+			int id = 0;
+			for (kit = GetPreprocessor()->mKdTree->mKdIntersectables.begin(); kit != kit_end; ++ kit, ++ id)
+			{
+				Intersectable *obj = *kit;
+				const AxisAlignedBox3 box = obj->GetBox();
+		  
+				obj->SetId(id);
+	
+				stream << "<BoundingBox" << " id=\"" << id << "\""
+					   << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
+					   << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\" />" << endl;
+			}
+		}
+		else
+		{
+			ObjectContainer::const_iterator oit, oit_end = objects.end();
+		
+			for (oit = objects.begin(); oit != oit_end; ++ oit)
+			{
+				const AxisAlignedBox3 box = (*oit)->GetBox();
+
+				////////////
+				//-- the bounding boxes
+
+				stream << "<BoundingBox" << " id=\"" << (*oit)->GetId() << "\""
+					   << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
+					   << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\" />" << endl;
+			}
+		}
+
+		stream << "</BoundingBoxes>" << endl;
+	}
+
+	
+	/////////////
+	//-- export the view cells and the pvs
+
+	const int numViewCells = mCurrentViewCellsStats.viewCells;
+	stream << "<ViewCells number=\"" << numViewCells << "\" >" << endl;
+
+	mViewCellsTree->Export(stream, exportPvs);
+
+	stream << "</ViewCells>" << endl;
+
+
+	//////////
+	//-- export the view space hierarchy
+	
+	stream << "<ViewSpaceHierarchy type=\"bsp\""
+		   << " min=\"" << mViewSpaceBox.Min().x << " " << mViewSpaceBox.Min().y << " " << mViewSpaceBox.Min().z << "\""
+		   << " max=\"" << mViewSpaceBox.Max().x << " " << mViewSpaceBox.Max().y << " " << mViewSpaceBox.Max().z << "\">" << endl;
+
+	mVspBspTree->Export(stream);
+	stream << "</ViewSpaceHierarchy>" << endl;
+
+	stream << "</VisibilitySolution>" << endl;
+
+	stream.close();
+	cout << "finished" << endl;
+
+	return true;
+}
+
 void ViewCellsManager::EvalViewCellHistogram(const string filename, 
 											 const int nViewCells)
@@ -757,148 +901,4 @@
 	outstream.close();
 }
-
-
-ViewCellsManager *ViewCellsManager::LoadViewCells(const string &filename, 
-												  ObjectContainer *objects,
-												  bool finalizeViewCells,
-												  BoundingBoxConverter *bconverter)
-												 
-{
-	ViewCellsParser parser;
-	ViewCellsManager *vm = NULL;
-
-	const long startTime = GetTime();
-	bool success = parser.ParseViewCellsFile(filename, &vm, objects, bconverter);
-
-	cout<<"viewcells parsed "<<endl<<flush;
-	
-	if (success)
-	{
-		//vm->ResetViewCells();
-		//hack
-		vm->mViewCells.clear();
-		ViewCellContainer leaves;
-		vm->mViewCellsTree->CollectLeaves(vm->mViewCellsTree->GetRoot(), leaves);
-
-		ViewCellContainer::const_iterator it, it_end = leaves.end();
-
-		for (it = leaves.begin(); it != it_end; ++ it)
-		{
-			vm->mViewCells.push_back(*it);
-		}
-		vm->mViewCellsFinished = true;
-		vm->mMaxPvsSize = (int)objects->size();
-
-		if (finalizeViewCells)
-		{
-			// create the meshes and compute volumes
-			vm->FinalizeViewCells(true);
-			// vm->mViewCellsTree->AssignRandomColors();
-		}
-
-		Debug << (int)vm->mViewCells.size() << " view cells loaded in "
-			  << TimeDiff(startTime, GetTime()) * 1e-3 << " secs" << endl;
-	}
-	else
-	{
-		Debug << "Error: loading view cells failed!" << endl;
-		DEL_PTR(vm);
-	}
-
-	return vm;
-}
-
-
-bool VspBspViewCellsManager::ExportViewCells(const string filename, 
-											 const bool exportPvs, 
-											 const ObjectContainer &objects)
-{
-	if (!ViewCellsConstructed() || !ViewCellsTreeConstructed())
-	{
-		return false;
-	}
-
-	cout << "exporting view cells to xml ... ";
-
-	OUT_STREAM stream(filename.c_str());
-
-	// for output we need unique ids for each view cell
-	CreateUniqueViewCellIds();
-
-	stream << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>"<<endl;
-	stream << "<VisibilitySolution>" << endl;
-
-	if (exportPvs) 
-	{
-		//-- export bounding boxes
-		stream << "<BoundingBoxes>" << endl;
-
-		if (mUseKdPvs)
-		{
-			vector<KdIntersectable *>::iterator kit, kit_end = GetPreprocessor()->mKdTree->mKdIntersectables.end();
-
-			int id = 0;
-			for (kit = GetPreprocessor()->mKdTree->mKdIntersectables.begin(); kit != kit_end; ++ kit, ++ id)
-			{
-				Intersectable *obj = *kit;
-				const AxisAlignedBox3 box = obj->GetBox();
-		  
-				obj->SetId(id);
-	
-				stream << "<BoundingBox" << " id=\"" << id << "\""
-					   << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
-					   << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\" />" << endl;
-			}
-		}
-		else
-		{
-			ObjectContainer::const_iterator oit, oit_end = objects.end();
-		
-			for (oit = objects.begin(); oit != oit_end; ++ oit)
-			{
-				const AxisAlignedBox3 box = (*oit)->GetBox();
-
-				////////////
-				//-- the bounding boxes
-
-				stream << "<BoundingBox" << " id=\"" << (*oit)->GetId() << "\""
-					   << " min=\"" << box.Min().x << " " << box.Min().y << " " << box.Min().z << "\""
-					   << " max=\"" << box.Max().x << " " << box.Max().y << " " << box.Max().z << "\" />" << endl;
-			}
-		}
-
-		stream << "</BoundingBoxes>" << endl;
-	}
-
-	
-	/////////////
-	//-- export the view cells and the pvs
-
-	const int numViewCells = mCurrentViewCellsStats.viewCells;
-	stream << "<ViewCells number=\"" << numViewCells << "\" >" << endl;
-
-	mViewCellsTree->Export(stream, exportPvs);
-
-	stream << "</ViewCells>" << endl;
-
-
-	//////////
-	//-- export the view space hierarchy
-	
-	stream << "<ViewSpaceHierarchy type=\"bsp\""
-		   << " min=\"" << mViewSpaceBox.Min().x << " " << mViewSpaceBox.Min().y << " " << mViewSpaceBox.Min().z << "\""
-		   << " max=\"" << mViewSpaceBox.Max().x << " " << mViewSpaceBox.Max().y << " " << mViewSpaceBox.Max().z << "\">" << endl;
-
-	mVspBspTree->Export(stream);
-	stream << "</ViewSpaceHierarchy>" << endl;
-
-	stream << "</VisibilitySolution>" << endl;
-
-	stream.close();
-	cout << "finished" << endl;
-
-	return true;
-}
-
 
 void ViewCellsManager::EvalViewCellHistogramForPvsSize(const string filename, 
@@ -1269,5 +1269,6 @@
 			filename = string(statsPrefix) + string(s);
 
-			EvalViewCellHistogramForPvsSize(filename, pass);
+			EvalViewCellHistogram(filename, pass);
+			//			EvalViewCellHistogramForPvsSize(filename, pass);
 		}
 	}
@@ -1992,5 +1993,7 @@
   stat.maxPvs = 0;
   stat.avgPvs = 0.0f;
+  stat.avgPvsEntries = 0.0f;
   stat.avgFilteredPvs = 0.0f;
+  stat.avgFilteredPvsEntries = 0.0f;
   stat.avgFilterContribution = 0.0f;
   stat.avgFilterRadius = 0;
@@ -2016,5 +2019,5 @@
 	  if (viewcell->GetValid()) {
 		const float pvsCost = mViewCellsTree->GetPvsCost(viewcell);
-		
+
 		if (pvsCost < stat.minPvs)
 		  stat.minPvs = pvsCost;
@@ -2024,4 +2027,7 @@
 		stat.avgPvs += pvsCost;
 
+		const float pvsEntries = mViewCellsTree->GetPvsEntries(viewcell);
+		stat.avgPvsEntries += pvsEntries;
+		
 
 		if (mPerViewCellStat[i].pvsSize > 0.0f)
@@ -2037,4 +2043,5 @@
 		if (evaluateFilter) {
 		  ObjectPvs filteredPvs;
+		  
 		  PvsFilterStatistics fstat = ApplyFilter2(viewcell,
 												   false,
@@ -2043,5 +2050,8 @@
 		  
 		  float filteredCost = filteredPvs.EvalPvsCost();
+
 		  stat.avgFilteredPvs += filteredCost;
+		  stat.avgFilteredPvsEntries += filteredPvs.GetSize();
+		  
 		  stat.avgFilterContribution += filteredCost - pvsCost;
 		  
@@ -2066,4 +2076,5 @@
   if (stat.viewcells) {
 	stat.avgPvs/=stat.viewcells;
+	stat.avgPvsEntries/=stat.viewcells;
 	stat.avgFilteredPvs/=stat.viewcells;
 	stat.avgFilterContribution/=stat.viewcells;
@@ -2089,5 +2100,7 @@
   GetPvsStatistics(pvsStat);
   s<<"#AVG_PVS\n"<<pvsStat.avgPvs<<endl;
+  s<<"#AVG_ENTRIES_PVS\n"<<pvsStat.avgPvsEntries<<endl;
   s<<"#AVG_FILTERED_PVS\n"<<pvsStat.avgFilteredPvs<<endl;
+  s<<"#AVG_FILTERED_ENTRIES_PVS\n"<<pvsStat.avgFilteredPvsEntries<<endl;
   s<<"#AVG_FILTER_CONTRIBUTION\n"<<pvsStat.avgFilterContribution<<endl;
   s<<"#AVG_FILTER_RADIUS\n"<<pvsStat.avgFilterRadius<<endl;
@@ -2991,5 +3004,10 @@
 	  objects.clear();
 	  // $$ warning collect objects takes only unmailed ones!
-	  CollectObjects(box, objects);
+	  if (mUseKdPvsAfterFiltering) {
+		float area = GetPreprocessor()->mKdTree->GetBox().SurfaceArea()*KD_PVS_AREA;
+		GetPreprocessor()->mKdTree->CollectKdObjects(box, objects, area);
+	  } else
+		CollectObjects(box, objects);
+	  
 	  //	cout<<"collected objects="<<objects.size()<<endl;
 	  ObjectContainer::const_iterator noi = objects.begin();
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Makefile
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 1967)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 1974)
@@ -1,5 +1,5 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (2.00a) (Qt 4.1.2) on: st 10. I 21:33:14 2007
+# Generated by qmake (2.00a) (Qt 4.1.2) on: ne 14. I 01:57:15 2007
 # Project:  preprocessor.pro
 # Template: app
@@ -63,5 +63,5 @@
 	$(MAKE) -f $(MAKEFILE).Debug uninstall
 
-Makefile: preprocessor.pro  C:/Qt/4.1.2/mkspecs/win32-msvc.net\qmake.conf C:/Qt/4.1.2/mkspecs/qconfig.pri \
+Makefile: preprocessor.pro  C:/Qt/4.1.2/mkspecs/win32-msvc2005\qmake.conf C:/Qt/4.1.2/mkspecs/qconfig.pri \
 		C:\Qt\4.1.2\mkspecs\features\qt_config.prf \
 		C:\Qt\4.1.2\mkspecs\features\exclusive_builds.prf \
Index: GTP/trunk/Lib/Vis/Preprocessing/src/KdTree.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/KdTree.cpp	(revision 1867)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/KdTree.cpp	(revision 1974)
@@ -826,4 +826,34 @@
 }
 
+void
+KdTree::CollectKdObjects(const AxisAlignedBox3 &box,
+						 ObjectContainer &objects,
+						 const float maxArea
+						 )
+{
+  stack<KdNode *> nodeStack;
+  
+  nodeStack.push(mRoot);
+
+  while (!nodeStack.empty()) {
+    KdNode *node = nodeStack.top();
+    nodeStack.pop();
+	if (node->IsLeaf() || (GetSurfaceArea(node) <= maxArea) ) {
+	  Intersectable *object = GetOrCreateKdIntersectable(node);
+	  if (!object->Mailed()) {
+		object->Mail();
+		objects.push_back(object);
+	  }
+	} else {
+      KdInterior *interior = (KdInterior *)node;
+	  
+	  if ( box.Max()[interior->mAxis] > interior->mPosition )
+		nodeStack.push(interior->mFront);
+	  
+	  if (box.Min()[interior->mAxis] < interior->mPosition)
+		nodeStack.push(interior->mBack);
+    }
+  }
+}
 
 void
Index: GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp	(revision 1966)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp	(revision 1974)
@@ -120,4 +120,40 @@
 }
 
+void
+RayCaster::SortRays(SimpleRayContainer &rays)
+{
+  _SortRays(rays, 0, rays.size()-1, 0);
+}
+					
+void
+RayCaster::_SortRays(SimpleRayContainer &rays,
+					 const int l,
+					 const int r,
+					 const int axis)
+{
+  // pick-up a pivot
+  int     i=l,j=r;
+  float x = rays[(l+r)/2].GetParam(axis);
+  do {
+	while(rays[i].GetParam(axis) < x)
+	  i++;
+	while(x < rays[j].GetParam(axis))
+	  j--;
+	
+	if (i <= j) {
+	  swap(rays[i], rays[j]);
+	  i++;
+	  j--;
+	}
+  } while (i<=j);
+  
+  if (l + 16 < j )
+	_SortRays(rays, l, j, (axis+1)%6);
+  
+  if (i + 16 < r)
+    _SortRays(rays, i, r, (axis+1)%6); 
+}
+  
+					
 
 int
Index: GTP/trunk/Lib/Vis/Preprocessing/src/KdTree.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/KdTree.h	(revision 1761)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/KdTree.h	(revision 1974)
@@ -345,4 +345,10 @@
 
   void
+  CollectKdObjects(const AxisAlignedBox3 &box,
+				   ObjectContainer &objects,
+				   const float maxArea
+				   );
+
+  void
   CollectObjects(KdNode *n, ObjectContainer &objects);
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro	(revision 1952)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro	(revision 1974)
@@ -111,5 +111,5 @@
 Trackball.cpp ObjExporter.cpp SubdivisionCandidate.cpp \
 Mailable.cpp \
-CombinedPreprocessor.cpp 
+CombinedPreprocessor.cpp Vector2.cpp
 
 SOURCES += BoostPreprocessorThread.cpp 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp	(revision 1968)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp	(revision 1974)
@@ -6,4 +6,6 @@
 #include "AxisAlignedBox3.h"
 #include "RssTree.h"
+#include "Vector2.h"
+#include "RndGauss.h"
 
 namespace GtpVisibilityPreprocessor {
@@ -566,6 +568,6 @@
 }
 
-#define RAY_CAST_TIME 0.9f
-#define VIEWCELL_CAST_TIME 0.1f
+#define RAY_CAST_TIME 0.7f
+#define VIEWCELL_CAST_TIME 0.3f
 
 void
@@ -589,28 +591,26 @@
   }
 
-  
-  
   if (sum == 0.0f)
 	sum = Limits::Small;
-
+  
   const float minratio = 0.01f;
-  float threshold = minratio*sum;
-
-  // recaluate the sum
+  
+  for (i=0; i < mDistributions.size(); i++) {
+	mDistributions[i]->mRatio /= sum;
+	if (mDistributions[i]->mRatio < minratio) 
+	  mDistributions[i]->mRatio = minratio;
+  }
+
+  // recaluate the sum after clip
   sum = 0.0f;
-  for (i=0; i < mDistributions.size(); i++) {
-	if (mDistributions[i]->mRatio < threshold) {
-	  mDistributions[i]->mRatio = threshold;
-	}
+  for (i=0; i < mDistributions.size(); i++)
 	sum += mDistributions[i]->mRatio;
-  }
-  
-  mDistributions[0]->mRatio /= sum;
+
+  for (i=0; i < mDistributions.size(); i++)
+	mDistributions[i]->mRatio /= sum;
   
   for (i=1; i < mDistributions.size(); i++) {
-	float r = mDistributions[i]->mRatio / sum;
-	mDistributions[i]->mRatio = mDistributions[i-1]->mRatio + r;
-  }
-  
+	mDistributions[i]->mRatio = mDistributions[i-1]->mRatio + mDistributions[i]->mRatio;
+  }
   
   cout<<"ratios: ";
@@ -692,8 +692,23 @@
   cerr<<"Muattion update..."<<endl;
   cerr<<"rays = "<<mRays.size()<<endl;
-  if (mRays.size())
-	cerr<<"Oversampling factor = "<<mRays[0].mSamples<<endl;
+  if (mRays.size()) {
+	cerr<<"Oversampling factors = "<<
+	  GetEntry(0).mSamples<<" "<<
+	  GetEntry(1).mSamples<<" "<<
+	  GetEntry(2).mSamples<<" "<<
+	  GetEntry(3).mSamples<<" "<<
+	  GetEntry(4).mSamples<<" "<<
+	  GetEntry(5).mSamples<<" ... "<<
+	  GetEntry(mRays.size()-6).mSamples<<" "<<
+	  GetEntry(mRays.size()-5).mSamples<<" "<<
+	  GetEntry(mRays.size()-4).mSamples<<" "<<
+	  GetEntry(mRays.size()-3).mSamples<<" "<<
+	  GetEntry(mRays.size()-2).mSamples<<" "<<
+	  GetEntry(mRays.size()-1).mSamples<<endl;
+  }
+  int contributingRays = 0;
   for (int i=0; i < vssRays.size(); i++) {
 	if (vssRays[i]->mPvsContribution) {
+	  contributingRays++;
 	  if (mRays.size() < mMaxRays) {
 		VssRay *newRay = new VssRay(*vssRays[i]);
@@ -712,60 +727,45 @@
 	}
   }
+  
+  float pContributingRays = contributingRays/(float)vssRays.size();
+  float importance = 1.0f/(pContributingRays + 1e-5);
+  // set this values for last contributingRays
+  int index = mBufferStart - 1;
+  int i;
+  for (i=0; i < contributingRays; i++, index--) {
+	if (index < 0)
+	  index = mRays.size()-1;
+	mRays[index].mImportance = importance;
+  }
+
+  
+  // compute cdf
+  mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mSamples+1);
+  for (i=1; i < mRays.size(); i++) 
+	mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mSamples+1);
+  
+  float scale = 1.0f/mRays[i-1].mCdf;
+  for (i=0; i < mRays.size(); i++) {
+	mRays[i].mCdf *= scale;
+  }
+
+  cout<<"Importance = "<<
+	GetEntry(0).mImportance<<" "<<
+	GetEntry(mRays.size()-1).mImportance<<endl;
+  
   cerr<<"Mutation update done."<<endl;
-
-}
-
-bool
-MutationBasedDistribution::GenerateSample(SimpleRay &sray)
-{
-  float r[5];
-
-  if (mRays.size() == 0) {
-	// use direction based distribution
-	Vector3 origin, direction;
-	static HaltonSequence halton;
-	
-	halton.GetNext(5, r);
-	mPreprocessor.mViewCellsManager->GetViewPoint(origin,
-												  Vector3(r[0], r[1], r[2]));
-	
-
-	direction = UniformRandomVector(r[3], r[4]);
-	
-	const float pdf = 1.0f;
-	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
-
-	return true;
-  } 
-  
-  //  int index = (int) (r[0]*mRays.size());
-  //  if (index >= mRays.size()) {
-  //	cerr<<"mutation: index out of bounds\n";
-  //	exit(1);
-  //  }
-
-  // get tail of the buffer
-  int index = (mLastIndex+1)%mRays.size();
-  if (mRays[index].mSamples > mRays[mLastIndex].mSamples) {
-	// search back for index where this is valid
-	index = (mLastIndex - 1 + mRays.size())%mRays.size();
-	for (int i=0; i < mRays.size(); i++) {
-	  if (mRays[index].mSamples > mRays[mLastIndex].mSamples)
-		break;
-	  index = (index - 1 + mRays.size())%mRays.size();
-	}
-	// go one step back
-	index = (index+1)%mRays.size();
-  }
-  
-  VssRay *ray = mRays[index].mRay;
-  mRays[index].mSamples++;
-  mLastIndex = index;
-
-  mRays[index].mHalton.GetNext(4, r);
-  
+}
+
+
+Vector3
+MutationBasedDistribution::ComputeOriginMutation(const VssRay &ray,
+												 const Vector3 &U,
+												 const Vector3 &V,
+												 const Vector2 vr2,
+												 const float radius
+												 )
+{
+#if 0
   Vector3 v;
-  // mutate the origin
-  Vector3 d = ray->GetDir();
   if (d.DrivingAxis() == 0)
 	v = Vector3(0, r[0]-0.5f, r[1]-0.5f);
@@ -775,9 +775,41 @@
 	else
 	  v = Vector3(r[0]-0.5f, r[1]-0.5f, 0);
-
-  v=v*mOriginMutationSize;
-
-  Vector3 origin = ray->mOrigin + v;
-  
+  return v*(2*radius);
+#endif
+#if 0
+  return (U*(r[0] - 0.5f) + V*(r[1] - 0.5f))*(2*radius);
+#endif
+
+
+  // Output random variable
+  Vector2 gaussvec2; 
+  
+  // Here we apply transform to gaussian, so 2D bivariate
+  // normal distribution
+  //  float sigma = ComputeSigmaFromRadius(radius);
+  float sigma = radius;
+  GaussianOn2D(vr2,
+			   sigma, // input
+			   gaussvec2);  // output
+
+	
+    // Here we tranform the point correctly to 3D space using base
+    // vectors of the 3D space defined by the direction
+    Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
+	
+	//	cout<<shift<<endl;
+	return shift;
+}
+
+Vector3
+MutationBasedDistribution::ComputeTerminationMutation(const VssRay &ray,
+													  const Vector3 &U,
+													  const Vector3 &V,
+													  const Vector2 vr2,
+													  const float radius
+													  )
+{
+#if 0
+  Vector3 v;
   // mutate the termination
   if (d.DrivingAxis() == 0)
@@ -788,29 +820,154 @@
 	else
 	  v = Vector3(r[2]-0.5f, r[3]-0.5f, 0);
+  
+  //   Vector3 nv;
+  
+  //   if (Magnitude(v) > Limits::Small)
+  // 	nv = Normalize(v);
+  //   else
+  // 	nv = v;
+  
+  //  v = nv*size + v*size;
+
+  return v*(4.0f*radius);
+#endif
+#if 0
+  return (U*(vr2.xx - 0.5f) + V*(vr2.yy - 0.5f))*(4.0f*radius);
+#endif
+  Vector2 gaussvec2; 
+#if 1
+  float sigma = radius;
+  GaussianOn2D(vr2,
+			   sigma, // input
+			   gaussvec2);  // output
+  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
+  //	cout<<shift<<endl;
+  return shift;
+#endif
+#if 0
+  // Here we estimate standard deviation (sigma) from radius
+  float sigma = 1.1f*ComputeSigmaFromRadius(radius);
+  Vector3 vr3(vr2.xx, vr2.yy, RandomValue(0,1));
+  PolarGaussianOnDisk(vr3,
+					  sigma,
+					  radius, // input
+					  gaussvec2); // output
+  
+  // Here we tranform the point correctly to 3D space using base
+  // vectors of the 3D space defined by the direction
+  Vector3 shift = gaussvec2.xx * U + gaussvec2.yy * V;
+  
+  //	cout<<shift<<endl;
+  return shift;
+#endif
+}
+
+
+bool
+MutationBasedDistribution::GenerateSample(SimpleRay &sray)
+{
+  float rr[5];
+
+  if (mRays.size() == 0) {
+	// use direction based distribution
+	Vector3 origin, direction;
+	static HaltonSequence halton;
+	
+	halton.GetNext(5, rr);
+	mPreprocessor.mViewCellsManager->GetViewPoint(origin,
+												  Vector3(rr[0], rr[1], rr[2]));
+	
+
+	direction = UniformRandomVector(rr[3], rr[4]);
+	
+	const float pdf = 1.0f;
+	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
+
+	return true;
+  } 
+
+  int index;
+#if 1
+  // get tail of the buffer
+  index = (mLastIndex+1)%mRays.size();
+  if (mRays[index].GetSamplingFactor() >
+	  mRays[mLastIndex].GetSamplingFactor()) {
+	// search back for index where this is valid
+	index = (mLastIndex - 1 + mRays.size())%mRays.size();
+	for (int i=0; i < mRays.size(); i++) {
+	  
+	  //	  if (mRays[index].mSamples > mRays[mLastIndex].mSamples)
+	  //		break;
+	  if (mRays[index].GetSamplingFactor() >
+		  mRays[mLastIndex].GetSamplingFactor() )
+		break;
+	  index = (index - 1 + mRays.size())%mRays.size();
+	}
+	// go one step back
+	index = (index+1)%mRays.size();
+  }
+#else
+  static HaltonSequence iHalton;
+  iHalton.GetNext(1, rr);
+  //rr[0] = RandomValue(0,1);
+  // use binary search to find index with this cdf
+  int l=0, r=mRays.size()-1;
+  while(l<r) {
+	int i = (l+r)/2;
+	if (rr[0] < mRays[i].mCdf )
+	  r = i;
+	else
+	  l = i+1;
+  }
+  index = l;
+  //  if (rr[0] >= mRays[r].mCdf)
+  //	index = r;
+  //  else
+  //	index = l;
+	
+	
+#endif
+  //  cout<<index<<" "<<rr[0]<<" "<<mRays[index].mCdf<<" "<<mRays[(index+1)%mRays.size()].mCdf<<endl;
+  
+  VssRay *ray = mRays[index].mRay;
+  mRays[index].mSamples++;
+  mLastIndex = index;
+
+  mRays[index].mHalton.GetNext(4, rr);
+
+  // mutate the origin
+  Vector3 d = ray->GetDir();
+
 
   AxisAlignedBox3 box = ray->mTerminationObject->GetBox();
-  float size = 2.0f*Magnitude(box.Diagonal());
-  if (size < Limits::Small)
+  float objectRadius = 0.5f*Magnitude(box.Diagonal());
+  //  cout<<objectRadius<<endl;
+  if (objectRadius < Limits::Small)
 	return false;
   
-//   Vector3 nv;
-  
-//   if (Magnitude(v) > Limits::Small)
-// 	nv = Normalize(v);
-//   else
-// 	nv = v;
-  
-//  v = nv*size + v*size;
-
-  v = v*size;
-
-  //  Vector3 termination = ray->mTermination + v;
-  Vector3 termination = box.Center() + v;
+  // Compute right handed coordinate system from direction
+  Vector3 U, V;
+  Vector3 nd = Normalize(d);
+  nd.RightHandedBase(U, V);
+
+  Vector3 origin = ray->mOrigin;
+  Vector3 termination = ray->mTermination; //box.Center(); //ray->mTermination; //box.Center();
+
+  float radiusExtension = 1.0f;
+  //  + mRays[index].mSamples/50.0f;
+  
+  //  origin += ComputeOriginMutation(*ray, U, V, Vector2(r[0], r[1]), 0.5f*mOriginMutationSize*radiusExtension);
+  origin += ComputeOriginMutation(*ray, U, V,
+								  Vector2(rr[0], rr[1]),
+								  objectRadius*radiusExtension);
+  termination += ComputeTerminationMutation(*ray, U, V,
+											Vector2(rr[2], rr[3]),
+											objectRadius*radiusExtension);
   
   Vector3 direction = termination - origin;
-
+  
   if (Magnitude(direction) < Limits::Small)
 	return false;
-
+  
   // shift the origin a little bit
   origin += direction*0.5f;
