Index: /GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp	(revision 1997)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp	(revision 2001)
@@ -9,4 +9,5 @@
 #include "RndGauss.h"
 #include "Mutation.h"
+#include "Exporter.h"
 
 #ifdef GTP_INTERNAL
@@ -20,6 +21,9 @@
 
 #define USE_SIL_TERMINATION_MUTATION 1
+#define MUTATE_ORIGIN 1
 
 #define EVALUATE_MUTATION_STATS 1
+
+#define Q_SEARCH_STEPS 2
 
 void
@@ -51,4 +55,6 @@
   int dummyNcMutations = 0;
   int dummyCMutations = 0;
+
+  int reverseCandidates = 0;
   
   for (int i=0; i < vssRays.size(); i++) {
@@ -86,16 +92,47 @@
 		*mRays[mBufferStart].mRay = *vssRays[i];
 		mRays[mBufferStart].mMutations = 0;
+		mRays[mBufferStart].ResetReverseMutation();
 		//		mRays[mBufferStart] = RayEntry(newRay);
 		mBufferStart++;
 		if (mBufferStart >= mMaxRays)
 		  mBufferStart = 0;
+		
+
 	  }
 	} else {
-#if EVALUATE_MUTATION_STATS 
 	  if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
 		  vssRays[i]->mGeneratorId != -1
 		  ) {
+		// check whether not to store a new backward mutation candidate
+		VssRay *oldRay = mRays[vssRays[i]->mGeneratorId].mRay;
+		VssRay *newRay = vssRays[i];
+
+#define DIST_THRESHOLD 3.0f
+
+		Intersectable *oldObject =
+		  mPreprocessor.mViewCellsManager->GetIntersectable(
+															*oldRay,
+															true);
+		
+
+		if (!mRays[newRay->mGeneratorId].HasReverseMutation()) {
+		  if (DotProd(oldRay->GetDir(), newRay->GetDir()) > 0.0f) {
+		  float oldDist = Magnitude(oldRay->mTermination - newRay->mOrigin);
+		  float newDist = Magnitude(newRay->mTermination - newRay->mOrigin);
+		  
+		  if (newDist < oldDist - oldObject->GetBox().Radius()*DIST_THRESHOLD) {
+			Vector3 origin, termination;
+			if (ComputeReverseMutation(*oldRay, *newRay, origin, termination)) {
+			  mRays[newRay->mGeneratorId].SetReverseMutation(origin, termination);
+			}
+			
+			reverseCandidates++;
+			//mReverseCandidates
+		  }
+		  }
+		}
+#if EVALUATE_MUTATION_STATS 
 		mutationRays++;
-
+		
 		Intersectable *newObject =
 		  mPreprocessor.mViewCellsManager->GetIntersectable(
@@ -103,16 +140,12 @@
 															true);
 
-		Intersectable *oldObject =
-		  mPreprocessor.mViewCellsManager->GetIntersectable(
-															*mRays[vssRays[i]->mGeneratorId].mRay,
-															true);
 
 		if (oldObject == newObject)
 		  dummyNcMutations++;
+#endif
 	  }
-#endif
 	}
   }
-
+  
   if (mutationRays) {
 	cout<<"Mutated rays:"<<mutationRays<<endl;
@@ -121,4 +154,5 @@
 	cout<<"Dummy NC mutations ratio:"<<100.0f*dummyNcMutations/(float)mutationRays<<"%"<<endl;
 	cout<<"Dummy C mutations ratio:"<<100.0f*dummyCMutations/(float)mutationRays<<"%"<<endl;
+	cout<<"Reverse candidates:"<<reverseCandidates<<endl;
   }
 
@@ -259,4 +293,129 @@
 }
 
+bool
+MutationBasedDistribution::ComputeReverseMutation(
+												  const VssRay &oldRay,
+												  const VssRay &newRay,
+												  Vector3 &origin,
+												  Vector3 &termination
+												  )
+{
+  // first reconstruct the termination point
+  Vector3 oldDir = Normalize(oldRay.GetDir());
+  Plane3 oldPlane(oldDir, oldRay.mTermination);
+  
+  termination = oldPlane.FindIntersection(newRay.mOrigin,
+										  newRay.mTermination);
+  
+  // now find the new origin of the ray by casting ray backward from the termination and termining
+  // silhouette point with respect to the occluding object (object containing the newRay termination)
+
+  Plane3 newPlane(oldDir, newRay.mTermination);
+
+  Vector3 oldPivot = newPlane.FindIntersection(oldRay.mOrigin,
+											   oldRay.mTermination);
+
+  Vector3 newPivot = newRay.mTermination;
+  Vector3 line = 2.0f*(oldPivot - newPivot);
+
+  Intersectable *occluder = mPreprocessor.mViewCellsManager->GetIntersectable(
+																			  newRay,
+																			  true);
+  
+  AxisAlignedBox3 box = occluder->GetBox();
+  box.Scale(2.0f);
+  
+  const int packetSize = 4;
+  static int hit_triangles[packetSize];
+  static float dist[packetSize];
+  static Vector3 dirs[packetSize];
+  static Vector3 shifts[packetSize];
+  // now find the silhouette along the line
+  int i;
+  float left = 0.0f;
+  float right = 1.0f;
+  // cast rays to find silhouette ray
+  for (int j=0; j < Q_SEARCH_STEPS; j++) {
+	for (i=0; i < packetSize; i++) {
+	  float r = left + (i+1)*(right-left)/(packetSize+1);
+	  shifts[i] = r*line;
+	  dirs[i] = Normalize(newPivot + shifts[i] - termination );
+	  mlrtaStoreRayASEye4(&termination.x,
+						  &dirs[i].x,
+						  i);
+	}
+	
+	mlrtaTraverseGroupASEye4(&box.Min().x,
+							 &box.Max().x,
+							 hit_triangles,
+							 dist);
+	
+	for (i=0; i < packetSize; i++) {
+	  if (hit_triangles[i] == -1) {
+		// break on first passing ray
+		break;
+	  }
+	}
+	float rr = left + (i+1)*(right-left)/(packetSize+1);
+	float rl = left + i*(right-left)/(packetSize+1);
+	left = rl;
+	right = rr;
+  }
+
+  float t = right;
+  if (right==1.0f)
+	return false;
+  
+  if (i == packetSize)
+	origin = newPivot + right*line;
+  else
+	origin = newPivot + shifts[i];
+
+  if (0) {
+	
+	static VssRayContainer rRays;
+	static int counter = 0;
+	char filename[256];
+
+	if (counter < 50) {
+	  sprintf(filename, "reverse_rays_%03d.x3d", counter++);
+	  
+	  VssRay tRay(origin, termination, NULL, NULL);
+	  rRays.push_back((VssRay *)&oldRay);
+	  rRays.push_back((VssRay *)&newRay);
+	  rRays.push_back(&tRay);
+	  
+	  Exporter *exporter = NULL;
+	  exporter = Exporter::GetExporter(filename);
+	  
+	  exporter->SetFilled();
+	  
+	  Intersectable *occludee = mPreprocessor.mViewCellsManager->GetIntersectable(
+																				  oldRay,
+																				  true);
+	  
+	  exporter->SetForcedMaterial(RgbColor(0,0,1));
+	  exporter->ExportIntersectable(occluder);
+	  exporter->SetForcedMaterial(RgbColor(0,1,0));
+	  exporter->ExportIntersectable(occludee);
+	  exporter->ResetForcedMaterial();
+	  
+	  exporter->SetWireframe();
+	  
+	  
+	  exporter->ExportRays(rRays, RgbColor(1, 0, 0));
+	  delete exporter;
+	  rRays.clear();
+	}
+  }
+
+
+  
+  return true;
+  
+  // now the origin and termination is swapped compred to the generator ray
+  // swap(origin, termination);???
+  // -> perhaps not neccessary as the reverse mutation wil only be used once!
+}
 
 Vector3
@@ -270,11 +429,12 @@
 																)
 {
-  const int packetSize = 16;
-
+  const int packetSize = 4;
   static int hit_triangles[packetSize];
   static float dist[packetSize];
   static Vector3 dirs[packetSize];
+  static Vector3 shifts[packetSize];
   // mutate the 
   float alpha = RandomValue(0.0f, 2.0f*M_PI);
+  //float alpha = vr2.x*2.0f*M_PI;
   
   // direction along which we will mutate the ray
@@ -283,23 +443,34 @@
   //  cout<<line<<endl;
   // create 16 rays along the selected dir
-  
+  int i;
+  float left = 0.0f;
+  float right = radius;
   // cast rays to find silhouette ray
-  for (int i=0; i < packetSize; i++) {
-	dirs[i] = Normalize(ray.mTermination + ((radius/(packetSize - i))*line) - origin );
-	mlrtaStoreRayAS16(&origin.x,
-					  &dirs[i].x,
-					  i);
-  }
-  
-  mlrtaTraverseGroupAS16(&box.Min().x,
-						 &box.Max().x,
-						 hit_triangles,
-						 dist);
-  
-  for (int i=0; i < packetSize; i++) {
-	if (hit_triangles[i] == -1 || !box.IsInside(origin + dist[i]*dirs[i])) {
-	  // break on first passing ray
-	  break;
+  for (int j=0; j < Q_SEARCH_STEPS; j++) {
+	for (i=0; i < packetSize; i++) {
+	  float r = left + (i+1)*(right-left)/(packetSize+1);
+	  shifts[i] = r*line;
+	  dirs[i] = Normalize(ray.mTermination + shifts[i] - origin );
+	  mlrtaStoreRayASEye4(&origin.x,
+						  &dirs[i].x,
+						  i);
 	}
+	
+	mlrtaTraverseGroupASEye4(&box.Min().x,
+							 &box.Max().x,
+							 hit_triangles,
+							 dist);
+	
+	for (i=0; i < packetSize; i++) {
+	  if (hit_triangles[i] == -1) {
+		//	  if (hit_triangles[i] == -1 || !box.IsInside(origin + dist[i]*dirs[i])) {
+		// break on first passing ray
+		break;
+	  }
+	}
+	float rr = left + (i+1)*(right-left)/(packetSize+1);
+	float rl = left + i*(right-left)/(packetSize+1);
+	left = rl;
+	right = rr;
   }
   
@@ -307,10 +478,11 @@
 	//	cerr<<"Warning: hit the same box here should never happen!"<<endl;
 	// shift the ray even a bit more
-	//	cout<<"W"<<i<<endl;
-	return ray.mTermination + (RandomValue(1.0f, 2.0f)*radius)*line;
-  }
-
+	//cout<<"W"<<i<<endl;
+	//	return (RandomValue(1.0f, 1.5f)*radius)*line;
+	return right*line;
+  }
+  
   //  cout<<i<<endl;
-  return dirs[i];
+  return shifts[i];
 }
 
@@ -386,4 +558,15 @@
   mLastIndex = index;
 
+  if (mRays[index].HasReverseMutation()) {
+	//cout<<"R "<<mRays[index].mutatedOrigin<<" "<<mRays[index].mutatedTermination<<endl;
+	sray = SimpleRay(mRays[index].mutatedOrigin,
+					 Normalize(mRays[index].mutatedTermination - mRays[index].mutatedOrigin),
+					 MUTATION_BASED_DISTRIBUTION, 1.0f);
+	sray.mGeneratorId = index;
+	mRays[index].ResetReverseMutation();
+	mRays[index].mMutations++;
+	return true;
+  }
+  
 #if USE_SILHOUETTE_MUTATIONS
   return GenerateSilhouetteMutation(index, sray);
@@ -391,5 +574,4 @@
   return GenerateMutation(index, sray);
 #endif
-  
 }
 
@@ -410,4 +592,8 @@
 
   mRays[index].mHalton.GetNext(4, rr);
+  //  rr[0] = RandomValue(0.0f,0.99999f);
+  //  rr[1] = RandomValue(0.0f,0.99999f);
+  //  rr[2] = RandomValue(0.0f,0.99999f);
+  //  rr[3] = RandomValue(0.0f,0.99999f);
   
   // mutate the origin
@@ -430,15 +616,17 @@
   // optimal for Vienna 0.5f
   
-  float radiusExtension = 0.05f;
+  float radiusExtension = 0.5f;
   //  + mRays[index].mMutations/50.0f;
-
+  
   float mutationRadius = objectRadius*radiusExtension;
-
+  
   // tmp for pompeii
-  mutationRadius = 0.22f;
-  
+  //  mutationRadius = 0.22f;
+
+#if MUTATE_ORIGIN
   origin += ComputeOriginMutation(*ray, U, V,
 								  Vector2(rr[0], rr[1]),
 								  mutationRadius);
+#endif
   
 #if USE_SIL_TERMINATION_MUTATION
@@ -448,5 +636,5 @@
 													  U, V,
 													  Vector2(rr[2], rr[3]),
-													  3.0f*objectRadius);
+													  2.0f*objectRadius);
 #else
   termination += ComputeTerminationMutation(*ray, U, V,
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/run_test2
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/run_test2	(revision 1974)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/run_test2	(revision 2001)
@@ -15,5 +15,5 @@
 
 #SCENE=../data/soda/soda.dat
-#VIEWCELLS=../data/soda/soda5-viewcells-1000.xml
+#VIEWCELLS=../data/soda/soda5-viewcells.xml
 #VIEWCELLS=../data/soda/soda-viewcells-5000.xml
 #VIEWCELLS=../data/soda/soda5-viewcells-single.xml
@@ -49,7 +49,13 @@
 $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
+ -preprocessor_visibility_file=$PREFIX-i-mixed-b1-n4g.xml \
+ -preprocessor_stats=$PREFIX-i-mixed-b1-n4g.log \
+ -preprocessor_histogram_file=$PREFIX-i-mixed-b1-n4g.hlog
+
+
+# e - no origin mutation, silh. tm Q=2, reverse samples
+# f - no origin, q=2, no reverse samples
+# g - gaussian origin, q=2, reverse samples
+
 
 # $COMMAND -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 1999)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 2001)
@@ -1294,5 +1294,6 @@
 Preprocessor::ExportRays(const char *filename,
 						 const VssRayContainer &vssRays,
-						 const int number
+						 const int number,
+						 const bool exportScene
 						 )
 {
@@ -1309,5 +1310,5 @@
   exporter->SetFilled();
   // $$JB temporarily do not export the scene
-  if (0)
+  if (exportScene)
 	exporter->ExportScene(mSceneGraph->GetRoot());
 
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.h
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.h	(revision 1996)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.h	(revision 2001)
@@ -123,8 +123,9 @@
 
 	virtual bool
-		ExportRays(const char *filename,
-		const VssRayContainer &vssRays,
-		const int number
-		);
+	ExportRays(const char *filename,
+			   const VssRayContainer &vssRays,
+			   const int number,
+			   const bool exportScene = false
+			   );
 
 	virtual int
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/Makefile
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 1997)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 2001)
@@ -1,5 +1,5 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (2.00a) (Qt 4.1.2) on: ?t 18. I 20:30:31 2007
+# Generated by qmake (2.00a) (Qt 4.1.2) on: pá 19. I 15:33:39 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/Mutation.h
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.h	(revision 1997)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.h	(revision 2001)
@@ -47,10 +47,13 @@
 	// halton sequence for generatin gmutations of this ray
 	VssRay *mRay;
-	int mMutations;
-	int mUnsuccessfulMutations;
+	short mMutations;
+	short mUnsuccessfulMutations;
 	HaltonSequence mHalton;
 	float mImportance;
 	float mCdf;
 
+	Vector3 mutatedOrigin;
+	Vector3 mutatedTermination;
+	
 	float GetSamplingFactor() const { return mMutations/mImportance; }
 	RayEntry() {}
@@ -59,7 +62,23 @@
 						mUnsuccessfulMutations(0),
 						mHalton(),
-						mImportance(1.0f) {}
+						mImportance(1.0f)
+	{
+	  ResetReverseMutation();
+	}
+
+	void ResetReverseMutation() {
+	  mutatedOrigin = mutatedTermination = Vector3(0,0,0);
+	}
+	bool HasReverseMutation() const {
+	  return !(mutatedOrigin == mutatedTermination);
+	}
+	
+	void SetReverseMutation(const Vector3 &a, const Vector3 &b) {
+	  mutatedOrigin = a;
+	  mutatedTermination = b;
+	}
+
   };
-
+  
 
   Vector3
@@ -90,4 +109,13 @@
 									   );
 
+
+  bool
+  ComputeReverseMutation(
+						 const VssRay &oldRay,
+						 const VssRay &newRay,
+						 Vector3 &origin,
+						 Vector3 &termination
+						 );
+
   RayEntry &GetEntry(const int index) {
 	return mRays[(mBufferStart+index)%mRays.size()];
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/common.cpp
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/common.cpp	(revision 1634)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/common.cpp	(revision 2001)
@@ -115,5 +115,8 @@
     QueryPerformanceCounter(&counter);
     // return in usec
-    return (long) (1000000*counter.QuadPart/(hrFreq.QuadPart));
+    //return (long) (1000000*counter.QuadPart/(hrFreq.QuadPart));
+	// $$
+	// tmp store time in ms
+    return (long) (1000*counter.QuadPart/(hrFreq.QuadPart));
   } else {
     static struct _timeb mtime;
@@ -152,5 +155,7 @@
 TimeDiff(long time1,long time2) // in ms
 {
-  const Real clk=1.0e-3f; // ticks per second
+  //  const Real clk=1.0e-3f; // ticks per second
+  // $$ tmp store time in ms
+  const Real clk=1.0f; // ticks per second
   long t=time2-time1;
   
Index: /GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp
===================================================================
--- /GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp	(revision 2000)
+++ /GTP/trunk/Lib/Vis/Preprocessing/src/SamplingStrategy.cpp	(revision 2001)
@@ -652,7 +652,7 @@
 				  } else
 					if (strcmp(curr, "mutation")==0) {
-						// temp matt: still no mutationstrategy!
-					  //mDistributions.push_back(new MutationBasedDistribution(mPreprocessor));
-					} 
+					  // temp matt: still no mutationstrategy!
+					  mDistributions.push_back(new MutationBasedDistribution(mPreprocessor));
+					}
 	
 	
