Index: GTP/trunk/Lib/Vis/Preprocessing/scripts/run_gvs_vps.sh
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/scripts/run_gvs_vps.sh	(revision 2727)
+++ GTP/trunk/Lib/Vis/Preprocessing/scripts/run_gvs_vps.sh	(revision 2729)
@@ -42,6 +42,6 @@
 -gvs_total_samples=1000000 \
 -gvs_samples_per_pass=1000000 \
--gvs_initial_samples=50 \
--gvs_max_viewcells=1 \
+-gvs_initial_samples=4 \
+-gvs_max_viewcells=2 \
 -gvs_min_contribution=500 \
 -gvs_per_viewcell=true \
Index: GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.h	(revision 2728)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.h	(revision 2729)
@@ -146,10 +146,9 @@
 	*/
 	inline bool HandleRay(VssRay *ray);
-
 	/** The adaptive border sampling step. It aims to find neighbouring 
 		triangles of the one hit by the current ray.
 	*/	
 	int AdaptiveBorderSampling(const VssRay &currentRay);
-	
+	int AdaptiveBorderSamplingOpt(const VssRay &currentRay);
 	/** The reverse sampling step. It is started once the cast
 		ray finds a discontinuity, i.e., a closer triangle.
@@ -157,16 +156,8 @@
 		triangle passing through a gap.
 	*/
-	VssRay *ReverseSampling(const VssRay &currentRay,
-							const Triangle3 &hitTriangle,
-							const VssRay &oldRay);
-
-	bool ReverseSampling2(const VssRay &currentRay,
-                          const Triangle3 &hitTriangle,
-						  const VssRay &oldRay,
-						  VssRayContainer &reverseRays);
-
-	virtual void CastRays4(SimpleRayContainer &rays,
-						   VssRayContainer &vssRays);
-
+	bool ComputeReverseRay(const VssRay &currentRay,
+                       		 const Triangle3 &hitTriangle,
+							 const VssRay &oldRay,
+							 SimpleRay &reverseRay);
 	/** Returns true if we sampled a closer triangle than with the previous ray. 
 		Does reverse sampling if gap found.
@@ -175,10 +166,8 @@
 						   const Triangle3 &hitTriangle,
 						   const VssRay &oldRay);
-
 	/** Adds new samples to the ray queue and classifies them
 		with respect to the previous ray.
 	*/
 	void EnqueueRays(VssRayContainer &samples);
-
 	/** Hepler function for adaptive border sampling. It finds 
 		new sample points around a triangle in a eps environment
@@ -264,4 +253,22 @@
 
 	inline void AddKdNodeToPvs(const Vector3 &termination);
+	/** For all rayys of the ray container, generates a ray bundle of 4 jittered rays 
+		and casts them using optimized 4 eye ray casting.
+	*/
+	void CastRayBundles4(const SimpleRayContainer &rays, VssRayContainer &vssRays);
+	
+	void CastRayBundles16(const SimpleRayContainer &rays, VssRayContainer &vssRays);
+
+	/** Generates a ray bundle of 4 jittered rays and casts them using optimized 4 eye ray casting.
+	*/
+	void CastRayBundle4(const SimpleRay &ray, VssRayContainer &vssRays, const AxisAlignedBox3 &box);
+	/**  Cast rays using the optimized 4 eye ray casting routine.
+	*/
+	void CastRays4(const SimpleRayContainer &rays, VssRayContainer &vssRays);
+	/** Wrapper for optized ray casting routine taking 4 eye rays and a bounding box for ray clipping.
+	*/
+	void Cast4Rays(float *dist, Vector3 *dirs, Vector3 *origs, VssRayContainer &vssRays, const AxisAlignedBox3 &box);
+
+	void CastRayBundle16(const SimpleRay &ray, VssRayContainer &vssRays);
 
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h	(revision 2660)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h	(revision 2729)
@@ -8,6 +8,5 @@
 #include "Vector3.h"
 #include "VssRay.h"
-
-using namespace std;
+#include "Timer/PerfTimer.h"
 
 //
@@ -124,56 +123,4 @@
   SortRays2(SimpleRayContainer &rays);
   
-  // pool of vss rays to be used in one pass of the sampling
-  struct VssRayPool 
-  {
-    VssRayPool(): mRays(NULL), mIndex(0), mNumber(0) 
-    {}
-    
-    ~VssRayPool() 
-    {
-      delete []mRays;
-    }
-    
-    void Reserve(const int number) 
-    {
-      DEL_PTR(mRays);
-      mRays = new VssRay[number];
-      mNumber = number;
-    }
-    
-    void Clear() 
-    {
-      mIndex = 0;
-    }
-    
-    VssRay *Alloc() 
-    {
-      // reset pool
-      if (mIndex == mNumber)
-	mIndex = 0;
-      return mRays + mIndex ++;
-    }
-  protected:
-    VssRay *mRays;
-    int mIndex;
-    int mNumber;
-  };
-
-
-  VssRayPool mVssRayPool;
-  
-  void ReserveVssRayPool(const int n) 
-  {
-    mVssRayPool.Reserve(n);
-  }
-  
-  void InitPass() 
-  {
-    mVssRayPool.Clear();
-  }
-  
-  
-protected:
-  
   VssRay *RequestRay(const Vector3 &origin, 
 		     const Vector3 &termination, 
@@ -182,4 +129,71 @@
 		     const int pass, 
 		     const float pdf);
+  
+  // pool of vss rays to be used in one pass of the sampling
+  struct VssRayPool 
+  {
+    VssRayPool(): mRays(NULL), mIndex(0), mNumber(0) 
+    {}
+    
+    ~VssRayPool() 
+    {
+      delete []mRays;
+    }
+    
+    void Reserve(const int number) 
+    {
+      DEL_PTR(mRays);
+      mRays = new VssRay[number];
+      mNumber = number;
+    }
+    
+    void Clear() 
+    {
+      mIndex = 0;
+    }
+    
+    VssRay *Alloc() 
+    {
+#if 1
+		// reset pool
+		if (mIndex == mNumber) 
+		{
+			std::cerr << "warning: ray pool too small! " << std::endl;
+			mIndex = 0;
+		}
+
+		// raypool larger index => enlarge ray pool
+#else
+		if (mNumber == mIndex)
+		{
+			cerr << "warning: ray pool too small! " << "reserving " << mNumber * 2 << " rays" << std::endl;
+			Reserve(mNumber * 2);
+		}
+#endif
+		return mRays + mIndex ++;
+    }
+  protected:
+	  VssRay *mRays;
+	  int mIndex;
+	  int mNumber;
+  };
+
+
+  VssRayPool mVssRayPool;
+  
+  void ReserveVssRayPool(const int n) 
+  {
+    mVssRayPool.Reserve(n);
+  }
+  
+  void InitPass() 
+  {
+    mVssRayPool.Clear();
+  }
+  
+  PerfTimer rawCastTimer;
+
+protected:
+  
   
   void _SortRays(SimpleRayContainer &rays,
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 2726)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 2729)
@@ -1181,12 +1181,16 @@
 	{
 		const int n = 4;
-		cout << "hack: setting ray pool size to multiple " << n << " of samples per pass" << endl;       
 		reserveRays = mSamplesPerPass * n;
+		cout << "setting ray pool size to " << reserveRays << endl;       
+
 	}
 
 	cout << "======================" << endl;
 	cout << "reserving " << reserveRays << " rays " << endl;
+
 	mRayCaster->ReserveVssRayPool(reserveRays);
+	
 	cout<<"done."<<endl<<flush;
+	
 	return true;
 }
Index: GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.cpp	(revision 2726)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/ViewCellsManager.cpp	(revision 2729)
@@ -3217,7 +3217,5 @@
 							   )
 {
-
-  
-  pvs.Reserve(viewCell->GetFilteredPvsSize());
+	pvs.Reserve(viewCell->GetFilteredPvsSize());
 	
 	PvsFilterStatistics stats;
Index: GTP/trunk/Lib/Vis/Preprocessing/src/GtpVisibility05.vcproj
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/GtpVisibility05.vcproj	(revision 2727)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/GtpVisibility05.vcproj	(revision 2729)
@@ -102,7 +102,15 @@
 			<Tool
 				Name="VCCLCompilerTool"
+				Optimization="3"
+				InlineFunctionExpansion="2"
+				EnableIntrinsicFunctions="true"
+				FavorSizeOrSpeed="1"
+				OmitFramePointers="true"
 				AdditionalIncludeDirectories="..\src;Timer;..\..\..\..\..\..\NonGTP\Devil\include;..\..\..\..\..\..\NonGTP\Zlib\include;&quot;..\..\..\..\..\..\NonGTP\Xerces\xerces-c_2_8_0\include&quot;;..\..\..\..\..\..\NonGTP\Boost;..\src\havran;..\MultiLevelRayTracing;&quot;$(CG_INC_PATH)&quot;"
 				PreprocessorDefinitions="WIN32;NDEBUG;_CONSOLE;USE_VERBOSE_PVS;USE_QT;GTP_INTERNAL;USE_CG"
+				StringPooling="true"
 				RuntimeLibrary="2"
+				BufferSecurityCheck="false"
+				RuntimeTypeInfo="false"
 				WarningLevel="3"
 				Detect64BitPortabilityProblems="true"
Index: GTP/trunk/Lib/Vis/Preprocessing/src/common.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/common.h	(revision 2726)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/common.h	(revision 2729)
@@ -531,5 +531,5 @@
 #endif
 
-#define DYN_ERROR_STRING do {std::cout << "error: dynamic objects not supported" << endl;} while (0)
+#define DYN_ERROR_STRING do {std::cout << "error: dynamic objects not supported" << std::endl;} while (0)
 
 /** view cell id belonging to empty view space.
Index: GTP/trunk/Lib/Vis/Preprocessing/src/TestPreprocessor05.vcproj
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/TestPreprocessor05.vcproj	(revision 2725)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/TestPreprocessor05.vcproj	(revision 2729)
@@ -104,4 +104,5 @@
 			InheritedPropertySheets="$(VCInstallDir)VCProjectDefaults\UpgradeFromVC71.vsprops"
 			CharacterSet="2"
+			WholeProgramOptimization="1"
 			>
 			<Tool
@@ -123,9 +124,9 @@
 				Name="VCCLCompilerTool"
 				Optimization="2"
-				InlineFunctionExpansion="0"
-				EnableIntrinsicFunctions="false"
+				InlineFunctionExpansion="2"
+				EnableIntrinsicFunctions="true"
 				FavorSizeOrSpeed="1"
-				OmitFramePointers="false"
-				EnableFiberSafeOptimizations="false"
+				OmitFramePointers="true"
+				EnableFiberSafeOptimizations="true"
 				WholeProgramOptimization="true"
 				AdditionalIncludeDirectories="..\src;..\include;..\..\..\..\..\..\NonGTP\Boost;..\..\..\..\..\..\NonGTP\Devil\include;..\..\..\..\..\..\NonGTP\Zlib\include;&quot;$(QTDIR)\include\QtOpenGl&quot;;&quot;$(QTDIR)\include\Qt&quot;;&quot;$(QTDIR)\include\QtCore&quot;;&quot;$(QTDIR)\include&quot;;QtInterface;..\src\sparsehash\src\;..\src\ootl\src\include\ootl;..\src\ootl\src\include\ootl\sandbox;$(NOINHERIT)"
@@ -135,4 +136,5 @@
 				ExceptionHandling="1"
 				RuntimeLibrary="2"
+				BufferSecurityCheck="false"
 				EnableEnhancedInstructionSet="0"
 				RuntimeTypeInfo="false"
Index: GTP/trunk/Lib/Vis/Preprocessing/src/IntelRayCaster.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/IntelRayCaster.cpp	(revision 2710)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/IntelRayCaster.cpp	(revision 2729)
@@ -80,5 +80,7 @@
   }
   
+  rawCastTimer.Entry();
   mlrtaTraverseGroupASEye4(&minBox.x, &maxBox.x, result4, dist4);
+  rawCastTimer.Exit();
 
   if (saveMutationRays) {
@@ -126,4 +128,6 @@
   float dist1, dist2;
   double normal[3];
+
+  rawCastTimer.Entry();
   
   hittriangle = mlrtaIntersectAS(
@@ -132,4 +136,6 @@
 				 normal,
 				 dist);
+
+  rawCastTimer.Exit();
   dist1 = dist;
   
@@ -152,9 +158,15 @@
   {
     Vector3 dir = -simpleRay.mDirection;
-    hittriangle = mlrtaIntersectAS(
+
+	rawCastTimer.Entry();
+
+	hittriangle = mlrtaIntersectAS(
 				   &simpleRay.mOrigin.x,
 				   &dir.x,
 				   normal,
 				   dist);
+	
+	rawCastTimer.Exit();
+
     dist2 = dist;
 
@@ -265,4 +277,5 @@
     }
 #endif
+
     mlrtaStoreRayAS16(&rays[k].mOrigin.x,
 		      &rays[k].mDirection.x,
@@ -273,9 +286,12 @@
   Debug<<"TA\n"<<flush;
 #endif
-  
+  rawCastTimer.Entry();
+
   mlrtaTraverseGroupAS16(&min.x,
 			 &max.x,
 			 forward_hit_triangles,
 			 forward_dist);
+  
+  rawCastTimer.Exit();
   
 #if DEBUG_RAYCAST
@@ -294,9 +310,11 @@
     Debug<<"TB\n"<<flush;
 #endif
-    
+    rawCastTimer.Entry();
     mlrtaTraverseGroupAS16(&min.x,
 			   &max.x,
 			   backward_hit_triangles,
 			   backward_dist);
+
+	rawCastTimer.Exit();
   }
 	
@@ -439,9 +457,11 @@
 			&rays[k].mDirection.x,
 			j);
-    
+    rawCastTimer.Entry();
+
     mlrtaTraverseGroupAS16(&min.x,
 			   &max.x,
 			   hit_triangles,
 			   dist);	
+	rawCastTimer.Exit();
   } // for
 
@@ -449,4 +469,5 @@
   for (; k < rays.size(); k++) {
 	double normal[3];
+	rawCastTimer.Entry();
 	hit_triangles[0] = mlrtaIntersectAS(
 					    &rays[k].mOrigin.x,
@@ -454,4 +475,5 @@
 					    normal,
 					    dist[0]);
+	rawCastTimer.Exit();
   }
   
Index: GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.cpp	(revision 2728)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.cpp	(revision 2729)
@@ -159,62 +159,44 @@
 {
 	// the predicted hitpoint: we expect to hit the same mesh again
-	const Vector3 predictedHit = CalcPredictedHitPoint(currentRay, hitTriangle, oldRay);
-
-	const float predictedLen = Magnitude(predictedHit - currentRay.mOrigin);
-	const float len = Magnitude(currentRay.mTermination - currentRay.mOrigin);
+	Vector3 predictedHit = CalcPredictedHitPoint(currentRay, hitTriangle, oldRay);
+
+	float predictedLen = Magnitude(predictedHit - currentRay.mOrigin);
+	float len = Magnitude(currentRay.mTermination - currentRay.mOrigin);
 	
 	// distance large => this is likely to be a discontinuity
-	if ((predictedLen - len) > mThreshold) 
+	if (!((predictedLen - len) > mThreshold))
 	// q: rather use relative distance?
 	//if ((predictedLen / len) > mThreshold)
-	{
-#if 0
-		// apply reverse sampling to find the gap
-		VssRay *newRay = ReverseSampling(currentRay, hitTriangle, oldRay);
-
-		// check if reverse sampling was successful
-		if (newRay)
-		{
-			// set flag for visualization
-			if (0) newRay->mFlags |= VssRay::ReverseSample;
-		
-			rayTimer.Entry();
-			// check if ray can be processed further
-			// note: ray deletion handled by ray pool)
-			HandleRay(newRay);
-			rayTimer.Exit();
-
-			return 1;
-		}
-#else
-		static VssRayContainer reverseRays;
-		reverseRays.clear();
-
-		// apply reverse sampling to find the gap
-		ReverseSampling2(currentRay, hitTriangle, oldRay, reverseRays);
-
-		rayTimer.Entry();
-
-		int castRays = 0;
-		// check if reverse sampling was successful
-		for (size_t i = 0; i < reverseRays.size(); ++ i)
-		{
-			// check if ray can be processed further
-			// note: ray deletion handled by ray pool)
-			if (reverseRays[i])
-			{
-				++ castRays;
-				HandleRay(reverseRays[i]);
-			}
-		}
-		
-		rayTimer.Exit();
-
-		return castRays;
-
-#endif
-	}
-
-	return 0;
+		return 0;
+
+	SimpleRay simpleRay;
+
+	// apply reverse sampling to find the gap
+	if (!ComputeReverseRay(currentRay, hitTriangle, oldRay, simpleRay))
+		return 0;
+
+	static VssRayContainer reverseRays;
+	reverseRays.clear();
+
+	if (0)
+	{
+		VssRay *reverseRay = 
+			mRayCaster->CastRay(simpleRay, mViewCellsManager->GetViewSpaceBox(), !mPerViewCell);
+
+		reverseRays.push_back(reverseRay);
+	}
+	else
+	{
+		if (0)
+			CastRayBundle4(simpleRay, reverseRays, mViewCellsManager->GetViewSpaceBox());
+		else
+			CastRayBundle16(simpleRay, reverseRays);
+	}
+
+	mGvsStats.mReverseSamples += (int)reverseRays.size();
+
+	EnqueueRays(reverseRays);
+
+	return (int)reverseRays.size();
 }
 
@@ -422,10 +404,13 @@
 	SimpleRay sray(oldRay.mOrigin, p - oldRay.mOrigin, SamplingStrategy::GVS, 1.0f);
 
-#if 1
+	VssRay *newRay;
+
+	// cast single ray
 	castTimer.Entry();
-	// cast single ray to check if a triangle was missed
-	// note: not efficient!!
-	VssRay *newRay = mRayCaster->CastRay(sray, mViewCellsManager->GetViewSpaceBox(), !mPerViewCell);
+
+	// cast single ray to check if a triangle was missed. note: not efficient!!
+	newRay = mRayCaster->CastRay(sray, mViewCellsManager->GetViewSpaceBox(), !mPerViewCell);
 	castTimer.Exit();
+
 	++ castRays;
 
@@ -440,54 +425,5 @@
 
 	rayTimer.Exit();
-
-#else
-
-	// cast ray into the new point
-	static SimpleRayContainer importanceRays;
-	importanceRays.clear();
-
-	generationTimer.Entry();
-
-	
-	importanceRays.push_back(sray);
-
-	const int numRandomRays = 15;
-
-	GenerateJitteredRays(importanceRays, sray, numRandomRays, 0);
-	GenerateRays(numRandomRays, *mDistribution, importanceRays, sInvalidSamples);
-
-	generationTimer.Exit();
-	
-	// for the original implementation we don't cast double rays
-	const bool castDoubleRays = !mPerViewCell;
-	// cannot prune invalid rays because we have to compare adjacent rays
-	const bool pruneInvalidRays = false;
-
-	static VssRayContainer vssRays;
-	vssRays.clear();
-
-	castTimer.Entry();
-	CastRays(importanceRays, vssRays, castDoubleRays, pruneInvalidRays);
-	castTimer.Exit();
-
-	castRays += (int)vssRays.size();
-
-	VssRay *newRay = vssRays.empty() ? NULL : vssRays[0];
-
-	// new triangle discovered? => add new ray to queue
-	rayTimer.Entry();
-	for (size_t i = 0; i < vssRays.size(); ++ i)
-	{
-		if (vssRays[i])
-			HandleRay(vssRays[i]);
-	}
-
-	rayTimer.Exit();
-
-	// this ray will not be further processed
-	// note: ray deletion is handled by ray pool
-	if (!newRay) return castRays;
-#endif
-
+	
 	//newRay->mFlags |= VssRay::BorderSample;
 
@@ -516,7 +452,9 @@
 
 	//cout << "type: " << Intersectable::GetTypeName(tObj) << endl;
-
-	VertexContainer enlargedTriangle;
-	
+	generationTimer.Entry();
+
+	static VertexContainer enlargedTriangle;
+	enlargedTriangle.clear();
+
 	/// create 3 new hit points for each vertex
 	EnlargeTriangle(enlargedTriangle, hitTriangle, currentRay);
@@ -525,5 +463,4 @@
 	static SimpleRayContainer simpleRays;
 	simpleRays.clear();
-	//simpleRays.reserve(9);
 
 	VertexContainer::const_iterator vit, vit_end = enlargedTriangle.end();
@@ -551,14 +488,19 @@
 
 	const int numRandomRays = 16 - (int)simpleRays.size();
-	SimpleRay mainRay = SimpleRay(currentRay.GetOrigin(), currentRay.GetNormalizedDir(), SamplingStrategy::GVS, 1.0f);
+	SimpleRay mainRay = SimpleRay(currentRay.GetOrigin(), 
+		                          currentRay.GetNormalizedDir(), 
+								  SamplingStrategy::GVS, 
+								  1.0f);
 
 	GenerateJitteredRays(simpleRays, mainRay, numRandomRays, 0);
 	//GenerateRays(numRandomRays, *mDistribution, simpleRays, sInvalidSamples);
 
+	generationTimer.Exit();
 
 	/////////////////////
 	//-- cast rays to vertices and determine visibility
 	
-	VssRayContainer vssRays;
+	static VssRayContainer vssRays;
+	vssRays.clear();
 
 	// for the original implementation we don't cast double rays
@@ -567,7 +509,11 @@
 	const bool pruneInvalidRays = false;
 
+	//gvsTimer.Entry();
 	castTimer.Entry();
+
 	CastRays(simpleRays, vssRays, castDoubleRays, pruneInvalidRays);
+	
 	castTimer.Exit();
+	//gvsTimer.Exit();
 
 	const int numBorderSamples = (int)vssRays.size() - numRandomRays;
@@ -576,12 +522,13 @@
 	// handle cast rays
 	EnqueueRays(vssRays);
-	
-#if 0
-	// set flags
-	VssRayContainer::const_iterator rit, rit_end = vssRays.end();
-	for (rit = vssRays.begin(); rit != rit_end; ++ rit, ++ i)
-		(*rit)->mFlags |= VssRay::BorderSample;
-#endif
-
+
+	if (0)
+	{
+		// set flags
+		VssRayContainer::const_iterator rit, rit_end = vssRays.end();
+	
+		for (rit = vssRays.begin(); rit != rit_end; ++ rit)
+			(*rit)->mFlags |= VssRay::BorderSample;
+	}
 
     // recursivly subdivide each edge
@@ -602,4 +549,121 @@
 
 
+int GvsPreprocessor::AdaptiveBorderSamplingOpt(const VssRay &currentRay)
+{
+	Intersectable *tObj = currentRay.mTerminationObject;
+
+	// question matt: can this be possible?
+	if (!tObj) return 0;
+
+	Triangle3 hitTriangle;
+
+	// other types not implemented yet
+	if (tObj->Type() == Intersectable::TRIANGLE_INTERSECTABLE)
+		hitTriangle = static_cast<TriangleIntersectable *>(tObj)->GetItem();
+	else
+		cout << "border sampling: " << Intersectable::GetTypeName(tObj) << " not yet implemented" << endl;
+
+	generationTimer.Entry();
+
+	static VertexContainer enlargedTriangle;
+	enlargedTriangle.clear();
+	
+	/// create 3 new hit points for each vertex
+	EnlargeTriangle(enlargedTriangle, hitTriangle, currentRay);
+
+	static VertexContainer subdivEnlargedTri;	
+	subdivEnlargedTri.clear();
+
+	for (size_t i = 0; i < enlargedTriangle.size(); ++ i)
+	{
+		const Vector3 p1 = enlargedTriangle[i];
+		const Vector3 p2 = enlargedTriangle[(i + 1) % enlargedTriangle.size()];
+
+		// the new subdivision point
+		const Vector3 p = (p1 + p2) * 0.5f;
+
+		subdivEnlargedTri.push_back(p1);
+		subdivEnlargedTri.push_back(p);
+	}
+
+	/// create rays from sample points and handle them
+	static SimpleRayContainer simpleRays;
+	simpleRays.clear();
+	
+	VertexContainer::const_iterator vit, vit_end = subdivEnlargedTri.end();
+
+	for (vit = subdivEnlargedTri.begin(); vit != vit_end; ++ vit)
+	{
+		const Vector3 rayDir = (*vit) - currentRay.GetOrigin();
+
+		SimpleRay sr(currentRay.GetOrigin(), rayDir, SamplingStrategy::GVS, 1.0f);
+		simpleRays.AddRay(sr);	
+	}
+
+	
+	//////////
+	//-- fill up simple rays with random rays so we can cast 16 rays at a time
+
+	const int numRandomRays = 16 - ((int)simpleRays.size() % 16);
+
+	SimpleRay mainRay = SimpleRay(currentRay.GetOrigin(), 
+								  currentRay.GetNormalizedDir(), 
+								  SamplingStrategy::GVS, 1.0f);
+
+	GenerateJitteredRays(simpleRays, mainRay, numRandomRays, 0);
+	//GenerateRays(numRandomRays, *mDistribution, simpleRays, sInvalidSamples);
+
+	generationTimer.Exit();
+
+
+	/////////////////////
+	//-- cast rays to vertices and determine visibility
+	
+	static VssRayContainer vssRays;
+	vssRays.clear();
+
+	// for the original implementation we don't cast double rays
+	const bool castDoubleRays = !mPerViewCell;
+	// cannot prune invalid rays because we have to compare adjacent rays
+	const bool pruneInvalidRays = false;
+
+	//if ((simpleRays.size() % 16) != 0) cerr << "ray size not aligned" << endl;
+	
+	//gvsTimer.Entry();
+	castTimer.Entry();
+
+	CastRays(simpleRays, vssRays, castDoubleRays, pruneInvalidRays);
+	
+	castTimer.Exit();
+	//gvsTimer.Exit();
+
+	const int numBorderSamples = (int)vssRays.size() - numRandomRays;
+	int castRays = (int)simpleRays.size();
+
+	// handle cast rays
+	EnqueueRays(vssRays);
+	
+#if 0
+	
+    // recursivly subdivide each edge
+	for (int i = 0; i < numBorderSamples; ++ i)
+	{
+		castRays += SubdivideEdge(hitTriangle,
+								  subdivEnlargedTri[i], 
+								  subdivEnlargedTri[(i + 1) % numBorderSamples], 
+								  *vssRays[i], 
+								  *vssRays[(i + 1) % numBorderSamples],
+								  currentRay);
+	}
+	
+#endif
+
+
+	mGvsStats.mBorderSamples += castRays;
+	
+	return castRays;
+}
+
+
 bool GvsPreprocessor::GetPassingPoint(const VssRay &currentRay,
 									  const Triangle3 &occluder,
@@ -607,5 +671,6 @@
 									  Vector3 &newPoint) const
 {
-	//-- The plane p = (xp, hit(x), hit(xold)) is intersected
+	/////////////////////////
+	//-- We interserct the plane p = (xp, hit(x), hit(xold))
 	//-- with the newly found occluder (xold is the previous ray from
 	//-- which x was generated). On the intersecting line, we select a point
@@ -619,11 +684,6 @@
 	Vector3 pt1, pt2;
 
-	const bool intersects = occluder.GetPlaneIntersection(plane, pt1, pt2);
-
-	if (!intersects)
-	{
-		//cerr << "big error!! no intersection " << pt1 << " " << pt2 << endl;
+	if (!occluder.GetPlaneIntersection(plane, pt1, pt2))
 		return false;
-	}
 
 	// get the intersection point on the old ray
@@ -635,5 +695,5 @@
 	// evaluate new hitpoint just outside the triangle
 	// the point is chosen to be on the side closer to the original ray
-	if (Distance(pt1, pt3) < Distance(pt2, pt3))
+	if (SqrDistance(pt1, pt3) < SqrDistance(pt2, pt3))
 		newPoint = pt1 + mEps * (pt1 - pt2);
 	else
@@ -645,7 +705,8 @@
 
 
-VssRay *GvsPreprocessor::ReverseSampling(const VssRay &currentRay,
-										 const Triangle3 &hitTriangle,
-										 const VssRay &oldRay)
+bool GvsPreprocessor::ComputeReverseRay(const VssRay &currentRay,
+										const Triangle3 &hitTriangle,
+										const VssRay &oldRay,
+										SimpleRay &reverseRay)
 {
 	// get triangle occluding the path to the hit mesh
@@ -655,5 +716,5 @@
 	// q: why can this happen?
 	if (!tObj)
-		return NULL;
+		return false;
 
 	// other types not implemented yet
@@ -668,81 +729,5 @@
 	// q: why is there sometimes no intersecton found?
 	if (!GetPassingPoint(currentRay, occluder, oldRay, newPoint))
-		return NULL;
-
-	const Vector3 predicted = CalcPredictedHitPoint(currentRay, hitTriangle, oldRay);
-
-	Vector3 newDir, newOrigin;
-
-	//////////////
-	//-- Construct the mutated ray with xnew,
-	//-- dir = predicted(x)- pnew as direction vector
-
-	newDir = predicted - newPoint;
-
-	// take xnew, p = intersect(viewcell, line(pnew, predicted(x)) as origin ?
-	// difficult to say!!
-	const float offset = 0.5f;
-	newOrigin = newPoint - newDir * offset;
-	
-
-	//////////////
-	//-- for per view cell sampling, we must check for intersection
-	//-- with the current view cell
-
-    if (mPerViewCell)
-	{
-		// send ray to view cell
-		static Ray ray;
-		ray.Clear();
-		ray.Init(newOrigin, -newDir, Ray::LOCAL_RAY);
-		
-		//cout << "z";
-		// check if ray intersects view cell
-		if (!mCurrentViewCell->CastRay(ray))
-			return NULL;
-
-		Ray::Intersection &hit = ray.intersections[0];
-	
-		//cout << "q";
-		// the ray starts from the view cell
-		newOrigin = ray.Extrap(hit.mT);
-	}
-
-	++ mGvsStats.mReverseSamples;
-
-	const SimpleRay simpleRay(newOrigin, newDir, SamplingStrategy::GVS, 1.0f);
-
-	VssRay *reverseRay = 
-		mRayCaster->CastRay(simpleRay, mViewCellsManager->GetViewSpaceBox(), !mPerViewCell);
-
-	return reverseRay;
-}
-
-
-bool GvsPreprocessor::ReverseSampling2(const VssRay &currentRay,
-									   const Triangle3 &hitTriangle,
-									   const VssRay &oldRay,
-									   VssRayContainer &reverseRays)
-{
-	// get triangle occluding the path to the hit mesh
-	Triangle3 occluder;
-	Intersectable *tObj = currentRay.mTerminationObject;
-
-	// q: why can this happen?
-	if (!tObj)
-		return NULL;
-
-	// other types not implemented yet
-	if (tObj->Type() == Intersectable::TRIANGLE_INTERSECTABLE)
-		occluder = static_cast<TriangleIntersectable *>(tObj)->GetItem();
-	else
-		cout << "reverse sampling: " << tObj->Type() << " not yet implemented" << endl;
-	
-	// get a point which is passing just outside of the occluder
-    Vector3 newPoint;
-
-	// q: why is there sometimes no intersecton found?
-	if (!GetPassingPoint(currentRay, occluder, oldRay, newPoint))
-		return NULL;
+		return false;
 
 	const Vector3 predicted = CalcPredictedHitPoint(currentRay, hitTriangle, oldRay);
@@ -784,17 +769,5 @@
 	}
 
-	static SimpleRayContainer reverseSimpleRays;
-	reverseSimpleRays.clear();
-
-	SimpleRay mainRay = SimpleRay(newOrigin, newDir, SamplingStrategy::GVS, 1.0f);
-
-	reverseSimpleRays.push_back(mainRay);
-	GenerateJitteredRays(reverseSimpleRays, mainRay, 15, 0);
-
-	castTimer.Entry();
-	CastRays(reverseSimpleRays, reverseRays, false, false);
-	castTimer.Exit();
-
-	mGvsStats.mReverseSamples += (int)reverseRays.size();
+	reverseRay = SimpleRay(newOrigin, newDir, SamplingStrategy::GVS, 1.0f);
 
 	return true;
@@ -818,46 +791,33 @@
 
 	// cast generated samples 
-	static VssRayContainer samples;
-	samples.clear();
-
-	const bool castDoubleRays = !mPerViewCell;
-	const bool pruneInvalidRays = true;
-	
-#if 0
+	static VssRayContainer vssRays;
+	vssRays.clear();
+
 	castTimer.Entry();
 
-	CastRays(simpleRays, samples, castDoubleRays, pruneInvalidRays);
-	
+	if (0)
+	{
+		const bool castDoubleRays = !mPerViewCell;
+		const bool pruneInvalidRays = false;
+		//const bool pruneInvalidRays = true;
+		CastRays(simpleRays, vssRays, castDoubleRays, pruneInvalidRays);
+	}
+	else
+	{
+		if (0)
+			// casting bundles of 4 rays generated from the simple rays
+			CastRayBundles4(simpleRays, vssRays);
+		else
+			CastRayBundles16(simpleRays, vssRays);
+	}
+
 	castTimer.Exit();
 
-#else
-	static SimpleRayContainer jitteredRays;
-
-	for (size_t i = 0; i < simpleRays.size(); ++ i)
-	{
-		SimpleRay mainRay = simpleRays[i];
-		
-		jitteredRays.clear();
-		jitteredRays.push_back(mainRay);
-
-		generationTimer.Entry();
-
-		GenerateJitteredRays(jitteredRays, mainRay, 15, 0);
-
-		generationTimer.Exit();
-		castTimer.Entry();
-
-		CastRays(jitteredRays, samples, castDoubleRays, pruneInvalidRays);
-
-		castTimer.Entry();
-	}
-#endif
-
 	// add to ray queue
-	EnqueueRays(samples);
+	EnqueueRays(vssRays);
 
 	initialTimer.Exit();
 
-	return (int)samples.size();
+	return (int)vssRays.size();
 }
 
@@ -894,5 +854,8 @@
 		mRayQueue.pop();
 		
-		castSamples += AdaptiveBorderSampling(*ray);
+		if (1)
+			castSamples += AdaptiveBorderSampling(*ray);
+		else
+			castSamples += AdaptiveBorderSamplingOpt(*ray);
 		// note: don't have to delete because handled by ray pool
 	}
@@ -976,5 +939,5 @@
 void GvsPreprocessor::VisualizeViewCells()
 {
-	char str[64]; sprintf(str, "tmp/pass%06d_%04d-", mProcessedViewCells, mPass);
+	char str[64]; sprintf_s(str, "tmp/pass%06d_%04d-", mProcessedViewCells, mPass);
 			
 	// visualization
@@ -1458,4 +1421,6 @@
 	// initialise
 	mGvsStats.mPerViewCellSamples = 0;
+	mGvsStats.mRandomSamples = 0;
+	mGvsStats.mGvsSamples = 0;
 
 	int oldContribution = 0;
@@ -1630,19 +1595,23 @@
 	mGvsStats.mTotalSamples += mGvsStats.mPerViewCellSamples;
 
-	cout << "id: " << mCurrentViewCell->GetId() << " pvs cost: " 
-	     << mGvsStats.mPvsCost << " pvs tri: " << mTrianglePvs.size() << endl;
+	cout << "id: " << mCurrentViewCell->GetId() << endl;
+	cout << "pvs cost "  << mGvsStats.mPvsCost / 1000000 << " M" << endl;
+	cout << "pvs tri: " << mTrianglePvs.size() << endl;
+	cout << "samples: " << mGvsStats.mTotalSamples / 1000000 << " M" << endl;
+	printf("contri: %f tri / K rays\n", 1000.0f * (float)mTrianglePvs.size() / mGvsStats.mTotalSamples);
 
 	//cout << "invalid samples ratio: " << (float)sInvalidSamples / mGvsStats.mPerViewCellSamples << endl;
 
-	float rayTime = rayTimer.TotalTime();
-	float kdPvsTime = kdPvsTimer.TotalTime();
-	float gvsTime = gvsTimer.TotalTime();
-	float initialTime = initialTimer.TotalTime();
-	float intersectionTime = intersectionTimer.TotalTime();
-	float preparationTime = preparationTimer.TotalTime();
-	float mainLoopTime = mainLoopTimer.TotalTime();
-	float contributionTime = contributionTimer.TotalTime();
-	float castTime = castTimer.TotalTime();
-	float generationTime = generationTimer.TotalTime();
+	double rayTime = rayTimer.TotalTime();
+	double kdPvsTime = kdPvsTimer.TotalTime();
+	double gvsTime = gvsTimer.TotalTime();
+	double initialTime = initialTimer.TotalTime();
+	double intersectionTime = intersectionTimer.TotalTime();
+	double preparationTime = preparationTimer.TotalTime();
+	double mainLoopTime = mainLoopTimer.TotalTime();
+	double contributionTime = contributionTimer.TotalTime();
+	double castTime = castTimer.TotalTime();
+	double generationTime = generationTimer.TotalTime();
+	double rawCastTime = mRayCaster->rawCastTimer.TotalTime();
 
 	cout << "handleRay              : " << rayTime << endl;
@@ -1655,17 +1624,18 @@
 	cout << "has contribution       : " << contributionTime << endl;
 	cout << "cast time              : " << castTime << endl;
-	cout << "generation time       : " << generationTime << endl;
+	cout << "generation time        : " << generationTime << endl;
+	cout << "raw cast time          : " << rawCastTime << endl;
 
 	float randomRaysPerSec = mGvsStats.mRandomSamples / (1e6f * initialTime);
 	float gvsRaysPerSec = mGvsStats.mGvsSamples / (1e6f * gvsTime);
+	float rawRaysPerSec = mGvsStats.mPerViewCellSamples / (1e6f * rawCastTime);
 
 	cout << "cast " << randomRaysPerSec << " M random rays/s: " << endl;
 	cout << "cast " << gvsRaysPerSec << " M gvs rays/s: " << endl;
-	
+	cout << "cast " << rawRaysPerSec << " M raw rays/s: " << endl;
+
 	mGvsStats.Stop();
 	mGvsStats.Print(mGvsStatsStream);
-
-}
-
+}
 
 
@@ -1707,6 +1677,6 @@
 
 		mTrianglePvs.push_back(triObj);
-		//mDummyBuffer.push_back(mGvsStats.mGvsRuns);
 		mGenericStats[0] = (int)mTrianglePvs.size();
+
 		return true;
 	}
@@ -1716,53 +1686,121 @@
 
 
-void Preprocessor::CastRays4(SimpleRayContainer &rays, VssRayContainer &vssRays)
-{
-	const int packetSize = 4;
-
-	static int hit_triangles[packetSize];
-	static float dist[packetSize];
-	static Vector3 dirs[packetSize];
-	static Vector3 shifts[packetSize];
-	static Vector3 origs[packetSize];
-
+void GvsPreprocessor::CastRayBundles4(const SimpleRayContainer &rays, VssRayContainer &vssRays)
+{
+	AxisAlignedBox3 box = mViewCellsManager->GetViewSpaceBox();
 
 	SimpleRayContainer::const_iterator it, it_end = rays.end();
 
 	for (it = rays.begin(); it != it_end; ++ it)
-	{
-
-		SimpleRay mainRay = *it;
-
-		origs[0] = mainRay.mOrigin;
-		dirs[0] = mainRay.mDirection;
-
-		for (i = 1; i < packetSize; i++) 
+		CastRayBundle4(*it, vssRays, box);
+}
+
+
+void GvsPreprocessor::CastRayBundle4(const SimpleRay &ray, 
+									 VssRayContainer &vssRays, 
+									 const AxisAlignedBox3 &box)
+{
+	const float pertubDir = 0.01f;
+	static Vector3 pertub;
+
+	static int hit_triangles[4];
+	static float dist[4];
+	static Vector3 dirs[4];
+	static Vector3 origs[4];
+
+	origs[0] = ray.mOrigin;
+	dirs[0] = ray.mDirection;
+
+	for (int i = 1; i < 4; ++ i) 
+	{
+		origs[i] = ray.mOrigin;
+
+		pertub.x = RandomValue(-pertubDir, pertubDir);
+		pertub.y = RandomValue(-pertubDir, pertubDir);
+		pertub.z = RandomValue(-pertubDir, pertubDir);
+
+		dirs[i] = ray.mDirection + pertub;
+		dirs[i] *= 1.0f / Magnitude(dirs[i]);
+	}
+
+	Cast4Rays(dist, dirs, origs, vssRays, box);
+}
+
+
+void GvsPreprocessor::CastRays4(const SimpleRayContainer &rays, VssRayContainer &vssRays)
+{
+	SimpleRayContainer::const_iterator it, it_end = rays.end();
+	
+	static float dist[4];
+	static Vector3 dirs[4];
+	static Vector3 origs[4];
+
+	AxisAlignedBox3 box = mViewCellsManager->GetViewSpaceBox();
+
+	for (it = rays.begin(); it != it_end; ++ it)
+	{
+		for (int i = 1; i < 4; ++ i) 
 		{
-			origs[i] = mainRay.mOrigin;
-
-			const float pertubDir = 0.01f;
-			Vector3 pertub;
-
-			pertub.x = RandomValue(-pertubDir, pertubDir);
-			pertub.y = RandomValue(-pertubDir, pertubDir);
-			pertub.z = RandomValue(-pertubDir, pertubDir);
-
-			dirs[i] = mainRay.mDirection + pertub;
-			const float c = Magnitude(newDir);
-			dirs{i] *= 1.0f / c;
+			origs[i] = rays[i].mOrigin;
+			dirs[i] = rays[i].mDirection;
 		}
-
-		AxisAlignedBox3 box;
-		preprocessor->mRayCaster->CastRaysPacket4(box.Min(),
-												  box.Max(),
-												  origs,
-												  dirs,
-												  hit_triangles,
-												  dist);       
-	}
-
-	DeterminePvsObjects(vssRays);
-}
-
-
-}
+	}
+
+	Cast4Rays(dist, dirs, origs, vssRays, box);
+}
+
+
+void GvsPreprocessor::Cast4Rays(float *dist,
+								Vector3 *dirs,
+								Vector3 *origs,
+								VssRayContainer &vssRays,
+								const AxisAlignedBox3 &box)
+{
+	static int hit_triangles[4];
+
+	mRayCaster->CastRaysPacket4(box.Min(), box.Max(), origs, dirs, hit_triangles,	dist);       
+
+	VssRay *ray;
+
+	for (int i = 0; i < 4; ++ i)
+	{
+		if (hit_triangles[i] != -1)
+		{
+			TriangleIntersectable *triObj = static_cast<TriangleIntersectable *>(mObjects[hit_triangles[i]]);
+
+			if (DotProd(triObj->GetItem().GetNormal(), dirs[i] >= 0))
+				ray = NULL;
+			else
+				ray = mRayCaster->RequestRay(origs[i], origs[i] + dirs[i] * dist[i], NULL, triObj, mPass, 1.0f);
+		}
+		else
+		{
+			ray = NULL;
+		}
+
+		vssRays.push_back(ray);
+	}
+}
+
+
+void GvsPreprocessor::CastRayBundle16(const SimpleRay &ray, VssRayContainer &vssRays)
+{
+	static SimpleRayContainer simpleRays;
+	simpleRays.clear();
+
+	simpleRays.push_back(ray);
+	GenerateJitteredRays(simpleRays, ray, 15, 0);
+
+	CastRays(simpleRays, vssRays, false, false);
+}
+
+
+void GvsPreprocessor::CastRayBundles16(const SimpleRayContainer &rays, VssRayContainer &vssRays)
+{
+	SimpleRayContainer::const_iterator it, it_end = rays.end();
+
+	for (it = rays.begin(); it != it_end; ++ it)
+		CastRayBundle16(*it, vssRays);
+}
+
+}
