Index: GTP/trunk/Lib/Vis/Preprocessing/src/OcclusionQuery.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/OcclusionQuery.cpp	(revision 2612)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/OcclusionQuery.cpp	(revision 2629)
@@ -12,5 +12,8 @@
 #endif
 
+#ifdef _WIN32
+// This Macro does not compile under LINUX
 #define _ARBGL
+#endif
 
 using namespace std;
Index: GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h	(revision 2621)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.h	(revision 2629)
@@ -25,5 +25,4 @@
 struct SimpleRay;
 class RayPacket2x2;
-class Matrix4x4;
 
 /** This class provides an interface for ray casting.
@@ -37,5 +36,6 @@
     INTERNAL_RAYCASTER = 0,
     INTEL_RAYCASTER = 1,
-    HAVRAN_RAYCASTER = 2
+    HAVRAN_RAYCASTER = 2,
+    HAVRAN_DYN_RAYCASTER = 3
   };
 
@@ -72,4 +72,9 @@
 			const bool pruneInvalidRays = true);
 
+  virtual void
+  CastSimpleForwardRays(SimpleRayContainer &rays,
+			const AxisAlignedBox3 &sbox
+			) { return;}
+  
   // Using packet of 4 rays supposing that these are coherent
   virtual void CastRaysPacket4(Vector3 origin4[],
@@ -82,9 +87,9 @@
   // ray shooting is computed !
   virtual void CastRaysPacket4(const Vector3 &minBox,
-							   const Vector3 &maxBox,
-							   const Vector3 origin4[],
-							   const Vector3 direction4[],
-							   int result4[],
-							   float dist4[]) {  }
+			       const Vector3 &maxBox,
+			       const Vector3 origin4[],
+			       const Vector3 direction4[],
+			       int result4[],
+			       float dist4[]) {  }
   
   // Just for testing concept
@@ -93,4 +98,14 @@
 				 const bool pruneInvalidRays = true)
   { }
+
+  virtual void AddDynamicObjecs(const ObjectContainer &objects, const Matrix4x4 &m)
+  { }
+
+  virtual void UpdateDynamicObjects(const Matrix4x4 &m)
+  { }
+
+  virtual void DeleteDynamicObjects()
+  { }
+
   
   /*virtual void CastRaysEye4(SimpleRayContainer &rays,
@@ -101,12 +116,12 @@
   */
 
+  // This sorts only rays by origin
   virtual void
   SortRays(SimpleRayContainer &rays);
 
-
-  virtual void AddDynamicObjecs(const ObjectContainer &objects, const Matrix4x4 &m) = 0;
-  virtual void UpdateDynamicObjects(const Matrix4x4 &m) = 0;
-  virtual void DeleteDynamicObjects() = 0;
-
+  // This sorts the ray by origin and direction
+  virtual void
+  SortRays2(SimpleRayContainer &rays);
+  
   // pool of vss rays to be used in one pass of the sampling
   struct VssRayPool 
@@ -173,4 +188,10 @@
 		 const int depth,
 		 float box[12]);
+
+  void _SortRays2(SimpleRayContainer &rays,
+		  const int l,
+		  const int r,
+		  const int depth,
+		  float box[12]);
   
   struct Intersection
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 2623)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp	(revision 2629)
@@ -1045,13 +1045,26 @@
 	GetIntValue("Preprocessor.rayCastMethod", rayCastMethod);
 
-  if (rayCastMethod == 2) {
-	//HavranRayCaster *hr =
-	//  dynamic_cast<HavranRayCaster*>(mRayCaster);
-	HavranRayCaster *hr =
-	  reinterpret_cast<HavranRayCaster*>(mRayCaster);
-	hr->Build(this->mObjects);
-	// $$ do not exit here, so that the "internal" kD-tree used 
-	// other puroposes can be loaded
-	//	return true;
+  if ((rayCastMethod == 2) || (rayCastMethod == 3)) {
+    HavranRayCaster *hr = 0;
+    if (rayCastMethod == 3)
+      hr = reinterpret_cast<HavranDynRayCaster*>(mRayCaster);
+    else
+      hr = reinterpret_cast<HavranRayCaster*>(mRayCaster);
+
+    string ktbFilename = internKdTree;
+    int l = ktbFilename.length();
+    ktbFilename[l-1] = 't';
+    ktbFilename[l-2] = 'b';
+    ktbFilename[l-3] = 'k';
+		
+    cout << "Trying to load tree from file " << ktbFilename << endl;
+    if (hr->ImportBinTree(ktbFilename, this->mObjects)) {
+      cout << "Loading failed - building kd-tree" << endl;
+      hr->Build(this->mObjects);
+      cout << "Exporting kd-tree to file " << ktbFilename << endl;
+      hr->ExportBinTree(ktbFilename);
+    }
+    else
+      cout << " done." << endl;
   }
   
@@ -1121,8 +1134,14 @@
 	if (rayCastMethod == 2) 
 	{
-		cout << "ray cast method: havran" << endl <<flush;
-		mRayCaster = new HavranRayCaster(*this);
-	}
-
+	  cout << "ray cast method: havran" << endl <<flush;
+	  mRayCaster = new GALIGN16 HavranRayCaster(*this);
+	}
+	if (rayCastMethod == 3) 
+	{
+	  cout << "ray cast method: havran - dyn" << endl <<flush;
+	  mRayCaster = new GALIGN16 HavranDynRayCaster(*this);
+	}
+
+	
 	/////////////////
 	//-- reserve constant block of rays
Index: GTP/trunk/Lib/Vis/Preprocessing/src/default.env
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/default.env	(revision 2616)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/default.env	(revision 2629)
@@ -4,6 +4,9 @@
 #############################################################################
 
-
 Scene {
+
+#filename ./data/arena-high-lods.obj
+filename ./data/vienna/vienna_cropped.obj
+#filename ./data/PowerPlantM.obj
 
 #filename ../data/Arena/arena-low-lods.obj
@@ -20,5 +23,4 @@
 #filename ../data/vienna/vienna.obj
 
-filename ../data/vienna/vienna_cropped.obj
 
 # filename ../data/vienna/viewcells-25-sel.x3d
@@ -41,5 +43,5 @@
 	totalSamples 250000000
 #	totalSamples 2000000
-	samplesPerPass 500000
+	samplesPerPass 100000
 	samplesPerEvaluation 5000000
 
@@ -54,5 +56,5 @@
 	useGlDebugger false
 # 0 = INTERNAL 	1 = MLRT 2 = HAVRAN
-	rayCastMethod 1
+#	rayCastMethod 2
 	
 #	type sampling
@@ -107,4 +109,5 @@
 #	distributions rss+spatial+object_direction
 #	distributions rss+object_direction+spatial+object+direction
+#	distributions mutation+spatial+object_direction+filter_based
 	distributions mutation+spatial+object_direction
 #	distributions rss+object_direction
@@ -159,7 +162,6 @@
 
 Limits {
-
        threshold	1e-6
-       small		1e-6
+       small		1e-5
        infinity		1e9
 }
@@ -183,7 +185,6 @@
 	splitMethod SAH
 	splitBorder 0.01
-	pvsArea 1e-5
-}
-
+	pvsArea 1e-3
+}
 
 
@@ -196,6 +197,6 @@
 	}
 
-	splitMethod spatialMedian
-#	splitMethod SAH
+#	splitMethod spatialMedian
+	splitMethod SAH
 	splitBorder 0.01
 }
@@ -204,5 +205,7 @@
 	importRandomViewCells false
 
-	useKdPvs true
+	# useKdPvs true
+	### VH - problem in ViewCellsManager.cpp - line 690
+	useKdPvs false
 	useKdPvsAfterFiltering true
 	# samples used for view cell construction
@@ -317,5 +320,5 @@
 #	filename ../data/vienna/vsposp-seq-viewCells.xml.gz
 
-	filename ../data/vienna/vienna_cropped-gradient-viewcells.xml.gz
+	filename ./data/vienna/vienna_cropped-gradient-viewcells.xml.gz
 
 #	filename ../data/vienna/vienna_cropped-2-sequential-30000-viewcells.xml.gz
@@ -506,6 +509,5 @@
 	}
 	
-
-
+	
 	Visualization {
 		# x3d visualization of the split planes
@@ -514,13 +516,11 @@
 }
 
-# settings for havran ray caster
-TestDoubleRays	true
+# Settings for havran ray caster
+
 
 BSP {
-  termCrit	  auto
-  # auto2 allows to set the number of primitives per leaf	
+  #termCrit	  auto
   #termCrit	  auto2
-  #adhoc allows to set maxdepth + number of primitives per leaf	
-  #termCrit	  adhoc
+  termCrit	  adhoc
   # For Arena
   #maxDepthAllowed 15
@@ -529,10 +529,38 @@
   #maxDepthAllowed 22
   #maxListLength    6
-  # For PowerPlant
-  maxDepthAllowed 18
-  maxListLength    8
-  # For Power plant
-  #maxDepthAllowed 8
-  #maxListLength    10
-}
-
+  # For Powerplant
+  maxDepthAllowed  27
+  maxListLength     4
+}
+
+Rays {
+  cnt	1000000
+  #cnt	1999000
+  #cnt	499000
+
+  #file	 data/fileRays_arena.txt
+  #file  data/rays-arena-4M-6M.txt
+
+  #file	 data/fileRays_vienna.txt
+  #file  data/fileRays_vienna2.txt
+  file  data/rays-vienna_cropped-4M-6M.txt
+  #file   data/fileRays_vienna_4.txt
+  #file   data/fileRays_vienna_test.txt
+}
+
+#TestDoubleRays	false
+TestDoubleRays	true
+
+Preprocessor {
+  # internal
+  #rayCastMethod 0
+  # Intel
+  #rayCastMethod 1
+  # Havran - 2 - static scene
+  rayCastMethod 2
+  # Havran - 3 - dynamic scene
+  rayCastMethod 3
+
+  loadMeshes false
+}
+
Index: GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp	(revision 2621)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp	(revision 2629)
@@ -26,20 +26,12 @@
 #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
+GALIGN16 RayPacket2x2
 HavranRayCaster::raypack;
 #endif // _USE_HAVRAN_SSE
-#endif // USE_HAVRAN_RAYCASTER
+
   
 HavranRayCaster::HavranRayCaster(const Preprocessor &preprocessor):
@@ -70,16 +62,17 @@
 
   if (!ImportBinTree(kdfile, objlist)) {
-	cout << "\nKd-tree for Havran ray caster imported."<<endl<<flush;
-  } else {
-	CTimer timer;
-	cout << "\nBuilding up kd-tree for Havran ray caster ..."<<endl<<flush;
-	
-	timer.Start();
-	mKtbtree->BuildUp(objlist);
-	timer.Stop();
-	cout << "\nBuilding up kd-tree is finished, user time = "
-		 << timer.UserTime() << " real time = " << timer.RealTime() <<
-	  endl <<flush;
-	ExportBinTree(kdfile);
+        cout << "\nKd-tree for Havran ray caster imported."<<endl<<flush;
+  }
+  else {
+    CTimer timer;
+    cout << "\nBuilding up kd-tree for Havran ray caster ..."<<endl<<flush;
+    
+    timer.Start();
+    mKtbtree->BuildUp(objlist);
+    timer.Stop();
+    cout << "\nBuilding up kd-tree is finished, user time = "
+	 << timer.UserTime() << " real time = " << timer.RealTime()
+	 << endl;
+    ExportBinTree(kdfile);
   }
 #endif
@@ -513,11 +506,10 @@
 #endif // USE_HAVRAN_RAYCASTER 
 }
+
+
+
+  
 #endif // _USE_HAVRAN_SSE
 
 
-    /*
-  virtual void HavranRayCaster::AddDynamicObjecs(const ObjectContainer &objects, const Matrix4x4 &m);
-  virtual void UpdateDynamicObjects(const Matrix4x4 &m);
-  virtual void DeleteDynamicObjects);
-*/
 } // the namespace
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Makefile
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 2622)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Makefile	(revision 2629)
@@ -1,5 +1,5 @@
 #############################################################################
 # Makefile for building: preprocessor
-# Generated by qmake (2.00a) (Qt 4.1.2) on: po 21. I 21:36:51 2008
+# Generated by qmake (2.00a) (Qt 4.1.2) on: út 22. I 23:58:58 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-msvc.net\qmake.conf C:/Qt/4.1.2/mkspecs/qconfig.pri \
+Makefile: preprocessor.pro  C:/Qt/4.1.2/mkspecs/win32-msvc2005\qmake.conf C:/Qt/4.1.2/mkspecs/qconfig.pri \
 		C:\Qt\4.1.2\mkspecs\features\qt_config.prf \
 		C:\Qt\4.1.2\mkspecs\features\exclusive_builds.prf \
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Camera.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Camera.h	(revision 2575)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Camera.h	(revision 2629)
@@ -4,4 +4,5 @@
 #include "Vector3.h"
 #include "AxisAlignedBox3.h"
+#include "ktbconf.h"
 
 namespace GtpVisibilityPreprocessor {
@@ -31,10 +32,9 @@
   }
 
-  Camera(int width, int height) {
+  Camera(int width, int height, float fieldOfView = 90.f) {
     mWidth = width;
     mHeight = height;
-    mFovy = 90.0f*(float)M_PI/180.0f;
+    mFovy = fieldOfView*(float)M_PI/180.0f;
   }
-
   
   void Precompute() {
@@ -87,4 +87,11 @@
 
   bool
+  SnapImage2(std::string filename,
+	     RayCaster *raycaster,
+	     AxisAlignedBox3 &bbox,
+	     SceneGraph *sceneGraph
+	     );
+  
+  bool
   SnapImagePacket(std::string filename,
 		  RayCaster *raycaster,
@@ -92,5 +99,12 @@
 		  SceneGraph *sceneGraph
 		  );
-  
+
+  bool
+  SnapImagePacket2(std::string filename,
+		  RayCaster *raycaster,
+		  AxisAlignedBox3 &bbox,
+		  SceneGraph *sceneGraph
+		  );
+
   void SetupRay(Ray &ray, const int x, const int y);
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp	(revision 2592)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp	(revision 2629)
@@ -5,4 +5,5 @@
 #include "ViewCellsManager.h"
 
+#include <cassert>
 
 namespace GtpVisibilityPreprocessor {
@@ -178,5 +179,4 @@
   //const int batchsize = 16384;
   const int batchsize = 8192;
-  // const int batchsize = 1024;
   //const int batchsize = 128;
 
@@ -238,5 +238,128 @@
 }
   
+
+void
+RayCaster::SortRays2(SimpleRayContainer &rays)
+{
+  AxisAlignedBox3 box =
+	mPreprocessor.mViewCellsManager->GetViewSpaceBox();
+
+  const float sizeBox = Magnitude(box.Diagonal());
+  // This is some size of the 
+  const float sizeDir = 0.2f * sizeBox;
+  
+  float b[12]={
+    box.Min().x,
+    box.Min().y,
+    box.Min().z,
+    -sizeDir,
+    -sizeDir,
+    -sizeDir,
+    box.Max().x,
+    box.Max().y,
+    box.Max().z,
+    sizeDir,
+    sizeDir,
+    sizeDir
+  };
+
+#if 0
+  static vector<SimpleRay *> pointerArray;
+
+  if (pointerArray.size()!=rays.size()) {
+	// realloc the pointerarray
+	pointerArray.resize(rays.size());
+  }
+
+  // init pointer array
+  SimpleRay *p = &pointerArray[0];
+  for (i=0; i < rays.size(); i++, p++)
+	pointerArray[i] = p;
+#endif
+  
+  _SortRays2(rays, 0, (int)rays.size()-1, 0, b);
+
+  return;
+}
 					
+void
+RayCaster::_SortRays2(SimpleRayContainer &rays,
+		      const int l,
+		      const int r,
+		      const int depth,
+		      float box[12])
+{
+  // pick-up a pivot
+  int axis;
+  
+  float maxDiff = -1.0f;
+  // get the largest axis
+  int offset = 0;
+  int i;
+
+  //const int batchsize = 16384;
+  //const int batchsize = 8192;
+  const int batchsize = 128;
+
+  //if (r - l < 16*batchsize)
+  //	offset = 3;
+	
+  //  if (depth%2==0)
+  //	offset = 3;
+  
+  for (i=offset; i < offset + 6; i++) {
+    float diff = box[i + 6] - box[i];
+    assert(diff >= 0.f);
+    if (diff > maxDiff) {
+      // Find the maximum
+      maxDiff = diff;
+      axis = i;
+    }
+  }
+
+  //  cout<<depth<<" "<<axis<<" "<<l<<" "<<r<<endl;
+  
+  i=l;
+  int j=r;
+
+  float x = (box[axis] + box[axis+6])*0.5f;
+  //  float x = rays[(l+r)/2].GetParam(axis);
+  do {
+	while(i<j && rays[i].GetParam(axis) < x)
+	  i++;
+	while(i<j && x < rays[j].GetParam(axis))
+	  j--;
+	
+	if (i <= j) {
+	  swap(rays[i], rays[j]);
+	  i++;
+	  j--;
+	}
+  } while (i<=j);
+
+  
+  if (l + batchsize < j ) {
+	// set new max
+	float save = box[axis+6];
+	box[axis+6] = x;
+	_SortRays2(rays, l, j, depth+1, box);
+	box[axis+6] = save;
+  } else {
+	//	for (int k=0; k < 6; k++)
+	//	  cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
+  }
+  
+  if (i + batchsize < r) {
+    // set new min
+    box[axis] = x;
+    _SortRays2(rays, i, r, depth+1, box);
+  } else {
+	//	for (int k=0; k < 6; k++)
+	//	  cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
+  }
+	
+}
+  
+  
 VssRay *RayCaster::RequestRay(const Vector3 &origin, 
 							  const Vector3 &termination, 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/main.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/main.cpp	(revision 2625)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/main.cpp	(revision 2629)
@@ -26,7 +26,8 @@
 #include "GlobalLinesRenderer.h"
 #include "RayCaster.h"
-//#include "vlastimil/testrt.h"
-
 #include "ViewCellsManager.h"
+
+#include "testrt.h"
+
 
 #ifdef USE_QT  
@@ -134,7 +135,8 @@
 #if 0
   // Test code by VH
-  //TestRTcamera(argc, argv);
+  TestRTcamera(argc, argv);
   //TestRTfromFile(argc, argv);
-  TestRTfromFilePackets(argc, argv);
+  //TestRTfromFilePackets(argc, argv);
+  //TestRT_4_fromFile(argc, argv);
   return 0;
 #endif
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h	(revision 2592)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h	(revision 2629)
@@ -10,5 +10,5 @@
 //
 // Initial coding by Vlastimil Havran, 2006. The data design is in fact
-// Jakko Biker layout as propose in the article on Intel Web Site in year 2005
+// Jakko Biker layout as proposed in the article on Intel Web Site in year 2005
 // http://www.intel.com/cd/ids/developer/asmo-na/eng/245711.htm?page=1
 
@@ -21,4 +21,5 @@
 
 #include "Vector3.h"
+#include "Matrix4x4.h"
 #include "ktbconf.h"
   
@@ -88,6 +89,5 @@
     // and also for shadow rays finish at the same object
     const Intersectable *_stopObject = NULL
-  )
-  {
+  ) {
     // location
     ox[0] = orf[0].x; ox[1] = orf[1].x; ox[2] = orf[2].x; ox[3] = orf[3].x;
@@ -182,4 +182,19 @@
   int GetType() const { return ttype; }
 
+  void ApplyTransform(const Matrix4x4 &tform) {
+    for (int i = 0; i < 4; i++) {
+      Vector3 o_orig(ox[i], oy[i], oz[i]);
+      Vector3 t_orig = tform * o_orig;
+      ox[i] = t_orig.x; oy[i] = t_orig.y; oz[i] = t_orig.z;
+      Vector3 o_dir(dx[i], dy[i], dz[i]);
+      Vector3 t_dir = RotateOnly(tform, o_dir);
+      dx[i] = t_dir.x; dy[i] = t_dir.y; dz[i] = t_dir.z;
+
+      // ?? note that normalization to the unit size of the direction
+      // ?? is NOT computed -- this is what we want.
+    }
+    Precompute();
+  }
+  
   // Reading and Setting origin of the ray and direction
   // Ray origin
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.h	(revision 2629)
@@ -20,4 +20,5 @@
 #include "ktb8b.h"
 #include "Containers.h"
+#include "IntersectableWrapper.h"
 
 namespace GtpVisibilityPreprocessor {
@@ -42,4 +43,22 @@
 #endif
 
+class AxisAlignedBox3Intersectable:
+    public IntersectableWrapper<AxisAlignedBox3>
+{
+public:
+  AxisAlignedBox3Intersectable(const AxisAlignedBox3 &item):
+    IntersectableWrapper<AxisAlignedBox3>(item) { }
+
+  AxisAlignedBox3 GetBox() const { return mItem;}
+
+  int Type() const
+  {
+    // This is not ture, but for our purposes it is OK
+    return Intersectable::TRIANGLE_INTERSECTABLE;
+  }
+
+};
+
+// ---------------------------------------------------------------  
 // The base class for KD-tree with irregular change of axes, where
 // the splitting plane can be positioned.
@@ -216,5 +235,5 @@
   void Check1List(SInputData *data, int axis, int countExpected);
     
-  //----------------------------------------------------------------------
+  //----------------------------------------------------------------
   // Termination criteria and fixing the splitting plane orientation
 
@@ -230,5 +249,5 @@
     CKTBAxes::Axes reqAxis;
 
-    // -------------- AUTOMATIC TERMINATION CRITERIA ---------------------
+    // -------------- AUTOMATIC TERMINATION CRITERIA ----------------
     // the ratio of improvement for the cost by subdivision and not-subdividing
     // for the previous subdivision
@@ -578,4 +597,8 @@
   int startEmptyCutDepth;
 
+  // Biasing the empty cuts (no objects are split). The cost is multiplied
+  // by the coefficient which is assumed to be 0.8-0.9
+  float  biasFreeCuts;
+  
   // ---------- Special improvements on the kd-tree construction --------
   // flag if to split bounding boxes during splitting
@@ -588,7 +611,4 @@
   int     minObjectsToCreateMinBox, minDepthDistanceBetweenMinBoxes;
   float   minSA2ratioMinBoxes;
-  // Biasing the empty cuts (no objects are split). The cost is multiplied
-  // by the coefficient which is assumed to be 0.8-0.9
-  float  biasFreeCuts;
   // Make min box here
   bool   makeMinBoxHere;
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbf2trv.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbf2trv.cpp	(revision 2608)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbf2trv.cpp	(revision 2629)
@@ -57,4 +57,5 @@
       results[i].intersectable;
     SimpleRay::IntersectionRes[i + copyOffset].tdist =
+    SimpleRay::IntersectionRes[i + copyOffset].maxt =
       results[i].tdist;
   } // for i
@@ -204,5 +205,5 @@
     FindNearestI(ray);
     rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-    rp.SetT(i, SimpleRay::IntersectionRes[0].tdist);
+    rp.SetT(i, SimpleRay::IntersectionRes[0].maxt);
     // SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
   } // for
@@ -403,5 +404,5 @@
     FindNearestI(ray);
     rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-    rp.SetT(i, SimpleRay::IntersectionRes[0].tdist);
+    rp.SetT(i, SimpleRay::IntersectionRes[0].maxt);
     // SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
   } // for
@@ -414,6 +415,6 @@
 void
 CKTBTraversal::FindNearestI(RayPacket2x2 &rp,
-							const Vector3 &boxmin,
-							const Vector3 &boxmax)
+			    const Vector3 &boxmin,
+			    const Vector3 &boxmax)
 {
   static AxisAlignedBox3 localbox;
@@ -606,5 +607,5 @@
     FindNearestI(ray, localbox);
     rp.SetObject(i, SimpleRay::IntersectionRes[0].intersectable);
-    rp.SetT(i, SimpleRay::IntersectionRes[0].tdist);
+    rp.SetT(i, SimpleRay::IntersectionRes[0].maxt);
     // SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
   } // for
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.cpp	(revision 2608)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.cpp	(revision 2629)
@@ -11,7 +11,7 @@
 
 // GOLEM library
-#include "ktbconf.h"
 #include "ktb.h"
 #include "ktbai.h"
+#include "ktbs.h"
 #include "ktball.h"
 #include "ktbtrav.h"
@@ -61,14 +61,16 @@
 {
   CKTBAllocManPredecessor *bc = 0;
-  bool useArray = true;
+  bool useArray = false;
   
   // this should be already initialised, but we want to be sure
   if (useArray) {
+    cout << "Using kbt-tree based on arrays" << endl;
     // the implementation based on arrays
     bc = new CKTBABuildUp(); // file ktbai.cpp
   }
   else {
-    // the implementation based on lists
-    // bc = new CKTBBuildUp(); // file ktbi.cpp
+    cout << "Using kbt-tree based on sampling" << endl;
+    // the implementation based on sampling
+    bc = new CKTBSBuildUp(); // file ktbs.cpp
   }
     
@@ -425,6 +427,6 @@
 void
 CKTB::FindNearestI(RayPacket2x2 &raypack,
-				   const Vector3 &boxmin,
-				   const Vector3 &boxmax)
+		   const Vector3 &boxmin,
+		   const Vector3 &boxmax)
 {
   if (!makeMinBoxes)
@@ -766,5 +768,20 @@
 }
 
-  
+
+static int
+cntInterior = 0;
+
+static int
+cntInteriorA[4] = {0, 0, 0, 0};
+
+static int
+cntLeaves = 0;
+
+static int
+cntEmptyLeaves = 0;
+
+static int
+objReferences = 0;
+
 SKTBNodeT*
 CKTB::ImportBinLeaf(IN_STREAM &stream, 
@@ -775,4 +792,6 @@
   int objId = leafId;
   int size;
+
+  cntLeaves++;
   
   stream.read(reinterpret_cast<char *>(&size), sizeof(int));
@@ -782,4 +801,5 @@
     SKTBNodeT* leaf = buildClass->AllocLeaf(0);
     *nodeToLink = buildClass->nodeToLink;
+    cntEmptyLeaves++;
     return leaf;
   }
@@ -790,4 +810,6 @@
 
   ObjectContainer *newobjlist = new ObjectContainer;
+
+  objReferences += size;
   
   // read object ids
@@ -826,5 +848,4 @@
 }
 
-
 SKTBNodeT *
 CKTB::ImportBinInterior(IN_STREAM  &stream, SKTBNodeT **nodeToLink)
@@ -839,5 +860,9 @@
     buildClass->AllocInteriorNode(axis, pos, 0, 0);
   *nodeToLink = buildClass->nodeToLink;
-	
+
+  cntInterior++;
+  assert( (axis >= 0) && (axis < 3));
+  cntInteriorA[axis]++;
+  
   return interiorNode;
 }
@@ -928,5 +953,5 @@
 
   if (!stream.is_open()) {
-    cerr << "Kd-tree description file (.kdh) cannot be opened for reading\n";
+    cerr << "Kd-tree description file (.kbt) cannot be opened for reading\n";
     return true; // error
   }
@@ -963,4 +988,7 @@
   bc->InitAux(0, CKTBNodeAbstract::MAX_HEIGHT - 1, maxItemsAtOnce);  
 
+  for(int i = 0; i < 4; i++)
+    cntInteriorA[i] = 0;
+  
   // Compute the box from all objects
   bbox.Initialize();
@@ -1040,4 +1068,12 @@
   builtUp = true;  
 
+  cout << "Importing ktb-tree cntInterior = "
+       << cntInterior << " cntLeaves = " << cntLeaves
+       << " cntEmptyLeaves = " << cntEmptyLeaves << endl;
+
+  cout << "X-splits = " << cntInteriorA[0] << " Y-splits = " << cntInteriorA[1] 
+       << " Z-splits = " << cntInteriorA[2] << "\nref_to_objs/full_leaf = "
+       << (float)objReferences/(cntLeaves - cntEmptyLeaves) << endl;
+  
   return false; // OK
 }
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.h	(revision 2608)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.h	(revision 2629)
@@ -16,5 +16,4 @@
 // GOLEM headers
 #include "configh.h"
-#include "ktbconf.h"
 //#include "ktbi.h"
 #include "ktbai.h"
@@ -155,5 +154,5 @@
   
   // allocate new build class and return it
-  static CKTBAllocManPredecessor* AllocBuildClass();
+  CKTBAllocManPredecessor* AllocBuildClass();
   
   CKTBNodeIteratorPredecessor* GetTraversalClass() {
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/configh.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/configh.h	(revision 2610)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/configh.h	(revision 2629)
@@ -13,13 +13,4 @@
 #ifndef __CONFIGH_H__
 #define __CONFIGH_H__
-
-
-#if defined(_MSC_VER)
-// use perftimer only on msvc
-// define __SSE__ macro as it is not defined under MSVC
-#define __SSE__
-// If we support the use of SSE instructions for ray shooting
-//#define _USE_HAVRAN_SSE
-#endif
 
 #ifndef USE_GOLEM_NAMESPACE
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbtrav.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbtrav.h	(revision 2602)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbtrav.h	(revision 2629)
@@ -150,5 +150,6 @@
   // of a packet to individual rays and traced individually.
   virtual void FindNearestI(RayPacket2x2 &raypack) { }
-  virtual void FindNearestI(RayPacket2x2 &raypack, const Vector3 &boxMin, const Vector3 &boxMax) { }
+  virtual void FindNearestI(RayPacket2x2 &raypack, const Vector3 &boxMin,
+			    const Vector3 &boxMax) { }
 #endif
   
@@ -328,5 +329,6 @@
   // of a packet to individual rays and traced individually.
   virtual void FindNearestI(RayPacket2x2 &raypack);
-  virtual void FindNearestI(RayPacket2x2 &raypack, const Vector3 &boxMin, const Vector3 &boxMax);
+  virtual void FindNearestI(RayPacket2x2 &raypack, const Vector3 &boxMin,
+			    const Vector3 &boxMax);
 #endif // __SSE__
 #endif // _USE_HAVRAN_SSE  
@@ -421,5 +423,6 @@
   // of a packet to individual rays and traced individually.
   virtual void FindNearestI(RayPacket2x2 &raypack);
-  virtual void FindNearestI(RayPacket2x2 &raypack, const Vector3 &boxMin, const Vector3 &boxMax);
+  virtual void FindNearestI(RayPacket2x2 &raypack, const Vector3 &boxMin,
+			    const Vector3 &boxMax);
 #endif // __SSE__
 #endif // _USE_HAVRAN_SSE  
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/sbbox.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/sbbox.h	(revision 2582)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/sbbox.h	(revision 2629)
@@ -146,5 +146,8 @@
 
   inline bool Equal(const SBBox &b, float eps = 0.f) const;
-    
+
+  // Returns the intersection of two axis-aligned boxes.
+  friend inline SBBox Intersect(const SBBox &x, const SBBox &y);
+  
   // Test if the box is really sensefull
   bool IsCorrect();
@@ -278,4 +281,19 @@
 DescribeXYZ(const SBBox &b, std::ostream &app, int indent);
 
+inline SBBox
+Intersect(const SBBox &x, const SBBox &y)
+{
+  SBBox ret = x;
+  if (OverlapS(ret, y)) {
+    Maximize(ret.pp[0], y.pp[0]);
+    Minimize(ret.pp[1], y.pp[1]);
+    return ret;
+  }
+
+  // Null intersection.
+  return SBBox(Vector3(0), Vector3(0));
+}
+
+
 #if 1
 // The implementation I, the first version implemented by Vlastimil
@@ -425,4 +443,5 @@
 #endif
 
+
 }
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbftrav.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbftrav.cpp	(revision 2602)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbftrav.cpp	(revision 2629)
@@ -166,5 +166,6 @@
 	    cout << "Full leaf HIT " << endl; 
 #endif 	
-	      
+	    SimpleRay::IntersectionRes[0].tdist =
+	      SimpleRay::IntersectionRes[0].maxt;
 #ifdef __TRAVERSAL_STATISTICS
 	    _allNodesTraversed += allNodesTraversed;
@@ -372,4 +373,6 @@
 	    cout << "Full leaf HIT " << endl; 
 #endif 	
+	    SimpleRay::IntersectionRes[0].tdist =
+	      SimpleRay::IntersectionRes[0].maxt;
 	      
 #ifdef __TRAVERSAL_STATISTICS
@@ -381,4 +384,5 @@
 	    // We have to add the distance from the original ray origin
 	    SimpleRay::IntersectionRes[0].tdist += tminOffset;
+	    SimpleRay::IntersectionRes[0].maxt = SimpleRay::IntersectionRes[0].tdist;
 	    
 	    // signed distance should be already set in TestFullLeaf
@@ -592,4 +596,7 @@
 	      tmax + Limits::Small;
 	    if (TestFullLeaf(rays[indexR+offset], currNode, indexR)) {
+	      // copy the result to tdist
+	      SimpleRay::IntersectionRes[indexR + rayOffset].tdist =
+		SimpleRay::IntersectionRes[indexR + rayOffset].maxt;
 
 	      // we remove the ray from the calculation
@@ -599,5 +606,5 @@
 #ifdef _DEBUGKTB
 	      cout << "Full leaf HIT " << endl; 
-#endif 	
+#endif
 	      
 #ifdef __TRAVERSAL_STATISTICS
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/testrt.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/testrt.cpp	(revision 2602)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/testrt.cpp	(revision 2629)
@@ -26,4 +26,6 @@
 #include "GlobalLinesRenderer.h"
 #include "RayCaster.h"
+#include "Triangle3.h"
+#include "IntersectableWrapper.h"
 #include "timer.h"
 #include "raypack.h"
@@ -64,34 +66,109 @@
 extern void Cleanup();
 
-static string ReplaceSuffix(const string &filename, const string &a, const string &b)
+void
+_SortRays2(SimpleRayContainer &rays,
+	   const int l,
+	   const int r,
+	   const int depth,
+	   float box[12])
 {
-	string result = filename;
-
-	int pos = (int)filename.rfind(a, (int)filename.size() - 1);
-	if (pos == filename.size() - a.size()) 
-	{
-		result.replace(pos, a.size(), b);
+  // pick-up a pivot
+  int axis;
+  
+  float maxDiff = -1.0f;
+  // get the largest axis
+  int offset = 0;
+  int i;
+
+  //const int batchsize = 16384;
+  //const int batchsize = 8192;
+  //const int batchsize = 128;
+  const int batchsize = 4;
+
+  //if (r - l < 16*batchsize)
+  //	offset = 3;
+	
+  //  if (depth%2==0)
+  //	offset = 3;
+  
+  for (i=offset; i < offset + 6; i++) {
+    float diff = box[i + 6] - box[i];
+    assert(diff >= 0.f);
+    if (diff > maxDiff) {
+      // Find the maximum
+      maxDiff = diff;
+      axis = i;
+    }
+  }
+
+  //  cout<<depth<<" "<<axis<<" "<<l<<" "<<r<<endl;
+  
+  i=l;
+  int j=r;
+
+  float x = (box[axis] + box[axis+6])*0.5f;
+  //  float x = rays[(l+r)/2].GetParam(axis);
+  do {
+	while(i<j && rays[i].GetParam(axis) < x)
+	  i++;
+	while(i<j && x < rays[j].GetParam(axis))
+	  j--;
+	
+	if (i <= j) {
+	  swap(rays[i], rays[j]);
+	  i++;
+	  j--;
 	}
-
-	return result;
+  } while (i<=j);
+
+  
+  if (l + batchsize < j ) {
+	// set new max
+	float save = box[axis+6];
+	box[axis+6] = x;
+	_SortRays2(rays, l, j, depth+1, box);
+	box[axis+6] = save;
+  } else {
+	//	for (int k=0; k < 6; k++)
+	//	  cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
+  }
+  
+  if (i + batchsize < r) {
+    // set new min
+    box[axis] = x;
+    _SortRays2(rays, i, r, depth+1, box);
+  } else {
+	//	for (int k=0; k < 6; k++)
+	//	  cout<<k<<" "<<box[k]<<" - "<<box[k+6]<<endl;
+  }	
 }
 
-
-static int SplitFilenames(const string &str, vector<string> &filenames)
+void
+SortRays2(SimpleRayContainer &rays, const AxisAlignedBox3 &box)
 {
-	int pos = 0;
-
-	while(1) {
-		int npos = (int)str.find(';', pos);
-		
-		if (npos < 0 || npos - pos < 1)
-			break;
-		filenames.push_back(string(str, pos, npos - pos));
-		pos = npos + 1;
-	}
-	
-	filenames.push_back(string(str, pos, str.size() - pos));
-	return (int)filenames.size();
-}
+  const float sizeBox = Magnitude(box.Diagonal());
+  // This is some size of the 
+  const float sizeDir = 0.1f * sizeBox;
+  
+  float b[12]={
+    box.Min().x,
+    box.Min().y,
+    box.Min().z,
+    -sizeDir,
+    -sizeDir,
+    -sizeDir,
+    box.Max().x,
+    box.Max().y,
+    box.Max().z,
+    sizeDir,
+    sizeDir,
+    sizeDir
+  };
+
+  
+  _SortRays2(rays, 0, (int)rays.size()-1, 0, b);
+
+  return;
+}					
 
 
@@ -254,11 +331,11 @@
   AxisAlignedBox3 bboxOrig = preprocessor->mSceneGraph->GetBox();
   AxisAlignedBox3 bbox = bboxOrig;
-  int sizeDiag = Magnitude(bbox.Diagonal());
+  float sizeDiag = Magnitude(bbox.Diagonal());
   bbox.Enlarge(sizeDiag * 1.5f);
 
   Vector3 origin, dir;
 
-//#define WIEN1
-#define WIEN2
+#define WIEN1
+//#define WIEN2
 //#define ARENA1
     
@@ -267,4 +344,5 @@
   origin = Vector3(1099.9f, 183.0f, -387.0f);
   dir = Vector3(-0.6f, 0.0001f, -0.8f);
+  //dir = -dir;
 #define DIREXISTS
 #endif  
@@ -292,25 +370,120 @@
   cout << "Computing image\n" << endl;
 
+  int id = preprocessor->mObjects.size() + 1;
+    
+  //We add here a few objects
+  ObjectContainer dynObjects;
+  Vector3 baseVec = origin + dir * 2.0f;
+  Triangle3 tr1(baseVec, baseVec + Vector3(1, 0, 0),
+		baseVec + Vector3(0, 1, 0));	       
+  TriangleIntersectable ti1(tr1);
+  dynObjects.push_back(&ti1);
+  ti1.mId = id+1;
+  
+  tr1.Init(baseVec, baseVec + Vector3(1, 0, 0),
+	   baseVec + Vector3(0, 0, 1));	       
+  TriangleIntersectable ti2(tr1);
+  dynObjects.push_back(&ti2);
+  ti2.mId = id+2;
+
+  tr1.Init(baseVec, baseVec + Vector3(0, 1, 0),
+	   baseVec + Vector3(0, 0, 1));	       
+  TriangleIntersectable ti3(tr1);
+  dynObjects.push_back(&ti3);
+  ti3.mId = id+3;
+
 #if 1
-  // ray by ray
-  cam.SnapImage("test-rays.tga",
-		preprocessor->mRayCaster,
-		bboxOrig,
-		preprocessor->mSceneGraph);
-#endif
+  // This is required for cam.SnapImagePacket2
+  // otherwise dynamic object cannot be identified!
+  preprocessor->mObjects.push_back(&ti1);
+  preprocessor->mObjects.push_back(&ti2);
+  preprocessor->mObjects.push_back(&ti3);
+#endif
+
+  if (0) {
+    Matrix4x4 mat;
+    mat = IdentityMatrix();
+    mat = mat * TranslationMatrix(Vector3(-1, 0, 0));
+    preprocessor->mRayCaster->AddDynamicObjecs(dynObjects, mat);
+  
 #if 1
-  // using ray packets
-  cam.SnapImagePacket("test-packet.tga",
-		      preprocessor->mRayCaster,
-		      bboxOrig,
-		      preprocessor->mSceneGraph);
+    // ray by ray
+    cam.SnapImage("test-rays.tga",
+		  preprocessor->mRayCaster,
+		  bboxOrig,
+		  preprocessor->mSceneGraph);
 #endif
 #if 1
-  // using ray packets
-  cam.SnapImage2("test-oneDir.tga",
-		 preprocessor->mRayCaster,
-		 bboxOrig,
-		 preprocessor->mSceneGraph);
-#endif
+    // using ray packets
+    cam.SnapImage2("test-oneDir.tga",
+		   preprocessor->mRayCaster,
+		   bboxOrig,
+		   preprocessor->mSceneGraph);
+#endif
+    
+#if 1
+#ifdef _USE_HAVRAN_SSE
+    // using ray packets
+    cam.SnapImagePacket("test-packet.tga",
+			preprocessor->mRayCaster,
+			bboxOrig,
+			preprocessor->mSceneGraph);
+#endif
+#endif    
+#if 1
+    // using ray packets
+    cam.SnapImagePacket2("test-packet4.tga",
+			 preprocessor->mRayCaster,
+			 bboxOrig,
+			 preprocessor->mSceneGraph);
+#endif
+  }
+  else {
+    cout << "Computing animation" << endl;    
+    for (int i = 0; i < 20; i++) {
+      preprocessor->mRayCaster->DeleteDynamicObjects();
+      Matrix4x4 mat;
+      mat = IdentityMatrix();
+      mat = mat * TranslationMatrix(Vector3(-0.1*(float)i, 0, 0));
+      preprocessor->mRayCaster->AddDynamicObjecs(dynObjects, mat);
+      char name[200];      
+#if 1
+    // ray by ray
+      sprintf(name, "test-rays-%03d.tga", i);
+      cam.SnapImage(name,
+		    preprocessor->mRayCaster,
+		    bboxOrig,
+		    preprocessor->mSceneGraph);
+#endif
+#if 1
+      // using ray packets
+      sprintf(name, "test-oneDir-%03d.tga", i);
+      cam.SnapImage2(name,
+		     preprocessor->mRayCaster,
+		     bboxOrig,
+		     preprocessor->mSceneGraph);
+#endif
+
+#if 1
+#ifdef _USE_HAVRAN_SSE  
+      // using ray packets
+      sprintf(name, "test-packet-%03d.tga", i);
+      cam.SnapImagePacket(name,
+			  preprocessor->mRayCaster,
+			  bboxOrig,
+			  preprocessor->mSceneGraph);
+#endif
+#endif
+#if 1
+      // using ray packets
+      sprintf(name, "test-packet4-%03d.tga", i);
+      cam.SnapImagePacket2(name,
+			   preprocessor->mRayCaster,
+			   bboxOrig,
+			   preprocessor->mSceneGraph);
+#endif
+    } // for
+  } // animation
+
   
   cout << "Done\n" << endl;
@@ -327,4 +500,5 @@
 };
 
+
 // This is for testing RT implementation
 void
@@ -427,5 +601,5 @@
   int cntMaxRays = 100000;
   Environment::GetSingleton()->GetIntValue("Rays.cnt", cntMaxRays);
-  vector<SimpleRay> rays;
+  SimpleRayContainer rays;
   SimpleRay rayTest;
   vector<RESult> results;
@@ -475,10 +649,34 @@
   if (castDoubleRays)
     mult = 2.0;
-  
+
+  cout << "Press a key to start ray shooting" << endl;
+  getchar();
+  cout << "Ray shooting " << cntRays << " rays started - "
+       << (castDoubleRays ? " double " : " single ")
+       << " dir " << endl;
+  
+  AxisAlignedBox3 bboxOrig = preprocessor->mSceneGraph->GetBox();
+  
+  cout << "Sorting rays " << endl;
+
+  CTimer timer;
+  timer.Reset();
+  timer.Start();
+  long t1 = GetTime();
+
+  //SortRays2(rays, bboxOrig);
+
+  timer.Stop();
+  long t2 = GetTime();
+  cout<<"\n#SORTING_TIME = ";
+  cout << TimeDiff(t1, t2)<<" [mikrosec]"
+       << " userTime = " << timer.UserTime() << " realTime = "
+       << timer.RealTime() << endl;
+
   cout << "Starting to shoot " << cntRays * mult << " rays" << endl;
+  timer.Start();
   
   RayCaster *raycaster = preprocessor->mRayCaster;
   VssRayContainer vssRays;
-  AxisAlignedBox3 bboxOrig = preprocessor->mSceneGraph->GetBox();
 
 //#define DEBUGRESULTS
@@ -494,17 +692,7 @@
 #endif
 
-  cout << "Press a key to start ray shooting" << endl;
-  getchar();
-  cout << "Ray shooting " << cntRays << " rays started - "
-       << (castDoubleRays ? " double " : " single ")
-       << " dir " << endl;
-
-  long t1 = GetTime();
-  CTimer timer;
-  timer.Reset();
-  timer.Start();
-  
   SimpleRayContainer raysToTest;
   for (int i = 0; i < cntRays - 16; i++) {
+
 #if 0
     int res = raycaster->CastRay(rays[i],
@@ -573,6 +761,6 @@
 
   timer.Stop();
-  long t2 = GetTime();
-  cout<<"\n#RAY_CAST_TIME = ";
+  t2 = GetTime();
+  cout<<"\n#SORTING + RAY_CAST_TIME = ";
   cout << TimeDiff(t1, t2)<<" [mikrosec]"
        << " userTime = " << timer.UserTime() << " realTime = "
@@ -695,5 +883,5 @@
   int cntMaxRays = 100000;
   Environment::GetSingleton()->GetIntValue("Rays.cnt", cntMaxRays);
-  vector<SimpleRay> rays;
+  SimpleRayContainer rays;
   SimpleRay rayTest;
   vector<RESult> results;
@@ -957,5 +1145,5 @@
   int cntMaxRays = 100000;
   Environment::GetSingleton()->GetIntValue("Rays.cnt", cntMaxRays);
-  vector<SimpleRay> rays;
+  SimpleRayContainer rays;
   SimpleRay rayTest;
   vector<RESult> results;
@@ -1038,4 +1226,6 @@
   timer.Start();
   Vector3 boxMin, boxMax;
+
+  bool printOut = false;
   
   SimpleRayContainer raysToTest;
@@ -1056,18 +1246,48 @@
 			       origin4, direction4,
 			       result4, dist4);
-    printf("I %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f\n",
-	   boxMin.x, boxMin.y, boxMin.z, boxMax.x, boxMax.y, boxMax.z);
+    if (printOut) {
+      printf("I %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f\n",
+	     boxMin.x, boxMin.y, boxMin.z, boxMax.x, boxMax.y, boxMax.z);
+    }
 
     for (int j = 0; j < 4; j++) {
-      printf("%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
-	     i+j,
-	     origin4[j].x,
-	     origin4[j].y,
-	     origin4[j].z,
-	     direction4[j].x,
-	     direction4[j].y,
-	     direction4[j].z,
-	     (result4[j] != -1) ? 1 : 0,
-	     (result4[j] != -1) ? dist4[j] : 0);      
+#if 0
+      if (result4[j] == 0) {
+	int res = raycaster->CastRay(rays[i+j],
+				     vssRays,
+				     bboxOrig,
+				     false,
+				     true);
+	if (res) {
+	  float tdist = SimpleRay::IntersectionRes[0].tdist;
+	  Vector3 point = rays[i+j].Extrap(tdist);
+	  AxisAlignedBox3 testbox(boxMin, boxMax);
+	  
+	  if (testbox.IsInside(point)) {
+	    cout << "Error in the algorithm - computed not in the box, but"
+		 << " it is later found in the box" << endl;
+	    cout << " j = " << j << endl;
+	    cout << " box = " << testbox << endl;
+	    cout << " point = " << point << endl;
+	    raycaster->CastRaysPacket4(boxMin, boxMax,
+			       origin4, direction4,
+			       result4, dist4);
+	  }
+	}
+      }
+#endif      
+      
+      if (printOut) {
+	printf("%d %4.7f %4.7f %4.7f %4.7f %4.7f %4.7f %d %4.7f\n",
+	       i+j,
+	       origin4[j].x,
+	       origin4[j].y,
+	       origin4[j].z,
+	       direction4[j].x,
+	       direction4[j].y,
+	       direction4[j].z,
+	       (result4[j] != -1) ? 1 : 0,
+	       (result4[j] != -1) ? dist4[j] : 0);
+      }
     } // for j    
   } // for i
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.cpp	(revision 2592)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.cpp	(revision 2629)
@@ -218,6 +218,6 @@
 
   // Maximum depth of the tree is set and stack is allocated
-  maxTreeDepth = 50;
-  stackDepth = 2*maxTreeDepth;
+  maxTreeDepth = MAX_HEIGHT;
+  stackDepth = maxTreeDepth + 2;
   stackID = new GALIGN16 SInputData[stackDepth];
   assert(stackID);
@@ -1149,5 +1149,5 @@
       int estHeight = (int)(log((float)initcnt)/log((float)2.0) + 0.9);
       // cout << "EstHeight=" << estHeight << endl;
-      maxDepthAllowed = (int)((float)estHeight * 1.2 + 2.0);
+      maxDepthAllowed = (int)((float)estHeight * 1.2f + 2.0f);
 
       // maximum number of trials to further subdivide
@@ -1205,4 +1205,8 @@
     }
   }
+  else {
+    if (maxDepthAllowed >= MAX_HEIGHT)
+      maxDepthAllowed = MAX_HEIGHT - 1;
+  }
 
   if (verbose)
@@ -1221,4 +1225,33 @@
   if (objlist.size() == 0)
     return 0; // nothing
+
+  // ---------------------------------------------------
+  // Rewriting the triangles to the just wrappers to save
+  // the memory during building!!!!
+  bool useWrappers = false;
+  if (useWrappers) {
+    cout << "WARNING: using only wrappers, not objects to save the memory!"
+	 << endl;
+    cout << "Size of(AxisAlignedBox3Intersectable) = " << sizeof(AxisAlignedBox3Intersectable) << endl;
+    cout << "Size of(TriangleIntersectable) = " << sizeof(TriangleIntersectable) << endl;
+    for (ObjectContainer::iterator it = objlist.begin();
+	 it != objlist.end(); it++) {
+      // take the triangle
+      Intersectable *i = *it;
+      // store the properties to new variables
+      AxisAlignedBox3 b = i->GetBox();
+      int mId = i->mId;
+      // delete the triangle
+      delete i;
+      // create the wrapper with the same box as triangle
+      AxisAlignedBox3Intersectable *a = new AxisAlignedBox3Intersectable(b);
+      a->mId = mId;
+      // and put it back to the list of objects
+      *it = a;
+    } // for    
+    cout << "Rewriting to wrappers is finished!" << endl;
+  } // if
+  // ---------------------------------------------------
+  
   // cerr<<"hh44"<<endl;
   // initialize the whole box of the kd-tree and sort
Index: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbconf.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbconf.h	(revision 2608)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbconf.h	(revision 2629)
@@ -23,4 +23,6 @@
 #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/SG08/run_test_arena
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/SG08/run_test_arena	(revision 2622)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/SG08/run_test_arena	(revision 2629)
@@ -4,14 +4,15 @@
 
 # NORMAL TEST
-#COMMAND="../scripts/preprocessor.sh -preprocessor_quit_on_finish+ -preprocessor_use_gl_renderer- -preprocessor_evaluate_filter- -preprocessor_detect_empty_viewspace+ -total_samples=1000000000  -samples_per_pass=1000000 -total_samples=1000000000 -view_cells_use_kd_pvs+ -af_use_kd_pvs+ -kd_pvs_area=1e-5"
+COMMAND="../scripts/preprocessor.sh -preprocessor_quit_on_finish+ -preprocessor_use_gl_renderer- -preprocessor_evaluate_filter- -preprocessor_detect_empty_viewspace+ -total_samples=1000000000  -samples_per_pass=1000000 -total_samples=1000000000 -view_cells_use_kd_pvs+ -af_use_kd_pvs+ -kd_pvs_area=1e-5"
 
 #TIME TERMINATION
-COMMAND="../scripts/preprocessor.sh -preprocessor_quit_on_finish+ -preprocessor_use_gl_renderer- -preprocessor_evaluate_filter- -preprocessor_detect_empty_viewspace- -total_samples=1000000000  -samples_per_pass=1000000 -total_samples=4000000000 -total_time=7200 -view_cells_use_kd_pvs+ -af_use_kd_pvs+ -kd_pvs_area=1e-5"
+#COMMAND="../scripts/preprocessor.sh -preprocessor_quit_on_finish+ -preprocessor_use_gl_renderer- -preprocessor_evaluate_filter- -preprocessor_detect_empty_viewspace- -total_samples=1000000000  -samples_per_pass=1000000 -total_samples=4000000000 -total_time=7200 -view_cells_use_kd_pvs+ -af_use_kd_pvs+ -kd_pvs_area=1e-5"
 
 
 SCENE=../data/Arena/arena-high-lods.obj
-VIEWCELLS=../data/Arena/viewcells-5000.xml.gz
+#VIEWCELLS=../data/Arena/viewcells-5000.xml.gz
+VIEWCELLS=../data/Arena/arena-high-lods-5000-viewcells.xml.gz
 
-PREFIX=../work/plots/osp-arena-SG08x-TIME
+PREFIX=../work/plots/osp-arena-SG08k
 
 # $COMMAND -preprocessor=combined -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
@@ -24,14 +25,14 @@
 
 
-$COMMAND -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
- -rss_distributions=direction -view_cells_filter_max_size=1 \
- -preprocessor_visibility_file=$PREFIX-r-reference-global.xml \
- -preprocessor_stats=$PREFIX-r-reference-global.log \
- -preprocessor_histogram_file=$PREFIX-r-reference-global.hlog
+# $COMMAND -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
+#  -rss_distributions=direction -view_cells_filter_max_size=1 \
+#  -preprocessor_visibility_file=$PREFIX-r-reference-global.xml \
+#  -preprocessor_stats=$PREFIX-r-reference-global.log \
+#  -preprocessor_histogram_file=$PREFIX-r-reference-global.hlog
 
-$COMMAND -preprocessor=sampling -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
--preprocessor_visibility_file=$PREFIX-r-reference.xml \
--view_cells_filter_max_size=1 -preprocessor_stats=$PREFIX-r-reference.log \
--preprocessor_histogram_file=$PREFIX-r-reference.hlog
+# $COMMAND -preprocessor=sampling -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
+# -preprocessor_visibility_file=$PREFIX-r-reference.xml \
+# -view_cells_filter_max_size=1 -preprocessor_stats=$PREFIX-r-reference.log \
+# -preprocessor_histogram_file=$PREFIX-r-reference.hlog
 
 # -total_samples=500000000000 -total_time=12208 \
Index: GTP/trunk/Lib/Vis/Preprocessing/src/SG08/run_test_powerplant
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/SG08/run_test_powerplant	(revision 2608)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/SG08/run_test_powerplant	(revision 2629)
@@ -2,5 +2,5 @@
 
 #COMMAND="./release/preprocessor.exe -preprocessor_quit_on_finish+"
-COMMAND="../scripts/preprocessor.sh -preprocessor_quit_on_finish+ -preprocessor_use_gl_renderer- -preprocessor_evaluate_filter- -preprocessor_detect_empty_viewspace+ -total_samples=500000000 -view_cells_use_kd_pvs+ -af_use_kd_pvs+"
+COMMAND="../scripts/preprocessor.sh -preprocessor_quit_on_finish+ -preprocessor_use_gl_renderer- -preprocessor_evaluate_filter- -preprocessor_detect_empty_viewspace+ -samples_per_evaluation=5000000 -samples_per_pass=1000000 -total_samples=250000000 -view_cells_use_kd_pvs+ -af_use_kd_pvs+ -preprocessor_ray_cast_method=2 -kd_pvs_area=1e-5"
 
 
@@ -8,15 +8,22 @@
 VIEWCELLS=../data/PowerPlant/powerplant-seq-viewcells.xml.gz
 
-PREFIX=../work/plots/osp-powerplant-SG08
+PREFIX=../work/plots/osp-powerplant-SG08c
 
 
-# $COMMAND -preprocessor=combined -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
-#  -rss_distributions=mutation+object_direction+spatial -view_cells_filter_max_size=1 \
-#  -view_cells_use_kd_pvs+ -af_use_kd_pvs+ \
-#  -preprocessor_visibility_file=$PREFIX-i-mixed-b1-n4a.xml \
-#  -preprocessor_stats=$PREFIX-i-mixed-b1-n4a.log \
-#  -preprocessor_histogram_file=$PREFIX-i-mixed-b1-n4a.hlog
+# n - no origin mutation, q=2, reverse samples
 
-# n - no origin mutation, q=2, reverse samples
+$COMMAND -preprocessor=combined -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
+ -rss_distributions=mutation+object_direction+spatial -view_cells_filter_max_size=1 \
+ -mutation_silhouette_prob=0.5 \
+ -mutation_reverse_samples_distance=1.0 \
+ -mutation_radius_origin=0.5 \
+ -mutation_radius_termination=0.2 \
+ -mutation_use_unsucc_count_importance- \
+ -mutation_use_pass_importance+ \
+ -mutation_buffer_size=2000000 \
+ -preprocessor_visibility_file=$PREFIX-i-mixed-b1-n4a.xml \
+ -preprocessor_stats=$PREFIX-i-mixed-b1-n4a.log \
+ -preprocessor_histogram_file=$PREFIX-i-mixed-b1-n4a.hlog
+
 
 $COMMAND -preprocessor=sampling -scene_filename=$SCENE -view_cells_filename=$VIEWCELLS \
@@ -31,3 +38,5 @@
  -preprocessor_histogram_file=$PREFIX-r-reference-global.hlog
 
+
+
 #mutation+object_direction+
Index: GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.h
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.h	(revision 2623)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.h	(revision 2629)
@@ -22,5 +22,8 @@
 struct SimpleRay;
 class CKTB;
-  
+
+// This macro should be undefined when testing ray tracing
+// by casting rays from file or using camera
+#define _PROCESS_RAY
 
 /** This class provides an interface for ray casting.
@@ -37,9 +40,5 @@
   
   int Type() const { return HAVRAN_RAYCASTER; }
-    
-  virtual void AddDynamicObjecs(const ObjectContainer &objects, const Matrix4x4 &m){};
-  virtual void UpdateDynamicObjects(const Matrix4x4 &m){};
-  virtual void DeleteDynamicObjects(){};
-
+  
   virtual int CastRay(
 		      const SimpleRay &simpleRay,
@@ -98,9 +97,114 @@
 protected:
   CKTB *mKtbtree;
-  static RayPacket2x2 raypack;
+#ifdef _USE_HAVRAN_SSE
+  static GALIGN16 RayPacket2x2 raypack;
+#endif
 };
 
 
+// --------------------------------------------------------------------
+// The implementation of ray caster with dynamic objects
+
+class HavranDynRayCaster: public HavranRayCaster
+{
+public:
+  /** Default constructor initialising e.g., KD tree
+   */
+  HavranDynRayCaster(const Preprocessor &preprocessor);
+  virtual ~HavranDynRayCaster();
+
+  int Type() const { return HAVRAN_DYN_RAYCASTER; }
+  
+  virtual int CastRay(const SimpleRay &simpleRay,
+		      VssRayContainer &vssRays,
+		      const AxisAlignedBox3 &box,
+		      const bool castDoubleRay,
+		      const bool pruneInvalidRays = true
+		      );
+  
+  virtual void CastRays16(SimpleRayContainer &rays, 
+			  VssRayContainer &vssRays,
+			  const AxisAlignedBox3 &sbox,
+			  const bool castDoubleRay,
+			  const bool pruneInvalidRays = true
+			  );
+
+  virtual void CastRays16(SimpleRayContainer &rays,
+			  int offset,
+			  VssRayContainer &vssRays,
+			  const AxisAlignedBox3 &sbox,
+			  const bool castDoubleRay,
+			  const bool pruneInvalidRays = true
+			  );
+  virtual void
+  CastSimpleForwardRays(SimpleRayContainer &rays,
+			const AxisAlignedBox3 &sbox);
+
+  virtual void CastRays(SimpleRayContainer &rays,
+			VssRayContainer &vssRays,
+			const AxisAlignedBox3 &sbox,
+			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(const Vector3 &minBox,
+			       const Vector3 &maxBox,
+			       const Vector3 origin4[],
+			       const 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  
+
+  virtual void AddDynamicObjecs(const ObjectContainer &objects, const Matrix4x4 &m);
+
+  virtual void UpdateDynamicObjects(const Matrix4x4 &m);
+
+  virtual void DeleteDynamicObjects();
+  
+protected:
+
+  // The kd-tree for dynamic objects
+  CKTB *mDynKtbtree;
+
+#ifdef _USE_HAVRAN_SSE
+  // This has to be aligned by 16 Bytes boundary - if HavranDynRayCaster
+  // is aligned, then also this data entity below !!!
+  static GALIGN16 RayPacket2x2 raypack_t;
+#endif
+
+  ObjectContainer *dynobjects;
+  bool dynamicFlag;
+
+  Matrix4x4 matTr, matTr_inv;
+
+  int     result4_t[4];
+  float   dist4_t[4];
+  SimpleRay sray_t;
+  
+  Vector3 orig[16];
+  Vector3 dirs[16];
+  float   tdist[32];
+  Intersectable* objI[32];
+  Vector3 normal[32];
+
+  // This transforms the ray
+  void ApplyTransform(SimpleRay &sray) {
+    sray.mOrigin = matTr_inv * sray.mOrigin;
+    sray.mDirection = RotateOnly(matTr_inv, sray.mDirection);
+    // note that normalization to the unit size of the direction
+    // is NOT computed -- this is what we want.
+  }
+};
+
+  
 } // namespace
 
-#endif
+#endif // _HavranRayCaster_H__
Index: GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro	(revision 2627)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro	(revision 2629)
@@ -34,5 +34,5 @@
   $$NONGTP/Zlib/include $$NONGTP/Boost $$NONGTP/Devil/include sparsehash/src
 
-win32:LIBPATH += GL $$XERCES/lib $$NONGTP/Devil/lib \
+win32:LIBPATH += GL $$NONGTP/Xerces/xerces/lib $$NONGTP/Devil/lib \
   $$NONGTP/glut $$NONGTP/Boost/lib $$NONGTP/Zlib/lib 
 
@@ -61,10 +61,8 @@
 CONFIG(release) {
 win32:LIBS += xerces-c_2.lib 
-#win32:LIBS += xerces-c_static_2.lib 
 }
 
 CONFIG(debug) {
 win32:LIBS += xerces-c_2D.lib 
-#win32:LIBS += xerces-c_static_2D.lib 
 }
 
@@ -127,5 +125,5 @@
 Intersectable.cpp TraversalTree.cpp ObjectPvs.cpp ObjectsParser.cpp \
 FilterBasedDistribution.cpp DifferenceSampling.cpp HavranRayCaster.cpp \
-OcclusionQuery.cpp RandomViewCellsHandler.cpp
+HavranDynRayCaster.cpp OcclusionQuery.cpp 
 
 
@@ -135,5 +133,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/ktbf2trv.cpp havran/timer.cpp
+havran/ktb8b.cpp havran/ktbftrav.cpp havran/timer.cpp 
 
 
Index: GTP/trunk/Lib/Vis/Preprocessing/src/Camera.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/Camera.cpp	(revision 2603)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/Camera.cpp	(revision 2629)
@@ -7,4 +7,5 @@
 #include <cassert>
 
+#include "common.h"
 #include "Camera.h"
 #include "Ray.h"
@@ -16,8 +17,8 @@
 #include "Preprocessor.h"
 #include "RayCaster.h"
-
-#ifdef USE_HAVRAN_RAYCASTER 
+#include "tgaimg.h"
+
+#ifdef USE_HAVRAN_RAYCASTER
 //#include "timer.h"
-#include "ktbconf.h"
 #include "raypack.h"
 #endif
@@ -54,11 +55,5 @@
 
   bool exportRays = false;
-
-  
-  InitDevIl();
-  int components = 4;
-  float *buffer = new float[components*mWidth*mHeight];
-  float *pbuffer = buffer;
-  //  - components*mWidth;
+  CTGAImage tgaimg; tgaimg.Init(mWidth, mHeight);
   
   vector<Ray *> rays;
@@ -69,5 +64,5 @@
 
   for (y = 0; y < mHeight; y++) {
-    cout<<"+";
+    cout<<"+" << flush;
     for (x = 0; x < mWidth; x++) {
       SetupRay(ray, mWidth - (x + 1), mHeight - (y + 1));
@@ -91,24 +86,13 @@
 		  if (mesh->GetMesh()->mMaterial)
 			color = mesh->GetMesh()->mMaterial->mDiffuseColor;
-		  
-		  pbuffer[0] = color.r;
-		  pbuffer[1] = color.g;
-		  pbuffer[2] = color.b;
-		  pbuffer[3] = 1.0f;
-		  
+
+		  tgaimg.SetPixel(x, y, color.r * 255.f, color.g * 255.f, color.b * 255.f);
 		}
 	  }
-      pbuffer+=components;
-      
-	  //      if (exportRays && ( 1||(x==222) && (y==97))) {
-	  //	  if (exportRays && ((x%4==0) && (y%4==0))) {
-
-	  //	  if (exportRays && ray.intersections.empty()) {
 	  if (debug) {
-		Ray *nray = new Ray(ray);
-		rays.push_back(nray);
+	    Ray *nray = new Ray(ray);
+	    rays.push_back(nray);
 	  }
     }
-    //    pbuffer-=2*components*mWidth;
   }
   
@@ -119,9 +103,5 @@
   
   cout<<"Saving image"<<endl;
-  
-  ilRegisterType(IL_FLOAT);
-  ilTexImage(mWidth, mHeight, 1, 4, IL_RGBA, IL_FLOAT, buffer);
-  ilSaveImage((char *const)filename.c_str());
-  delete buffer;
+  tgaimg.SaveToFile(filename.c_str());  
 
   cout<<"done."<<endl<<flush;
@@ -183,6 +163,5 @@
   Vector3 xv = mRight*((x - mWidth/2)/(float)mWidth);
   Vector3 yv = mUp*((y - mHeight/2)/(float)mHeight);
-  Vector3 target = xv + yv + mDirection;
-  Vector3 dir = target - mPosition;
+  Vector3 dir = mDirection - xv + yv;
   dir.Normalize();
   
@@ -202,12 +181,6 @@
 
   bool exportRays = false;
-
-  
-  InitDevIl();
-  int components = 4;
-  float *buffer = new float[components*mWidth*mHeight];
-  assert(buffer);
-  float *pbuffer = buffer;
-  //  - components*mWidth;
+  CTGAImage tgaimg; tgaimg.Init(mWidth, mHeight);
+  tgaimg.FillInImage(0.0, 0.0f, 200.f);
   
   vector<Ray *> rays;
@@ -218,42 +191,48 @@
   VssRayContainer vssRays;
 
+  bool doubleRays = true;
+  bool exploitDoubleRays = true;
+  if (!doubleRays)
+    exploitDoubleRays = false;
+  
   //CTimer timer;
   //timer.Start();
   
-  for (y = 0; y < mHeight; y++) {
-    cout<<"+";
+  //int ystart = 92;
+  int xstart = 0;
+  int ymax = mHeight;
+  for (y = 0; y < ymax; y++) {
+    cout<<"+" << flush;
     for (x = 0; x < mWidth; x++) {
-      SetupRay(ray, mWidth - (x + 1), mHeight - (y + 1));
-
-      bool debug = true;
+      SetupRay(ray, x, y);
+      if (exploitDoubleRays)
+        ray.mDirection = - ray.mDirection;
+
+      bool debug = false;
       //	  (y == mHeight/2) && (x== mWidth/3);
       //	  MeshDebug = debug;
-      
 
       int res = raycaster->CastRay(ray,
 				   vssRays,
 				   bbox,
-				   false, // castDoubleRay,
+				   doubleRays, // castDoubleRay,
 				   false); // pruneInvalidRays
+
+      // cout << "End ------------------------ "
+      // << SimpleRay::IntersectionRes[0].intersectable
+      // << " t = " << SimpleRay::IntersectionRes[0].tdist
+      // << endl;
       
       if (res) {
-	Vector3 normal = raycaster->intersect.mNormal;
+	Vector3 normal = SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
 	float v =
 	  ray.mDirection.x * normal.x +
 	  ray.mDirection.y * normal.y +
 	  ray.mDirection.z * normal.z;
-	v *= 1.0f;
-	pbuffer[0] = v;
-	pbuffer[1] = v;
-	pbuffer[2] = v;
-	pbuffer[3] = 1.0f;		  
+	v = fabs(v);
+	v *= 200.0f;
+	if (v < 50.f) v = 50.f;
+	tgaimg.SetPixel(x, y, v, v, v);
       }
-      else {
-	pbuffer[0] = 0.0;
-	pbuffer[1] = 0.0;
-	pbuffer[2] = 0.0;
-	pbuffer[3] = 1.0f;
-      }
-      pbuffer+=components;
       
       if (debug) {
@@ -263,5 +242,4 @@
       }
     }
-    //    pbuffer-=2*components*mWidth;
   } // for y
 
@@ -273,9 +251,11 @@
   
   cout<<"Saving image"<<endl;
-  
-  ilRegisterType(IL_FLOAT);
-  ilTexImage(mWidth, mHeight, 1, 4, IL_RGBA, IL_FLOAT, buffer);
-  ilSaveImage((char *const)filename.c_str());
-  delete buffer;
+
+  tgaimg.SaveToFile(filename.c_str());  
+  
+  //ilRegisterType(IL_FLOAT);
+  //ilTexImage(mWidth, mHeight, 1, 4, IL_RGBA, IL_FLOAT, buffer);
+  //ilSaveImage((char *const)filename.c_str());
+  //delete buffer;
 
   cout<<"done."<<endl<<flush;
@@ -284,4 +264,123 @@
 }
 
+
+bool
+Camera::SnapImage2(string filename,
+		   RayCaster *raycaster,
+		   AxisAlignedBox3 &bbox,
+		   SceneGraph *sceneGraph)
+{
+  int x;
+  int y;
+
+  bool exportRays = false;
+  CTGAImage tgaimg; tgaimg.Init(mWidth, mHeight);
+  tgaimg.FillInImage(0.0, 0.0f, 200.f);
+  
+  vector<Ray *> raysDebug;
+
+  long t1 = GetTime();
+
+  SimpleRay ray;
+  SimpleRayContainer rays;
+  VssRayContainer vssRays;
+
+  bool doubleRays = false;
+  bool exploitDoubleRays = false;
+  int batchSize = 16;
+  if (!doubleRays)
+    exploitDoubleRays = false;
+  int shiftIndex = 0;
+  if (exploitDoubleRays)
+    shiftIndex = batchSize;
+    
+  //CTimer timer;
+  //timer.Start();
+  
+  //int ystart = 92;
+  int xstart = 0;
+  int ymax = mHeight;
+  for (y = 0; y < ymax; y++) {
+    if (y == 250) {
+      cout << "Debug y = " << y << endl;
+    }
+    
+    cout<<"+" << flush;
+    for (x = 0; x < mWidth; x += batchSize ) {
+      rays.erase(rays.begin(), rays.end());
+      int xi;
+      for (xi = 0; (xi < batchSize) && (x+xi < mWidth); xi++ ) {
+	SetupRay(ray, x + xi, y);
+	if (exploitDoubleRays)
+	  ray.mDirection = -ray.mDirection;
+	rays.push_back(ray);
+      } // for
+
+      bool debug = false;
+      //	  (y == mHeight/2) && (x== mWidth/3);
+      //	  MeshDebug = debug;
+
+#if 1
+      assert(batchSize == 16);
+      raycaster->CastRays16(rays,
+			    vssRays,
+			    bbox,
+			    doubleRays, // castDoubleRay,
+			    false); // pruneInvalidRays
+#else
+      raycaster->CastSimpleForwardRays(rays, bbox);
+#endif
+      
+      // cout << "End ------------------------ "
+      // << SimpleRay::IntersectionRes[0].intersectable
+      // << " t = " << SimpleRay::IntersectionRes[0].tdist
+      // << endl;
+      
+      for (xi = 0; (xi < batchSize) && (x+xi < mWidth); xi++ ) {
+	Intersectable* res =
+	  SimpleRay::IntersectionRes[xi + shiftIndex].intersectable;
+
+	if (res) {
+	  Vector3 normal = res->GetNormal(0);
+	  float v =
+	    rays[xi].mDirection.x * normal.x +
+	    rays[xi].mDirection.y * normal.y +
+	    rays[xi].mDirection.z * normal.z;
+	  v = fabs(v);
+	  v *= 200.0f;
+	  if (v < 50.f) v = 50.f;
+	  tgaimg.SetPixel(x + xi, y, v, v, v);
+	}
+      }
+      
+      if (debug) {
+	Ray *nray = new Ray(ray.mOrigin, ray.mDirection,
+			    Ray::LOCAL_RAY);
+	raysDebug.push_back(nray);
+      }
+    } // for x
+  } // for y
+
+  //timer.Stop();
+  
+  long t2 = GetTime();
+  cout<<"\n#RAY_CAST_TIME = ";
+  cout << TimeDiff(t1, t2)<<" [mikrosec]\n";  
+  
+  cout<<"Saving image"<<endl;
+
+  tgaimg.SaveToFile(filename.c_str());  
+  
+  //ilRegisterType(IL_FLOAT);
+  //ilTexImage(mWidth, mHeight, 1, 4, IL_RGBA, IL_FLOAT, buffer);
+  //ilSaveImage((char *const)filename.c_str());
+  //delete buffer;
+
+  cout<<"done."<<endl<<flush;
+    
+  return true;
+}
+
+  
 bool
 Camera::SnapImagePacket(string filename,
@@ -291,18 +390,15 @@
 			)
 {
-#if defined(USE_HAVRAN_RAYCASTER) && defined(_USE_HAVRAN_SSE)
+#ifdef USE_HAVRAN_RAYCASTER 
+#ifdef _USE_HAVRAN_SSE
   
   int x;
   int y;
 
-  bool exportRays = false;
-
-  
-  InitDevIl();
-  int components = 4;
-  float *buffer = new float[components*mWidth*mHeight];
-  assert(buffer);
-  float *pbuffer = buffer;
+  bool exportRays = false;  
+
   //  - components*mWidth;
+  CTGAImage tgaimg; tgaimg.Init(mWidth, mHeight);
+  tgaimg.FillInImage(0.0, 0.0f, 200.f); // blue color image
   
   vector<Ray *> rays;
@@ -315,49 +411,77 @@
   //CTimer timer;
   //timer.Start();
-  
-  RayPacket2x2 rp;
+
+  bool doubleRays = false;
+  bool exploitDoubleRays = false;
+  if (!doubleRays)
+    exploitDoubleRays = false;
+  
+  GALIGN16 RayPacket2x2 rp;
+  //int ystart = 92;
+  //int xstart = 738;
+  int pixels = 0;;
   for (y = 0; y < mHeight-1; y+=2) {
-    float *pppbuffer = pbuffer;
-    cout<<"+";
+  //for (y = ystart-2; y > 0; y-=2) {
     for (x = 0; x < mWidth-1; x+=2) {
       int i = 0;
       for (int yi = 0; yi < 2; yi++) {
 	for (int xi = 0; xi < 2; xi++) {
-	  SetupRay(ray, mWidth - (x+xi + 1), mHeight - (y+yi + 1));
+	  SetupRay(ray, x+xi, y+yi);
 	  rp.SetLoc(i, ray.mOrigin);
 	  rp.SetDir(i, ray.mDirection);
+	  if (exploitDoubleRays)
+	    rp.SetDir(i, -ray.mDirection);
 	  i++;
 	} // for xi
       } // for yi      
-
-      raycaster->CastRaysPacket2x2(rp, false, false);
+      pixels += 4;
+      if (pixels > mWidth) {
+	cout << "+" << flush;
+	pixels = 0;
+      }
+      
+      raycaster->CastRaysPacket2x2(rp, doubleRays, false);
 
       i = 0;
-      float *ppbuffer = pbuffer;
       for (int yi = 0; yi < 2; yi++) {
 	for (int xi = 0; xi < 2; xi++) {
 	  Intersectable* res = rp.GetObject(i);
 	  if (res) {
-	    pbuffer[0] = 1.0f;
-	    pbuffer[1] = 1.0f;
-	    pbuffer[2] = 1.0f;
-	    pbuffer[3] = 1.0f;		  
-	  }
-	  else {
-	    pbuffer[0] = 0.0;
-	    pbuffer[1] = 0.0;
-	    pbuffer[2] = 0.0;
-	    pbuffer[3] = 1.0f;
+#if 0
+	    // This is debugging code to find a bug
+	    for (int j = 0; j < 4; j++) {
+	      cout << " j = " << j << " obj = "
+		   << rp.GetObject(j)
+		   << " t = " << rp.GetT(j);
+	      ray.Set(rp.GetLoc(j), rp.GetDir(j), 0, 1.0f, 0);
+	      int res = raycaster->CastRay(ray,
+					   vssRays,
+					   bbox,
+					   doubleRays, // castDoubleRay,
+					   false); // pruneInvalidRays
+	      cout << " correct result obj = "
+		   << SimpleRay::IntersectionRes[0].intersectable
+		   << " t = "
+		   << SimpleRay::IntersectionRes[0].tdist
+		   << endl;		
+	    } // for
+
+	    raycaster->CastRaysPacket2x2(rp, false, false);
+#endif
+	    
+	    Vector3 normal = res->GetNormal(0);
+	    float v =
+	      ray.mDirection.x * normal.x +
+	      ray.mDirection.y * normal.y +
+	      ray.mDirection.z * normal.z;
+	    v = fabs(v);
+	    v *= 200.0f;
+	    if (v < 50.f) v = 50.f;
+	    tgaimg.SetPixel(x+xi, y+yi, v, v, v);
 	  }
 	  i++;
-	  ppbuffer += components;
 	} // xi
-	ppbuffer += components * (mWidth-2);
       } // yi
-
-      pbuffer = ppbuffer + 2 * components;
     } // for x
-    pbuffer = pppbuffer + mWidth * 2 * components;
-    pppbuffer = pbuffer; // for the next time;
   } // for y
 
@@ -369,10 +493,8 @@
   
   cout<<"Saving image"<<endl;
-  
-  ilRegisterType(IL_FLOAT);
-  ilTexImage(mWidth, mHeight, 1, 4, IL_RGBA, IL_FLOAT, buffer);
-  ilSaveImage((char *const)filename.c_str());
-  delete buffer;
-#endif
+  tgaimg.SaveToFile(filename.c_str());  
+    
+#endif // _USE_HAVRAN_SSE
+#endif // USE_HAVRAN_RAYCASTER 
   cout<<"done."<<endl<<flush;
     
@@ -380,4 +502,131 @@
 }
   
-
-}
+inline static bool ilt(Intersectable *obj1, Intersectable *obj2)
+{
+  return obj1->mId < obj2->mId;
+}
+
+
+bool
+Camera::SnapImagePacket2(string filename,
+			RayCaster *raycaster,
+			AxisAlignedBox3 &bbox,
+			SceneGraph *sceneGraph
+			)
+{
+#ifdef USE_HAVRAN_RAYCASTER 
+#ifdef _USE_HAVRAN_SSE
+  
+  int x;
+  int y;
+
+  bool exportRays = false;  
+
+  //  - components*mWidth;
+  CTGAImage tgaimg; tgaimg.Init(mWidth, mHeight);
+  tgaimg.FillInImage(0.0, 0.0f, 200.f); // blue color image
+  
+  vector<Ray *> rays;
+
+  long t1 = GetTime();
+
+  SimpleRay ray;
+  VssRayContainer vssRays;
+
+  //CTimer timer;
+  //timer.Start();
+
+  bool doubleRays = false;
+  bool exploitDoubleRays = false;
+  if (!doubleRays)
+    exploitDoubleRays = false;
+  
+  Vector3 origin4[4];
+  Vector3 direction4[4];
+  int     result4[4];
+  float   dist4[4];
+  
+  GALIGN16 RayPacket2x2 rp;
+  //int ystart = 92;
+  //int xstart = 738;
+  int pixels = 0;;
+  for (y = 0; y < mHeight-1; y+=2) {
+  //for (y = ystart-2; y > 0; y-=2) {
+    for (x = 0; x < mWidth-1; x+=2) {
+      int i = 0;
+      for (int yi = 0; yi < 2; yi++) {
+	for (int xi = 0; xi < 2; xi++) {
+	  SetupRay(ray, x+xi, y+yi);
+	  origin4[i] = ray.mOrigin;
+	  direction4[i] = ray.mDirection;
+	  if (exploitDoubleRays)
+	    direction4[i] = -ray.mDirection;
+	  i++;
+	} // for xi
+      } // for yi      
+      pixels += 4;
+      if (pixels > mWidth) {
+	cout << "+" << flush;
+	pixels = 0;
+      }
+      
+      raycaster->CastRaysPacket4(bbox.Min(), bbox.Max(),
+				 origin4,
+				 direction4,
+				 result4,
+				 dist4);
+
+      i = 0;
+      for (int yi = 0; yi < 2; yi++) {
+	for (int xi = 0; xi < 2; xi++) {
+	  int objId = result4[i];
+	  if (objId != -1) {
+	    MeshInstance dummyInst(NULL);
+	    dummyInst.SetId(objId);
+	    ObjectContainer &objects = preprocessor->mObjects;
+	    ObjectContainer::const_iterator oit =
+	      lower_bound(objects.begin(), objects.end(),
+			  (Intersectable *)&dummyInst, ilt);
+      
+	    if ((oit != objects.end()) && ((*oit)->GetId() == objId)) {
+	      Intersectable *res = *oit;
+	    
+	      Vector3 normal = res->GetNormal(0);
+	      float v =
+		ray.mDirection.x * normal.x +
+		ray.mDirection.y * normal.y +
+		ray.mDirection.z * normal.z;
+	      v = fabs(v);
+	      v *= 200.0f;
+	      if (v < 50.f) v = 50.f;
+	      tgaimg.SetPixel(x+xi, y+yi, v, v, v);
+	    }
+	    else {
+	      cout <<"*" << endl;
+	    }
+	    
+	  } // object intersected
+	  i++;
+	} // xi
+      } // yi
+    } // for x
+  } // for y
+
+  //timer.Stop();
+  
+  long t2 = GetTime();
+  cout<<"\n#RAY_CAST_TIME = ";
+  cout << TimeDiff(t1, t2)<<" [mikrosec]\n";  
+  
+  cout<<"Saving image"<<endl;
+  tgaimg.SaveToFile(filename.c_str());  
+    
+#endif // _USE_HAVRAN_SSE
+#endif // USE_HAVRAN_RAYCASTER 
+  cout<<"done."<<endl<<flush;
+    
+  return true;
+}
+
+}
+
Index: GTP/trunk/Lib/Vis/Preprocessing/src/MeshKdTree.cpp
===================================================================
--- GTP/trunk/Lib/Vis/Preprocessing/src/MeshKdTree.cpp	(revision 2614)
+++ GTP/trunk/Lib/Vis/Preprocessing/src/MeshKdTree.cpp	(revision 2629)
@@ -418,5 +418,5 @@
   mRoot = Subdivide(TraversalData(leaf, NULL, GetBox(), 0));
   
-  cout <<  "created " << 
+  cout <<  "created " << endl;
   // remove the allocated array
   delete mSubdivisionCandidates;
