WebSVN – gelsvn – Diff – /trunk/src/Geometry/BSPTree.cpp

 #include "CGLA/Vec3d.h"
 #include "Geometry/TriMesh.h"
 #include "BSPTree.h"
 using namespace std;
 using namespace CGLA;
 namespace Geometry
 {
   int BSPTree::node_calls;
   int BSPTree::tri_calls;
   void create_tri_accel(Vec3f A, Vec3f B, Vec3f C, TriAccel &tri_accel)
   {
     Vec3f N = normalize(cross(B-A, C-A));
     // Find projection dir
     int k,u,v;
     if (fabs(N[0]) > fabs(N[1]))
       if (fabs(N[0]) > fabs(N[2]))
-	k = 0; /* X */
+    k = 0; /* X */
       else
-	k=2;   /* Z */
+    k=2;   /* Z */
     else
       if (fabs(N[1]) > fabs(N[2]))
-	k = 1; /* Y */
+    k = 1; /* Y */
       else
-	k=2;   /* Z */
+    k=2;   /* Z */
     u = (k+1)% 3;
     v = (k+2)% 3;
     Vec3f b = C - A;
     Vec3f c = B - A;
     double div = (b[u]*c[v]-b[v]*c[u]);
     tri_accel.n_u = N[u]/N[k];
     tri_accel.n_v = N[v]/N[k];
     tri_accel.n_d = dot(A, N/N[k]);
     tri_accel.k = k;
     tri_accel.b_nu = -b[v]/div;
     tri_accel.b_nv = b[u]/div;
     tri_accel.b_d = (b[v]*A[u]-b[u]*A[v])/div;
     tri_accel.c_nu = c[v]/div;
     tri_accel.c_nv = -c[u]/div;
     tri_accel.c_d = (c[u]*A[v]-c[v]*A[u])/div;
   }
   // Most of this file is a direct copy from Henrik's notes
   BSPTree::BSPTree()
   {
     root=0;
     b_is_build = false;
   }
   BSPTree::~BSPTree()
   {
     clear();
   }
   void BSPTree::clear()
   {
     if (root!=0)
     {
       delete_node(root);
       root=0;
       b_is_build=false;
     }
     isecttris.clear();
     all_objects.clear();
     all_triaccel.clear();
   }
   void BSPTree::delete_node(BSPNode *node)
   {
     if (node->left!=0)
       delete_node(node->left);
     if (node->right!=0)
       delete_node(node->right);
     delete node;
   }
   void BSPTree::subdivide_node(BSPNode &node, BBox &bbox,
-			       unsigned int level,
+                   unsigned int level,
-			       vector<ISectTri*>& objects,
+                   vector<ISectTri*>& objects,
-			       vector<TriAccel*>& tri_objects)
+                   vector<TriAccel*>& tri_objects)
   {
     const int TESTS = 2;
     if (objects.size()<=max_objects || level==max_level)
     {
       node.axis_leaf = 4; // Means that this is a leaf
       node.plane = 0;
       node.left = 0;
       node.right = 0;
       node.id = all_objects.size();
       node.count = objects.size();
       for(unsigned int i = 0; i < objects.size(); ++i)
       {
-	all_objects.push_back(objects[i]);
+    all_objects.push_back(objects[i]);
-	all_triaccel.push_back(tri_objects[i]);
+    all_triaccel.push_back(tri_objects[i]);
       }
     }
     else
     {
       bool right_zero=false;
       bool left_zero=false;
       unsigned int i;
       BSPNode* left_node  = new BSPNode();
       BSPNode* right_node = new BSPNode();
       vector<ISectTri*> left_objects;
       vector<ISectTri*> right_objects;
       vector<TriAccel*> tri_left_objects;
       vector<TriAccel*> tri_right_objects;
       node.left = left_node;
       node.right = right_node;
       int new_axis=-1;
       double min_cost=CGLA::BIG;
       int new_pos = -1;
-	  for(i=0;i<3;i++)
+      for(i=0;i<3;i++)
       {
-	for(int k=1;k<TESTS;k++)
+    for(int k=1;k<TESTS;k++)
-	{
+    {
-	  BBox left_bbox = bbox;
+      BBox left_bbox = bbox;
-	  BBox right_bbox = bbox;
+      BBox right_bbox = bbox;
-	  double center = (bbox.max_corner[i]- bbox.min_corner[i])*(double)k/(double)TESTS + bbox.min_corner[i];
+      double center = (bbox.max_corner[i]- bbox.min_corner[i])*(double)k/(double)TESTS + bbox.min_corner[i];
-	  node.plane = center;
+      node.plane = center;
-	  left_bbox.max_corner[i] = center;
+      left_bbox.max_corner[i] = center;
-	  right_bbox.min_corner[i] = center;
+      right_bbox.min_corner[i] = center;
-	  // Try putting the triangles in the left and right boxes
+      // Try putting the triangles in the left and right boxes
-	  int left_count = 0;
+      int left_count = 0;
-	  int right_count = 0;
+      int right_count = 0;
-	  for(unsigned int j=0;j<objects.size();j++)
+      for(unsigned int j=0;j<objects.size();j++)
-	  {
+      {
-	    ISectTri* tri = objects[j];
+        ISectTri* tri = objects[j];
-	    left_count += left_bbox.intersect_triangle(*tri);
+        left_count += left_bbox.intersect_triangle(*tri);
-	    right_count += right_bbox.intersect_triangle(*tri);
+        right_count += right_bbox.intersect_triangle(*tri);
-	  }
+      }
-	  //double len = bbox.max_corner[i] - bbox.min_corner[i];
+      //double len = bbox.max_corner[i] - bbox.min_corner[i];
-	  double cost = left_count*left_bbox.area() + right_count*right_bbox.area(); // - len*len;
+      double cost = left_count*left_bbox.area() + right_count*right_bbox.area(); // - len*len;
-	  if(cost < min_cost)
+      if(cost < min_cost)
-	  {
+      {
-	    min_cost = cost;
+        min_cost = cost;
-	    new_axis = i;
+        new_axis = i;
-	    new_pos = k;
+        new_pos = k;
-	    right_zero = (right_count==0);
+        right_zero = (right_count==0);
-	    left_zero = (left_count==0);
+        left_zero = (left_count==0);
-	  }
+      }
-	}
+    }
       }
       node.axis_leaf = new_axis;
       left_node->axis_leaf = static_cast<unsigned char>(-1);
       right_node->axis_leaf = static_cast<unsigned char>(-1);
       // Now chose the right splitting plane
       BBox left_bbox = bbox;
       BBox right_bbox = bbox;
       double size = bbox.max_corner[node.axis_leaf]- bbox.min_corner[node.axis_leaf];
       double center = size*(double)new_pos/(double)TESTS + bbox.min_corner[node.axis_leaf];
       double diff = f_eps < size/8.0 ? f_eps : size/8.0;
       if (left_zero)
       {
-	// Find min position of all triangle vertices and place the center there
+    // Find min position of all triangle vertices and place the center there
-	center = bbox.max_corner[node.axis_leaf];
+    center = bbox.max_corner[node.axis_leaf];
-	for(unsigned int j=0;j<objects.size();j++)
+    for(unsigned int j=0;j<objects.size();j++)
-	{
+    {
-	  ISectTri* tri = objects[j];
+      ISectTri* tri = objects[j];
-	  if (tri->point0[node.axis_leaf]<center)
+      if (tri->point0[node.axis_leaf]<center)
-	    center=tri->point0[node.axis_leaf];
+        center=tri->point0[node.axis_leaf];
-	  if (tri->point1[node.axis_leaf]<center)
+      if (tri->point1[node.axis_leaf]<center)
-	    center=tri->point1[node.axis_leaf];
+        center=tri->point1[node.axis_leaf];
-	  if (tri->point2[node.axis_leaf]<center)
+      if (tri->point2[node.axis_leaf]<center)
-	    center=tri->point2[node.axis_leaf];
+        center=tri->point2[node.axis_leaf];
-	}
+    }
-	center -= diff;
+    center -= diff;
       }
       if (right_zero)
       {
-	// Find max position of all triangle vertices and place the center there
+    // Find max position of all triangle vertices and place the center there
-	center = bbox.min_corner[node.axis_leaf];
+    center = bbox.min_corner[node.axis_leaf];
-	for(unsigned int j=0;j<objects.size();j++)
+    for(unsigned int j=0;j<objects.size();j++)
-	{
+    {
-	  ISectTri* tri = objects[j];
+      ISectTri* tri = objects[j];
-	  if (tri->point0[node.axis_leaf]>center)
+      if (tri->point0[node.axis_leaf]>center)
-	    center=tri->point0[node.axis_leaf];
+        center=tri->point0[node.axis_leaf];
-	  if (tri->point1[node.axis_leaf]>center)
+      if (tri->point1[node.axis_leaf]>center)
-	    center=tri->point1[node.axis_leaf];
+        center=tri->point1[node.axis_leaf];
-	  if (tri->point2[node.axis_leaf]>center)
+      if (tri->point2[node.axis_leaf]>center)
-	    center=tri->point2[node.axis_leaf];
+        center=tri->point2[node.axis_leaf];
-	}
+    }
-	center += diff;
+    center += diff;
       }
       node.plane = center;
       left_bbox.max_corner[node.axis_leaf] = center;
       right_bbox.min_corner[node.axis_leaf] = center;
       // Now put the triangles in the right and left node
       for(i=0;i<objects.size();i++)
       {
-	ISectTri* tri = objects[i];
+    ISectTri* tri = objects[i];
-	TriAccel *tri_accel = tri_objects[i];
+    TriAccel *tri_accel = tri_objects[i];
-	if (left_bbox.intersect_triangle(*tri))
+    if (left_bbox.intersect_triangle(*tri))
-	{
+    {
-	  left_objects.push_back(tri);
+      left_objects.push_back(tri);
-	  tri_left_objects.push_back(tri_accel);
+      tri_left_objects.push_back(tri_accel);
-	}
+    }
-	if (right_bbox.intersect_triangle(*tri))
+    if (right_bbox.intersect_triangle(*tri))
-	{
+    {
-	  right_objects.push_back(tri);
+      right_objects.push_back(tri);
-	  tri_right_objects.push_back(tri_accel);
+      tri_right_objects.push_back(tri_accel);
-	}
+    }
       }
     //if (left_zero||right_zero)
-	//  cout << left_objects.size() << "," << right_objects.size() << "," << level << endl;
+    //  cout << left_objects.size() << "," << right_objects.size() << "," << level << endl;
       objects.clear();
       subdivide_node(*left_node , left_bbox , level+1, left_objects, tri_left_objects);
       subdivide_node(*right_node, right_bbox, level+1, right_objects, tri_right_objects);
     }
   }
   void BSPTree::init()
   {
     root = new BSPNode();
     bbox.compute_bbox(isecttris);
     bbox.min_corner-=Vec3f(1.0);
     bbox.max_corner+=Vec3f(1.0);
   }
   void BSPTree::init(vector<const TriMesh*>& _trimesh,
-		     vector<Mat4x4f>& _transforms,
+             vector<Mat4x4f>& _transforms,
-		     int _max_objects, int _max_level)
+             int _max_objects, int _max_level)
   {
     trimesh = _trimesh;
     transforms = _transforms;
     for(unsigned int i=0;i<trimesh.size();i++)
     {
       const TriMesh *mesh = trimesh[i];
       // Loop through all triangles and add them to intersection structure
       for(int j=0;j<mesh->geometry.no_faces();j++)
       {
-	Vec3i face = mesh->geometry.face(j);
+    Vec3i face = mesh->geometry.face(j);
-	ISectTri new_tri;
+    ISectTri new_tri;
-	new_tri.point0 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[0]));
+    new_tri.point0 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[0]));
-	new_tri.point1 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[1]));
+    new_tri.point1 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[1]));
-	new_tri.point2 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[2]));
+    new_tri.point2 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[2]));
-	new_tri.edge0 = new_tri.point1 - new_tri.point0;
+    new_tri.edge0 = new_tri.point1 - new_tri.point0;
-	new_tri.edge1 = new_tri.point2 - new_tri.point0;
+    new_tri.edge1 = new_tri.point2 - new_tri.point0;
-	new_tri.mesh_id = i;
+    new_tri.mesh_id = i;
-	new_tri.tri_id = j;
+    new_tri.tri_id = j;
-	isecttris.push_back(new_tri);
+    isecttris.push_back(new_tri);
-	TriAccel ta;
+    TriAccel ta;
-	create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);
+    create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);
-	ta.mesh_id = i;
+    ta.mesh_id = i;
-	ta.tri_id = j;
+    ta.tri_id = j;
-	triaccel.push_back(ta);
+    triaccel.push_back(ta);
       }
     }
     max_objects = _max_objects;
     max_level = _max_level;
     init();
   }
   void BSPTree::init(const TriMesh* mesh, Mat4x4f transform,
-		     vector<int> &trilist,
+             vector<int> &trilist,
-		     int _max_objects, int _max_level)
+             int _max_objects, int _max_level)
   {
     trimesh.push_back(mesh);
     transforms.push_back(transform);
     // Loop through all triangles and add them to intersection structure
     for(unsigned int j=0;j<trilist.size();j++)
     {
       Vec3i face = mesh->geometry.face(trilist[j]);
       ISectTri new_tri;
       new_tri.point0 = transform.mul_3D_point(mesh->geometry.vertex(face[0]));
       new_tri.point1 = transform.mul_3D_point(mesh->geometry.vertex(face[1]));
       new_tri.point2 = transform.mul_3D_point(mesh->geometry.vertex(face[2]));
       new_tri.edge0 = new_tri.point1 - new_tri.point0;
       new_tri.edge1 = new_tri.point2 - new_tri.point0;
       new_tri.mesh_id = 0;
       new_tri.tri_id = trilist[j];
       isecttris.push_back(new_tri);
       TriAccel ta;
       create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);
       ta.mesh_id = 0;
       ta.tri_id = trilist[j];
       triaccel.push_back(ta);
     }
     max_objects = _max_objects;
     max_level = _max_level;
     init();
   }
   void BSPTree::build()
   {
     if (!b_is_build)
     {
       vector<ISectTri*> objects;
       vector<TriAccel*> tri_objects;
       for(unsigned int i=0;i<isecttris.size();i++)
       {
-	ISectTri* tri = &isecttris[i];
+    ISectTri* tri = &isecttris[i];
-	TriAccel* tri_accel = &triaccel[i];
+    TriAccel* tri_accel = &triaccel[i];
-	objects.push_back(tri);
+    objects.push_back(tri);
-	tri_objects.push_back(tri_accel);
+    tri_objects.push_back(tri_accel);
       }
       subdivide_node(*root, bbox, 0, objects, tri_objects);
       b_is_build = true;
     }
     isecttris.clear();
     all_objects.clear();
     make_fast_tree(root);
   }
   bool BSPTree::is_build()
   {
     return b_is_build;
   }
   void BSPTree::print(BSPNode *node, int depth)
   {
     if (node==0)
       return;
     for(int i=0;i<depth;i++)
       cout << " ";
-//	cout << "axis:" << node->axis_leaf << ", count:" << node->objects.size() << ", plane:" << node->plane << ", " << endl;
+//  cout << "axis:" << node->axis_leaf << ", count:" << node->objects.size() << ", plane:" << node->plane << ", " << endl;
     print(node->left, depth+1);
     print(node->right, depth+1);
   }
   int BSPTree::size(BSPNode *node)
   {
     if (node==0)
       return 0;
     int s = sizeof(BSPNode);
     s+= node->count * sizeof(ISectTri);
     s+=size(node->left);
     s+=size(node->right);
     return s;
   }
   int BSPTree::size()
   {
     return size(root);
   }
 /*__declspec(align(16))*/ static const unsigned int modulo[] = {0,1,2,0,1};
   inline bool intersect2(Ray &ray, const TriAccel &acc, double t_max)
   {
-//	cout << "Acc: " << (int)&acc << " ";
+//  cout << "Acc: " << (int)&acc << " ";
 //inline bool Intersect(TriAccel &acc,Ray &ray)
 #define ku modulo[acc.k+1]
 #define kv modulo[acc.k+2]
     // dont prefetch here, assume data has already been prefetched
     // start high-latency division as early as possible
-    const double nd = 1.0/((double)ray.direction[acc.k]	+ (double)acc.n_u * (double)ray.direction[ku] + (double)acc.n_v * (double)ray.direction[kv]);
+    const double nd = 1.0/((double)ray.direction[acc.k] + (double)acc.n_u * (double)ray.direction[ku] + (double)acc.n_v * (double)ray.direction[kv]);
-    const double f = ((double)acc.n_d - (double)ray.origin[acc.k]	- (double)acc.n_u * (double)ray.origin[ku] - (double)acc.n_v * (double)ray.origin[kv]) * nd;
+    const double f = ((double)acc.n_d - (double)ray.origin[acc.k]   - (double)acc.n_u * (double)ray.origin[ku] - (double)acc.n_v * (double)ray.origin[kv]) * nd;
     // check for valid distance.
     if (!(t_max > f && f > 0.001)||ray.dist<f) return false;
     // compute hitpoint positions on uv plane
     const double hu = (ray.origin[ku] + f * ray.direction[ku]);
     const double hv = (ray.origin[kv] + f * ray.direction[kv]);
     // check first barycentric coordinate
     const double lambda = (hu * (double)acc.b_nu + hv * (double)acc.b_nv + (double)acc.b_d);
     if (lambda < 0.0) return false;
     // check second barycentric coordinate
     const double mue = (hu * (double)acc.c_nu + hv * (double)acc.c_nv + (double)acc.c_d);
     if (mue < 0.0) return false;
     // check third barycentric coordinate
     if (lambda+mue > 1.0) return false;
     // have a valid hitpoint here. store it.
     ray.dist = f;
     ray.u = lambda;
     ray.v = mue;
     ray.hit_object = (TriMesh*)acc.mesh_id;
     ray.hit_face_id = acc.tri_id;
     ray.has_hit=true;
     return true;
   }
   bool BSPTree::intersect_node(Ray &ray, const BSPNode &node, double t_min, double t_max) const
   {
 //    cout << node.plane << endl;
     node_calls++;
     static bool found;
     static int i;
     if (node.axis_leaf==4)
     {
       found = false;
       for(i=0;i<node.count;i++)
       {
-//	const TriAccel* tri2 = all_triaccel[node.id+i];
+//  const TriAccel* tri2 = all_triaccel[node.id+i];
-	const ISectTri* tri = all_objects[node.id+i];
+    const ISectTri* tri = all_objects[node.id+i];
-//			if (intersect2(ray, *tri2, t_max))
+//          if (intersect2(ray, *tri2, t_max))
-//				found=true;
+//              found=true;
-	if (intersect(ray, *tri, t_max))
+    if (intersect(ray, *tri, t_max))
-	  found=true;
+      found=true;
       }
       if (found)
-	return true;
+    return true;
       else
-	return false;
+    return false;
     }
     else
     {
       BSPNode *near_node;
       BSPNode *far_node;
       if (ray.direction[node.axis_leaf]>=0)
       {
-	near_node = node.left;
+    near_node = node.left;
-	far_node = node.right;
+    far_node = node.right;
       }
       else
       {
-	near_node = node.right;
+    near_node = node.right;
-	far_node = node.left;
+    far_node = node.left;
       }
       // In order to avoid instability
       double t;
       if (fabs(ray.direction[node.axis_leaf])<d_eps)
-	t = (node.plane - ray.origin[node.axis_leaf])/d_eps;// intersect node plane;
+    t = (node.plane - ray.origin[node.axis_leaf])/d_eps;// intersect node plane;
       else
-	t = (node.plane - ray.origin[node.axis_leaf])/ray.direction[node.axis_leaf];// intersect node plane;
+    t = (node.plane - ray.origin[node.axis_leaf])/ray.direction[node.axis_leaf];// intersect node plane;
       if (t>t_max)
-      	return intersect_node(ray, *near_node, t_min, t_max);
+        return intersect_node(ray, *near_node, t_min, t_max);
       else if (t<t_min)
-	return intersect_node(ray, *far_node, t_min, t_max);
+    return intersect_node(ray, *far_node, t_min, t_max);
       else
       {
-	if (intersect_node(ray, *near_node, t_min, t))
+    if (intersect_node(ray, *near_node, t_min, t))
-	  return true;
+      return true;
-	else
+    else
-	  return intersect_node(ray, *far_node, t, t_max);
+      return intersect_node(ray, *far_node, t, t_max);
       }
     }
   }
   bool BSPTree::intersect(Ray &ray) const
   {
     double t_min, t_max;
     bbox.intersect_min_max(ray, t_min, t_max);
     if (t_min>t_max)
       return false;
-//    if (!intersect_node(ray, *root, t_min, t_max))
-//      return false;
+    printf("%.2f %.2f %.2f\n", ray.dist, t_min, t_max);
-//	cout << "____" << endl;
+//    if (!intersect_node(ray, *root, t_min, t_max))
-//	ray.reset();
+//      return false;
-//	cout << "Here " << endl;
+//  cout << "____" << endl;
-    intersect_fast_node(ray, &fast_tree[0], t_min, t_max);
+//  ray.reset();
-    if (!ray.has_hit)
+//  cout << "Here " << endl;
-      return false;
+    intersect_fast_node(ray, &fast_tree[0], t_min, t_max);
-    else
+    if (!ray.has_hit)
-    {
+      return false;
-      // Calculate the normal at the intersection
+    else
-		ray.id = (int)ray.hit_object;
+    {
-		ray.hit_object = trimesh[ray.id];
+        printf("HIT\n");
-		Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);
-		Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);
-		Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);
+      // Calculate the normal at the intersection
-		Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);
+        ray.id = (int)ray.hit_object;
-		ray.hit_normal = transforms[ray.id].mul_3D_vector(
+        ray.hit_object = trimesh[ray.id];
-			normalize(normal0*(1 - ray.u - ray.v)
-			+normal1*ray.u
+        Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);
-			+normal2*ray.v));
+        Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);
+        Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);
-		ray.hit_pos = ray.origin + ray.direction*ray.dist;
+        Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);
-/*
+        ray.hit_normal = transforms[ray.id].mul_3D_vector(
-      Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);
+            normalize(normal0*(1 - ray.u - ray.v)
-      Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);
+            +normal1*ray.u
-      Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);
+            +normal2*ray.v));
-      Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);
-      ray.hit_normal = normalize(normal0*(1 - ray.u - ray.v)
+        ray.hit_pos = ray.origin + ray.direction*ray.dist;
-				 +normal1*ray.u
+/*
-				 +normal2*ray.v);
+      Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);
-      ray.hit_pos = ray.origin + ray.direction*ray.dist;
+      Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);
-*/
+      Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);
+      Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);
-	  return true;
+      ray.hit_normal = normalize(normal0*(1 - ray.u - ray.v)
-    }
+                 +normal1*ray.u
-  }
+                 +normal2*ray.v);
+      ray.hit_pos = ray.origin + ray.direction*ray.dist;
-  const int MAX_DEPTH=25;
+*/
-  void BSPTree::make_fast_tree(BSPNode *node)
+      return true;
-  {
+    }
-    fast_tree.resize(1000000); //
+  }
-    FastBSPNode f_node;
-    fast_tree.push_back(f_node);
+  const int MAX_DEPTH=25;
-    push_fast_bsp_node(node,0);
-  }
+  void BSPTree::make_fast_tree(BSPNode *node)
+  {
-  void BSPTree::push_fast_bsp_node(BSPNode *node, int id)
+    fast_tree.resize(1000000); //
-  {
+    FastBSPNode f_node;
-    if (node->axis_leaf==4)  // It is a leaf
+    fast_tree.push_back(f_node);
-    {
+    push_fast_bsp_node(node,0);
-      fast_tree[id].leaf.flagAndOffset = (unsigned int)1<<31 | (unsigned int)(&all_triaccel[node->id]);
+  }
-      fast_tree[id].leaf.count = node->count;
-    }
+  void BSPTree::push_fast_bsp_node(BSPNode *node, int id)
-    else // It is an inner node
+  {
-    {
+    if (node->axis_leaf==4)  // It is a leaf
-      FastBSPNode fnode;
+    {
-      int p_l = fast_tree.size();
+      fast_tree[id].leaf.flagAndOffset = (unsigned int)1<<31 | (unsigned int)(&all_triaccel[node->id]);
-      fast_tree.push_back(fnode); // left
+      fast_tree[id].leaf.count = node->count;
-      fast_tree.push_back(fnode); // right
+    }
-      push_fast_bsp_node(node->left, p_l);
+    else // It is an inner node
-      push_fast_bsp_node(node->right, p_l+1);
+    {
-      fast_tree[id].inner.flagAndOffset = (unsigned int) &fast_tree[p_l] | node->axis_leaf;
+      FastBSPNode fnode;
-      fast_tree[id].inner.splitCoordinate = node->plane;
+      int p_l = fast_tree.size();
-      node->ref = fast_tree[id].inner.flagAndOffset;
+      fast_tree.push_back(fnode); // left
-    }
+      fast_tree.push_back(fnode); // right
-  }
+      push_fast_bsp_node(node->left, p_l);
+      push_fast_bsp_node(node->right, p_l+1);
-#define ABSP_ISLEAF(n)       (n->inner.flagAndOffset & (unsigned int)1<<31)
+      fast_tree[id].inner.flagAndOffset = (unsigned int) &fast_tree[p_l] | node->axis_leaf;
-#define ABSP_DIMENSION(n)    (n->inner.flagAndOffset & 0x3)
+      fast_tree[id].inner.splitCoordinate = node->plane;
-#define ABSP_OFFSET(n)       (n->inner.flagAndOffset & (0x7FFFFFFC))
+      node->ref = fast_tree[id].inner.flagAndOffset;
-#define ABSP_NEARNODE(n)     (FastBSPNode*)(ray.direction[dimension]>=0?ABSP_OFFSET(node):ABSP_OFFSET(node)+sizeof(*node))
+    }
-#define ABSP_FARNODE(n)      (FastBSPNode*)(ray.direction[dimension]>=0?ABSP_OFFSET(node)+sizeof(*node):ABSP_OFFSET(node))
+  }
-#define ABSP_TRIANGLENODE(n) (vector<TriAccel*>::iterator)(n->flagAndOffset & (0x7FFFFFFF))
+#define ABSP_ISLEAF(n)       (n->inner.flagAndOffset & (unsigned int)1<<31)
-  struct Stack
+#define ABSP_DIMENSION(n)    (n->inner.flagAndOffset & 0x3)
-  {
+#define ABSP_OFFSET(n)       (n->inner.flagAndOffset & (0x7FFFFFFC))
-    FastBSPNode *node;
+#define ABSP_NEARNODE(n)     (FastBSPNode*)(ray.direction[dimension]>=0?ABSP_OFFSET(node):ABSP_OFFSET(node)+sizeof(*node))
-    double t_min;
+#define ABSP_FARNODE(n)      (FastBSPNode*)(ray.direction[dimension]>=0?ABSP_OFFSET(node)+sizeof(*node):ABSP_OFFSET(node))
-    double t_max;
+#define ABSP_TRIANGLENODE(n) (vector<TriAccel*>::iterator)(n->flagAndOffset & (0x7FFFFFFF))
-  };
+  struct Stack
-  inline void IntersectAlltrianglesInLeaf(const BSPLeaf* leaf, Ray &ray, double t_max) {
+  {
-    TriAccel** tri_acc_ptr = reinterpret_cast<TriAccel**>(leaf->flagAndOffset & (0x7FFFFFFF));
+    FastBSPNode *node;
-    vector<TriAccel*>::iterator acc = vector<TriAccel*>::iterator(tri_acc_ptr);
+    double t_min;
-//	vector<TriAccel*>::iterator acc = ABSP_TRIANGLENODE(leaf);
+    double t_max;
-    for(unsigned int i=0;i<leaf->count;++i)
+  };
-      intersect2(ray, *(*acc++), t_max);
-  }
+  inline void IntersectAlltrianglesInLeaf(const BSPLeaf* leaf, Ray &ray, double t_max) {
+    TriAccel** tri_acc_ptr = reinterpret_cast<TriAccel**>(leaf->flagAndOffset & (0x7FFFFFFF));
-  void BSPTree::intersect_fast_node(Ray &ray, const FastBSPNode *node, double t_min, double t_max) const
+    vector<TriAccel*>::iterator acc = vector<TriAccel*>::iterator(tri_acc_ptr);
-  {
+//  vector<TriAccel*>::iterator acc = ABSP_TRIANGLENODE(leaf);
-    Stack stack[MAX_DEPTH];
+    for(unsigned int i=0;i<leaf->count;++i)
-    int stack_id=0;
+      intersect2(ray, *(*acc++), t_max);
-    double t;
+  }
-    // Precalculate one over dir
-    double one_over_dir[3];
+  void BSPTree::intersect_fast_node(Ray &ray, const FastBSPNode *node, double t_min, double t_max) const
-    for(int i=0;i<3;i++)
+  {
-    {
+    Stack stack[MAX_DEPTH];
-      if (ray.direction[i]!=0)
+    int stack_id=0;
-	one_over_dir[i]=1.0/ray.direction[i];
+    double t;
-      else
+    // Precalculate one over dir
-	one_over_dir[i]=1.0/d_eps;
+    double one_over_dir[3];
-    }
+    for(int i=0;i<3;i++)
+    {
-    int dimension;
+      if (ray.direction[i]!=0)
-    while(1)
+    one_over_dir[i]=1.0/ray.direction[i];
-    {
+      else
-      while(!ABSP_ISLEAF(node))
+    one_over_dir[i]=1.0/d_eps;
-      {
+    }
-	dimension = ABSP_DIMENSION(node);
-	t = (node->inner.splitCoordinate - ray.origin[dimension])*one_over_dir[dimension];
+    int dimension;
-	if (t>=t_max)
+    while(1)
-	  node = ABSP_NEARNODE(node);
+    {
-	else if (t<=t_min)
+      while(!ABSP_ISLEAF(node))
-	  node = ABSP_FARNODE(node);
+      {
-	else
+    dimension = ABSP_DIMENSION(node);
-	{
+    t = (node->inner.splitCoordinate - ray.origin[dimension])*one_over_dir[dimension];
-	  // Stack push
+    if (t>=t_max)
-	  stack[stack_id].node = ABSP_FARNODE(node);
+      node = ABSP_NEARNODE(node);
-	  stack[stack_id].t_min = t;
+    else if (t<=t_min)
-	  stack[stack_id++].t_max = t_max;
+      node = ABSP_FARNODE(node);
-	  // Set current node to near side
+    else
-	  node = ABSP_NEARNODE(node);
+    {
-	  t_max = t;
+      // Stack push
-	}
+      stack[stack_id].node = ABSP_FARNODE(node);
-      }
+      stack[stack_id].t_min = t;
+      stack[stack_id++].t_max = t_max;
-      IntersectAlltrianglesInLeaf(&node->leaf, ray, t_max);
+      // Set current node to near side
-      if (ray.dist<t_max)
+      node = ABSP_NEARNODE(node);
-	return;
+      t_max = t;
-      if (stack_id==0)
+    }
-	return;
+      }
-      // Stack pop
+      IntersectAlltrianglesInLeaf(&node->leaf, ray, t_max);
-      node = stack[--stack_id].node;
+      if (ray.dist<t_max)
-      t_min = stack[stack_id].t_min;
+    return;
-      t_max = stack[stack_id].t_max;
+      if (stack_id==0)
-    }
+    return;
-  }
+      // Stack pop
-  bool BSPTree::intersect(Ray &ray, const ISectTri &isecttri, double t_max) const
+      node = stack[--stack_id].node;
-  {
+      t_min = stack[stack_id].t_min;
-    tri_calls++;
+      t_max = stack[stack_id].t_max;
-    // This is the Möller-Trumbore method
+    }
-    static Vec3d direction; // = Vec3d(ray.direction);
+  }
-    static Vec3d origin; // = Vec3d(ray.direction);
-    static Vec3d edge1; // = Vec3d(ray.direction);
+  bool BSPTree::intersect(Ray &ray, const ISectTri &isecttri, double t_max) const
-    static Vec3d edge0; // = Vec3d(ray.direction);
+  {
-    static Vec3d point0; // = Vec3d(ray.direction);
+    tri_calls++;
-    static Vec3d p;
+    // This is the Möller-Trumbore method
-    static Vec3d q;
+    static Vec3d direction; // = Vec3d(ray.direction);
-    static Vec3d s;
+    static Vec3d origin; // = Vec3d(ray.direction);
-    static Vec3f point;
+    static Vec3d edge1; // = Vec3d(ray.direction);
-    static double a, f, u, v, t;
+    static Vec3d edge0; // = Vec3d(ray.direction);
-    static double ray_dist_sq;
+    static Vec3d point0; // = Vec3d(ray.direction);
+    static Vec3d p;
-    static double dist_sq;
+    static Vec3d q;
+    static Vec3d s;
-    static Vec3f found_point;
+    static Vec3f point;
+    static double a, f, u, v, t;
-    direction.set((double)ray.direction[0], (double)ray.direction[1], (double)ray.direction[2]);
+    static double ray_dist_sq;
-    edge0.set((double)isecttri.edge0[0], (double)isecttri.edge0[1], (double)isecttri.edge0[2]);
-    edge1.set((double)isecttri.edge1[0], (double)isecttri.edge1[1], (double)isecttri.edge1[2]);
+    static double dist_sq;
-// Ray-tri intersection
+    static Vec3f found_point;
-//		const double eps = 0.001f;
-/********* Complain!!!!!!!!!!!!!!!! *****************/
+    direction.set((double)ray.direction[0], (double)ray.direction[1], (double)ray.direction[2]);
-	//p = cross(direction, edge1);
+    edge0.set((double)isecttri.edge0[0], (double)isecttri.edge0[1], (double)isecttri.edge0[2]);
-	// Why the &%¤/ is this so much faster????? - Because of MS Compiler - the intel compiler does it right!!
+    edge1.set((double)isecttri.edge1[0], (double)isecttri.edge1[1], (double)isecttri.edge1[2]);
-    p.set(direction[1] * edge1[2] - direction[2] * edge1[1],
-	  direction[2] * edge1[0] - direction[0] * edge1[2],
+// Ray-tri intersection
-	  direction[0] * edge1[1] - direction[1] * edge1[0]);
+//      const double eps = 0.001f;
-/****************************************************/
+/********* Complain!!!!!!!!!!!!!!!! *****************/
-    a = dot(edge0,p);
+    //p = cross(direction, edge1);
-    if (a>-d_eps && a<d_eps)
+    // Why the &%¤/ is this so much faster????? - Because of MS Compiler - the intel compiler does it right!!
-      return false;
+    p.set(direction[1] * edge1[2] - direction[2] * edge1[1],
-  // Just delay these
+      direction[2] * edge1[0] - direction[0] * edge1[2],
-    origin.set((double)ray.origin[0], (double)ray.origin[1], (double)ray.origin[2]);
+      direction[0] * edge1[1] - direction[1] * edge1[0]);
-    point0.set((double)isecttri.point0[0], (double)isecttri.point0[1], (double)isecttri.point0[2]);
+/****************************************************/
+    a = dot(edge0,p);
-    f = 1.0/a;
+    if (a>-d_eps && a<d_eps)
-    s = origin - point0;
+      return false;
-    u = f*dot(s,p);
+  // Just delay these
-    if (u<0.0 || u>1.0)
+    origin.set((double)ray.origin[0], (double)ray.origin[1], (double)ray.origin[2]);
-      return false;
+    point0.set((double)isecttri.point0[0], (double)isecttri.point0[1], (double)isecttri.point0[2]);
-/********* Complain!!!!!!!!!!!!!!!! *****************/
-  //q = cross(s, edge0);
+    f = 1.0/a;
-  // Why the &%¤/ is this so much faster?????
+    s = origin - point0;
-    q.set(s[1] * edge0[2] - s[2] * edge0[1],
+    u = f*dot(s,p);
-	  s[2] * edge0[0] - s[0] * edge0[2],
+    if (u<0.0 || u>1.0)
-	  s[0] * edge0[1] - s[1] * edge0[0]);
+      return false;
-/****************************************************/
+/********* Complain!!!!!!!!!!!!!!!! *****************/
+  //q = cross(s, edge0);
-    v = f * dot(direction, q);
+  // Why the &%¤/ is this so much faster?????
-    if (v<0.0 || u+v>1.0)
+    q.set(s[1] * edge0[2] - s[2] * edge0[1],
-      return false;
+      s[2] * edge0[0] - s[0] * edge0[2],
-    t = f*dot(edge1, q);
+      s[0] * edge0[1] - s[1] * edge0[0]);
-    if (t<0)
+/****************************************************/
-      return false;
-    if (fabs(t)<d_eps)
+    v = f * dot(direction, q);
-      return false;
+    if (v<0.0 || u+v>1.0)
-    if (t_max<t)
+      return false;
-      return false;
+    t = f*dot(edge1, q);
-    point = ray.origin + ray.direction*t;
+    if (t<0)
+      return false;
-    dist_sq = dot(point-ray.origin, point-ray.origin);
+    if (fabs(t)<d_eps)
-    ray_dist_sq = ray.dist * ray.dist;
+      return false;
+    if (t_max<t)
-    if (dist_sq<f_eps)
+      return false;
-      return false;
+    point = ray.origin + ray.direction*t;
-    if (dist_sq>ray_dist_sq)
-      return false;
+    dist_sq = dot(point-ray.origin, point-ray.origin);
+    ray_dist_sq = ray.dist * ray.dist;
-    ray.dist = sqrt(dist_sq);
-    ray.u = u;
+    if (dist_sq<f_eps)
-    ray.v = v;
+      return false;
-    ray.hit_object = (TriMesh*)isecttri.mesh_id;
+    if (dist_sq>ray_dist_sq)
-    ray.hit_face_id = isecttri.tri_id;
+      return false;
-    ray.has_hit=true;
-    return true;
+    ray.dist = sqrt(dist_sq);
-  }
+    ray.u = u;
-}
+    ray.v = v;
+    ray.hit_object = (TriMesh*)isecttri.mesh_id;
+    ray.hit_face_id = isecttri.tri_id;
+    ray.has_hit=true;
+    return true;
+  }
+}

Subversion Repositories gelsvn

(root)/trunk/src/Geometry/BSPTree.cpp @ 430 – Rev 319 → 349