Index: GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp	(revision 2590)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp	(revision 2592)
@@ -94,5 +94,6 @@
 
 		Intersectable *oldObject = mRays[vssRays[i]->mGeneratorId].mRay->mTerminationObject;
-		
+		// the ray generated a contribution although it hits the same object
+		// mark this using a counter
 		if (oldObject == newObject)
 		  dummyCMutations++;
@@ -106,5 +107,5 @@
 		mRays.push_back(RayEntry(newRay));
 	  } else {
-		// unref the old ray
+		// unref the old ray and add the new ray into the mutation buffer
 		*mRays[mBufferStart].mRay = *vssRays[i];
 		mRays[mBufferStart].mMutations = 0;
@@ -117,15 +118,17 @@
 	  }
 	} else {
+	  // the ray did not generate any contribution
 	  if (vssRays[i]->mDistribution == MUTATION_BASED_DISTRIBUTION &&
 		  vssRays[i]->mGeneratorId != -1
 		  ) {
 		// check whether not to store a new backward mutation candidate
+		// this happens if the new ray is occluded significantly closer than
+		// the generator ray
 		VssRay *oldRay = mRays[vssRays[i]->mGeneratorId].mRay;
 		VssRay *newRay = vssRays[i];
 
-
 		Intersectable *oldObject = oldRay->mTerminationObject;
-		
-
+
+		// only allow one reverse mutation per generator ray
 		if (!mRays[newRay->mGeneratorId].HasReverseMutation()) {
 		  if (DotProd(oldRay->GetDir(), newRay->GetDir()) > 0.0f) {
@@ -149,5 +152,4 @@
 		Intersectable *newObject = vssRays[i]->mTerminationObject;
 
-
 		if (oldObject == newObject)
 		  dummyNcMutations++;
@@ -163,10 +165,10 @@
 	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;
+	cout<<"Reverse candidates:"<<100.0f*reverseCandidates/(float)mutationRays<<endl;
   }
   
   float pContributingRays = contributingRays/(float)vssRays.size();
   
-  cout<<"Percentage of contributing rays:"<<pContributingRays<<endl;
+  cout<<"Ratio of contributing rays:"<<pContributingRays<<endl;
   
   if (mUseUnsuccCountImportance) {
@@ -233,15 +235,4 @@
 												 )
 {
-#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);
@@ -339,4 +330,6 @@
   
   AxisAlignedBox3 box = occluder->GetBox();
+  // consider slightly larger neighborhood of the occluder
+  // in search for unocclude reverse ray
   box.Scale(2.0f);
   
@@ -430,9 +423,10 @@
   }
 
+  return true;
 #else
   cerr << "warning: reverse mutation not supported!" << endl;
-#endif
-  
-  return true;
+  return false;
+#endif
+  
   
   // now the origin and termination is swapped compred to the generator ray
@@ -510,6 +504,6 @@
 
 #else
-	cerr << "warning: shiluette mutation not supported" << endl;
-	return Vector3(0, 0, 0);
+  //cerr << "warning: shiluette mutation not supported" << endl;
+  return Vector3(0, 0, 0);
 #endif
   
@@ -523,5 +517,5 @@
   if (mRays.size() == 0) {
 	float rr[5];
-	// use direction based distribution
+	// use direction based distribution until we have some mutation candidates
 	Vector3 origin, direction;
 
@@ -532,5 +526,5 @@
 												  Vector3(rr[0], rr[1], rr[2]));
 	
-
+	
 	direction = UniformRandomVector(rr[3], rr[4]);
 	
@@ -538,7 +532,7 @@
 	sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
 	sray.mGeneratorId = -1;
-
+	
 	return true;
-  } 
+  }
 
   int index;
@@ -552,5 +546,5 @@
 	if (
 		mRays[index].GetSamplingFactor() >= mRays[mLastIndex].GetSamplingFactor()) {
-	  // make another round
+	  // make another round to equalize the oversampling factor
 
 	  //	  cout<<"R2"<<endl;
@@ -666,5 +660,4 @@
   
   if (mUseSilhouetteSamples && a < mSilhouetteProb) {
-	
 	termination += ComputeSilhouetteTerminationMutation(*ray,
 														origin,
@@ -692,5 +685,10 @@
   
   // shift the origin a little bit
+#if 1
   origin += direction*0.5f;
+#else
+  // $$JB - 14.1. 2008 increase shift to test HavranRayCaster
+  origin = 0.5f*(origin + termination);
+#endif
   
   direction.Normalize();
@@ -718,5 +716,4 @@
     mRays[index].mMutations++;
 	mRays[index].mUnsuccessfulMutations++;
-
 	return true;
   }
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 2590)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 2592)
@@ -1230,8 +1230,8 @@
 	if (
 	    rays.size() > 10000)
-	{
+	  {
 
 	  mRayCaster->SortRays(rays);
-	  //	cout<<"Rays sorted in "<<TimeDiff(t1, GetTime())<<" ms."<<endl;
+	  cout<<"Rays sorted in "<<TimeDiff(t1, GetTime())<<" ms."<<endl;
 	}
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/default.env
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/default.env	(revision 2588)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/default.env	(revision 2592)
@@ -53,5 +53,5 @@
 	useGlDebugger false
 # 0 = INTERNAL 	1 = MLRT 2 = HAVRAN
-	rayCastMethod 1
+	rayCastMethod 2
 	
 #	type sampling
Index: GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp	(revision 2588)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp	(revision 2592)
@@ -6,5 +6,4 @@
 
 
-#include "HavranRayCaster.h"
 #include "VssRay.h"
 #include "KdTree.h"
@@ -13,19 +12,31 @@
 
 #ifdef USE_HAVRAN_RAYCASTER
+#include "ktbconf.h"
 #include "timer.h"
 #include "ktball.h"
 #include "ktb.h"
+#ifdef _USE_HAVRAN_SSE
+
 #include "raypack.h"
+#endif // _USE_HAVRAN_SSE
 #endif // USE_HAVRAN_RAYCASTER
+#include "HavranRayCaster.h"
 
 #define DEBUG_RAYCAST 0
 
+// This macro should be undefined when testing ray tracing
+// by casting rays from file or using camera
+#define _PROCESS_RAY
+
 namespace GtpVisibilityPreprocessor {
 
 
 #ifdef USE_HAVRAN_RAYCASTER 
+#ifdef _USE_HAVRAN_SSE
+  
 // static rays
 RayPacket2x2
 HavranRayCaster::raypack;
+#endif // _USE_HAVRAN_SSE
 #endif // USE_HAVRAN_RAYCASTER
   
@@ -65,15 +76,28 @@
 }
 
-
+bool
+HavranRayCaster::ExportBinTree(const string &filename)
+{
+  return mKtbtree->ExportBinTree(filename);
+}
+
+bool
+HavranRayCaster::ImportBinTree(const string &filename, ObjectContainer &objects)
+{
+  return mKtbtree->ImportBinTree(filename, objects);
+}
+
+  
 // Using packet of 4 rays supposing that these are coherent
 // We give a box to which each ray is clipped to before the
 // ray shooting is computed !
 void HavranRayCaster::CastRaysPacket4(Vector3 &boxmax,
-									  Vector3 &boxmin,
-									  Vector3 origin4[],
-									  Vector3 direction4[],
-									  int     result4[],
-									  float   dist4[])
-{
+				      Vector3 &boxmin,
+				      Vector3 origin4[],
+				      Vector3 direction4[],
+				      int     result4[],
+				      float   dist4[])
+{
+#ifdef _USE_HAVRAN_SSE  
 #ifdef USE_HAVRAN_RAYCASTER 
   for (int i = 0; i < 4; i++) {
@@ -92,12 +116,33 @@
     // this object
     Intersectable* intersectable;
-    if ( (intersectable = raypack.GetObject(i)) )
-	  // $$JB this is very hacky - > should be replaced with fetching the index of the triangle
-	  result4[i] = (int)intersectable;
+    if ( (intersectable = raypack.GetObject(i)) ) {
+      // $$JB this is very hacky - > should be replaced with fetching the index of the triangle
+      result4[i] = intersectable->mId;
+      dist4[i] = raypack.GetT(i);
+    }
   }
   
   return;
-#endif
-}
+#endif // USE_HAVRAN_RAYCASTER 
+#else // _USE_HAVRAN_SSE
+  // Compute the result ray by ray
+  SimpleRay sray;
+  for (int i = 0; i < 4; i++) {
+    sray.mOrigin = origin4[i];
+    sray.mDirection = direction4[i];
+    mKtbtree->FindNearestI(sray, boxmin, boxmax);
+    if (SimpleRay::IntersectionRes[0].intersectable) {
+      // This is object ID - is this the triangle index ???
+      result4[i] = SimpleRay::IntersectionRes[0].intersectable->mId;
+      dist4[i] = SimpleRay::IntersectionRes[0].tdist;
+    }
+    else {
+      result4[i] = -1; // no intersection
+    }
+  } // for i;
+  
+#endif // _USE_HAVRAN_SSE  
+}
+
 
 int HavranRayCaster::CastRay(const SimpleRay &simpleRay,
@@ -124,6 +169,6 @@
   
   // ray.mFlags &= ~Ray::CULL_BACKFACES;
-
-  if (mKtbtree->FindNearestI(simpleRay))  {
+  bool result;
+  if ((result = mKtbtree->FindNearestI(simpleRay)))  {
     hitA.mObject = SimpleRay::IntersectionRes[0].intersectable;
     float tdist = SimpleRay::IntersectionRes[0].tdist;
@@ -151,4 +196,5 @@
   }
 
+#ifdef _PROCESS_RAY  
   // This code is also in IntelRayCaster.cpp
   return ProcessRay(
@@ -161,7 +207,11 @@
 		    pruneInvalidRays
 		    );
-#else // USE_HAVRAN_RAYCASTER 
+#else // _PROCESS_RAY
+  return result;
+#endif // _PROCESS_RAY  
+
+#else
   return 0;
-#endif
+#endif // USE_HAVRAN_RAYCASTER   
 }
 
@@ -190,7 +240,7 @@
 
 #if 0
-  SimpleRayContainer::const_iterator sit = rays.begin() + offset;
-  SimpleRayContainer::const_iterator sit_end = sit + 16;
-  
+  // This is shooting ray by ray
+  SimpleRayContainer::const_iterator sit, sit_end = rays.end();
+
   // no acceleration for ray bundles implemented right now
   for (sit = rays.begin(); sit != sit_end; ++ sit) 
@@ -199,30 +249,18 @@
   }
 #else
+  // Use special algorithm for 16 rays at once
   if (castDoubleRays) {
-#if 0
+
+#if 0 // ------------------------------------------------
     // Special routine to cast double sided rays
     mKtbtree->SetOffset(0);
+#ifdef _USE_HAVRAN_SSE  
     mKtbtree->FindNearestI_16twoDir(rays);
 #else
+    mKtbtree->FindNearestI_16twoDirNoSSE(rays);    
+#endif
+#else // -------------------------------------------------
     // Here we decompose shooting into two phases
-#if 1    
-    // Here we shoot first forward rays and then backward ones
-    mKtbtree->SetOffset(0);
-    mKtbtree->FindNearestI_16oneDir(rays, offset);
-    SimpleRayContainer::iterator sit = rays.begin() + offset; 
-    SimpleRayContainer::const_iterator sit_end = rays.begin() + offset + 16;
-    for (; sit != sit_end; ++ sit) 
-    {
-      (*sit).mDirection = - (*sit).mDirection;
-    }
-    // We store the results at different place
-    mKtbtree->SetOffset(16);
-    mKtbtree->FindNearestI_16oneDir(rays, offset);
-    sit = rays.begin() + offset; 
-    for ( ; sit != sit_end; ++ sit) 
-    {
-      (*sit).mDirection = - (*sit).mDirection;
-    }
-#else
+
     // Here we shoot first backward rays and forward ones
     SimpleRayContainer::iterator sit = rays.begin() + offset; 
@@ -232,7 +270,12 @@
       (*sit).mDirection = - (*sit).mDirection;
     }
-    // backward rays to be shot
+    // backward rays to be shot - saving with offset 16
+#ifdef _USE_HAVRAN_SSE  
+    mKtbtree->SetOffset(0);
+    mKtbtree->FindNearestI_16oneDir(rays, offset, 16);
+#else
     mKtbtree->SetOffset(16);
-    mKtbtree->FindNearestI_16oneDir(rays, offset);
+    mKtbtree->FindNearestI_16oneDirNoSSE(rays, offset);
+#endif // _USE_HAVRAN_SSE
     sit = rays.begin() + offset; 
     for ( ; sit != sit_end; ++ sit) 
@@ -241,16 +284,25 @@
     }
     // forward rays to be shot
+#ifdef _USE_HAVRAN_SSE
     mKtbtree->SetOffset(0);
-    mKtbtree->FindNearestI_16oneDir(rays, offset);
-#endif    
-#endif    
-  }
+    mKtbtree->FindNearestI_16oneDir(rays, offset, 0);
+#else
+    mKtbtree->SetOffset(0);
+    mKtbtree->FindNearestI_16oneDirNoSSE(rays, offset);
+#endif // _USE_HAVRAN_SSE  
+#endif // ------------------------------------------------
+  } // cast double rays
   else {
+    // ONLY single rays
     // Shoot all 16 rays  at the same time using a special algorithm
     mKtbtree->SetOffset(0);
-    mKtbtree->FindNearestI_16oneDir(rays, offset);    
-  }
-#endif
-
+#ifdef _USE_HAVRAN_SSE  
+    mKtbtree->FindNearestI_16oneDir(rays, offset, 0);    
+#else
+    mKtbtree->FindNearestI_16oneDirNoSSE(rays, offset);    
+#endif // _USE_HAVRAN_SSE  
+  }
+#endif
+  
   
   for (int i=0, k=offset; i < 16; i++, k++) 
@@ -282,13 +334,12 @@
     {
       Intersectable *intersect =
-		SimpleRay::IntersectionRes[i+16].intersectable;
+	SimpleRay::IntersectionRes[i+16].intersectable;
 
       if (intersect)
       { 
-		hitB.mObject = intersect;
-		hitB.mNormal = intersect->GetNormal(0);
-		// $$ JB : i -> i+16
-		float tdist = SimpleRay::IntersectionRes[i+16].tdist;
-		hitB.mPoint = rays[k].Extrap(-tdist);
+	hitB.mObject = intersect;
+	hitB.mNormal = intersect->GetNormal(0);;
+	float tdist = SimpleRay::IntersectionRes[16+i].tdist;
+	hitB.mPoint = rays[k].Extrap(-tdist);
       }
     }
@@ -298,13 +349,13 @@
 #endif
 
-#if 1
+#ifdef _PROCESS_RAY  
     ProcessRay(rays[k],
-			   hitA,
-			   hitB,
-			   vssRays,
-			   sbox,
-			   castDoubleRays,
-			   pruneInvalidRays
-			   );
+	       hitA,
+	       hitB,
+	       vssRays,
+	       sbox,
+	       castDoubleRays,
+	       pruneInvalidRays
+	       );
 #endif
   } // for
@@ -340,5 +391,5 @@
   for (i=0; i < packets; i++) {
     int offset = i * 16;
-    mKtbtree->FindNearestI_16oneDir(rays, offset);
+    mKtbtree->FindNearestI_16oneDir(rays, offset, 0);
     // ??? What to do with the results ? These are
     // not used at the moment also in IntelRayCaster.cpp
@@ -372,5 +423,6 @@
   cout<<1e-3*2*rays.size()/TimeDiff(time, GetTime())<<" Mrays/sec"<<endl;
 #endif
-  
+
+  // Cast only by 16 rays at once
   for (int i=0; i < buckets; i++, offset+=16) {
     CastRays16(rays, offset, vssRays, sbox,
@@ -381,9 +433,12 @@
   }
 
+  // Cast the rest of the rays
   for (; offset < (int)rays.size(); offset++)
     CastRay(rays[offset], vssRays, sbox, castDoubleRay, pruneInvalidRays);
-}
-
-#if 0
+
+  return;
+}
+
+#ifdef _USE_HAVRAN_SSE
 // BUG1 41579  196948 1064111
 // BUG2 254    1672   10869
@@ -392,8 +447,9 @@
 void
 HavranRayCaster::CastRaysPacket2x2(RayPacket2x2 &raysPack,
-								   bool castDoubleRay,
-								   const bool pruneInvalidRays)
-{
-#ifdef USE_HAVRAN_RAYCASTER 
+				   bool castDoubleRay,
+				   const bool pruneInvalidRays)
+{
+#ifdef USE_HAVRAN_RAYCASTER 
+#ifdef _USE_HAVRAN_SSE
 
   if (castDoubleRay) {
@@ -443,7 +499,8 @@
 
   return;
+#endif // _USE_HAVRAN_SSE
 #endif // USE_HAVRAN_RAYCASTER 
 }
-#endif
+#endif // _USE_HAVRAN_SSE
 
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Makefile
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 2589)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 2592)
@@ -1,5 +1,5 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (2.00a) (Qt 4.1.2) on: ?t 10. I 13:25:38 2008
+# Generated by qmake (2.00a) (Qt 4.1.2) on: po 14. I 15:10:11 2008
 # Project:  preprocessor.pro
 # Template: app
@@ -63,5 +63,5 @@
 	$(MAKE) -f $(MAKEFILE).Debug uninstall
 
-Makefile: preprocessor.pro  C:/Qt/4.1.2/mkspecs/win32-msvc2005\qmake.conf C:/Qt/4.1.2/mkspecs/qconfig.pri \
+Makefile: preprocessor.pro  C:/Qt/4.1.2/mkspecs/win32-msvc.net\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/RayCaster.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp	(revision 2583)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp	(revision 2592)
@@ -178,4 +178,5 @@
   //const int batchsize = 16384;
   const int batchsize = 8192;
+  // const int batchsize = 1024;
   //const int batchsize = 128;
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h	(revision 2583)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h	(revision 2592)
@@ -21,4 +21,7 @@
 
 #include "Vector3.h"
+#include "ktbconf.h"
+  
+#ifdef _USE_HAVRAN_SSE
   
 #ifdef __SSE__
@@ -337,20 +340,20 @@
   // initialize inverted dir ?
 #ifdef _USE_INVDIR_RP
-  // inverted dir - maybe an overkill for SSE
-  // to be checked !
-  idx[0] = 1.0f / dx[0];
-  idx[1] = 1.0f / dx[1];
-  idx[2] = 1.0f / dx[2];
-  idx[3] = 1.0f / dx[3];
-  
-  idy[0] = 1.0f / dy[0];
-  idy[1] = 1.0f / dy[1];
-  idy[2] = 1.0f / dy[2];
-  idy[3] = 1.0f / dy[3];
-
-  idz[0] = 1.0f / dz[0];
-  idz[1] = 1.0f / dz[1];
-  idz[2] = 1.0f / dz[2];
-  idz[3] = 1.0f / dz[3];
+  //
+  const float epsAdd = -1e-25;
+  idx[0] = 1.0f / (dx[0] + epsAdd);
+  idx[1] = 1.0f / (dx[1] + epsAdd);
+  idx[2] = 1.0f / (dx[2] + epsAdd);
+  idx[3] = 1.0f / (dx[3] + epsAdd);
+  
+  idy[0] = 1.0f / (dy[0] + epsAdd);
+  idy[1] = 1.0f / (dy[1] + epsAdd);
+  idy[2] = 1.0f / (dy[2] + epsAdd);
+  idy[3] = 1.0f / (dy[3] + epsAdd);
+  
+  idz[0] = 1.0f / (dz[0] + epsAdd);
+  idz[1] = 1.0f / (dz[1] + epsAdd);
+  idz[2] = 1.0f / (dz[2] + epsAdd);
+  idz[3] = 1.0f / (dz[3] + epsAdd);
 #endif
 }
@@ -430,11 +433,8 @@
   tmax[i] = tmaxnew;
 }
-#else // __SSE__
-
-// ? What to do here
-//#error "AAA"
-
 #endif // __SSE__
 
+#endif // _USE_HAVRAN_SSE
+  
 } // namespace
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.cpp	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.cpp	(revision 2592)
@@ -22,4 +22,5 @@
 #include <fstream>
 #include <iomanip>
+#include <stack>
 
 namespace GtpVisibilityPreprocessor {
@@ -96,12 +97,12 @@
   const int size1D = 150;
 
-  cerr<<"th1"<<endl;
+  // cerr<<"th1"<<endl;
   // $$JB tmp
   size2D = 1;
-  cerr<<"size2d="<<size2D<<endl;
+  // cerr<<"size2d="<<size2D<<endl;
   // initial index
   index2D = 0;
   AllocBuffer(size2D, size1D);
-  cerr<<"th2"<<endl;
+  //  cerr<<"th2"<<endl;
 
 #if 0  
@@ -117,5 +118,5 @@
   }
 #else
-  cerr<<"th3"<<endl;
+  // cerr<<"th3"<<endl;
 
   traversalClass = new CKTBTraversal; 
@@ -308,5 +309,5 @@
 #endif // _DEBUG
 
-  cerr<<"hh2"<<endl;
+  // cerr<<"hh2"<<endl;
   // possibly remove the old building class
   // and allocate the class for CKTB building up
@@ -315,9 +316,9 @@
   buildClass = AllocBuildClass();
 
-  cerr<<"hh3"<<endl;
+  // cerr<<"hh3"<<endl;
 
   // and for traversal
   AllocTraversalClass();
-  cerr<<"hh4"<<endl;
+  // cerr<<"hh4"<<endl;
 
   if (boxToBoundObjects) {
@@ -330,13 +331,13 @@
     bbox.Include(*boxToInclude);  
 
-  cerr<<"hh5"<<endl;
+  // cerr<<"hh5"<<endl;
 
   // build up the CKTB tree and make setup for traversal
   SKTBNodeT *rootNode =
     buildClass->BuildUp(objlist, bbox, verbose);
-  cerr<<"hh6"<<endl;
+  // cerr<<"hh6"<<endl;
 
   traversalClass->Setup(bbox, rootNode);
-  cerr<<"hh7"<<endl;
+  // cerr<<"hh7"<<endl;
 
 #ifdef _DEBUG
@@ -400,4 +401,14 @@
 }
 
+int
+CKTB::FindNearestI(const SimpleRay &ray, Vector3 &boxmin, Vector3 &boxmax)
+{
+  const AxisAlignedBox3 localbox(boxmin, boxmax);
+  return traversalClass->FindNearestI(ray, localbox);      
+}
+
+  
+#ifdef _USE_HAVRAN_SSE  
+
 #ifdef __SSE__
 // --------------------------------------------------------------
@@ -420,4 +431,7 @@
 
 #endif // __SSE__
+
+#endif // _USE_HAVRAN_SSE  
+
   
 void
@@ -752,5 +766,5 @@
 SKTBNodeT*
 CKTB::ImportBinLeaf(IN_STREAM &stream, 
-		    SKTBNodeT *parent,
+		    SKTBNodeT **nodeToLink,
 		    const ObjectContainer &objects)
 {
@@ -764,4 +778,5 @@
     // since for 
     SKTBNodeT* leaf = buildClass->AllocLeaf(0);
+    *nodeToLink = buildClass->nodeToLink;
     return leaf;
   }
@@ -769,4 +784,5 @@
   MeshInstance dummyInst(NULL);
   SKTBNodeT* leaf = buildClass->AllocLeaf(size);
+  *nodeToLink = buildClass->nodeToLink;
 
   ObjectContainer *newobjlist = new ObjectContainer;
@@ -809,5 +825,5 @@
 
 SKTBNodeT *
-CKTB::ImportBinInterior(IN_STREAM  &stream, SKTBNodeT *parent)
+CKTB::ImportBinInterior(IN_STREAM  &stream, SKTBNodeT **nodeToLink)
 {
   int axis;
@@ -819,4 +835,5 @@
   SKTBNodeT *interiorNode =
     buildClass->AllocInteriorNode(axis, pos, 0, 0);
+  *nodeToLink = buildClass->nodeToLink;
 	
   return interiorNode;
@@ -830,231 +847,196 @@
   if (!stream.is_open()) return false;
 
-  // DEBUG << "CKTB::GatherStats called" << endl;
-  struct EDir {
-    SKTBNodeT *node;
-    int depth;
-    AxisAlignedBox3 bbox;
-  };
-  EDir stack[CKTBTraversal::MAX_TRAVERSAL_HEIGHT];
-  int stackTop = 0; // pointer to the stack
-  stack[stackTop].node = 0;
-  stack[stackTop].depth = 0;
-  stack[stackTop].bbox = bbox;
-  int tmpdir;
+  int startMark = 0xf0f0;
+  stream.write(reinterpret_cast<char *>(&startMark), sizeof(int));  
   
   // export binary version of mesh
-  queue<SKTBNodeT *> tStack;
-  // start from the root node
-  SKTBNodeT *currNode = buildClass->GetRootNode();
-
-  // copy the box over the scene objects
-  AxisAlignedBox3 bwbox = this->bbox;
-  int currDepth = 0;  
-
-  // traverse through whole CKTB tree
-  do {
-    CKTBAxes::Axes nodeType = buildClass->GetNodeType(currNode);
-    if (nodeType == CKTBAxes::EE_Link) {
-      // just relink
-      currNode = buildClass->GetLinkNode(currNode);
-      // cout << "Link node was accessed" << endl;
-      continue;	
+  std::stack<SKTBNodeT *> tStack;
+  tStack.push(buildClass->GetRootNode());
+  int cntSaves = 0, cntSavesLeaf = 0, cntSavesIN = 0;
+
+  while(!tStack.empty())
+  {
+    SKTBNodeT *node = tStack.top();
+    tStack.pop();    
+
+    int nodeType = buildClass->GetNodeType(node);
+    
+    if (nodeType == CKTBAxes::EE_Link)
+    {
+      tStack.push(buildClass->GetLinkNode(node));
     }
-      
-    if ( (nodeType == CKTBAxes::EE_X_axis)||
-	 (nodeType == CKTBAxes::EE_Y_axis)||
-	 (nodeType == CKTBAxes::EE_Z_axis)||
-	 (nodeType == CKTBAxes::EE_X_axisBox)||
-	 (nodeType == CKTBAxes::EE_Y_axisBox)||
-	 (nodeType == CKTBAxes::EE_Z_axisBox)
-       ) {
-      
-      // push onto the stack
-      AxisAlignedBox3 rightBox = bwbox;
-      float value = currNode->splitPlane;
-      switch(nodeType) {
-        case CKTBAxes:: EE_X_axisBox: {
-	  rightBox.SetMin(0, value);
-	  bwbox.SetMax(0, value);
-	  break;	  
-	}
-        case CKTBAxes:: EE_X_axis: {
-	  rightBox.SetMin(0, value);
-	  bwbox.SetMax(0, value);
-	  break;	  
-	}
-        case CKTBAxes:: EE_Y_axisBox: {
-	  rightBox.SetMin(0, value);
-	  bwbox.SetMax(0, value);
-	  break;	  
-	}
-        case CKTBAxes:: EE_Y_axis: {
-	  rightBox.SetMin(1, value);
-	  bwbox.SetMax(1, value);
-	  break;	  
-	}
-        case CKTBAxes:: EE_Z_axisBox: {
-	  rightBox.SetMin(0, value);
-	  bwbox.SetMax(0, value);
-	  break;	  
-	}
-        case CKTBAxes:: EE_Z_axis: {
-	  rightBox.SetMin(2, value);
-	  bwbox.SetMax(2, value);
-	  break;	  
-	}
-      } // switch
-
-      // Exporting interior node
-      ExportBinInterior(stream, currNode);
-      
-      stack[stackTop].node = buildClass->GetRight(currNode);
-      stack[stackTop].depth = currDepth + 1;
-      stack[stackTop].bbox = rightBox;
-      stackTop++;
-      // where to go
-      if ( (nodeType == CKTBAxes:: EE_X_axisBox)||
-	   (nodeType == CKTBAxes:: EE_Y_axisBox)||
-	   (nodeType == CKTBAxes:: EE_Z_axisBox) )
-	currNode = buildClass->GetLeftMinBox(currNode);
+    else
+    if (nodeType == CKTBAxes::EE_Leaf)
+    {
+      //Debug << "l";
+      ExportBinLeaf(stream, node);
+      cntSavesLeaf++; cntSaves++;      
+    }
+    else
+    { // interior node
+      //Debug << "i";
+      ExportBinInterior(stream, node);
+      cntSavesIN++; cntSaves++;      
+
+      tStack.push(buildClass->GetRight(node));
+      if (nodeType >= CKTBAxes::EE_X_axisBox)
+	tStack.push(buildClass->GetLeftLong(node));
       else
-	currNode = buildClass->GetLeft(currNode);
-      // the box was already set
-      currDepth++;
-      
-      continue;
-    } // interior nodes
-
-    assert(nodeType == CKTBAxes::EE_Leaf);
-      
-    while ( (nodeType == CKTBAxes::EE_Leaf) && currDepth) {
-      ObjectContainer *olist = buildClass->GetObjList(currNode);
-      // Exporting leaf node
-      ExportBinLeaf(stream, static_cast<SKTBNodeT *>(currNode));
-	
-      // traverse up the tree
-      stackTop--;
-      if (stackTop < 0)
-	goto FINISH;
-      // pop the node from the stack
-      currNode = stack[stackTop].node;
-      currDepth = stack[stackTop].depth;
-      bwbox = stack[stackTop].bbox;
-      
-      nodeType = buildClass->GetNodeType(currNode);
-	
-    } // while is a leaf
-  }
-  while (stackTop >= 0);
-
-FINISH:
-  // close the stream
+	tStack.push(buildClass->GetLeft(node));
+    }
+  } // while
+
+  int endMark = 0xf0ff;
+  stream.write(reinterpret_cast<char *>(&endMark), sizeof(int));  
+
+  cout << "Writes " << cntSaves << " cntLeafs "
+       << cntSavesLeaf << " cntIN " << cntSavesIN << endl;
+  
   stream.close();
-  
-  return true;
-}
-
-#if 0
-KdNode *KdTree::ImportNextNode(IN_STREAM &stream, 
-			       KdInterior *parent,
-			       const ObjectContainer &objects)
-{
-	int nodeType;
-	stream.read(reinterpret_cast<char *>(&nodeType), sizeof(int));
-
-	if (nodeType == TYPE_LEAF)
-		return ImportBinLeaf(stream, static_cast<KdInterior *>(parent), objects);
-
-	if (nodeType == TYPE_INTERIOR)
-		return ImportBinInterior(stream, static_cast<KdInterior *>(parent));
-
-	Debug << "error! loading failed!" << endl;
-	return NULL;
-}
-
-
-bool KdTree::ImportBinTree(const string &filename, ObjectContainer &objects)
-{
-  // export binary version of mesh
-  queue<TraversalData> tStack;
+  return true; // saved
+}
+
+
+SKTBNodeT *
+CKTB::ImportNextNode(IN_STREAM &stream, 
+		     SKTBNodeT **nodeToLink,
+		     const ObjectContainer &objects)
+{
+  int nodeType;
+  stream.read(reinterpret_cast<char *>(&nodeType), sizeof(int));
+
+  if (nodeType == TYPE_LEAF)
+    return ImportBinLeaf(stream, nodeToLink, objects);
+
+  if (nodeType == TYPE_INTERIOR)
+    return ImportBinInterior(stream, nodeToLink);
+
+  Debug << "error! loading failed!" << endl;
+  return NULL;
+}
+
+
+// creates the ASDS according to the file description
+bool
+CKTB::ImportBinTree(const string &filename, ObjectContainer &objects)
+{
+  // open the stream in text mode
   IN_STREAM stream(filename.c_str(), IN_BIN_MODE);
-  
-  if (!stream.is_open()) return false;
-  
-  // sort objects by their id
+
+  if (!stream.is_open()) {
+    cerr << "Kd-tree description file (.kbt) cannot be opened for reading\n";
+    return true; // error
+  }
+
+  int mark;
+  stream.read(reinterpret_cast<char *>(&mark), sizeof(int));
+  if (mark != 0xf0f0) {
+    cout << "Something wrong with the tree - heading" << endl;
+    abort();
+  }
+  
+  // sort objects by their id to allow list creation in kd-tree leaves
   //	if (!is_sorted(objects.begin(), objects.end(), ilt))
   sort(objects.begin(), objects.end(), ilt);
   
-  mBox.Initialize();
-  ObjectContainer::const_iterator oit, oit_end = objects.end();
-  
-  ///////////////////////////
-  //-- compute bounding box of object space
-  
-  for (oit = objects.begin(); oit != oit_end; ++ oit)
-  {
-    const AxisAlignedBox3 box = (*oit)->GetBox();
-    mBox.Include(box);
-  }
-  
-  // hack: we make a new root
-  DEL_PTR(mRoot);
-  
-  mRoot = ImportNextNode(stream, NULL, objects);
-  
-  tStack.push(TraversalData(mRoot, mBox, 0));
-  mStat.Reset();
-  mStat.nodes = 1;
-  
+  // remove all the data in the current data structure
+  CKTBAllocManPredecessor *bc = GetBuildClass();
+  if (!bc)
+    bc = buildClass = AllocBuildClass();
+  else
+    bc->Remove();
+
+  // How many items can be allocated at once
+  int maxItemsAtOnce = 1;
+  if (makeMinBoxes) {
+#ifdef _KTB8Bytes    
+    // We need to allocate for boxes the memory in a row
+    maxItemsAtOnce = 5; // 8x5=40 = 16+24;
+#else
+    maxItemsAtOnce = 4; // 12x4=48 = 24+24;
+#endif // _KTB8Bytes
+  }
+  
+  bc->InitAux(0, CKTBNodeAbstract::MAX_HEIGHT - 1, maxItemsAtOnce);  
+
+  // Compute the box from all objects
+  bbox.Initialize();
+  for (ObjectContainer::iterator sc = objects.begin();
+       sc != objects.end(); sc++) {
+    AxisAlignedBox3 abox = (*sc)->GetBox();
+    bbox.Include(abox);
+  } // for
+   
+  // initialize the root node
+  SKTBNodeT *node;
+
+  bc->root = ImportNextNode(stream, &node, objects);
+  assert(bc->root == node);
+
+  std::stack<STraversalData> tStack;
+  tStack.push(STraversalData(bc->root, 1, 0));
+  tStack.push(STraversalData(bc->root, 0, 0));
+  int depth, lastDepth;
+  
+  //int cntReads = 1, cntReadsLeaf = 0, cntReadsIN = 1;
   while(!tStack.empty())
-  {
-    TraversalData tData = tStack.front();
+  {    
+    STraversalData tData = tStack.top();
     tStack.pop();
     
-    KdNode *node = tData.mNode;
+    node = tData.mNode;
+    lastDepth = depth = tData.depth;
+    //int nodeType = bc->GetNodeType(node);
+    //cout << "nodeType = " << nodeType << " adr = " << (void*)node << endl;
+
+    SKTBNodeT *childNodeLink = 0;
+    SKTBNodeT *childNode = ImportNextNode(stream, &childNodeLink, objects);
+    //    cout << "childNode = " << (void*)childNode
+    // << " linkNode = " << (void*)childNodeLink << endl;
+
+    if (tData.dir)
+      bc->SetInteriorNodeRight(node, childNodeLink);
+    else
+      bc->SetInteriorNodeLeft(node, childNodeLink);
     
-    if (!node->IsLeaf())
+    //cntReads++;
+    if (!bc->IsLeaf_(childNode))
     {
-      mStat.nodes += 2;
-      
-      //Debug << "i" ;
-      KdInterior *interior = static_cast<KdInterior *>(node);
-      interior->mBox = tData.mBox;
-      
-      KdNode *front = ImportNextNode(stream, interior, objects);
-      KdNode *back = ImportNextNode(stream, interior, objects);
-      
-      interior->SetupChildLinks(back, front);
-      
-      ++ mStat.splits[interior->mAxis];
-      
-      // compute new bounding box
-      AxisAlignedBox3 frontBox, backBox;
-      
-      tData.mBox.Split(interior->mAxis, 
-		       interior->mPosition, 
-		       frontBox, 
-		       backBox);
-      
-      tStack.push(TraversalData(front, frontBox, tData.mDepth + 1));			
-      tStack.push(TraversalData(back, backBox, tData.mDepth + 1));
+      //cntReadsIN++;
+      tStack.push(STraversalData(childNode, 1, depth+1));
+      tStack.push(STraversalData(childNode, 0, depth+1));
     }
-    else
-    {
-      EvaluateLeafStats(tData);
-      //cout << "l";
-    }
-  }
-
-  float area = GetBox().SurfaceArea()*mKdPvsArea;
-  
-  SetPvsTerminationNodes(area);
-  
-  Debug << mStat << endl;
-  
-  return true;
-}
-#endif
-
+    //else {
+    //  cntReadsLeaf++;
+    //  //cout << "leaf" << endl;
+    // }
+  } // while
+
+  //cout << "Last depth = " << lastDepth << endl;
+  //cout << "Reads " << cntReads << " cntLeafs "
+  //   << cntReadsLeaf << " cntIN " << cntReadsIN << endl;
+  
+  stream.read(reinterpret_cast<char *>(&mark), sizeof(int));
+  if (mark != 0xf0ff) {
+    cout << "Something wrong with the kd-tree in the file"
+	 << endl;
+    abort();
+  }
+  
+  if (!traversalClass)
+    AllocTraversalClass();
+    
+  // Make setup for traversal
+  traversalClass->Setup(bbox, bc->GetRootNode());
+  // Handles unbounded objects in traversal class
+  // traversalClass->Setup2(0);
+  
+  stream.close();
+
+  // now the tree is surely well constructed
+  builtUp = true;  
+
+  return false; // OK
+}
 
 void
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.cpp	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.cpp	(revision 2592)
@@ -63,5 +63,5 @@
   float tclosest = Limits::Infinity;
   int intersected = 0;
-  // indexRay += rayOffset;
+  indexRay += rayOffset;
   // iterate the whole list and find out the nearest intersection
   for (ObjectContainer::iterator sc = list->begin(); sc != sc_end; sc++) {
@@ -76,5 +76,4 @@
   return intersected;
 }
-
   
 // test all objects in the leaf for intersection with ray
@@ -244,5 +243,6 @@
 CKTB8BAllocMan::AllocInteriorNodeWithBox(int axis, float position,
 					 int cntLeft, int cntRight,
-					 const SBBox &tsbox, SKTBNodeT* prevMinBoxNode,
+					 const SBBox &tsbox,
+					 SKTBNodeT* prevMinBoxNode,
 					 int depthStore)
 {
@@ -354,4 +354,12 @@
   leftChild = leftChild; // to satisfy the compiler
 
+#if 1
+  if (! ((node+1 == leftChild) || (node+2 == leftChild) ||
+	 (node+5 == leftChild))) {
+
+    cout << "Problem with left node linking " << endl;
+  }
+
+#endif      
   // Check on correctness of DFS order
   assert( (node+1 == leftChild) || (node+2 == leftChild) || (node+5 == leftChild) );
@@ -390,4 +398,49 @@
 }
 
+// Set the pointers to children for the interior node
+void
+CKTB8BAllocMan::SetInteriorNodeLeft(SKTBNodeT *node,
+				    SKTBNodeT *leftChild)
+{
+  leftChild = leftChild; // to satisfy the compiler
+
+#if 1
+  if (! ((node+1 == leftChild) || (node+2 == leftChild) ||
+	 (node+5 == leftChild))) {
+
+    cout << "Problem with left node linking " << endl;
+    abort();
+  }
+
+#endif      
+  // Check on correctness of DFS order
+  assert( (node+1 == leftChild) || (node+2 == leftChild) || (node+5 == leftChild) );
+
+  // Nothing to do - left link is linked automatically
+  return;
+}
+
+// Set the pointers to children for the interior node
+void
+CKTB8BAllocMan::SetInteriorNodeRight(SKTBNodeT *node,
+				     SKTBNodeT *rightChild)
+{
+  assert(node->offset <= CKTBAxes::EE_Z_axisBox);
+  //assert(node->offset >= 0);
+  //(unsigned long)node->offset |= (unsigned long)rightChild;
+  node->offset = node->offset | (unsigned int)rightChild;
+#ifdef _DEBUG
+  // Compute right child
+  SKTBNodeT *p = GetRight(node);
+  if (p != rightChild) {
+    cerr << "Problem with implementation of setting right child" << endl;
+    abort();;
+  } 
+#endif
+
+  return;
+}
+
+  
 // Create the representation of the leaf. Note that possibly there
 // can be special cases, such as 0, 1, 2, or 3 objects, or in general
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.h	(revision 2592)
@@ -174,10 +174,11 @@
   
   int FindNearestI(const SimpleRay &ray);
+  int FindNearestI(const SimpleRay &ray, Vector3 &boxmin, Vector3 &boxmax);
   void SetOffset(int offset) { traversalClass->SetOffset(offset); }
-  int FindNearestI_16oneDir(SimpleRayContainer &rays) {
-    return traversalClass->FindNearestI_16oneDir(rays);
-  }
-  int FindNearestI_16oneDir(SimpleRayContainer &rays, int offset) {
-    return traversalClass->FindNearestI_16oneDir(rays, offset);
+  int FindNearestI_16oneDir(SimpleRayContainer &rays, int offset, int copyOffset) {
+    return traversalClass->FindNearestI_16oneDir(rays, offset, copyOffset);
+  }
+  int FindNearestI_16oneDirNoSSE(SimpleRayContainer &rays, int offset) {
+    return traversalClass->FindNearestI_16oneDirNoSSE(rays, offset);
   }
   int FindNearestI_16twoDir(SimpleRayContainer &rays) {
@@ -186,4 +187,5 @@
   
 #ifdef __SSE__
+#ifdef _USE_HAVRAN_SSE  
   // The same operations for packets of rays, if implemented by
   // a particular ASDS, otherwise it is emulated by decomposition
@@ -191,4 +193,5 @@
   void FindNearestI(RayPacket2x2 &raypack);
   void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax);
+#endif
 #else
   void FindNearestI(RayPacket2x2 &raypack) { }
@@ -199,10 +202,18 @@
   void ExportBinLeaf(OUT_STREAM &stream, SKTBNodeT *leaf);
   SKTBNodeT* ImportBinLeaf(IN_STREAM &stream, 
-			   SKTBNodeT *parent,
+			   SKTBNodeT **nodeToLink,
 			   const ObjectContainer &objects);
+
   void ExportBinInterior(OUT_STREAM &stream, SKTBNodeT *interior);
-  SKTBNodeT* ImportBinInterior(IN_STREAM  &stream,
-			       SKTBNodeT *parent);  
+  SKTBNodeT *ImportBinInterior(IN_STREAM  &stream,
+			       SKTBNodeT **nodeToLink);  
+
+  SKTBNodeT * ImportNextNode(IN_STREAM &stream, 
+			     SKTBNodeT **nodeToLink,
+			     const ObjectContainer &objects);
+
   bool ExportBinTree(const string &filename);
+  bool ImportBinTree(const string &filename, ObjectContainer &objects);
+
   
 protected:  
@@ -211,4 +222,18 @@
   int sizeBuffer, size1D;
   void AllocBuffer(int size2D, int size1Dv);
+
+  struct STraversalData
+  {  
+    SKTBNodeT *mNode;
+    int        dir;
+    int        depth;
+    STraversalData() {}
+
+    STraversalData(SKTBNodeT *n): mNode(n)
+    { }
+    
+    STraversalData(SKTBNodeT *n, int ndir, int ndepth):
+      mNode(n), dir(ndir), depth(ndepth) { }
+  };
 };
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.h	(revision 2592)
@@ -161,7 +161,4 @@
 #endif
   
-  // init the stack of auxiliary variables from min to max
-  void InitAux(int min, int max, int maxItemsAtOnce);
-
   // the direction of traversal through CKTB tree
   enum EDirection { LEFT = 0, RIGHT = 1, NOWHERE = 2 };
@@ -185,4 +182,7 @@
   ~CKTB8BAllocMan() { }
 
+  // init the stack of auxiliary variables from min to max
+  void InitAux(int min, int max, int maxItemsAtOnce);
+  
   // forget the content that is created by previous kd-tree construction
   // or just init for the first use.
@@ -232,4 +232,9 @@
   void SetInteriorNodeLinks(SKTBNodeT *node,
 			    SKTBNodeT *leftChild,
+			    SKTBNodeT *rightChild);
+
+  void SetInteriorNodeLeft(SKTBNodeT *node,
+			   SKTBNodeT *leftChild);
+  void SetInteriorNodeRight(SKTBNodeT *node,
 			    SKTBNodeT *rightChild);
 
@@ -275,4 +280,5 @@
   void DecDepth() {_currDepth--;}
 
+  
   // --------------------------------------------------------------
   // Statistics
@@ -300,5 +306,5 @@
   float _sumSurfaceAreaInteriorNodes;
 #endif
-
+  
   // General statistics on the depth of current node
   // the depth of the currently accessed node during building up
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.cpp	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.cpp	(revision 2592)
@@ -61,4 +61,5 @@
   float tclosest = Limits::Infinity;
   int intersected = 0;
+  // rayIndex += rayOffset;
   // iterate the whole list and find out the nearest intersection
   for (ObjectContainer::iterator sc = list->begin(); sc != sc_end; sc++) {
@@ -352,4 +353,24 @@
 }
 
+// Set the pointers to children for the interior node
+void
+CKTBAllocMan::SetInteriorNodeLeft(SKTBNodeT *node,
+				  SKTBNodeT *leftChild)
+{
+  leftChild = leftChild; // to satisfy the compiler
+
+  // Check on correctness of DFS order
+  assert( (node+1 == leftChild) || (node+2 == leftChild) || (node+4 == leftChild) );
+}
+
+// Set the pointers to children for the interior node
+void
+CKTBAllocMan::SetInteriorNodeRight(SKTBNodeT *node,
+				   SKTBNodeT *rightChild)
+{
+  node->right = rightChild;
+}
+
+  
 // Create the representation of the leaf. Note that possibly there
 // can be special cases, such as 0, 1, 2, or 3 objects, or in general
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/configh.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/configh.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/configh.h	(revision 2592)
@@ -59,4 +59,5 @@
 
 #ifdef _MSC_VER
+#define __SSE__
 // disable some Microsoft Visual Compiler warnings if necessary
 #pragma warning( disable:4305 )
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbtrav.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbtrav.h	(revision 2583)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbtrav.h	(revision 2592)
@@ -30,6 +30,4 @@
 
 namespace GtpVisibilityPreprocessor {
-
-class RayPacket2x2;
   
 //#define TRV000
@@ -136,8 +134,10 @@
 public:
   virtual int FindNearestI(const SimpleRay &ray) { return 0; }
+  virtual int FindNearestI(const SimpleRay &ray, const AxisAlignedBox3 &localbox) { return 0;}
   void SetOffset(int offset) { rayOffset = offset; }
   virtual int FindNearestI_16oneDir(SimpleRayContainer &rays) { return 0; }
-  virtual int FindNearestI_16oneDir(SimpleRayContainer &rays, int offset)
-  { return 0; } 
+  virtual int FindNearestI_16oneDir(SimpleRayContainer &rays,
+				    int offset, int copyOffset) { return 0; } 
+  virtual int FindNearestI_16oneDirNoSSE(SimpleRayContainer &rays, int offset) { return 0; }
   virtual int FindNearestI_16twoDir(SimpleRayContainer &rays) { return 0; }
 
@@ -145,4 +145,5 @@
   enum { MAX_NESTING = 2};
   
+#ifdef _USE_HAVRAN_SSE  
   // The same operations for packets of rays, if implemented by
   // a particular ASDS, otherwise it is emulated by decomposition
@@ -150,5 +151,6 @@
   virtual void FindNearestI(RayPacket2x2 &raypack) { }
   virtual void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax) { }
-
+#endif
+  
   virtual void PrintStatsTR(ostream &) { }
   virtual void DynamicStatsReset() { }
@@ -306,10 +308,13 @@
   virtual int FindNearestI(const SimpleRay &ray);
   virtual int FindNearestI(const SimpleRay &ray, const AxisAlignedBox3 &localbox);
-
-  virtual int FindNearestI_16oneDir(SimpleRayContainer &rays) {
-    return FindNearestI_16oneDir(rays, 0);
-  }
-  // $$JB correction
-  virtual int FindNearestI_16oneDir(SimpleRayContainer &rays, int offset);
+  virtual int FindNearestI_16oneDir(SimpleRayContainer &rays,
+				    int offset, int copyOffset)
+#ifdef _USE_HAVRAN_SSE  
+    ; // this is implemented with SSE in ktbftrav.cpp !
+#else
+  { return 0;}
+#endif
+    
+  virtual int FindNearestI_16oneDirNoSSE(SimpleRayContainer &rays, int offset);
 
   int PrecomputeData(SimpleRayContainer &rays);
@@ -317,4 +322,5 @@
   virtual int FindNearestI_16twoDir(SimpleRayContainer &rays);
 
+#ifdef _USE_HAVRAN_SSE  
 #ifdef __SSE__
   // The same operations for packets of rays, if implemented by
@@ -324,5 +330,6 @@
   virtual void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax);
 #endif // __SSE__
-
+#endif // _USE_HAVRAN_SSE  
+  
   virtual void PrintStatsTR(ostream &) { }
 
@@ -408,4 +415,5 @@
   virtual int FindNearestI(SimpleRay &ray);
 
+#ifdef _USE_HAVRAN_SSE  
 #ifdef __SSE__
   // The same operations for packets of rays, if implemented by
@@ -415,4 +423,5 @@
   virtual void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax);
 #endif // __SSE__
+#endif // _USE_HAVRAN_SSE  
 
   virtual void DynamicStatsReset();
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbftrav.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbftrav.cpp	(revision 2583)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbftrav.cpp	(revision 2592)
@@ -20,130 +20,6 @@
 #ifdef TRV00F
 
-#if 0
-// the depth of traversal when kd-trees are nested. The global (the highest
-// level) kd-tree is at the depth 0.
-int
-CKTBTraversal::traversalDepth = 0;
-
-// sets the stack pointers that are used for the traversal.
-void
-CKTBTraversal::GetStackPointers(struct SStackElem *&entryPointer,
-				struct SStackElem *&exitPointer)
-{
-  assert(traversalDepth < MAX_NESTING);
-  
-  int index = traversalDepth * MAX_HEIGHT;
-  entryPointer = &(stack[index]);
-  exitPointer = &(stack[index + 1]);
-
-  traversalDepth++;  
-  return;
-}
-
-// sets the stack pointers that are used for the traversal.
-void
-CKTBTraversal::RestoreStackPointers()
-{
-  assert(traversalDepth >= 0);
-  traversalDepth--;  
-  return;
-}
-
-
-// default constructor
-CKTBTraversal::CKTBTraversal()
-{
-  bbox = AxisAlignedBox3(Vector3(MAXFLOAT),
-			 Vector3(-MAXFLOAT));
-  root = 0;
-  epsilon = 0;
-
-#ifdef __TRAVERSAL_STATISTICS
-  _allNodesTraversed = 0L;
-  _fullLeavesTraversed = 0L;
-  _emptyLeavesTraversed = 0L;
-#endif
-}
-  
-// Here we find the node (preferably minbox node) containing
-// the point
-const SKTBNodeT*
-CKTBTraversal::Locate(const Vector3 &position)
-{
-  const SKTBNodeT *current, *next = root;
-  const SKTBNodeT *returnNode = 0;
-
-  do {
-    current = next;
-    switch (GetNodeType(current)) {
-      // -------------------------------------------------
-      case CKTBAxes::EE_X_axis: {
-	if (position[CKTBAxes::EE_X_axis] < current->splitPlane)
-	  next =  GetLeft(current);
-	else
-	  next = GetRight(current);
-	break;
-      }
-      case CKTBAxes::EE_Y_axis: {
-	if (position[CKTBAxes::EE_Y_axis] < current->splitPlane)
-	  next =  GetLeft(current);
-	else
-	  next = GetRight(current);
-	break;
-      }
-      case CKTBAxes::EE_Z_axis: {
-	if (position[CKTBAxes::EE_Z_axis] < current->splitPlane)
-	  next =  GetLeft(current);
-	else
-	  next = GetRight(current);
-	break;
-      }
-      case CKTBAxes::EE_X_axisBox: {
-	if (position[CKTBAxes::EE_X_axis] < current->splitPlane)
-	  next =  GetLeft(current);
-	else
-	  next = GetRight(current);
-	returnNode = current;
-	break;
-      }
-      case CKTBAxes::EE_Y_axisBox: {
-	if (position[CKTBAxes::EE_Y_axis] < current->splitPlane)
-	  next =  GetLeft(current);
-	else
-	  next = GetRight(current);
-	returnNode = current;
-	break;
-      }
-      case CKTBAxes::EE_Z_axisBox: {
-	if (position[CKTBAxes::EE_Z_axis] < current->splitPlane)
-	  next =  GetLeft(current);
-	else
-	  next = GetRight(current);
-	returnNode = current;
-	break;
-      }
-      case CKTBAxes::EE_Leaf: {
-	next = 0; // finishing
-	break;
-      }
-      case CKTBAxes::EE_Link:{
-	next = GetLinkNode(current);
-	// cout << "Link node was accessed" << endl;
-	break;
-      }
-    } // switch
-  }
-  while (next != NULL);
-
-  if (returnNode)
-    return returnNode;
-
-  // return the last (leaf node) visited on the path
-  return current;
-}
-#endif
-  
-  
-#if 1
+// --------------------------------------------------------------
+// Shooting a single ray without SSE
 int
 CKTBTraversal::FindNearestI(const SimpleRay &ray)
@@ -162,11 +38,11 @@
   // passing through parameters
   float tmin, tmax;
+  SimpleRay::IntersectionRes[0].intersectable = 0;
   
   // test if the whole CKTB tree is missed by the input ray
   if ( (!root) ||
        (!bbox.ComputeMinMaxT(ray.mOrigin, ray.mDirection, &tmin, &tmax)) ||
-       (tmax <= tmin) ||
+       (tmax < tmin) ||
        (tmax <= 0.f) ) {
-    SimpleRay::IntersectionRes[0].intersectable = 0;
     return 0; // no object can be intersected
   }
@@ -185,7 +61,7 @@
 
   Vector3 invertedDir;
-  invertedDir.x = 1.0f / ray.mDirection.x;
-  invertedDir.y = 1.0f / ray.mDirection.y;
-  invertedDir.z = 1.0f / ray.mDirection.z;
+  invertedDir.x = 1.0f / (ray.mDirection.x - 1e-25f);
+  invertedDir.y = 1.0f / (ray.mDirection.y - 1e-25f);
+  invertedDir.z = 1.0f / (ray.mDirection.z - 1e-25f);
   
   // start from the root node
@@ -216,12 +92,16 @@
     }
 #endif
+
+CONTINUE_LINK:    
+
+    assert(tmin <= tmax);
+#ifdef __TRAVERSAL_STATISTICS
+    allNodesTraversed++;
+#endif // __TRAVERSAL_STATISTICS
+
+    register const int nodeType = GetNodeType(currNode);
+
+    // cout << " tmin = " << tmin << " tmax = " << tmax << " nodeType = " << (int)nodeType << endl;
     
-    assert(tmin <= tmax);
-    
-#ifdef __TRAVERSAL_STATISTICS
-    allNodesTraversed++;
-#endif // __TRAVERSAL_STATISTICS
-    register const int nodeType = GetNodeType(currNode);
-	
     if (nodeType < CKTBAxes::EE_Leaf) {
       float tval = (GetSplitValue(currNode) - ray.mOrigin[nodeType]);
@@ -265,5 +145,6 @@
 	  // test the objects in the full leaf against the ray
 	  
-	  SimpleRay::IntersectionRes[0].maxt = stack3[index].tmax;
+	  SimpleRay::IntersectionRes[0].maxt =
+	    stack3[index].tmax + Limits::Small;
 #if 0
 	  // using subroutine
@@ -320,7 +201,7 @@
 	cout << "Link " << endl;
 #endif
-	stack3[index].nodep = GetLinkNode(currNode);
 	// cout << "Link node was accessed" << endl;
-	continue;
+	currNode = GetLinkNode(currNode);
+	goto CONTINUE_LINK;
       }
     }
@@ -336,7 +217,8 @@
   return 0;
 } // FindNearestI - single ray
-#endif
-
-#if 1
+
+
+// Shooting a single ray without SSE with prespecified box of the scene, assuming
+// to be contained in the scene box!!!
 int
 CKTBTraversal::FindNearestI(const SimpleRay &ray, const AxisAlignedBox3 &localbox)
@@ -355,4 +237,5 @@
   // passing through parameters
   float tmin, tmax;
+  SimpleRay::IntersectionRes[0].intersectable = 0;
   
   // test if the whole CKTB tree is missed by the input ray
@@ -361,5 +244,4 @@
        (tmax <= tmin) ||
        (tmax <= 0.f) ) {
-    SimpleRay::IntersectionRes[0].intersectable = 0;
     return 0; // no object can be intersected
   }
@@ -378,7 +260,7 @@
 
   Vector3 invertedDir;
-  invertedDir.x = 1.0f / ray.mDirection.x;
-  invertedDir.y = 1.0f / ray.mDirection.y;
-  invertedDir.z = 1.0f / ray.mDirection.z;
+  invertedDir.x = 1.0f / (ray.mDirection.x - 1e-25f);
+  invertedDir.y = 1.0f / (ray.mDirection.y - 1e-25f);
+  invertedDir.z = 1.0f / (ray.mDirection.z - 1e-25f);
   
   // start from the root node
@@ -409,10 +291,12 @@
     }
 #endif
-    
+        
+CONTINUE_LINK:
     assert(tmin <= tmax);
-    
+
 #ifdef __TRAVERSAL_STATISTICS
     allNodesTraversed++;
 #endif // __TRAVERSAL_STATISTICS
+
     register const int nodeType = GetNodeType(currNode);
     if (nodeType < CKTBAxes::EE_Leaf) {
@@ -457,5 +341,6 @@
 	  // test the objects in the full leaf against the ray
 	  
-	  SimpleRay::IntersectionRes[0].maxt = stack3[index].tmax;
+	  SimpleRay::IntersectionRes[0].maxt =
+	    stack3[index].tmax + Limits::Small;
 #if 0
 	  // using subroutine
@@ -512,7 +397,7 @@
 	cout << "Link " << endl;
 #endif
-	stack3[index].nodep = GetLinkNode(currNode);
 	// cout << "Link node was accessed" << endl;
-	continue;
+	currNode = GetLinkNode(currNode);
+	goto CONTINUE_LINK;
       }
     }
@@ -528,196 +413,11 @@
   return 0;
 } // FindNearestI - single ray
-#endif
-  
-
-// ---------------------------------------------------------------------------
-// This is an attempt using SSE instructions - not successfull 31/12/2007
-#if 0
-int
-CKTBTraversal::FindNearestI(const SimpleRay &ray)
-{
-#if 0
-  static int counter = 0;
-  counter++;
-  bool debug = false;
-  if (counter == 530) {
-    debug = true;
-    cout << "COUNTER = " << counter << endl;
-    cout << "DEBUG starts" << endl;    
-  }
-#endif
-  
-  // passing through parameters
-  float tmin, tmax;
-  
-  // test if the whole CKTB tree is missed by the input ray
-  if ( (!root) ||
-       (!bbox.ComputeMinMaxT(ray.mOrigin, ray.mDirection, &tmin, &tmax)) ||
-       (tmax <= tmin) ||
-       (tmax <= 0.f) ) {
-    SimpleRay::IntersectionRes[0].intersectable = 0;
-    return 0; // no object can be intersected
-  }
-
-//#define _DEBUGKTB
-#ifdef _DEBUGKTB
-  int ib = 0;
-  int depth = 0; 
-#endif
-  
-#ifdef __TRAVERSAL_STATISTICS
-  int allNodesTraversed = 0L;
-  int fullLeavesTraversed = 0L;
-  int emptyLeavesTraversed = 0L;
-#endif // __TRAVERSAL_STATISTICS
-
-  Vector3 invertedDir;
-  invertedDir.x = 1.0f / ray.mDirection.x;
-  invertedDir.y = 1.0f / ray.mDirection.y;
-  invertedDir.z = 1.0f / ray.mDirection.z;
-  
-  // start from the root node
-  if (tmin < 0.f)
-    tmin = 0.f;
-
-  int index = 1;
-  stack3[1].nodep = root;
-  stack3[1].tmax = tmax;
-  tmax = tmin;
-  SKTBNodeT * childNodes[2];
-  int RayDirs[3];
-  RayDirs[0] = ray.mDirection.x < 0.f ? 1 : 0;
-  RayDirs[1] = ray.mDirection.y < 0.f ? 1 : 0;
-  RayDirs[2] = ray.mDirection.z < 0.f ? 1 : 0;
-  
-  // we have to check the node
-  // current node is not the leaf, empty leaves are NULL pointers
-  while (index) {
-    SKTBNodeT *currNode = stack3[index].nodep;
-    tmin = tmax;
-    tmax = stack3[index].tmax;
-#if 0
-    if (debug) {
-      cout << "node = " << (void*)currNode
-	   << " tmin = " << tmin << " tmax = " << tmax
-	   << endl;
-    }
-#endif
-    
-    assert(tmin <= tmax);
-    
-#ifdef __TRAVERSAL_STATISTICS
-    allNodesTraversed++;
-#endif // __TRAVERSAL_STATISTICS
-    const unsigned int nodeType = GetNodeType(currNode);
-    if (nodeType < CKTBAxes::EE_Leaf) {
-      float tval = (GetSplitValue(currNode) - ray.mOrigin[nodeType]);
-      tval *= invertedDir[nodeType];
-      SKTBNodeT *near, *far;
-      childNodes[0] = GetLeft(currNode);
-      childNodes[1] = GetRight(currNode);
-      int rayDir = RayDirs[nodeType];
-      near = childNodes[rayDir];
-      far = childNodes[rayDir ^ 0x1];
-      // This code is slower than above !!!!!
-      stack3[index].nodep = far;
-      // stack3[index].tmax = tmax; // this is already there, not necessary!
-      __m128 tsim = _mm_set_ss(tval);
-      // index += tval < tmax ? 1 : 0;
-      __m128 tother = _mm_set_ss(tmax);
-      index += _mm_ucomilt_ss(tsim, tother);
-      stack3[index].nodep = near;
-      // stack3[index].tmax = Min(tval, tmax);
-      _mm_store_ss(&(stack3[index].tmax), _mm_min_ps(tsim, tother));
-      tmax = tmin;
-      // index += tval < tmin ? -1 : 0;
-      __m128 tother2 = _mm_set_ss(tmin);
-      index -= _mm_ucomilt_ss(tsim, tother2);
-    }
-    else {
-      if (nodeType == CKTBAxes::EE_Leaf) {
-	// test objects for intersection
-#ifdef _DEBUGKTB
-	cout << "Leaf " << endl;
-	depth++;
-#endif
-#ifdef _DEBUGKTB
-	DEBUG << "currNode = " << currNode << " entp.t = " << entp->t
-	      << " extp.t = " << extp->t << endl;
-#endif
-	if (!IsEmptyLeaf_(currNode)) {
-#ifdef _DEBUGKTB
-	  cout << "Full leaf at depth= " << depth << endl; 
-#endif      
-
-#ifdef __TRAVERSAL_STATISTICS
-	  fullLeavesTraversed++;
-#endif // __TRAVERSAL_STATISTICS
-	  // test the objects in the full leaf against the ray
-	  
-	  SimpleRay::IntersectionRes[0].maxt = stack3[index].tmax;
-	  if (TestFullLeaf(ray, currNode)) {
-#ifdef _DEBUGKTB
-	    cout << "Full leaf HIT " << endl; 
-#endif 	
-	      
-#ifdef __TRAVERSAL_STATISTICS
-	    _allNodesTraversed += allNodesTraversed;
-	    _fullLeavesTraversed += fullLeavesTraversed;
-	    _emptyLeavesTraversed += emptyLeavesTraversed;
-#endif // __TRAVERSAL_STATISTICS
-	    
-	    // signed distance should be already set in TestFullLeaf
-	    // the first object intersected was found	
-	    return 1;
-	  }
-	} // full leaf
-#ifdef __TRAVERSAL_STATISTICS
-	else {
-#ifdef _DEBUGKTB
-	  cout << "Empty leaf at depth= " << depth << endl; 
-#endif      
-	  emptyLeavesTraversed++;
-	}
-#endif // __TRAVERSAL_STATISTICS
-
-#ifdef _DEBUGKTB
-	cout << "Pop the node" << endl;
-#endif      
-    
-	// pop farChild from the stack
-	// restore the current values
-	index--;
-	continue;
-      }
-      else {
-	assert(nodeType == CKTBAxes::EE_Link);
-#ifdef _DEBUGKTB
-	cout << "Link " << endl;
-#endif
-	stack3[index].nodep = GetLinkNode(currNode);
-	// cout << "Link node was accessed" << endl;
-	continue;
-      }
-    }
-  } // while current node is not the leaf
-
-#ifdef __TRAVERSAL_STATISTICS
-  _allNodesTraversed += allNodesTraversed;
-  _fullLeavesTraversed += fullLeavesTraversed;
-  _emptyLeavesTraversed += emptyLeavesTraversed;
-#endif // __TRAVERSAL_STATISTICS
-
-  // no objects found along the ray path
-  return 0;
-} // FindNearestI - single ray
-#endif
-
-#if 0
+
+
 // Reasonably fast - about 101,500 rays per second for single dir!
-// It allows fast switching context from one ray to the next ray so it it
+// It allows fast switching context from one ray to the next ray so it is
 // virtually independent of memory latency !
 int
-CKTBTraversal::FindNearestI_16oneDir(SimpleRayContainer &rays, int offset)
+CKTBTraversal::FindNearestI_16oneDirNoSSE(SimpleRayContainer &rays, int offset)
 {
   // passing through parameters
@@ -742,5 +442,7 @@
   float tmin, tmax;
   for (int i = 0; i < cntMaxRays; i++) {
-  // test if the whole CKTB tree is missed by the input ray
+    // Setting zero intersection as original result
+    SimpleRay::IntersectionRes[i+rayOffset].intersectable = 0;
+    // test if the whole CKTB tree is missed by the input ray
     if ((!bbox.ComputeMinMaxT(rays[i+offset].mOrigin,
 			      rays[i+offset].mDirection,
@@ -748,5 +450,4 @@
 	(tmax <= tmin) ||
 	(tmax <= 0.f) ) {
-      SimpleRay::IntersectionRes[i].intersectable = 0;
     }
     else {
@@ -756,7 +457,7 @@
       rayOrig[indexR + 2] = rays[i+offset].mOrigin.z;
       //rayOrig[indexR + 3] = 0.f;
-      invertedDir[indexR + 0] = 1.0f / (rays[i+offset].mDirection.x);
-      invertedDir[indexR + 1] = 1.0f / (rays[i+offset].mDirection.y);
-      invertedDir[indexR + 2] = 1.0f / (rays[i+offset].mDirection.z);
+      invertedDir[indexR + 0] = 1.0f / (rays[i+offset].mDirection.x - 1e-25f);
+      invertedDir[indexR + 1] = 1.0f / (rays[i+offset].mDirection.y - 1e-25f);
+      invertedDir[indexR + 2] = 1.0f / (rays[i+offset].mDirection.z - 1e-25f);
       //invertedDir[indexR + 2] = 0.f;
       rayDirs[indexR + 0] = rays[i+offset].mDirection.x < 0.f ? 1 : 0;
@@ -819,5 +520,4 @@
       }
 #endif      
-      assert(tmin <= tmax);
     
 #ifdef __TRAVERSAL_STATISTICS
@@ -876,5 +576,6 @@
 	    // which ray is processed
 	    int indexR = indexRay[indexA];
-	    SimpleRay::IntersectionRes[indexR].maxt = tmax;
+	    SimpleRay::IntersectionRes[indexR + rayOffset].maxt =
+	      tmax + Limits::Small;
 	    if (TestFullLeaf(rays[indexR+offset], currNode, indexR)) {
 
@@ -934,4 +635,5 @@
 #endif
 	  stackA[indexSA].nodep = GetLinkNode(currNode);
+	  GPREFETCH(stackA[indexSA].nodep, PREF_DEFAULT);
 	  // cout << "Link node was accessed" << endl;
 	  continue;
@@ -950,588 +652,4 @@
   return cntHits;
 }
-#endif
-
-#ifdef __SSE__
-
-#if 1
-// Even faster - about 125,500 rays per second for single dir and 164 rps
-// for double dir !
-int
-CKTBTraversal::FindNearestI_16oneDir(SimpleRayContainer &rays, int offset)
-{
-  static RayPacket2x2 raypack;
-  struct SResultI {
-    Intersectable *intersectable;
-    float          tdist;    
-  };
-  static SResultI results[16];
-    
-  for (int i = 0; i < 4; i++) {
-    int k = i * 4 + offset;
-    for (int j = 0; j < 4; j++, k++) {
-      raypack.SetLoc(j, rays[k].mOrigin);
-      raypack.SetDir(j, rays[k].mDirection);      
-    }
-    // Here either use ray packet traversal or
-    // casting individual rays
-    FindNearestI(raypack);
-    k = i * 4;
-    for (int j = 0; j < 4; j++, k++) {
-      results[k].intersectable = raypack.GetObject(j);
-      results[k].tdist = raypack.GetT(j);
-    } // for j
-  } // for i
-
-  // Copy the results to the output array
-  for (int i = 0; i < 16; i++) {
-    SimpleRay::IntersectionRes[i].intersectable =
-      results[i].intersectable;
-    SimpleRay::IntersectionRes[i].tdist =
-      results[i].tdist;
-  } // for i
-  return 0;
-}
-#endif
-  
-#if 0
-// This code works well 1/1/2008 - 11:00
-// The same operations for packets of rays for the same signs,
-// otherwise it is emulated by decomposition
-// of a packet to individual rays and traced individually.
-void
-CKTBTraversal::FindNearestI(RayPacket2x2 &rp)
-{
-  int m1 = _mm_movemask_ps(rp.dx4);
-  if ((m1 == 0)||(m1 == 15)) {
-  m1 = _mm_movemask_ps(rp.dy4);
-  if ((m1 == 0)||(m1 == 15)) {
-  m1 = _mm_movemask_ps(rp.dz4);
-  if ((m1 == 0)||(m1 == 15)) {
-    rp.Init();
-    // all the signs for 4 rays are the same, use
-    // ray packet traversal
-    // Compute min and max distances
-    GALIGN16 union { float tmin4[4]; __m128 tmin_4; };
-    GALIGN16 union { float tmax4[4]; __m128 tmax_4; };
-    SimpleRay sray[4];
-    int maxIntersections = 4;
-    GALIGN16 union { int inters[4]; __m128 inters_4; };
-    inters[0] = inters[1] = inters[2] = inters[3] = 1;
-    unsigned int inters32 = 0xf;
-    for (int i = 0; i < 4; i++) {
-      bbox.ComputeMinMaxT(rp.GetLoc(i), rp.GetDir(i), &(tmin4[i]), &(tmax4[i]));
-      if ( (tmin4[i] >= tmax4[i]) ||
-	   (tmax4[i] < 0.f) ) {
-	inters[i] = 0; // finished
-	inters32 &= ~(1 << i); // bit zero when ray is invalid
-	maxIntersections--;
-      }      
-      if (tmin4[i] < 0.f)
-	tmin4[i] = 0.f;
-      sray[i].mOrigin = rp.GetLoc(i);
-      sray[i].mDirection = rp.GetDir(i);
-    } // for i
-    if (maxIntersections == 0)
-      return;
-
-    SKTBNodeT * childNodes[2];
-    int RayDirs[3];
-    RayDirs[0] = (rp.dx[0] > 0.f) ? 1 : 0;
-    RayDirs[1] = (rp.dy[0] > 0.f) ? 1 : 0;
-    RayDirs[2] = (rp.dz[0] > 0.f) ? 1 : 0;
-    //int activeMask=_mm_movemask_ps(_mm_cmplt_ps( tmin_4, tmax_4 ))&inters32;
-    int activeMask = inters32;
-    int indexStack = 0;
-    SKTBNodeT *currNode = root;
-    unsigned int k = GetNodeType(currNode);
-    for (;;) {      
-      while (k < CKTBAxes::EE_Leaf) {
-	// the 3 operations below can be brought down to 3 simple float
-	// calculations by precomputing min/max of the inverse dir
-	const __m128 node_split = _mm_set_ps1(GetSplitValue(currNode));
-	const __m128 t4 =
-	  _mm_mul_ps(_mm_sub_ps(node_split, rp.orig[k]), rp.idir[k]);
-	childNodes[0] = GetLeft(currNode);
-	childNodes[1] = GetRight(currNode);
-	int rayDir = RayDirs[k];
-        SKTBNodeT *far = childNodes[rayDir];
-	if (!(_mm_movemask_ps(_mm_cmpgt_ps(t4, tmin_4)) & activeMask))
-	{
-	  currNode = far;
-	  k = GetNodeType(currNode);
-	  continue;
-	}
-	currNode = childNodes[rayDir ^ 0x1]; // this is near node
-	k = GetNodeType(currNode);	
-	if (! (_mm_movemask_ps(_mm_cmplt_ps( t4, tmax_4)) & activeMask))
-	  continue;
-
-	// pop far node to the stack
-	stack4[indexStack].nodep = far;
-	stack4[indexStack].tmax_4 = tmax_4;
-	stack4[indexStack].tmin_4 = _mm_max_ps(t4, tmin_4);
-	// stack4[indexStack].mask = activeMask;
-	indexStack++;
-	
-	tmax_4 = _mm_min_ps(t4, tmax_4); 
-	activeMask &= _mm_movemask_ps(_mm_cmplt_ps( tmin_4, tmax_4 ));
-      } // while this is an interior node
-
-      // either a leaf or a link
-      if (k == CKTBAxes::EE_Leaf) {
-	// test objects for intersection
-	if (!IsEmptyLeaf_(currNode)) {
-	  // cout << "Full leaf" << endl;
-	  
-	  // test the objects in the full leaf against the ray	  
-	  for (int i = 0; i < 4; i++) {
-	    if (inters[i] ) {
-	      // no intersection so far !
-	      SimpleRay::IntersectionRes[i].maxt = tmax4[i];
-	      // Test only rays that were not finished
-	      if (TestFullLeaf(sray[i], currNode, i)) {
-		// intersection for this ray found
-		inters[i] = 0;
-		inters32 &= ~(1 << i);
-		rp.SetT(i, SimpleRay::IntersectionRes[0].maxt);
-		rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);		
-		// signed distance should be already set in TestFullLeaf
-		// the first object intersected was found
-		if (--maxIntersections == 0)
-		  return;
-	      }
-	    } // if this ray did not hit the triangle so far
-	  } // for all 4 rays
-	} // full leaf
-	// pop farChild from the stack
-	// restore the current values
-	// update the minimum distance since we traverse to the next one
-
-	if (indexStack == 0)
-	  return;
-	indexStack--;
-	currNode = stack4[indexStack].nodep;
-	k = GetNodeType(currNode);
-	tmin_4 = stack4[indexStack].tmin_4;
-	tmax_4 = stack4[indexStack].tmax_4;
-	activeMask = _mm_movemask_ps(_mm_cmple_ps( tmin_4, tmax_4 )) & inters32;
-	continue;
-      }
-      // cout << "Link node was accessed" << endl;
-      assert(k == CKTBAxes::EE_Link);
-      currNode = GetLinkNode(currNode);
-      k = GetNodeType(currNode);
-    } // for
-    return;
-  }}}
-
-  // Trace ray by ray
-  SimpleRay ray;
-  for (int i = 0; i < 4; i++) {
-    ray.mOrigin = rp.GetLoc(i);
-    ray.mDirection = rp.GetDir(i);
-    FindNearestI(ray);
-    rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-    rp.SetT(i, SimpleRay::IntersectionRes[0].tdist);
-    // SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
-  } // for
-}
-#endif
-
-
-#if 1
-// This code also works well 1/1/2008 - 14:00
-// Using mask of 128-bits width - the code works as well, only a bit
-// faster than the code above
-void
-CKTBTraversal::FindNearestI(RayPacket2x2 &rp)
-{
-  int m1 = _mm_movemask_ps(rp.dx4);
-  if ((m1 == 0)||(m1 == 15)) {
-  m1 = _mm_movemask_ps(rp.dy4);
-  if ((m1 == 0)||(m1 == 15)) {
-  m1 = _mm_movemask_ps(rp.dz4);
-  if ((m1 == 0)||(m1 == 15)) {
-    rp.Init();
-    // all the signs for 4 rays are the same, use
-    // ray packet traversal
-    // Compute min and max distances
-    GALIGN16 union { float tmin4[4]; __m128 tmin_4; };
-    GALIGN16 union { float tmax4[4]; __m128 tmax_4; };
-    GALIGN16 union { float activeMask[4]; __m128 activeMask_4; };
-    GALIGN16 union { float liveMask[4]; __m128 liveMask_4; };
-    liveMask[0] = liveMask[1] = liveMask[2] = liveMask[3] = 0xffffffff;
-
-    GALIGN16 SimpleRay sray[4];
-    int maxIntersections = 4;
-    // unsigned int inters32 = 0xf;
-    for (int i = 0; i < 4; i++) {
-      rp.SetObject(i, 0);
-      bbox.ComputeMinMaxT(rp.GetLoc(i), rp.GetDir(i), &(tmin4[i]), &(tmax4[i]));
-      if ( (tmin4[i] >= tmax4[i]) ||
-	   (tmax4[i] < 0.f) ) {
-	liveMask[i] = 0; // finished
-	// inters32 &= ~(1 << i); // bit zero when ray is invalid
-	maxIntersections--;
-      }
-      if (tmin4[i] < 0.f)
-	tmin4[i] = 0.f;
-      sray[i].mOrigin = rp.GetLoc(i);
-      sray[i].mDirection = rp.GetDir(i);
-    } // for i
-    if (maxIntersections == 0)
-      return;
-
-    // This is the mask 128 bits witdth
-    //activeMask_4 =
-    //  _mm_and_ps(_mm_cmple_ps(tmin_4, tmax_4),
-    //		 _mm_cmplt_ps(tmax_4, _mm_setzero_ps()));
-    activeMask_4 = liveMask_4;
-
-    SKTBNodeT * childNodes[2];
-    int RayDirs[4];
-    RayDirs[0] = (rp.dx[0] > 0.f) ? 1 : 0;
-    RayDirs[1] = (rp.dy[0] > 0.f) ? 1 : 0;
-    RayDirs[2] = (rp.dz[0] > 0.f) ? 1 : 0;
-    int indexStack = 0;
-    SKTBNodeT *currNode = root;
-    unsigned int k = GetNodeType(currNode);
-    for (;;) {
-      // traverse until we find a leaf
-      while (k < CKTBAxes::EE_Leaf) {
-	// the 3 operations below can be brought down to 3 simple float
-	// calculations by precomputing min/max of the inverse dir
-	// const __m128 node_split = ;
-	const __m128 t4 =
-	  _mm_mul_ps(_mm_sub_ps(_mm_set_ps1(GetSplitValue(currNode)),
-				rp.orig[k]), rp.idir[k]);
-	childNodes[0] = GetLeft(currNode);
-	childNodes[1] = GetRight(currNode);
-	int rayDir = RayDirs[k];
-        SKTBNodeT *far = childNodes[rayDir];
-	if (_mm_movemask_ps(_mm_and_ps(_mm_cmpge_ps(t4, tmin_4),
-				       activeMask_4))) {	  
-	  currNode = far;
-	  k = GetNodeType(currNode);
-	  continue;
-	}
-
-	currNode = childNodes[rayDir ^ 0x1]; // this is near node
-	k = GetNodeType(currNode);	
-	if (_mm_movemask_ps(_mm_and_ps(_mm_cmple_ps(t4, tmax_4),
-				       activeMask_4)))
-	  continue;
-
-	// pop far node to the stack
-	stack4[indexStack].nodep = far;
-	stack4[indexStack].tmax_4 = tmax_4;
-
-// Uncomenting this macro is unsafe!
-// Not convinced if for packet of 4 rays we can say that since when
-// one ray is different than the others, it could bring to wrong state
-// It is surely true for one ray when tmin < t < tmax, but for a packet
-// of rays this condition can be true only for a single ray
-//   tmin4 = max(t4, tmin4) = min(t4, tmax4)
-//#define _NOT_STORE_MINT
-
-#ifdef _NOT_STORE_MINT	
-#else
-	// store mint onto the stack
-	stack4[indexStack].tmin_4 = _mm_max_ps(t4, tmin_4);
-#endif	
-	// stack4[indexStack].mask = activeMask;
-	indexStack++;
-	
-	tmax_4 = _mm_min_ps(t4, tmax_4); 
-	activeMask_4 = _mm_cmplt_ps( tmin_4, tmax_4 );
-      } // while this is an interior node
-
-      // either a leaf or a link
-      if (k == CKTBAxes::EE_Leaf) {
-	// test objects for intersection
-	if (!IsEmptyLeaf_(currNode)) {
-	  // cout << "Full leaf" << endl;
-	  
-	  // test the objects in the full leaf against the ray	  
-	  for (int i = 0; i < 4; i++) {
-	    if (liveMask[i] ) {
-	      // no intersection so far !
-	      SimpleRay::IntersectionRes[i].maxt = tmax4[i];
-#if 0
-	      // Using subroutine
-	      // Test only rays that were not finished
-	      if (TestFullLeaf(sray[i], currNode, i))
-#else
-	      // avoiding one call
-	      const ObjectContainer * const list = GetObjList(currNode); 
-	      int intersected = 0;
-	      // iterate the whole list and find out the nearest intersection
-	      ObjectContainer::const_iterator sc_end = list->end();
-	      for (ObjectContainer::const_iterator sc = list->begin(); sc != sc_end; sc++) {
-		// if the intersection realy lies in the node        
-		intersected |= ((*sc)->CastSimpleRay(sray[i], i));
-	      } // for all objects
-	      if (intersected)
-#endif
-	      {
-		rp.SetT(i, SimpleRay::IntersectionRes[0].maxt);
-		rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-		// signed distance should be already set in TestFullLeaf
-		// the first object intersected was found
-		if (--maxIntersections == 0)
-		  return;
-		// inters32 &= ~(1 << i);
-		liveMask[i] = 0;
-	      }
-	    } // if this ray did not hit the triangle so far
-	  } // for all 4 rays
-	} // full leaf
-
-	// pop farChild from the stack
-	// restore the current values
-	// update the minimum distance since we traverse to the next one
-	do {
-	  if (indexStack == 0)
-	    return;
-	  indexStack--;
-	  currNode = stack4[indexStack].nodep;
-	  k = GetNodeType(currNode);
-#ifdef _NOT_STORE_MINT
-	  // this is an attempt !
-	  tmin_4 = tmax_4;
-#else
- 	  // This surrely works
-	  tmin_4 = stack4[indexStack].tmin_4;
-#endif	  
-	  tmax_4 = stack4[indexStack].tmax_4;
-	  activeMask_4 = _mm_and_ps(_mm_cmple_ps( tmin_4, tmax_4 ), liveMask_4);
-	}
-	while (_mm_movemask_ps(activeMask_4) == 0);
-      }
-      else {
-	// cout << "Link node was accessed" << endl;
-	assert(k == CKTBAxes::EE_Link);
-	currNode = GetLinkNode(currNode);
-	k = GetNodeType(currNode);
-      }
-    } // for(;;)
-    return;
-  }}}
-
-  // Trace ray by ray
-  SimpleRay ray;
-  for (int i = 0; i < 4; i++) {
-    ray.mOrigin = rp.GetLoc(i);
-    ray.mDirection = rp.GetDir(i);
-    FindNearestI(ray);
-    rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-    rp.SetT(i, SimpleRay::IntersectionRes[0].tdist);
-    // SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
-  } // for
-}
-#endif
-
-#if 1
-// This code also works well 1/1/2008 - 14:00
-// Using mask of 128-bits width - the code works as well, only a bit
-// faster than the code above
-void
-CKTBTraversal::FindNearestI(RayPacket2x2 &rp, Vector3 &boxmin, Vector3 &boxmax)
-{
-
-  static AxisAlignedBox3 localbox;
-  localbox.SetMin(boxmin);
-  localbox.SetMax(boxmax);
-
-#define USE_SIMPLE_VERSION 1
-
-#if !USE_SIMPLE_VERSION
-  int m1 = _mm_movemask_ps(rp.dx4);
-  if ((m1 == 0)||(m1 == 15)) {
-  m1 = _mm_movemask_ps(rp.dy4);
-  if ((m1 == 0)||(m1 == 15)) {
-  m1 = _mm_movemask_ps(rp.dz4);
-  if ((m1 == 0)||(m1 == 15)) {
-    rp.Init();
-    
-    // all the signs for 4 rays are the same, use
-    // ray packet traversal
-    // Compute min and max distances
-    GALIGN16 union { float tmin4[4]; __m128 tmin_4; };
-    GALIGN16 union { float tmax4[4]; __m128 tmax_4; };
-    GALIGN16 union { float activeMask[4]; __m128 activeMask_4; };
-    GALIGN16 union { float liveMask[4]; __m128 liveMask_4; };
-    liveMask[0] = liveMask[1] = liveMask[2] = liveMask[3] = 0xffffffff;
-
-    GALIGN16 SimpleRay sray[4];
-    int maxIntersections = 4;
-    // unsigned int inters32 = 0xf;
-    for (int i = 0; i < 4; i++) {
-      rp.SetObject(i, 0);
-      localbox.ComputeMinMaxT(rp.GetLoc(i), rp.GetDir(i), &(tmin4[i]), &(tmax4[i]));
-      if ( (tmin4[i] >= tmax4[i]) ||
-	   (tmax4[i] < 0.f) ) {
-	liveMask[i] = 0; // finished
-	// inters32 &= ~(1 << i); // bit zero when ray is invalid
-	maxIntersections--;
-      }
-      if (tmin4[i] < 0.f)
-	tmin4[i] = 0.f;
-      sray[i].mOrigin = rp.GetLoc(i);
-      sray[i].mDirection = rp.GetDir(i);
-    } // for i
-    if (maxIntersections == 0)
-      return;
-
-    // This is the mask 128 bits witdth
-    //activeMask_4 =
-    //  _mm_and_ps(_mm_cmple_ps(tmin_4, tmax_4),
-    //		 _mm_cmplt_ps(tmax_4, _mm_setzero_ps()));
-    activeMask_4 = liveMask_4;
-
-    SKTBNodeT * childNodes[2];
-    int RayDirs[4];
-    RayDirs[0] = (rp.dx[0] > 0.f) ? 1 : 0;
-    RayDirs[1] = (rp.dy[0] > 0.f) ? 1 : 0;
-    RayDirs[2] = (rp.dz[0] > 0.f) ? 1 : 0;
-    int indexStack = 0;
-    SKTBNodeT *currNode = root;
-    unsigned int k = GetNodeType(currNode);
-    for (;;) {
-      // traverse until we find a leaf
-      while (k < CKTBAxes::EE_Leaf) {
-	// the 3 operations below can be brought down to 3 simple float
-	// calculations by precomputing min/max of the inverse dir
-	// const __m128 node_split = ;
-	const __m128 t4 =
-	  _mm_mul_ps(_mm_sub_ps(_mm_set_ps1(GetSplitValue(currNode)),
-				rp.orig[k]), rp.idir[k]);
-	childNodes[0] = GetLeft(currNode);
-	childNodes[1] = GetRight(currNode);
-	int rayDir = RayDirs[k];
-        SKTBNodeT *far = childNodes[rayDir];
-	if (_mm_movemask_ps(_mm_and_ps(_mm_cmpge_ps(t4, tmin_4),
-				       activeMask_4))) {	  
-	  currNode = far;
-	  k = GetNodeType(currNode);
-	  continue;
-	}
-
-	currNode = childNodes[rayDir ^ 0x1]; // this is near node
-	k = GetNodeType(currNode);	
-	if (_mm_movemask_ps(_mm_and_ps(_mm_cmple_ps(t4, tmax_4),
-				       activeMask_4)))
-	  continue;
-
-	// pop far node to the stack
-	stack4[indexStack].nodep = far;
-	stack4[indexStack].tmax_4 = tmax_4;
-
-// Uncomenting this macro is unsafe!
-// Not convinced if for packet of 4 rays we can say that since when
-// one ray is different than the others, it could bring to wrong state
-// It is surely true for one ray when tmin < t < tmax, but for a packet
-// of rays this condition can be true only for a single ray
-//   tmin4 = max(t4, tmin4) = min(t4, tmax4)
-//#define _NOT_STORE_MINT
-
-#ifdef _NOT_STORE_MINT	
-#else
-	// store mint onto the stack
-	stack4[indexStack].tmin_4 = _mm_max_ps(t4, tmin_4);
-#endif	
-	// stack4[indexStack].mask = activeMask;
-	indexStack++;
-	
-	tmax_4 = _mm_min_ps(t4, tmax_4); 
-	activeMask_4 = _mm_cmplt_ps( tmin_4, tmax_4 );
-      } // while this is an interior node
-
-      // either a leaf or a link
-      if (k == CKTBAxes::EE_Leaf) {
-	// test objects for intersection
-	if (!IsEmptyLeaf_(currNode)) {
-	  // cout << "Full leaf" << endl;
-	  
-	  // test the objects in the full leaf against the ray	  
-	  for (int i = 0; i < 4; i++) {
-	    if (liveMask[i] ) {
-	      // no intersection so far !
-	      SimpleRay::IntersectionRes[i].maxt = tmax4[i];
-#if 0
-	      // Using subroutine
-	      // Test only rays that were not finished
-	      if (TestFullLeaf(sray[i], currNode, i))
-#else
-	      // avoiding one call
-	      const ObjectContainer * const list = GetObjList(currNode); 
-	      int intersected = 0;
-	      // iterate the whole list and find out the nearest intersection
-	      ObjectContainer::const_iterator sc_end = list->end();
-	      for (ObjectContainer::const_iterator sc = list->begin(); sc != sc_end; sc++) {
-		// if the intersection realy lies in the node        
-		intersected |= ((*sc)->CastSimpleRay(sray[i], i));
-	      } // for all objects
-	      if (intersected)
-#endif
-	      {
-		rp.SetT(i, SimpleRay::IntersectionRes[0].maxt);
-		rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-		// signed distance should be already set in TestFullLeaf
-		// the first object intersected was found
-		if (--maxIntersections == 0)
-		  return;
-		// inters32 &= ~(1 << i);
-		liveMask[i] = 0;
-	      }
-	    } // if this ray did not hit the triangle so far
-	  } // for all 4 rays
-	} // full leaf
-
-	// pop farChild from the stack
-	// restore the current values
-	// update the minimum distance since we traverse to the next one
-	do {
-	  if (indexStack == 0)
-	    return;
-	  indexStack--;
-	  currNode = stack4[indexStack].nodep;
-	  k = GetNodeType(currNode);
-#ifdef _NOT_STORE_MINT
-	  // this is an attempt !
-	  tmin_4 = tmax_4;
-#else
- 	  // This surrely works
-	  tmin_4 = stack4[indexStack].tmin_4;
-#endif	  
-	  tmax_4 = stack4[indexStack].tmax_4;
-	  activeMask_4 = _mm_and_ps(_mm_cmple_ps( tmin_4, tmax_4 ), liveMask_4);
-	}
-	while (_mm_movemask_ps(activeMask_4) == 0);
-      }
-      else {
-	// cout << "Link node was accessed" << endl;
-	assert(k == CKTBAxes::EE_Link);
-	currNode = GetLinkNode(currNode);
-	k = GetNodeType(currNode);
-      }
-    } // for(;;)
-    return;
-  }}}
-#endif
-  // Trace ray by ray
-  SimpleRay ray;
-  for (int i = 0; i < 4; i++) {
-    ray.mOrigin = rp.GetLoc(i);
-    ray.mDirection = rp.GetDir(i);
-    FindNearestI(ray, localbox);
-    rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-    rp.SetT(i, SimpleRay::IntersectionRes[0].tdist);
-    // SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
-  } // for
-}
-#endif
-
-#endif // __SSE__
 
 #endif //  TRV00F
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.h	(revision 2592)
@@ -156,7 +156,4 @@
 #endif
   
-  // init the stack of auxiliary variables from min to max
-  virtual void InitAux(int min, int max, int maxItemsAtOnce);
-
   // the direction of traversal through CKTB tree
   enum EDirection { LEFT = 0, RIGHT = 1, NOWHERE = 2 };
@@ -180,4 +177,7 @@
   virtual ~CKTBAllocMan() { }
 
+  // init the stack of auxiliary variables from min to max
+  virtual void InitAux(int min, int max, int maxItemsAtOnce);
+  
   // forget the content that is created by previous kd-tree construction
   // or just init for the first use.
@@ -224,4 +224,8 @@
   void SetInteriorNodeLinks(SKTBNodeT *node,
 			    SKTBNodeT *leftChild,
+			    SKTBNodeT *rightChild);
+  void SetInteriorNodeLeft(SKTBNodeT *node,
+			   SKTBNodeT *leftChild);
+  void SetInteriorNodeRight(SKTBNodeT *node,
 			    SKTBNodeT *rightChild);
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.cpp	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.cpp	(revision 2592)
@@ -1221,10 +1221,10 @@
   if (objlist.size() == 0)
     return 0; // nothing
-  cerr<<"hh44"<<endl;
+  // cerr<<"hh44"<<endl;
   // initialize the whole box of the kd-tree and sort
   // the boundary entries
   SInputData *d = Init(&objlist, box);
 
-  cerr<<"hh45"<<endl;
+  //cerr<<"hh45"<<endl;
 
   // copy verbose to be used consistenly further on
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbconf.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbconf.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbconf.h	(revision 2592)
@@ -22,4 +22,7 @@
 #include <xmmintrin.h>
 #endif
+
+// If we support the use of SSE instructions for ray shooting
+#define _USE_HAVRAN_SSE
 
 namespace GtpVisibilityPreprocessor {
Index: GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.h	(revision 2592)
@@ -5,6 +5,6 @@
 #include "Containers.h"
 #include <string>
-#include "havran/raypack.h"
-
+#include "ktbconf.h"
+#include "raypack.h"
 
 namespace GtpVisibilityPreprocessor {
@@ -15,5 +15,4 @@
 class KdTree;
 class Ray;
-class RayPacket2x2;
 class SimpleRayContainer;
 class AxisAlignedBox3;
@@ -38,14 +37,4 @@
   
   int Type() const { return HAVRAN_RAYCASTER; }
-
-  // Using packet of 4 rays supposing that these are coherent
-  // We give a box to which each ray is clipped to before the
-  // ray shooting is computed !
-  virtual void CastRaysPacket4(Vector3 &minBox,
-			       Vector3 &maxBox,
-			       Vector3 origin4[],
-			       Vector3 direction4[],
-			       int     result4[],
-			       float   dist4[]);
   
   virtual int CastRay(
@@ -82,4 +71,25 @@
 			const bool castDoubleRay,
 			const bool pruneInvalidRays = true);
+
+  // Using packet of 4 rays supposing that these are coherent
+  // We give a box to which each ray is clipped to before the
+  // ray shooting is computed !
+  virtual void CastRaysPacket4(Vector3 &minBox,
+			       Vector3 &maxBox,
+			       Vector3 origin4[],
+			       Vector3 direction4[],
+			       int     result4[],
+			       float   dist4[]);
+
+#ifdef _USE_HAVRAN_SSE  
+  // Just for testing concept
+  virtual void CastRaysPacket2x2(RayPacket2x2 &raysPack,
+				 bool castDoubleRay,
+				 const bool pruneInvalidRays = true);
+#endif  
+
+  bool ExportBinTree(const string &filename);
+  bool ImportBinTree(const string &filename, ObjectContainer &objects);
+  
 protected:
   CKTB *mKtbtree;
Index: GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro	(revision 2583)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro	(revision 2592)
@@ -129,5 +129,5 @@
 SOURCES += havran/allocgo2.cpp havran/ktbai.cpp havran/ktbtrav.cpp \
 havran/ktb.cpp havran/ktball.cpp havran/sbbox.cpp \
-havran/ktb8b.cpp havran/ktbftrav.cpp havran/timer.cpp
+havran/ktb8b.cpp havran/ktbftrav.cpp havran/ktbf2trv.cpp havran/timer.cpp
 
 
