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