Index: GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.cpp	(revision 1928)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/GvsPreprocessor.cpp	(revision 1932)
@@ -46,5 +46,5 @@
 	Debug << "eps: " << mEps << endl;
 
-	mStats.open("gvspreprocessor.log");
+	mGvsStats.open("gvspreprocessor.log");
 }
 
@@ -54,19 +54,20 @@
 										 const VssRay &oldRay)
 {
-	const float dist = Magnitude(oldRay.GetDir());
-	const float newDist = Magnitude(currentRay.GetDir());
- 
-#if 0
-	if ((dist - newDist) > mThresHold)
-#else 
-	// rather take relative distance
-	if ((dist / newDist) > mThreshold)
-#endif
-	{
+	// 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);
+	
+	// distance large => this is likely to be a discontinuity
+	if ((predictedLen - len) > mThreshold)
+	{
+		cout << "d";
 		VssRay *newRay = ReverseSampling(currentRay, hitTriangle, oldRay);
+		
 		// set flag for visualization
 		newRay->mFlags |= VssRay::ReverseSample;
 		
-		// ray is not pushed into the queue => can delete ray
+		// if ray is not further processed => can delete ray
 		if (!HandleRay(newRay))
 			delete newRay;
@@ -81,7 +82,10 @@
 bool GvsPreprocessor::HandleRay(VssRay *vssRay)
 {
+	// compute the contribution to the view cells
 	const bool storeRaysForViz = true;
-	mViewCellsManager->ComputeSampleContribution(*vssRay, true, storeRaysForViz);
-		
+	mViewCellsManager->ComputeSampleContribution(*vssRay, 
+												 true, 
+												 storeRaysForViz);
+
 	// some pvs contribution for this ray?
 	if (vssRay->mPvsContribution > 0)
@@ -92,11 +96,15 @@
 		if (storeRaysForViz)
 		{
+			VssRay *nray = new VssRay(*vssRay);
+			nray->mFlags = vssRay->mFlags;
+
+			// store ray in contributing view cell
 			ViewCellContainer::const_iterator vit, vit_end = vssRay->mViewCells.end();
 			for (vit = vssRay->mViewCells.begin(); vit != vit_end; ++ vit)
-			{
-				VssRay *nray = new VssRay(*vssRay);
+			{			
 				(*vit)->GetOrCreateRays()->push_back(nray);				
 			}
 		}
+
 		//mVssRays.push_back(new VssRay(*vssRay));
         ++ mSampleContriPerPass;
@@ -111,9 +119,8 @@
 /** Creates 3 new vertices for triangle vertex with specified index.
 */
-static void CreateNewVertices(VertexContainer &vertices,
-							  const Triangle3 &hitTriangle,
-							  const VssRay &ray, 
-							  const int index,
-							  const float eps)
+void GvsPreprocessor::CreateDisplacedVertices(VertexContainer &vertices,
+											  const Triangle3 &hitTriangle,
+											  const VssRay &ray, 
+											  const int index) const
 {
 	const int indexU = (index + 1) % 3;
@@ -125,5 +132,5 @@
 	
 	const float len = Magnitude(a);
-
+	
 	const Vector3 dir1 = Normalize(CrossProd(a, b)); //N((pi-xp)×(pi+1- pi));
 	const Vector3 dir2 = Normalize(CrossProd(a, c)); // N((pi-xp)×(pi- pi-1))
@@ -133,9 +140,9 @@
 	// compute the new three hit points
 	// pi, i + 1:  pi+ e·|pi-xp|·di, j
-	const Vector3 pt1 = hitTriangle.mVertices[index] + eps * len * dir1;
+	const Vector3 pt1 = hitTriangle.mVertices[index] + mEps * len * dir1;
 	// pi, i - 1:  pi+ e·|pi-xp|·di, j
-    const Vector3 pt2 = hitTriangle.mVertices[index] + eps * len * dir2;
+    const Vector3 pt2 = hitTriangle.mVertices[index] + mEps * len * dir2;
 	// pi, i:  pi+ e·|pi-xp|·di, j
-	const Vector3 pt3 = hitTriangle.mVertices[index] + eps * len * dir3;
+	const Vector3 pt3 = hitTriangle.mVertices[index] + mEps * len * dir3;
 	
 	vertices.push_back(pt2);
@@ -147,16 +154,18 @@
 void GvsPreprocessor::EnlargeTriangle(VertexContainer &vertices,
 									  const Triangle3 &hitTriangle,
-									  const VssRay &ray)
-{
-	CreateNewVertices(vertices, hitTriangle, ray, 0, mEps);
-	CreateNewVertices(vertices, hitTriangle, ray, 1, mEps);
-	CreateNewVertices(vertices, hitTriangle, ray, 2, mEps);
-}
-
-
-static Vector3 CalcPredictedHitPoint(const VssRay &newRay, 
-									 const Triangle3 &hitTriangle,
-									 const VssRay &oldRay)
-{
+									  const VssRay &ray) const
+{
+	CreateDisplacedVertices(vertices, hitTriangle, ray, 0);
+	CreateDisplacedVertices(vertices, hitTriangle, ray, 1);
+	CreateDisplacedVertices(vertices, hitTriangle, ray, 2);
+}
+
+
+Vector3 GvsPreprocessor::CalcPredictedHitPoint(const VssRay &newRay, 
+											   const Triangle3 &hitTriangle,
+											   const VssRay &oldRay) const
+{
+	// find the intersection of the plane induced by the 
+	// hit triangle with the new ray
 	Plane3 plane(hitTriangle.GetNormal(), hitTriangle.mVertices[0]);
 
@@ -177,16 +186,12 @@
 								   const Vector3 &p1, 
 								   const Vector3 &p2, 
-								   const VssRay &x, 
-								   const VssRay &y,
+								   const VssRay &ray1, 
+								   const VssRay &ray2,
 								   const VssRay &oldRay)
 {
-	// the predicted hitpoint expects to hit the same mesh again
-	const Vector3 predictedHitX = CalcPredictedHitPoint(x, hitTriangle, oldRay);
-	const Vector3 predictedHitY = CalcPredictedHitPoint(y, hitTriangle, oldRay);
-
-	CheckDiscontinuity(x, hitTriangle, oldRay);
-	CheckDiscontinuity(y, hitTriangle, oldRay);
-
-	if (EqualVisibility(x, y))
+	CheckDiscontinuity(ray1, hitTriangle, oldRay);
+	CheckDiscontinuity(ray2, hitTriangle, oldRay);
+
+	if (EqualVisibility(ray1, ray2) || (Magnitude(p1 - p2) <= mEps))
 	{
 		return 0;
@@ -194,11 +199,15 @@
 	else
 	{
-		cout << "s";
+		// the new subdivision point
 		const Vector3 p = (p1 + p2) * 0.5f;
+	
+		//cout << "tobj " << ray1.mTerminationObject << " " << ray2.mTerminationObject << " " << p1 << " " << p2 << endl;
+		//cout << "term " << ray1.mTermination << " " << ray2.mTermination << endl;
+
+		// cast ray into the new point
 		SimpleRay sray(oldRay.mOrigin, p - oldRay.mOrigin, SamplingStrategy::GVS, 1.0f);
 	
-		// cast ray into the new subdivision point
-		VssRay *newRay = mRayCaster->CastRay(sray, mViewCellsManager->GetViewSpaceBox(), false);
-		
+		VssRay *newRay = mRayCaster->CastRay(sray, mViewCellsManager->GetViewSpaceBox());
+
 		if (!newRay) return 0;
 
@@ -209,11 +218,12 @@
 		
 		// subdivide further
-		const int s1 = SubdivideEdge(hitTriangle, p1, p, x, *newRay, oldRay);
-		const int s2 = SubdivideEdge(hitTriangle, p, p2, *newRay, y, oldRay);
-				
+		const int samples1 = SubdivideEdge(hitTriangle, p1, p, ray1, *newRay, oldRay);
+		const int samples2 = SubdivideEdge(hitTriangle, p, p2, *newRay, ray2, oldRay);
+			
+		// this ray will not be further processed
 		if (!enqueued)
 			delete newRay;
 		
-		return s1 + s2 + 1;
+		return samples1 + samples2 + 1;
 	}
 }
@@ -222,5 +232,4 @@
 int GvsPreprocessor::AdaptiveBorderSampling(const VssRay &currentRay)
 {
-	cout << "a";
 	Intersectable *tObj = currentRay.mTerminationObject;
 	Triangle3 hitTriangle;
@@ -245,4 +254,6 @@
 	simpleRays.reserve(9);
 
+	//cout << "currentRay: " << currentRay.mOrigin << " dir: " << currentRay.GetDir() << endl;
+
 	VertexContainer::const_iterator vit, vit_end = enlargedTriangle.end();
 
@@ -252,13 +263,14 @@
 		SimpleRay sr(currentRay.GetOrigin(), rayDir, SamplingStrategy::GVS, 1.0f);
 		simpleRays.AddRay(sr);
-	}
-
-	// visualize enlarged triangles
+
+		//cout << "new: " << sr.mOrigin << " dist: " << sr.mDirection << endl;
+	}
+
 	if (0)
 	{
+		// visualize enlarged triangles
 		VizStruct dummy;
 		dummy.enlargedTriangle = new Polygon3(enlargedTriangle);
 		dummy.originalTriangle = hitTriangle;
-		//dummy.ray = new VssRay(currentRay);
 		vizContainer.push_back(dummy);
 	}
@@ -266,5 +278,12 @@
 	// cast rays to triangle vertices and determine visibility
 	VssRayContainer vssRays;
-	CastRays(simpleRays, vssRays, false, false);
+
+	// don't cast double rays as we need only the forward rays
+	const bool castDoubleRays = false;
+	// cannot prune invalid rays because we have to 
+	// compare adjacent  rays.
+	const bool pruneInvalidRays = false;
+
+	CastRays(simpleRays, vssRays, castDoubleRays, pruneInvalidRays);
 
 	// set flags
@@ -275,5 +294,5 @@
 	}
 
-	// add to ray queue
+	// handle rays
 	EnqueueRays(vssRays);
 	
@@ -282,25 +301,25 @@
 
     // recursivly subdivide each edge
-	for (int i = 0; 1 && (i < n); ++ i)
-	{
-		castRays += SubdivideEdge(
-			hitTriangle,
-			enlargedTriangle[i], 
-			enlargedTriangle[(i + 1) % n], 
-			*vssRays[i], 
-			*vssRays[(i + 1) % n],
-			currentRay);
+	for (int i = 0; i < n; ++ i)
+	{
+		castRays += SubdivideEdge(hitTriangle,
+								  enlargedTriangle[i], 
+								  enlargedTriangle[(i + 1) % n], 
+								  *vssRays[i], 
+								  *vssRays[(i + 1) % n],
+								  currentRay);
 	}
 
 	mBorderSamples += castRays;
+
 	return castRays;
 }
 
 
-static Vector3 GetPassingPoint(const VssRay &currentRay,
-							   const Triangle3 &hitTriangle,
-							   const VssRay &oldRay)
-{
-	// intersect triangle plane with plane spanned by current samples
+Vector3 GvsPreprocessor::GetPassingPoint(const VssRay &currentRay,
+										 const Triangle3 &hitTriangle,
+										 const VssRay &oldRay) const
+{
+	//-- intersect triangle plane with plane spanned by current samples
 	Plane3 plane(currentRay.GetOrigin(), currentRay.GetTermination(), oldRay.GetTermination());
 	Plane3 triPlane(hitTriangle.GetNormal(), hitTriangle.mVertices[0]);
@@ -310,5 +329,5 @@
 	// Evaluate new hitpoint just outside the triangle
 	const float factor = 0.95f;
-	float t = triPlane.FindT(intersectLine);
+	const float t = triPlane.FindT(intersectLine);
 	const Vector3 newPoint = intersectLine.mOrigin + t * factor * intersectLine.mDirection;
 
@@ -329,7 +348,8 @@
 	const Vector3 predicted = CalcPredictedHitPoint(currentRay, hitTriangle, oldRay);
 
-	//-- Construct the mutated ray with xnew,dir = predicted(x)- pnew 
-	//-- as direction vector
+	//-- Construct the mutated ray with xnew,
+	//-- dir = predicted(x)- pnew as direction vector
 	const Vector3 newDir = predicted - newPoint ;
+
 	// take xnew,p = intersect(viewcell, line(pnew, predicted(x)) as origin ?
 	// difficult to say!!
@@ -438,5 +458,4 @@
 	else if (0)
 	{	
-		//-- load view cells from file
 		//-- test successful view cells loading by exporting them again
 		VssRayContainer dummies;
@@ -451,8 +470,8 @@
 		////////
 		//-- stats
+
 		cout << "\nPass " << mPass << " #samples: " << castSamples << " of " << mTotalSamples << endl;
-
-		//mVssRays.PrintStatistics(mStats);
-		mStats 
+		//mVssRays.PrintStatistics(mGvsStats);
+		mGvsStats 
 			<< "#Pass\n" << mPass << endl
 			<< "#Time\n" << TimeDiff(startTime, GetTime())*1e-3 << endl
@@ -464,5 +483,5 @@
 			<< "#GvsRuns\n" << mGvsPass << endl;
 
-		mViewCellsManager->PrintPvsStatistics(mStats);
+		mViewCellsManager->PrintPvsStatistics(mGvsStats);
 
 		char str[64]; sprintf(str, "tmp/pass%04d-", mPass);
@@ -489,5 +508,5 @@
 	}
 
-	cout << 2 * castSamples / (1e3f * TimeDiff(startTime, GetTime())) << "M rays/s" << endl;
+	cout << "cast " << 2 * castSamples / (1e3f * TimeDiff(startTime, GetTime())) << "M rays/s" << endl;
 	Visualize();
 
