Changeset 2592 for GTP/trunk/Lib/Vis/Preprocessing/src

Timestamp:

01/14/08 15:54:15 (17 years ago)

Author:

bittner

Message:

havran ray caster update

Location:

GTP/trunk/Lib/Vis/Preprocessing/src

Files:

• HavranRayCaster.cpp (modified) (18 diffs)
• HavranRayCaster.h (modified) (4 diffs)
• Makefile (modified) (2 diffs)
• Mutation.cpp (modified) (16 diffs)
• Preprocessor.cpp (modified) (1 diff)
• RayCaster.cpp (modified) (1 diff)
• default.env (modified) (1 diff)
• havran/configh.h (modified) (1 diff)
• havran/ktb.cpp (modified) (2 diffs)
• havran/ktb.h (modified) (3 diffs)
• havran/ktb8b.cpp (modified) (5 diffs)
• havran/ktb8b.h (modified) (5 diffs)
• havran/ktbai.cpp (modified) (1 diff)
• havran/ktball.cpp (modified) (14 diffs)
• havran/ktball.h (modified) (5 diffs)
• havran/ktbconf.h (modified) (1 diff)
• havran/ktbftrav.cpp (modified) (21 diffs)
• havran/ktbtrav.h (modified) (9 diffs)
• havran/raypack.h (modified) (3 diffs)
• preprocessor.pro (modified) (1 diff)

GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.cpp

@@ -6 +6 @@
 
 
-#include "HavranRayCaster.h"
 #include "VssRay.h"
 #include "KdTree.h"
@@ -13 +12 @@
 
 #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 +76 @@
 }
 
-
+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 +116 @@
     // 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 +169 @@
   
   // 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 +196 @@
   }
 
+#ifdef _PROCESS_RAY  
   // This code is also in IntelRayCaster.cpp
   return ProcessRay(
@@ -161 +207 @@
                     pruneInvalidRays
                     );
-#else // USE_HAVRAN_RAYCASTER 
+#else // _PROCESS_RAY
+  return result;
+#endif // _PROCESS_RAY  
+
+#else
   return 0;
-#endif
+#endif // USE_HAVRAN_RAYCASTER   
 }
 
@@ -190 +240 @@
 
 #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 +249 @@
   }
 #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 +270 @@
       (*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 +284 @@
     }
     // 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 +334 @@
     {
       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 +349 @@
 #endif
 
-#if 1
+#ifdef _PROCESS_RAY  
     ProcessRay(rays[k],
-                           hitA,
-                           hitB,
-                           vssRays,
-                           sbox,
-                           castDoubleRays,
-                           pruneInvalidRays
-                           );
+               hitA,
+               hitB,
+               vssRays,
+               sbox,
+               castDoubleRays,
+               pruneInvalidRays
+               );
 #endif
   } // for
@@ -340 +391 @@
   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 +423 @@
   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 +433 @@
   }
 
+  // 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 +447 @@
 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 +499 @@
 
   return;
+#endif // _USE_HAVRAN_SSE
 #endif // USE_HAVRAN_RAYCASTER 
 }
-#endif
+#endif // _USE_HAVRAN_SSE
 
 

GTP/trunk/Lib/Vis/Preprocessing/src/HavranRayCaster.h

@@ -5 +5 @@
 #include "Containers.h"
 #include <string>
-#include "havran/raypack.h"
-
+#include "ktbconf.h"
+#include "raypack.h"
 
 namespace GtpVisibilityPreprocessor {
@@ -15 +15 @@
 class KdTree;
 class Ray;
-class RayPacket2x2;
 class SimpleRayContainer;
 class AxisAlignedBox3;
@@ -38 +37 @@
   
   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 +71 @@
                         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;

GTP/trunk/Lib/Vis/Preprocessing/src/Makefile

@@ -1 +1 @@
 #############################################################################
 # 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 +63 @@
         $(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 \

GTP/trunk/Lib/Vis/Preprocessing/src/Mutation.cpp

@@ -94 +94 @@
 
                 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 +107 @@
                 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 +118 @@
           }
         } 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 +152 @@
                 Intersectable *newObject = vssRays[i]->mTerminationObject;
 
-
                 if (oldObject == newObject)
                   dummyNcMutations++;
@@ -163 +165 @@
         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 +235 @@
                                                                                                  )
 {
-#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 +330 @@
   
   AxisAlignedBox3 box = occluder->GetBox();
+  // consider slightly larger neighborhood of the occluder
+  // in search for unocclude reverse ray
   box.Scale(2.0f);
   
@@ -430 +423 @@
   }
 
+  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 +504 @@
 
 #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 +517 @@
   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 +526 @@
                                                                                                   Vector3(rr[0], rr[1], rr[2]));
         
-
+        
         direction = UniformRandomVector(rr[3], rr[4]);
         
@@ -538 +532 @@
         sray = SimpleRay(origin, direction, MUTATION_BASED_DISTRIBUTION, pdf);
         sray.mGeneratorId = -1;
-
+        
         return true;
-  } 
+  }
 
   int index;
@@ -552 +546 @@
         if (
                 mRays[index].GetSamplingFactor() >= mRays[mLastIndex].GetSamplingFactor()) {
-          // make another round
+          // make another round to equalize the oversampling factor
 
           //      cout<<"R2"<<endl;
@@ -666 +660 @@
   
   if (mUseSilhouetteSamples && a < mSilhouetteProb) {
-        
         termination += ComputeSilhouetteTerminationMutation(*ray,
                                                                                                                 origin,
@@ -692 +685 @@
   
   // 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 +716 @@
     mRays[index].mMutations++;
         mRays[index].mUnsuccessfulMutations++;
-
         return true;
   }

GTP/trunk/Lib/Vis/Preprocessing/src/Preprocessor.cpp

@@ -1230 +1230 @@
         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;
         }
 

GTP/trunk/Lib/Vis/Preprocessing/src/RayCaster.cpp

@@ -178 +178 @@
   //const int batchsize = 16384;
   const int batchsize = 8192;
+  // const int batchsize = 1024;
   //const int batchsize = 128;
 

GTP/trunk/Lib/Vis/Preprocessing/src/default.env

@@ -53 +53 @@
         useGlDebugger false
 # 0 = INTERNAL  1 = MLRT 2 = HAVRAN
-        rayCastMethod 1
+        rayCastMethod 2
         
 #       type sampling

GTP/trunk/Lib/Vis/Preprocessing/src/havran/configh.h

@@ -59 +59 @@
 
 #ifdef _MSC_VER
+#define __SSE__
 // disable some Microsoft Visual Compiler warnings if necessary
 #pragma warning( disable:4305 )

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.cpp

@@ -61 +61 @@
   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 +353 @@
 }
 
+// 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

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb.h

@@ -156 +156 @@
 #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 +177 @@
   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 +224 @@
   void SetInteriorNodeLinks(SKTBNodeT *node,
                             SKTBNodeT *leftChild,
+                            SKTBNodeT *rightChild);
+  void SetInteriorNodeLeft(SKTBNodeT *node,
+                           SKTBNodeT *leftChild);
+  void SetInteriorNodeRight(SKTBNodeT *node,
                             SKTBNodeT *rightChild);
 

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.cpp

@@ -63 +63 @@
   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 +76 @@
   return intersected;
 }
-
   
 // test all objects in the leaf for intersection with ray
@@ -244 +243 @@
 CKTB8BAllocMan::AllocInteriorNodeWithBox(int axis, float position,
                                          int cntLeft, int cntRight,
-                                         const SBBox &tsbox, SKTBNodeT* prevMinBoxNode,
+                                         const SBBox &tsbox,
+                                         SKTBNodeT* prevMinBoxNode,
                                          int depthStore)
 {
@@ -354 +354 @@
   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 +398 @@
 }
 
+// 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

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktb8b.h

@@ -161 +161 @@
 #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 +182 @@
   ~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 +232 @@
   void SetInteriorNodeLinks(SKTBNodeT *node,
                             SKTBNodeT *leftChild,
+                            SKTBNodeT *rightChild);
+
+  void SetInteriorNodeLeft(SKTBNodeT *node,
+                           SKTBNodeT *leftChild);
+  void SetInteriorNodeRight(SKTBNodeT *node,
                             SKTBNodeT *rightChild);
 
@@ -275 +280 @@
   void DecDepth() {_currDepth--;}
 
+  
   // --------------------------------------------------------------
   // Statistics
@@ -300 +306 @@
   float _sumSurfaceAreaInteriorNodes;
 #endif
-
+  
   // General statistics on the depth of current node
   // the depth of the currently accessed node during building up

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbai.cpp

@@ -1221 +1221 @@
   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

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.cpp

@@ -22 +22 @@
 #include <fstream>
 #include <iomanip>
+#include <stack>
 
 namespace GtpVisibilityPreprocessor {
@@ -96 +97 @@
   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 +118 @@
   }
 #else
-  cerr<<"th3"<<endl;
+  // cerr<<"th3"<<endl;
 
   traversalClass = new CKTBTraversal; 
@@ -308 +309 @@
 #endif // _DEBUG
 
-  cerr<<"hh2"<<endl;
+  // cerr<<"hh2"<<endl;
   // possibly remove the old building class
   // and allocate the class for CKTB building up
@@ -315 +316 @@
   buildClass = AllocBuildClass();
 
-  cerr<<"hh3"<<endl;
+  // cerr<<"hh3"<<endl;
 
   // and for traversal
   AllocTraversalClass();
-  cerr<<"hh4"<<endl;
+  // cerr<<"hh4"<<endl;
 
   if (boxToBoundObjects) {
@@ -330 +331 @@
     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 +401 @@
 }
 
+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 +431 @@
 
 #endif // __SSE__
+
+#endif // _USE_HAVRAN_SSE  
+
   
 void
@@ -752 +766 @@
 SKTBNodeT*
 CKTB::ImportBinLeaf(IN_STREAM &stream, 
-                    SKTBNodeT *parent,
+                    SKTBNodeT **nodeToLink,
                     const ObjectContainer &objects)
 {
@@ -764 +778 @@
     // since for 
     SKTBNodeT* leaf = buildClass->AllocLeaf(0);
+    *nodeToLink = buildClass->nodeToLink;
     return leaf;
   }
@@ -769 +784 @@
   MeshInstance dummyInst(NULL);
   SKTBNodeT* leaf = buildClass->AllocLeaf(size);
+  *nodeToLink = buildClass->nodeToLink;
 
   ObjectContainer *newobjlist = new ObjectContainer;
@@ -809 +825 @@
 
 SKTBNodeT *
-CKTB::ImportBinInterior(IN_STREAM  &stream, SKTBNodeT *parent)
+CKTB::ImportBinInterior(IN_STREAM  &stream, SKTBNodeT **nodeToLink)
 {
   int axis;
@@ -819 +835 @@
   SKTBNodeT *interiorNode =
     buildClass->AllocInteriorNode(axis, pos, 0, 0);
+  *nodeToLink = buildClass->nodeToLink;
         
   return interiorNode;
@@ -830 +847 @@
   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

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktball.h

@@ -174 +174 @@
   
   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 +187 @@
   
 #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 +193 @@
   void FindNearestI(RayPacket2x2 &raypack);
   void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax);
+#endif
 #else
   void FindNearestI(RayPacket2x2 &raypack) { }
@@ -199 +202 @@
   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 +222 @@
   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) { }
+  };
 };
 

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbconf.h

@@ -22 +22 @@
 #include <xmmintrin.h>
 #endif
+
+// If we support the use of SSE instructions for ray shooting
+#define _USE_HAVRAN_SSE
 
 namespace GtpVisibilityPreprocessor {

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbftrav.cpp

@@ -20 +20 @@
 #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 +38 @@
   // 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 +61 @@
 
   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 +92 @@
     }
 #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 +145 @@
           // 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 +201 @@
         cout << "Link " << endl;
 #endif
-        stack3[index].nodep = GetLinkNode(currNode);
         // cout << "Link node was accessed" << endl;
-        continue;
+        currNode = GetLinkNode(currNode);
+        goto CONTINUE_LINK;
       }
     }
@@ -336 +217 @@
   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 +237 @@
   // passing through parameters
   float tmin, tmax;
+  SimpleRay::IntersectionRes[0].intersectable = 0;
   
   // test if the whole CKTB tree is missed by the input ray
@@ -361 +244 @@
        (tmax <= tmin) ||
        (tmax <= 0.f) ) {
-    SimpleRay::IntersectionRes[0].intersectable = 0;
     return 0; // no object can be intersected
   }
@@ -378 +260 @@
 
   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 +291 @@
     }
 #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 +341 @@
           // 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 +397 @@
         cout << "Link " << endl;
 #endif
-        stack3[index].nodep = GetLinkNode(currNode);
         // cout << "Link node was accessed" << endl;
-        continue;
+        currNode = GetLinkNode(currNode);
+        goto CONTINUE_LINK;
       }
     }
@@ -528 +413 @@
   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 +442 @@
   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 +450 @@
         (tmax <= tmin) ||
         (tmax <= 0.f) ) {
-      SimpleRay::IntersectionRes[i].intersectable = 0;
     }
     else {
@@ -756 +457 @@
       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 +520 @@
       }
 #endif      
-      assert(tmin <= tmax);
     
 #ifdef __TRAVERSAL_STATISTICS
@@ -876 +576 @@
             // 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 +635 @@
 #endif
           stackA[indexSA].nodep = GetLinkNode(currNode);
+          GPREFETCH(stackA[indexSA].nodep, PREF_DEFAULT);
           // cout << "Link node was accessed" << endl;
           continue;
@@ -950 +652 @@
   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

GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbtrav.h

@@ -30 +30 @@
 
 namespace GtpVisibilityPreprocessor {
-
-class RayPacket2x2;
   
 //#define TRV000
@@ -136 +134 @@
 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 +145 @@
   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 +151 @@
   virtual void FindNearestI(RayPacket2x2 &raypack) { }
   virtual void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax) { }
-
+#endif
+  
   virtual void PrintStatsTR(ostream &) { }
   virtual void DynamicStatsReset() { }
@@ -306 +308 @@
   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 +322 @@
   virtual int FindNearestI_16twoDir(SimpleRayContainer &rays);
 
+#ifdef _USE_HAVRAN_SSE  
 #ifdef __SSE__
   // The same operations for packets of rays, if implemented by
@@ -324 +330 @@
   virtual void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax);
 #endif // __SSE__
-
+#endif // _USE_HAVRAN_SSE  
+  
   virtual void PrintStatsTR(ostream &) { }
 
@@ -408 +415 @@
   virtual int FindNearestI(SimpleRay &ray);
 
+#ifdef _USE_HAVRAN_SSE  
 #ifdef __SSE__
   // The same operations for packets of rays, if implemented by
@@ -415 +423 @@
   virtual void FindNearestI(RayPacket2x2 &raypack, Vector3 &boxmin, Vector3 &boxmax);
 #endif // __SSE__
+#endif // _USE_HAVRAN_SSE  
 
   virtual void DynamicStatsReset();

GTP/trunk/Lib/Vis/Preprocessing/src/havran/raypack.h

@@ -21 +21 @@
 
 #include "Vector3.h"
+#include "ktbconf.h"
+  
+#ifdef _USE_HAVRAN_SSE
   
 #ifdef __SSE__
@@ -337 +340 @@
   // 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 +433 @@
   tmax[i] = tmaxnew;
 }
-#else // __SSE__
-
-// ? What to do here
-//#error "AAA"
-
 #endif // __SSE__
 
+#endif // _USE_HAVRAN_SSE
+  
 } // namespace
 

GTP/trunk/Lib/Vis/Preprocessing/src/preprocessor.pro

@@ -129 +129 @@
 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