Index: GTP/trunk/Lib/Vis/Preprocessing/src/DifferenceSampling.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/DifferenceSampling.h	(revision 2546)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/DifferenceSampling.h	(revision 2546)
@@ -0,0 +1,45 @@
+#ifndef _DifferenceSampling_H__
+#define _DifferenceSampling_H__
+
+#include <vector>
+//
+
+#include "common.h"
+#include "Halton.h"
+
+namespace GtpVisibilityPreprocessor {
+
+
+
+/** This sampling strategy utilizes the spatial filter to find 
+	areas that are likely to be visible.
+
+	First the filter is applied. The difference set of objects (filtered pvs - unfiltered pvs) gives you the objects 
+	which are not found visible yet but have a good probability that they are in fact visible, based on the current 
+	state of the sampling and the size of the filter. In fact the difference set is a probably visible set.
+    
+	Idea:
+	We apply the filter for n times, each time with a higher filter radius. The difference set gives you the probable visible set. 
+	Then we use a certain #random samples for sampling the view cell border + the object border. 
+	We can use different filter sizes to estimate the probablility distribution. The objects shouldreceived samples 
+	roughly according to their probability of becoming invisible.
+*/	
+class DifferenceSampling: public SamplingStrategy
+{
+public:
+	DifferenceSampling(Preprocessor &preprocessor);
+	/** Compute the difference set between filtered and unfiltered pvs.
+	*/
+	void ComputeDifferenceSet(ViewCell *vc, ObjectContainer &differenceSet);
+
+private:
+
+  virtual bool GenerateSample(SimpleRay &ray);
+  
+};
+
+
+}
+
+#endif
+
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.vcproj
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.vcproj	(revision 2538)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.vcproj	(revision 2546)
@@ -339,4 +339,10 @@
 			<File
 				RelativePath=".\DefaultPvs.h">
+			</File>
+			<File
+				RelativePath=".\DifferenceSampling.cpp">
+			</File>
+			<File
+				RelativePath=".\DifferenceSampling.h">
 			</File>
 			<File
@@ -717,4 +723,10 @@
 			</File>
 			<File
+				RelativePath=".\ReverseGvs.cpp">
+			</File>
+			<File
+				RelativePath=".\ReverseGvs.h">
+			</File>
+			<File
 				RelativePath="..\src\RssPreprocessor.cpp">
 			</File>
@@ -1041,8 +1053,8 @@
 			</File>
 			<File
+				RelativePath=".\sparsehash\src\windows\hash_fun.h">
+			</File>
+			<File
 				RelativePath=".\sparsehash\src\google\sparsehash\hash_fun.h">
-			</File>
-			<File
-				RelativePath=".\sparsehash\src\windows\hash_fun.h">
 			</File>
 			<File
Index: GTP/trunk/Lib/Vis/Preprocessing/src/DifferenceSampling.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/DifferenceSampling.cpp	(revision 2546)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/DifferenceSampling.cpp	(revision 2546)
@@ -0,0 +1,25 @@
+#include "SamplingStrategy.h"
+#include "Intersectable.h"
+#include "DifferenceSampling.h"
+#include "ViewCell.h"
+
+
+#ifdef GTP_INTERNAL
+#include "ArchModeler2MLRT.hxx"
+#endif
+
+
+
+namespace GtpVisibilityPreprocessor {
+
+
+DifferenceSampling::DifferenceSampling(Preprocessor &preprocessor):
+SamplingStrategy(preprocessor)
+{}
+
+
+void DifferenceSampling::ComputeDifferenceSet(ViewCell *vc, ObjectContainer &differenceSet)
+{
+}
+
+}
Index: GTP/trunk/Lib/Vis/Preprocessing/src/VspTree.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/VspTree.cpp	(revision 2544)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/VspTree.cpp	(revision 2546)
@@ -3989,7 +3989,4 @@
 	queue<MyTraversalData> tStack;
 
-	// sort objects by their id
-	//sort(objects.begin(), objects.end(), ilt);
-
 	// hack: we make a new root
 	DEL_PTR(mRoot);
Index: GTP/trunk/Lib/Vis/Preprocessing/src/ReverseGvs.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/ReverseGvs.h	(revision 2546)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/ReverseGvs.h	(revision 2546)
@@ -0,0 +1,171 @@
+#ifndef _ReverserGvs_H__
+#define _ReverseGvs_H__
+
+#include <vector>
+//
+
+#include "common.h"
+#include "Halton.h"
+
+namespace GtpVisibilityPreprocessor {
+
+class Vector2;
+class Vector3;
+class VssRay;
+class Preprocessor;
+struct SimpleRay;
+class SimpleRayContainer;
+
+struct VssRayContainer;
+
+
+/** This algorithm can be best described as Reverse Guided Visibility Sampling. We exchange the roles of 
+	objects (triangles in gvs) and view cells, in order to make Gvs feasible for
+	global sampling (where all view cells are updated in parallel).
+	
+	Algorithm idea:
+
+	We loop through the objects (instead of the view cells). For each object, there s a directional 
+	random sampling step, which will find some view cells. This is again dual to Gvs.
+	
+	The view cells found in the first step are used as a seed in order to grow visibility regions towards 
+	neighbouring view cells. 
+	We use guided samples from the object towards these neighbouring view cells. Subsequently the algorithm 
+	would explore the boundaries of the view cells where the object can still be seen. This is dual to
+	the border sampling step of Gvs.
+	
+	There will also be a reverse sampling step. We take a new point on the object, abd connect it to a view cell 
+	found invisible by the border sampling. This would only make sense if the objects are reasonable large.
+	
+	The new algorithm would fully explore the global nature of the sampling. The progressive quality of the 
+	algorithm could be more emphasised by looping through the objects in some kind of importance fashion, 
+	where the importance comes from an initial estimate of the #view cells the object can be seen from.
+*/
+class ReverseGvs: public SamplingStrategy
+{
+public:
+  ReverseGvs(Preprocessor &preprocessor);
+  virtual void Update(VssRayContainer &vssRays);
+
+  virtual bool RequiresRays() { return true; }
+
+private:
+  virtual bool GenerateSample(SimpleRay &ray);
+
+  bool
+  GenerateMutation(const int index, SimpleRay &ray);
+
+  bool
+  GenerateMutationCandidate(const int index,
+							SimpleRay &ray,
+							Intersectable *object,
+							const AxisAlignedBox3 &box
+							);
+
+
+  
+  struct RayEntry {
+	// halton sequence for generatin gmutations of this ray
+	VssRay *mRay;
+	short mMutations;
+	short mUnsuccessfulMutations;
+	//	Halton<4> mHalton;
+	HaltonSequence mHalton;
+	float mImportance;
+	float mCdf;
+
+	Vector3 mutatedOrigin;
+	Vector3 mutatedTermination;
+	
+	float GetSamplingFactor() const { return mMutations/mImportance; }
+	
+	friend bool operator<(const RayEntry &a, const RayEntry &b) {
+	  return a.GetSamplingFactor() > b.GetSamplingFactor();
+	}
+	
+	RayEntry() {}
+	RayEntry(VssRay *r):mRay(r),
+						mMutations(0),
+						mUnsuccessfulMutations(0),
+						mHalton(),
+						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
+  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
+							 );
+
+
+  Vector3
+  ComputeSilhouetteTerminationMutation(const VssRay &ray,
+									   const Vector3 &origin,
+									   const AxisAlignedBox3 &box,
+									   const Vector3 &U,
+									   const Vector3 &V,
+									   const float radius
+									   );
+
+
+  bool
+  ComputeReverseMutation(
+						 const VssRay &oldRay,
+						 const VssRay &newRay,
+						 Vector3 &origin,
+						 Vector3 &termination
+						 );
+
+  RayEntry &GetEntry(const int index);
+
+  
+  vector<RayEntry> mRays;
+  int mBufferStart;
+  int mLastIndex;
+
+
+  int mMaxRays;
+  float mMutationRadiusOrigin;
+  float mMutationRadiusTermination;
+  bool mUseReverseSamples;
+  float mReverseSamplesDistance;
+  bool mUseSilhouetteSamples;
+  int mSilhouetteSearchSteps;
+  float mSilhouetteProb;
+  bool mUsePassImportance;
+  bool mUseUnsuccCountImportance;
+
+};
+
+
+}
+
+#endif
+
Index: GTP/trunk/Lib/Vis/Preprocessing/src/ReverseGvs.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/ReverseGvs.cpp	(revision 2546)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/ReverseGvs.cpp	(revision 2546)
@@ -0,0 +1,807 @@
+#include "SamplingStrategy.h"
+#include "Ray.h"
+#include "Intersectable.h"
+#include "Preprocessor.h"
+#include "ViewCellsManager.h"
+#include "AxisAlignedBox3.h"
+#include "RssTree.h"
+#include "Vector2.h"
+#include "RndGauss.h"
+#include "Exporter.h"
+#include "Environment.h"
+#include "RayCaster.h"
+#include "ReverseGvs.h"
+
+
+#ifdef GTP_INTERNAL
+#include "ArchModeler2MLRT.hxx"
+#endif
+
+namespace GtpVisibilityPreprocessor {
+
+#define MUTATION_USE_CDF 0
+#define EVALUATE_MUTATION_STATS 1
+#define SORT_RAY_ENTRIES 1
+
+// use avg ray contribution as importance
+// if 0 the importance is evaluated from the succ of mutations
+#define USE_AVG_CONTRIBUTION 1
+
+ReverseGvs::RayEntry &
+ReverseGvs::GetEntry(const int index)
+{
+#if SORT_RAY_ENTRIES
+  return mRays[index];
+#else
+  return mRays[(mBufferStart+index)%mRays.size()];
+#endif
+}
+
+void
+ReverseGvs::Update(VssRayContainer &vssRays)
+{
+  //  for (int i=0; i < mRays.size(); i++)
+  //	cout<<mRays[i].mMutations<<" ";
+  //  cout<<endl;
+  cerr<<"Mutation update..."<<endl;
+  cerr<<"rays = "<<mRays.size()<<endl;
+  if ((int)mRays.size()) {
+	cerr<<"Oversampling factors = "<<
+	  GetEntry(0).mMutations<<" "<<
+	  GetEntry(1).mMutations<<" "<<
+	  GetEntry(2).mMutations<<" "<<
+	  GetEntry(3).mMutations<<" "<<
+	  GetEntry(4).mMutations<<" "<<
+	  GetEntry(5).mMutations<<" ... "<<
+	  GetEntry((int)mRays.size()-6).mMutations<<" "<<
+	  GetEntry((int)mRays.size()-5).mMutations<<" "<<
+	  GetEntry((int)mRays.size()-4).mMutations<<" "<<
+	  GetEntry((int)mRays.size()-3).mMutations<<" "<<
+	  GetEntry((int)mRays.size()-2).mMutations<<" "<<
+	  GetEntry((int)mRays.size()-1).mMutations<<endl;
+  }
+  int contributingRays = 0;
+
+  int mutationRays = 0;
+  int dummyNcMutations = 0;
+  int dummyCMutations = 0;
+
+  int reverseCandidates = 0;
+
+#if 0
+  sort(mRays.begin(), mRays.end());
+  // reset the start of the buffer
+  mBufferStart = 0;
+#endif
+
+  for (int i=0; i < vssRays.size(); i++) {
+	if (vssRays[i]->mPvsContribution) {
+	  // reset the counter of unsuccsseful mutation for a generating ray (if it exists)
+	  if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
+		  vssRays[i]->mGeneratorId != -1
+		  ) {
+		mRays[vssRays[i]->mGeneratorId].mUnsuccessfulMutations = 0;
+#if EVALUATE_MUTATION_STATS 
+		mutationRays++;
+		
+		Intersectable *newObject = vssRays[i]->mTerminationObject;
+
+		Intersectable *oldObject = mRays[vssRays[i]->mGeneratorId].mRay->mTerminationObject;
+		
+		if (oldObject == newObject)
+		  dummyCMutations++;
+#endif
+	  }
+	  contributingRays++;
+	  if (mRays.size() < mMaxRays) {
+		VssRay *newRay = new VssRay(*vssRays[i]);
+		// add this ray
+		newRay->Ref();
+		mRays.push_back(RayEntry(newRay));
+	  } else {
+		// unref the old ray
+		*mRays[mBufferStart].mRay = *vssRays[i];
+		mRays[mBufferStart].mMutations = 0;
+		mRays[mBufferStart].mUnsuccessfulMutations = 0;
+		mRays[mBufferStart].ResetReverseMutation();
+		//		mRays[mBufferStart] = RayEntry(newRay);
+		mBufferStart++;
+		if (mBufferStart >= mMaxRays)
+		  mBufferStart = 0;
+	  }
+	} else {
+	  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];
+
+
+		Intersectable *oldObject = oldRay->mTerminationObject;
+		
+
+		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()*mReverseSamplesDistance) {
+			Vector3 origin, termination;
+			if (ComputeReverseMutation(*oldRay, *newRay, origin, termination)) {
+			  mRays[newRay->mGeneratorId].SetReverseMutation(origin, termination);
+			}
+			
+			reverseCandidates++;
+			//mReverseCandidates
+		  }
+		  }
+		}
+#if EVALUATE_MUTATION_STATS 
+		mutationRays++;
+		
+		Intersectable *newObject = vssRays[i]->mTerminationObject;
+
+
+		if (oldObject == newObject)
+		  dummyNcMutations++;
+#endif
+	  }
+	}
+  }
+  
+  if (mutationRays) {
+	cout<<"Mutated rays:"<<mutationRays<<endl;
+	cout<<"Dummy mutations ratio:"<<100.0f*(dummyCMutations + dummyNcMutations)/
+	  (float)mutationRays<<"%"<<endl;
+	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;
+  }
+  
+  float pContributingRays = contributingRays/(float)vssRays.size();
+  
+  cout<<"Percentage of contributing rays:"<<pContributingRays<<endl;
+  
+  if (mUseUnsuccCountImportance) {
+	// use unsucc mutation samples as feedback on importance
+	for (int i=0; i < mRays.size(); i++) {
+	  const float  minImportance = 0.1f;
+	  const int minImportanceSamples = 20;
+	  mRays[i].mImportance = minImportance +
+		(1-minImportance)*exp(-3.0f*mRays[i].mUnsuccessfulMutations/minImportanceSamples);
+	}
+  } else {
+	float importance = 1.0f;
+	if (mUsePassImportance) 
+	  importance = 1.0f/(pContributingRays + 1e-5);
+	// set this values for last contributingRays
+	int index = mBufferStart - 1;
+	
+	for (int i=0; i < contributingRays; i++, index--) {
+	  if (index < 0)
+		index = (int)mRays.size()-1;
+	  mRays[index].mImportance = importance;
+	}
+  }
+  
+#if SORT_RAY_ENTRIES
+  long t1 = GetTime();
+  sort(mRays.begin(), mRays.end());
+  // reset the start of the buffer
+  mBufferStart = 0;
+  mLastIndex = (int)mRays.size();
+  cout<<"Mutation candidates sorted in "<<TimeDiff(t1, GetTime())<<" ms."<<endl;
+#endif
+  
+#if MUTATION_USE_CDF
+  // compute cdf
+  mRays[0].mCdf = mRays[0].mImportance/(mRays[0].mMutations+1);
+  for (int i=1; i < mRays.size(); i++) 
+	mRays[i].mCdf = mRays[i-1].mCdf + mRays[i].mImportance/(mRays[i].mMutations+1);
+  
+  float scale = 1.0f/mRays[i-1].mCdf;
+  for (i=0; i < mRays.size(); i++) {
+	mRays[i].mCdf *= scale;
+  }
+#endif
+  
+  cout<<"Importance = "<<
+	GetEntry(0).mImportance<<" "<<
+	GetEntry((int)mRays.size()-1).mImportance<<endl;
+
+  cout<<"Sampling factor = "<<
+	GetEntry(0).GetSamplingFactor()<<" "<<
+	GetEntry((int)mRays.size()-1).GetSamplingFactor()<<endl;
+
+  cerr<<"Mutation update done."<<endl;
+}
+
+
+Vector3
+ReverseGvs::ComputeOriginMutation(const VssRay &ray,
+												 const Vector3 &U,
+												 const Vector3 &V,
+												 const Vector2 vr2,
+												 const float radius
+												 )
+{
+#if 0
+  Vector3 v;
+  if (d.DrivingAxis() == 0)
+	v = Vector3(0, r[0]-0.5f, r[1]-0.5f);
+  else
+	if (d.DrivingAxis() == 1)
+	  v = Vector3(r[0]-0.5f, 0, r[1]-0.5f);
+	else
+	  v = Vector3(r[0]-0.5f, r[1]-0.5f, 0);
+  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
+ReverseGvs::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)
+	v = Vector3(0, r[2]-0.5f, r[3]-0.5f);
+  else
+	if (d.DrivingAxis() == 1)
+	  v = Vector3(r[2]-0.5f, 0, r[3]-0.5f);
+	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
+ReverseGvs::ComputeReverseMutation(
+												  const VssRay &oldRay,
+												  const VssRay &newRay,
+												  Vector3 &origin,
+												  Vector3 &termination
+												  )
+{
+#ifdef GTP_INTERNAL
+
+  // 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 = newRay.mTerminationObject;
+  
+  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 < mSilhouetteSearchSteps; 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 = 
+		oldRay.mTerminationObject;
+	  
+	  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();
+	}
+  }
+
+#else
+cerr << "warning: reverse mutation not supported!" << endl;
+#endif
+  
+  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
+ReverseGvs::ComputeSilhouetteTerminationMutation(const VssRay &ray,
+																const Vector3 &origin,
+																const AxisAlignedBox3 &box,
+																const Vector3 &U,
+																const Vector3 &V,
+																const float radius
+																)
+{
+#ifdef GTP_INTERNAL
+
+  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, Real(2.0*M_PI));
+  //float alpha = vr2.x*2.0f*M_PI;
+  
+  // direction along which we will mutate the ray
+  Vector3 line = sin(alpha)*U + cos(alpha)*V;
+  
+  //  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 j=0; j < mSilhouetteSearchSteps; 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;
+  }
+  
+  if (i == packetSize) {
+	//	cerr<<"Warning: hit the same box here should never happen!"<<endl;
+	// shift the ray even a bit more
+	//cout<<"W"<<i<<endl;
+	//	return (RandomValue(1.0f, 1.5f)*radius)*line;
+	return right*line;
+  }
+  
+  //  cout<<i<<endl;
+  return shifts[i];
+
+#else
+	cerr << "warning: shiluette mutation not supported" << endl;
+	return Vector3(0, 0, 0);
+#endif
+  
+}
+
+
+bool
+ReverseGvs::GenerateSample(SimpleRay &sray)
+{
+
+  if (mRays.size() == 0) {
+	float rr[5];
+	// use direction based distribution
+	Vector3 origin, direction;
+
+	for (int i=0; i < 5; i++)
+	  rr[i] = RandomValue(0,1);
+	
+	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);
+	sray.mGeneratorId = -1;
+
+	return true;
+  } 
+
+  int index;
+
+#if !MUTATION_USE_CDF
+#if SORT_RAY_ENTRIES
+  index = mLastIndex - 1;
+  if (index < 0 || index >= mRays.size()-1) {
+	index = (int)mRays.size() - 1;
+  } else
+	if (
+		mRays[index].GetSamplingFactor() >= mRays[mLastIndex].GetSamplingFactor()) {
+	  // make another round
+
+	  //	  cout<<"R2"<<endl;
+	  //	  cout<<mLastIndex<<endl;
+	  //	  cout<<index<<endl;
+	  index = (int)mRays.size() - 1;
+	}
+#else
+  // 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].mMutations > mRays[mLastIndex].mMutations)
+	  //		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();
+  }
+#endif
+#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=(int)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;
+
+  mLastIndex = index;
+  //  Debug<<index<<" "<<mRays[index].GetSamplingFactor()<<endl;
+  
+  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++;
+	mRays[index].mUnsuccessfulMutations++;
+
+	return true;
+  }
+  
+  return GenerateMutation(index, sray);
+}
+
+
+
+
+
+bool
+ReverseGvs::GenerateMutationCandidate(const int index,
+													 SimpleRay &sray,
+													 Intersectable *object,
+													 const AxisAlignedBox3 &box
+													 )
+{
+  float rr[4];
+
+  VssRay *ray = mRays[index].mRay;
+
+  //  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
+  Vector3 d = ray->GetDir();
+  
+  float objectRadius = box.Radius();
+  //  cout<<objectRadius<<endl;
+  if (objectRadius < Limits::Small)
+	return false;
+  
+  // 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();
+
+  
+  // use probabilitistic approach to decide for the type of mutation
+  float a = RandomValue(0.0f,1.0f);
+  bool bidirectional = true;
+  
+  if (mUseSilhouetteSamples && a < mSilhouetteProb) {
+	
+	termination += ComputeSilhouetteTerminationMutation(*ray,
+														origin,
+														box,
+														U, V,
+														2.0f*objectRadius);
+	bidirectional = false;
+  } else {
+	mRays[index].mHalton.GetNext(4, rr);
+
+	// fuzzy random mutation
+	origin += ComputeOriginMutation(*ray, U, V,
+									Vector2(rr[0], rr[1]),
+									mMutationRadiusOrigin*objectRadius);
+	
+	termination += ComputeTerminationMutation(*ray, U, V,
+											  Vector2(rr[2], rr[3]),
+											  mMutationRadiusTermination*objectRadius);
+  }
+
+  Vector3 direction = termination - origin;
+  
+  if (Magnitude(direction) < Limits::Small)
+	return false;
+  
+  // shift the origin a little bit
+  origin += direction*0.5f;
+  
+  direction.Normalize();
+  
+  // $$ jb the pdf is yet not correct for all sampling methods!
+  const float pdf = 1.0f;
+  
+  sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
+  sray.mGeneratorId = index;
+  sray.SetBidirectional(bidirectional);
+  
+  return true;
+}
+
+bool
+ReverseGvs::GenerateMutation(const int index, SimpleRay &sray)
+{
+  VssRay *ray = mRays[index].mRay;
+
+  Intersectable *object = ray->mTerminationObject;
+  
+  AxisAlignedBox3 box = object->GetBox();
+  
+  if (GenerateMutationCandidate(index, sray, object, box)) {
+    mRays[index].mMutations++;
+	mRays[index].mUnsuccessfulMutations++;
+
+	return true;
+  }
+  return false;
+}
+  
+
+
+ReverseGvs::ReverseGvs(Preprocessor &preprocessor
+													 ) :
+  SamplingStrategy(preprocessor)
+{
+  mType = MUTATION_BASED_DISTRIBUTION;
+  mBufferStart = 0;
+  mLastIndex = 0;
+
+  Environment::GetSingleton()->GetIntValue("Mutation.bufferSize",
+										   mMaxRays);
+  
+  mRays.reserve(mMaxRays);
+  
+  Environment::GetSingleton()->GetFloatValue("Mutation.radiusOrigin",
+											 mMutationRadiusOrigin);
+
+  Environment::GetSingleton()->GetFloatValue("Mutation.radiusTermination",
+											 mMutationRadiusTermination);
+  
+  bool useEnhancedFeatures;
+
+#ifdef GTP_INTERNAL
+  int rayCaster;
+  Environment::GetSingleton()->GetIntValue("Preprocessor.rayCastMethod",
+										   rayCaster);
+  useEnhancedFeatures = (rayCaster ==
+						 RayCaster::INTEL_RAYCASTER);
+#else
+  useEnhancedFeatures = false;
+#endif
+  
+  if (useEnhancedFeatures) {
+	Environment::GetSingleton()->GetBoolValue("Mutation.useReverseSamples",
+											  mUseReverseSamples);
+	
+	Environment::GetSingleton()->GetFloatValue("Mutation.reverseSamplesDistance",
+
+											   mReverseSamplesDistance);
+
+	Environment::GetSingleton()->GetFloatValue("Mutation.silhouetteProb",
+											   mSilhouetteProb);
+
+  } else {
+	mUseReverseSamples = false;
+	mReverseSamplesDistance = 1e20f;
+	mSilhouetteProb = 0.0f;
+  }
+  
+  Environment::GetSingleton()->GetBoolValue("Mutation.useSilhouetteSamples",
+											mUseSilhouetteSamples);
+
+  Environment::GetSingleton()->GetIntValue("Mutation.silhouetteSearchSteps",
+										   mSilhouetteSearchSteps);
+  
+  
+  Environment::GetSingleton()->GetBoolValue("Mutation.usePassImportance",
+											mUsePassImportance);
+  
+
+  Environment::GetSingleton()->GetBoolValue("Mutation.useUnsuccCountImportance",
+											mUseUnsuccCountImportance);
+
+  
+  
+}
+
+
+
+}
