WebSVN – gelsvn – Diff – /branches/SPD/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;
      double min_cost=CGLA::BIG;
      int new_pos = -1;      
/*
       float max_size=0;
       for(i=0;i<3;i++) {
         if (bbox.max_corner[i]-bbox.min_corner[i]>max_size) {
	   max_size = bbox.max_corner[i]-bbox.min_corner[i];
	   new_axis=i;
	   new_pos = 1;
	 }
       }
*/
      unsigned int j;		
 
      for(i=0;i<3;i++) 
      {
	for(int k=1;k<TESTS;k++) 
	{
	  BBox left_bbox = bbox;
	  BBox right_bbox = bbox;
				
	  double center = (bbox.max_corner[i]- bbox.min_corner[i])*(double)k/(double)TESTS+bbox.min_corner[i];
	  node.plane = center;
	  
	  left_bbox.max_corner[i] = center; 
	  right_bbox.min_corner[i] = center; 
	  // Try putting the triangles in the left and right boxes
	  int left_count = 0;
	  int right_count = 0;
	  for(j=0;j<objects.size();j++) 
	  {
	    ISectTri* tri = objects[j];
	    if (left_bbox.intersect_triangle(*tri)) 
	      left_count++;
	    if (right_bbox.intersect_triangle(*tri)) 
	      right_count++;
	  }
	  
	  if ((double)left_count*left_bbox.area()+(double)right_count*right_bbox.area()<min_cost) 
	  {
	    min_cost = (double)left_count*left_bbox.area()+(double)right_count*right_bbox.area();
	    new_axis=i;
	    new_pos = k;
	    if (right_count==0)
	      right_zero=true;
	    else
	      right_zero=false;
	    if (left_count==0)
	      left_zero=true;
	    else
	      left_zero=false;	    
	  }
	}
      }
      node.axis_leaf = new_axis;
      //	node.axis_leaf = level %3;
      left_node->axis_leaf = -1; //(node.axis + 1)%3;
      right_node->axis_leaf = -1; //(node.axis + 1)%3;
      // Now chose the right splitting plane
      BBox left_bbox = bbox;
      BBox right_bbox = bbox;
 
      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];
 
      // This doesn't help that much
      if (left_zero) 
      {
//			cout << "Old left center: " << center;
	// Find min position af alle trekanter og placer center der
	center = bbox.max_corner[node.axis_leaf];
	for(j=0;j<objects.size();j++) 
	{
	  ISectTri* tri = objects[j];
	  if (tri->point0[node.axis_leaf]<center)
	    center=tri->point0[node.axis_leaf];
	  if (tri->point1[node.axis_leaf]<center)
	    center=tri->point1[node.axis_leaf];
	  if (tri->point2[node.axis_leaf]<center)
	    center=tri->point2[node.axis_leaf];
	}
//			cout << "new left center: " << center << endl;
      }
      if (right_zero) 
      {
//			cout << "Old right center: " << center;
	// Find min position af alle trekanter og placer center der
	center = left_bbox.min_corner[node.axis_leaf];
	for(j=0;j<objects.size();j++) 
	{
	  ISectTri* tri = objects[j];
	  if (tri->point0[node.axis_leaf]>center)
	    center=tri->point0[node.axis_leaf];
	  if (tri->point1[node.axis_leaf]>center)
	    center=tri->point1[node.axis_leaf];
	  if (tri->point2[node.axis_leaf]>center)
	    center=tri->point2[node.axis_leaf];
	}
//			cout << "new right center: " << center << endl;
      }
 
      node.plane = center;
      left_bbox.max_corner[node.axis_leaf] = center+f_eps; 
      right_bbox.min_corner[node.axis_leaf] = center-f_eps; 
      
 
		// Now put the triangles in the right and left node
      for(i=0;i<objects.size();i++) 
      {
	ISectTri* tri = objects[i];
	TriAccel *tri_accel = tri_objects[i];
	if (left_bbox.intersect_triangle(*tri)) 
	{
	  left_objects.push_back(tri);
	  tri_left_objects.push_back(tri_accel);
	}
	if (right_bbox.intersect_triangle(*tri)) 
	{
	  right_objects.push_back(tri);
	  tri_right_objects.push_back(tri_accel);
	}
      }
/*	if (left_zero||right_zero) {
			cout << left_objects.size() << "," << right_objects.size() << 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(f_eps*100,f_eps*100,f_eps*100);
    bbox.max_corner+=Vec3f(f_eps*100,f_eps*100,f_eps*100);
  }
 
  void BSPTree::init(vector<TriMesh*>& _trimesh, 
		     vector<Mat4x4f>& _transforms, 
		     int _max_objects, int _max_level) 
  {
    trimesh = _trimesh;
    transforms = _transforms;
    for(unsigned int i=0;i<trimesh.size();i++) 
    {
      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);
	ISectTri new_tri;
	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.point2 = transforms[i].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 = i;
	new_tri.tri_id = j;
	isecttris.push_back(new_tri);
	TriAccel ta;
	create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);
	ta.mesh_id = i;
	ta.tri_id = j;
	triaccel.push_back(ta);
      }
    }
 
    max_objects = _max_objects;
    max_level = _max_level;
    init();
  }
 
  void BSPTree::init(TriMesh* mesh, Mat4x4f transform, 
		     vector<int> &trilist, 
		     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];
	TriAccel& tri_accel = triaccel[i];
	objects.push_back(&tri);
	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;
    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 << " ";
//inline bool Intersect(TriAccel &acc,Ray &ray)
#define ku modulo[acc.k+1]
#define kv modulo[acc.k+2]
    // dont 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 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 > f_eps)||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 ISectTri* tri = all_objects[node.id+i];
//			if (intersect2(ray, *tri2, t_max))  
//				found=true;
	if (intersect(ray, *tri, t_max))  
	  found=true;
      }
      if (found)
	return true;
      else 
	return false;
    } 
    else 
    {
      BSPNode *near_node;
      BSPNode *far_node;
      if (ray.direction[node.axis_leaf]>=0) 
      {
	near_node = node.left;
	far_node = node.right;
      } 
      else 
      {
	near_node = node.right;
	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;
      else
	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);
      
      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))
	  return true;
	else 
	  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;
//	cout << "____" << endl;
//	ray.reset();
//	cout << "Here " << endl;
//	intersect_fast_node(ray, &fast_tree[0], t_min, t_max);
//	if (!ray.has_hit)
//		return false;
    else 
    {
      // Calculate the normal at the intersection
      ray.hit_object = trimesh[(int)ray.hit_object];
      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]);
      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);
      ray.hit_pos = ray.origin + ray.direction*ray.dist;
      return true;
    }
  }
 
  const int MAX_DEPTH=25;
 
  void BSPTree::make_fast_tree(BSPNode *node) 
  {
    fast_tree.resize(1000000); // 
    FastBSPNode f_node;
    fast_tree.push_back(f_node);
    push_fast_bsp_node(node,0);
  }
 
  void BSPTree::push_fast_bsp_node(BSPNode *node, int id) 
  {
    if (node->axis_leaf==4)  // It is a leaf
    {
      fast_tree[id].leaf.flagAndOffset = (unsigned int)1<<31 | (unsigned int)(&all_triaccel[node->id]);
      fast_tree[id].leaf.count = node->count;
    } 
    else // It is an inner node
    { 
      FastBSPNode fnode;
      int p_l = fast_tree.size();
      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);
      fast_tree[id].inner.flagAndOffset = (unsigned int) &fast_tree[p_l] | node->axis_leaf;
      fast_tree[id].inner.splitCoordinate = node->plane;
      node->ref = fast_tree[id].inner.flagAndOffset;
    }
  }
 
#define ABSP_ISLEAF(n)       (n->inner.flagAndOffset & (unsigned int)1<<31)
#define ABSP_DIMENSION(n)    (n->inner.flagAndOffset & 0x3)
#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))
  
  struct Stack 
  {
    FastBSPNode *node;
    double t_min;
    double 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 = vector<TriAccel*>::iterator(tri_acc_ptr);
//	vector<TriAccel*>::iterator acc = ABSP_TRIANGLENODE(leaf);
    for(unsigned int i=0;i<leaf->count;++i)
      intersect2(ray, *(*acc++), t_max);
  }
 
  void BSPTree::intersect_fast_node(Ray &ray, const FastBSPNode *node, double t_min, double t_max) const 
  {
    Stack stack[MAX_DEPTH];
    int stack_id=0;
    double t;
    // Precalculate one over dir
    double one_over_dir[3];
    for(int i=0;i<3;i++) 
    {
      if (ray.direction[i]!=0)
	one_over_dir[i]=1.0/ray.direction[i];
      else
	one_over_dir[i]=1.0/d_eps;
    }
 
    int dimension;
    while(1) 
    {
      while(!ABSP_ISLEAF(node)) 
      {
	dimension = ABSP_DIMENSION(node);
	t = (node->inner.splitCoordinate - ray.origin[dimension])*one_over_dir[dimension];
	if (t>=t_max) 
	  node = ABSP_NEARNODE(node);
	else if (t<=t_min)
	  node = ABSP_FARNODE(node);
	else 
	{
	  // Stack push
	  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
	  node = ABSP_NEARNODE(node);
	  t_max = t;
	}
      }
      
      IntersectAlltrianglesInLeaf(&node->leaf, ray, t_max);
      if (ray.dist<t_max)
	return;
      if (stack_id==0)
	return;
      // Stack pop
      
      node = stack[--stack_id].node;
      t_min = stack[stack_id].t_min;
      t_max = stack[stack_id].t_max;
    }
  }
 
  bool BSPTree::intersect(Ray &ray, const ISectTri &isecttri, double t_max) const 
  {
    tri_calls++;
    // 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);
    static Vec3d edge0; // = Vec3d(ray.direction);
    static Vec3d point0; // = Vec3d(ray.direction);
    static Vec3d p;
    static Vec3d q;
    static Vec3d s;
    static Vec3f point;
    static double a, f, u, v, t;
    static double ray_dist_sq;
    
    static double dist_sq;
    
    static Vec3f found_point;
    
    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]);
    edge1.set((double)isecttri.edge1[0], (double)isecttri.edge1[1], (double)isecttri.edge1[2]);
 
// Ray-tri intersection
//		const double eps = 0.001f;
/********* Complain!!!!!!!!!!!!!!!! *****************/		
	//p = cross(direction, edge1);
	// 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], 
	  direction[2] * edge1[0] - direction[0] * edge1[2], 
	  direction[0] * edge1[1] - direction[1] * edge1[0]);
/****************************************************/
    a = dot(edge0,p);
    if (a>-d_eps && a<d_eps)
      return false;
  // Just delay these 
    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]);
    
    f = 1.0/a;
    s = origin - point0;
    u = f*dot(s,p);
    if (u<0.0 || u>1.0)
      return false;
/********* Complain!!!!!!!!!!!!!!!! *****************/		
  //q = cross(s, edge0);
  // Why the &%¤/ is this so much faster?????
    q.set(s[1] * edge0[2] - s[2] * edge0[1], 
	  s[2] * edge0[0] - s[0] * edge0[2], 
	  s[0] * edge0[1] - s[1] * edge0[0]);
/****************************************************/
	
    v = f * dot(direction, q);  
    if (v<0.0 || u+v>1.0)
      return false;
    t = f*dot(edge1, q);
    if (t<0)
      return false;
    if (fabs(t)<d_eps)
      return false;
    if (t_max<t)
      return false;
    point = ray.origin + ray.direction*t;
    
    dist_sq = dot(point-ray.origin, point-ray.origin);
    ray_dist_sq = ray.dist * ray.dist;
    
    if (dist_sq<f_eps)
      return false;
    if (dist_sq>ray_dist_sq)
      return false;
  
    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; 
  }
}
 

Rev 208	Rev 290
1	`#include "CGLA/Vec3d.h"`	1	`#include "CGLA/Vec3d.h"`
2	`#include "Geometry/TriMesh.h"`	2	`#include "Geometry/TriMesh.h"`
3		3
4	`#include "BSPTree.h"`	4	`#include "BSPTree.h"`
5		5
6	`using namespace std;`	6	`using namespace std;`
7	`using namespace CGLA;`	7	`using namespace CGLA;`
8	`using namespace Geometry;`	-
9	`using namespace BRender;`	-
10		8
11	`namespace TriMeshTrace`	9	`namespace Geometry`
12	`{`	10	`{`
13	`int BSPTree::node_calls;`	11	`int BSPTree::node_calls;`
14	`int BSPTree::tri_calls;`	12	`int BSPTree::tri_calls;`
15		13
16	`void create_tri_accel(Vec3f A, Vec3f B, Vec3f C, TriAccel &tri_accel)`	14	`void create_tri_accel(Vec3f A, Vec3f B, Vec3f C, TriAccel &tri_accel)`
17	`{`	15	`{`
18	`Vec3f N = normalize(cross(B-A, C-A));`	16	`Vec3f N = normalize(cross(B-A, C-A));`
19	`// Find projection dir`	17	`// Find projection dir`
20	`int k,u,v;`	18	`int k,u,v;`
21	`if (fabs(N[0]) > fabs(N[1]))`	19	`if (fabs(N[0]) > fabs(N[1]))`
22	`if (fabs(N[0]) > fabs(N[2]))`	20	`if (fabs(N[0]) > fabs(N[2]))`
23	`k = 0; /* X */`	21	`k = 0; /* X */`
24	`else`	22	`else`
25	`k=2; /* Z */`	23	`k=2; /* Z */`
26	`else`	24	`else`
27	`if (fabs(N[1]) > fabs(N[2]))`	25	`if (fabs(N[1]) > fabs(N[2]))`
28	`k = 1; /* Y */`	26	`k = 1; /* Y */`
29	`else`	27	`else`
30	`k=2; /* Z */`	28	`k=2; /* Z */`
31		29
32	`u = (k+1)% 3;`	30	`u = (k+1)% 3;`
33	`v = (k+2)% 3;`	31	`v = (k+2)% 3;`
34		32
35	`Vec3f b = C - A;`	33	`Vec3f b = C - A;`
36	`Vec3f c = B - A;`	34	`Vec3f c = B - A;`
37		35
38	`double div = (b[u]c[v]-b[v]c[u]);`	36	`double div = (b[u]c[v]-b[v]c[u]);`
39		37
40	`tri_accel.n_u = N[u]/N[k];`	38	`tri_accel.n_u = N[u]/N[k];`
41	`tri_accel.n_v = N[v]/N[k];`	39	`tri_accel.n_v = N[v]/N[k];`
42	`tri_accel.n_d = dot(A, N/N[k]);`	40	`tri_accel.n_d = dot(A, N/N[k]);`
43	`tri_accel.k = k;`	41	`tri_accel.k = k;`
44		42
45	`tri_accel.b_nu = -b[v]/div;`	43	`tri_accel.b_nu = -b[v]/div;`
46	`tri_accel.b_nv = b[u]/div;`	44	`tri_accel.b_nv = b[u]/div;`
47	`tri_accel.b_d = (b[v]A[u]-b[u]A[v])/div;`	45	`tri_accel.b_d = (b[v]A[u]-b[u]A[v])/div;`
48		46
49	`tri_accel.c_nu = c[v]/div;`	47	`tri_accel.c_nu = c[v]/div;`
50	`tri_accel.c_nv = -c[u]/div;`	48	`tri_accel.c_nv = -c[u]/div;`
51	`tri_accel.c_d = (c[u]A[v]-c[v]A[u])/div;`	49	`tri_accel.c_d = (c[u]A[v]-c[v]A[u])/div;`
52	`}`	50	`}`
53		51
54	`// Most of this file is a direct copy from Henrik's notes`	52	`// Most of this file is a direct copy from Henrik's notes`
55	`BSPTree::BSPTree()`	53	`BSPTree::BSPTree()`
56	`{`	54	`{`
57	`root=0;`	55	`root=0;`
58	`b_is_build = false;`	56	`b_is_build = false;`
59	`}`	57	`}`
60		58
61	`BSPTree::~BSPTree()`	59	`BSPTree::~BSPTree()`
62	`{`	60	`{`
63	`clear();`	61	`clear();`
64	`}`	62	`}`
65		63
66	`void BSPTree::clear()`	64	`void BSPTree::clear()`
67	`{`	65	`{`
68	`if (root!=0)`	66	`if (root!=0)`
69	`{`	67	`{`
70	`delete_node(root);`	68	`delete_node(root);`
71	`root=0;`	69	`root=0;`
72	`b_is_build=false;`	70	`b_is_build=false;`
73	`}`	71	`}`
74	`isecttris.clear();`	72	`isecttris.clear();`
75	`all_objects.clear();`	73	`all_objects.clear();`
76	`all_triaccel.clear();`	74	`all_triaccel.clear();`
77	`}`	75	`}`
78		76
79	`void BSPTree::delete_node(BSPNode *node)`	77	`void BSPTree::delete_node(BSPNode *node)`
80	`{`	78	`{`
81	`if (node->left!=0)`	79	`if (node->left!=0)`
82	`delete_node(node->left);`	80	`delete_node(node->left);`
83	`if (node->right!=0)`	81	`if (node->right!=0)`
84	`delete_node(node->right);`	82	`delete_node(node->right);`
85	`delete node;`	83	`delete node;`
86	`}`	84	`}`
87		85
88	`void BSPTree::subdivide_node(BSPNode &node, BBox &bbox,`	86	`void BSPTree::subdivide_node(BSPNode &node, BBox &bbox,`
89	`unsigned int level,`	87	`unsigned int level,`
90	`vector<ISectTri*>& objects,`	88	`vector<ISectTri*>& objects,`
91	`vector<TriAccel*>& tri_objects)`	89	`vector<TriAccel*>& tri_objects)`
92	`{`	90	`{`
93	`const int TESTS = 2;`	91	`const int TESTS = 2;`
94		92
95	`if (objects.size()<=max_objects \|\| level==max_level)`	93	`if (objects.size()<=max_objects \|\| level==max_level)`
96	`{`	94	`{`
97	`node.axis_leaf = 4; // Means that this is a leaf`	95	`node.axis_leaf = 4; // Means that this is a leaf`
98	`node.plane = 0;`	96	`node.plane = 0;`
99	`node.left = 0;`	97	`node.left = 0;`
100	`node.right = 0;`	98	`node.right = 0;`
101	`node.id = all_objects.size();`	99	`node.id = all_objects.size();`
102	`node.count = objects.size();`	100	`node.count = objects.size();`
103		101
104	`for(unsigned int i = 0; i < objects.size(); ++i)`	102	`for(unsigned int i = 0; i < objects.size(); ++i)`
105	`{`	103	`{`
106	`all_objects.push_back(objects[i]);`	104	`all_objects.push_back(objects[i]);`
107	`all_triaccel.push_back(tri_objects[i]);`	105	`all_triaccel.push_back(tri_objects[i]);`
108	`}`	106	`}`
109	`}`	107	`}`
110	`else`	108	`else`
111	`{`	109	`{`
112	`bool right_zero=false;`	110	`bool right_zero=false;`
113	`bool left_zero=false;`	111	`bool left_zero=false;`
114	`unsigned int i;`	112	`unsigned int i;`
115	`BSPNode* left_node = new BSPNode();`	113	`BSPNode* left_node = new BSPNode();`
116	`BSPNode* right_node = new BSPNode();`	114	`BSPNode* right_node = new BSPNode();`
117	`vector<ISectTri*> left_objects;`	115	`vector<ISectTri*> left_objects;`
118	`vector<ISectTri*> right_objects;`	116	`vector<ISectTri*> right_objects;`
119	`vector<TriAccel*> tri_left_objects;`	117	`vector<TriAccel*> tri_left_objects;`
120	`vector<TriAccel*> tri_right_objects;`	118	`vector<TriAccel*> tri_right_objects;`
121		119
122	`node.left = left_node;`	120	`node.left = left_node;`
123	`node.right = right_node;`	121	`node.right = right_node;`
124	`//node.is_leaf = false;`	122	`//node.is_leaf = false;`
125		123
126	`// Make test on all three axis`	124	`// Make test on all three axis`
127	`int new_axis=-1;`	125	`int new_axis=-1;`
128	`double min_cost=CGLA::BIG;`	126	`double min_cost=CGLA::BIG;`
129	`int new_pos = -1;`	127	`int new_pos = -1;`
130	`/*`	128	`/*`
131	`float max_size=0;`	129	`float max_size=0;`
132	`for(i=0;i<3;i++) {`	130	`for(i=0;i<3;i++) {`
133	`if (bbox.max_corner[i]-bbox.min_corner[i]>max_size) {`	131	`if (bbox.max_corner[i]-bbox.min_corner[i]>max_size) {`
134	`max_size = bbox.max_corner[i]-bbox.min_corner[i];`	132	`max_size = bbox.max_corner[i]-bbox.min_corner[i];`
135	`new_axis=i;`	133	`new_axis=i;`
136	`new_pos = 1;`	134	`new_pos = 1;`
137	`}`	135	`}`
138	`}`	136	`}`
139	`*/`	137	`*/`
140	`unsigned int j;`	138	`unsigned int j;`
141		139
142	`for(i=0;i<3;i++)`	140	`for(i=0;i<3;i++)`
143	`{`	141	`{`
144	`for(int k=1;k<TESTS;k++)`	142	`for(int k=1;k<TESTS;k++)`
145	`{`	143	`{`
146	`BBox left_bbox = bbox;`	144	`BBox left_bbox = bbox;`
147	`BBox right_bbox = bbox;`	145	`BBox right_bbox = bbox;`
148		146
149	`double center = (bbox.max_corner[i]- bbox.min_corner[i])*(double)k/(double)TESTS+bbox.min_corner[i];`	147	`double center = (bbox.max_corner[i]- bbox.min_corner[i])*(double)k/(double)TESTS+bbox.min_corner[i];`
150	`node.plane = center;`	148	`node.plane = center;`
151		149
152	`left_bbox.max_corner[i] = center;`	150	`left_bbox.max_corner[i] = center;`
153	`right_bbox.min_corner[i] = center;`	151	`right_bbox.min_corner[i] = center;`
154	`// Try putting the triangles in the left and right boxes`	152	`// Try putting the triangles in the left and right boxes`
155	`int left_count = 0;`	153	`int left_count = 0;`
156	`int right_count = 0;`	154	`int right_count = 0;`
157	`for(j=0;j<objects.size();j++)`	155	`for(j=0;j<objects.size();j++)`
158	`{`	156	`{`
159	`ISectTri* tri = objects[j];`	157	`ISectTri* tri = objects[j];`
160	`if (left_bbox.intersect_triangle(*tri))`	158	`if (left_bbox.intersect_triangle(*tri))`
161	`left_count++;`	159	`left_count++;`
162	`if (right_bbox.intersect_triangle(*tri))`	160	`if (right_bbox.intersect_triangle(*tri))`
163	`right_count++;`	161	`right_count++;`
164	`}`	162	`}`
165		163
166	`if ((double)left_countleft_bbox.area()+(double)right_countright_bbox.area()<min_cost)`	164	`if ((double)left_countleft_bbox.area()+(double)right_countright_bbox.area()<min_cost)`
167	`{`	165	`{`
168	`min_cost = (double)left_countleft_bbox.area()+(double)right_countright_bbox.area();`	166	`min_cost = (double)left_countleft_bbox.area()+(double)right_countright_bbox.area();`
169	`new_axis=i;`	167	`new_axis=i;`
170	`new_pos = k;`	168	`new_pos = k;`
171	`if (right_count==0)`	169	`if (right_count==0)`
172	`right_zero=true;`	170	`right_zero=true;`
173	`else`	171	`else`
174	`right_zero=false;`	172	`right_zero=false;`
175	`if (left_count==0)`	173	`if (left_count==0)`
176	`left_zero=true;`	174	`left_zero=true;`
177	`else`	175	`else`
178	`left_zero=false;`	176	`left_zero=false;`
179	`}`	177	`}`
180	`}`	178	`}`
181	`}`	179	`}`
182	`node.axis_leaf = new_axis;`	180	`node.axis_leaf = new_axis;`
183	`// node.axis_leaf = level %3;`	181	`// node.axis_leaf = level %3;`
184	`left_node->axis_leaf = -1; //(node.axis + 1)%3;`	182	`left_node->axis_leaf = -1; //(node.axis + 1)%3;`
185	`right_node->axis_leaf = -1; //(node.axis + 1)%3;`	183	`right_node->axis_leaf = -1; //(node.axis + 1)%3;`
186	`// Now chose the right splitting plane`	184	`// Now chose the right splitting plane`
187	`BBox left_bbox = bbox;`	185	`BBox left_bbox = bbox;`
188	`BBox right_bbox = bbox;`	186	`BBox right_bbox = bbox;`
189		187
190	`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];`	188	`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];`
191		189
192	`// This doesn't help that much`	190	`// This doesn't help that much`
193	`if (left_zero)`	191	`if (left_zero)`
194	`{`	192	`{`
195	`// cout << "Old left center: " << center;`	193	`// cout << "Old left center: " << center;`
196	`// Find min position af alle trekanter og placer center der`	194	`// Find min position af alle trekanter og placer center der`
197	`center = bbox.max_corner[node.axis_leaf];`	195	`center = bbox.max_corner[node.axis_leaf];`
198	`for(j=0;j<objects.size();j++)`	196	`for(j=0;j<objects.size();j++)`
199	`{`	197	`{`
200	`ISectTri* tri = objects[j];`	198	`ISectTri* tri = objects[j];`
201	`if (tri->point0[node.axis_leaf]<center)`	199	`if (tri->point0[node.axis_leaf]<center)`
202	`center=tri->point0[node.axis_leaf];`	200	`center=tri->point0[node.axis_leaf];`
203	`if (tri->point1[node.axis_leaf]<center)`	201	`if (tri->point1[node.axis_leaf]<center)`
204	`center=tri->point1[node.axis_leaf];`	202	`center=tri->point1[node.axis_leaf];`
205	`if (tri->point2[node.axis_leaf]<center)`	203	`if (tri->point2[node.axis_leaf]<center)`
206	`center=tri->point2[node.axis_leaf];`	204	`center=tri->point2[node.axis_leaf];`
207	`}`	205	`}`
208	`// cout << "new left center: " << center << endl;`	206	`// cout << "new left center: " << center << endl;`
209	`}`	207	`}`
210	`if (right_zero)`	208	`if (right_zero)`
211	`{`	209	`{`
212	`// cout << "Old right center: " << center;`	210	`// cout << "Old right center: " << center;`
213	`// Find min position af alle trekanter og placer center der`	211	`// Find min position af alle trekanter og placer center der`
214	`center = left_bbox.min_corner[node.axis_leaf];`	212	`center = left_bbox.min_corner[node.axis_leaf];`
215	`for(j=0;j<objects.size();j++)`	213	`for(j=0;j<objects.size();j++)`
216	`{`	214	`{`
217	`ISectTri* tri = objects[j];`	215	`ISectTri* tri = objects[j];`
218	`if (tri->point0[node.axis_leaf]>center)`	216	`if (tri->point0[node.axis_leaf]>center)`
219	`center=tri->point0[node.axis_leaf];`	217	`center=tri->point0[node.axis_leaf];`
220	`if (tri->point1[node.axis_leaf]>center)`	218	`if (tri->point1[node.axis_leaf]>center)`
221	`center=tri->point1[node.axis_leaf];`	219	`center=tri->point1[node.axis_leaf];`
222	`if (tri->point2[node.axis_leaf]>center)`	220	`if (tri->point2[node.axis_leaf]>center)`
223	`center=tri->point2[node.axis_leaf];`	221	`center=tri->point2[node.axis_leaf];`
224	`}`	222	`}`
225	`// cout << "new right center: " << center << endl;`	223	`// cout << "new right center: " << center << endl;`
226	`}`	224	`}`
227		225
228	`node.plane = center;`	226	`node.plane = center;`
229	`left_bbox.max_corner[node.axis_leaf] = center+f_eps;`	227	`left_bbox.max_corner[node.axis_leaf] = center+f_eps;`
230	`right_bbox.min_corner[node.axis_leaf] = center-f_eps;`	228	`right_bbox.min_corner[node.axis_leaf] = center-f_eps;`
231		229
232		230
233	`// Now put the triangles in the right and left node`	231	`// Now put the triangles in the right and left node`
234	`for(i=0;i<objects.size();i++)`	232	`for(i=0;i<objects.size();i++)`
235	`{`	233	`{`
236	`ISectTri* tri = objects[i];`	234	`ISectTri* tri = objects[i];`
237	`TriAccel *tri_accel = tri_objects[i];`	235	`TriAccel *tri_accel = tri_objects[i];`
238	`if (left_bbox.intersect_triangle(*tri))`	236	`if (left_bbox.intersect_triangle(*tri))`
239	`{`	237	`{`
240	`left_objects.push_back(tri);`	238	`left_objects.push_back(tri);`
241	`tri_left_objects.push_back(tri_accel);`	239	`tri_left_objects.push_back(tri_accel);`
242	`}`	240	`}`
243	`if (right_bbox.intersect_triangle(*tri))`	241	`if (right_bbox.intersect_triangle(*tri))`
244	`{`	242	`{`
245	`right_objects.push_back(tri);`	243	`right_objects.push_back(tri);`
246	`tri_right_objects.push_back(tri_accel);`	244	`tri_right_objects.push_back(tri_accel);`
247	`}`	245	`}`
248	`}`	246	`}`
249	`/* if (left_zero\|\|right_zero) {`	247	`/* if (left_zero\|\|right_zero) {`
250	`cout << left_objects.size() << "," << right_objects.size() << endl;`	248	`cout << left_objects.size() << "," << right_objects.size() << endl;`
251	`}`	249	`}`
252	`*/`	250	`*/`
253	`objects.clear();`	251	`objects.clear();`
254	`subdivide_node(*left_node , left_bbox , level+1, left_objects, tri_left_objects);`	252	`subdivide_node(*left_node , left_bbox , level+1, left_objects, tri_left_objects);`
255	`subdivide_node(*right_node, right_bbox, level+1, right_objects, tri_right_objects);`	253	`subdivide_node(*right_node, right_bbox, level+1, right_objects, tri_right_objects);`
256	`}`	254	`}`
257	`}`	255	`}`
258		256
259	`void BSPTree::init()`	257	`void BSPTree::init()`
260	`{`	258	`{`
261	`root = new BSPNode();`	259	`root = new BSPNode();`
262	`bbox.compute_bbox(isecttris);`	260	`bbox.compute_bbox(isecttris);`
263	`bbox.min_corner-=Vec3f(f_eps100,f_eps100,f_eps*100);`	261	`bbox.min_corner-=Vec3f(f_eps100,f_eps100,f_eps*100);`
264	`bbox.max_corner+=Vec3f(f_eps100,f_eps100,f_eps*100);`	262	`bbox.max_corner+=Vec3f(f_eps100,f_eps100,f_eps*100);`
265	`}`	263	`}`
266		264
267	`void BSPTree::init(vector<TriMesh*>& _trimesh,`	265	`void BSPTree::init(vector<TriMesh*>& _trimesh,`
268	`vector<Mat4x4f>& _transforms,`	266	`vector<Mat4x4f>& _transforms,`
269	`int _max_objects, int _max_level)`	267	`int _max_objects, int _max_level)`
270	`{`	268	`{`
271	`trimesh = _trimesh;`	269	`trimesh = _trimesh;`
272	`transforms = _transforms;`	270	`transforms = _transforms;`
273	`for(unsigned int i=0;i<trimesh.size();i++)`	271	`for(unsigned int i=0;i<trimesh.size();i++)`
274	`{`	272	`{`
275	`TriMesh *mesh = trimesh[i];`	273	`TriMesh *mesh = trimesh[i];`
276	`// Loop through all triangles and add them to intersection structure`	274	`// Loop through all triangles and add them to intersection structure`
277	`for(int j=0;j<mesh->geometry.no_faces();j++)`	275	`for(int j=0;j<mesh->geometry.no_faces();j++)`
278	`{`	276	`{`
279	`Vec3i face = mesh->geometry.face(j);`	277	`Vec3i face = mesh->geometry.face(j);`
280	`ISectTri new_tri;`	278	`ISectTri new_tri;`
281	`new_tri.point0 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[0]));`	279	`new_tri.point0 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[0]));`
282	`new_tri.point1 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[1]));`	280	`new_tri.point1 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[1]));`
283	`new_tri.point2 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[2]));`	281	`new_tri.point2 = transforms[i].mul_3D_point(mesh->geometry.vertex(face[2]));`
284	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`	282	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`
285	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`	283	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`
286	`new_tri.mesh_id = i;`	284	`new_tri.mesh_id = i;`
287	`new_tri.tri_id = j;`	285	`new_tri.tri_id = j;`
288	`isecttris.push_back(new_tri);`	286	`isecttris.push_back(new_tri);`
289	`TriAccel ta;`	287	`TriAccel ta;`
290	`create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);`	288	`create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);`
291	`ta.mesh_id = i;`	289	`ta.mesh_id = i;`
292	`ta.tri_id = j;`	290	`ta.tri_id = j;`
293	`triaccel.push_back(ta);`	291	`triaccel.push_back(ta);`
294	`}`	292	`}`
295	`}`	293	`}`
296		294
297	`max_objects = _max_objects;`	295	`max_objects = _max_objects;`
298	`max_level = _max_level;`	296	`max_level = _max_level;`
299	`init();`	297	`init();`
300	`}`	298	`}`
301		299
302	`void BSPTree::init(TriMesh* mesh, Mat4x4f transform,`	300	`void BSPTree::init(TriMesh* mesh, Mat4x4f transform,`
303	`vector<int> &trilist,`	301	`vector<int> &trilist,`
304	`int _max_objects, int _max_level)`	302	`int _max_objects, int _max_level)`
305	`{`	303	`{`
306	`trimesh.push_back(mesh);`	304	`trimesh.push_back(mesh);`
307	`transforms.push_back(transform);`	305	`transforms.push_back(transform);`
308	`// Loop through all triangles and add them to intersection structure`	306	`// Loop through all triangles and add them to intersection structure`
309	`for(unsigned int j=0;j<trilist.size();j++)`	307	`for(unsigned int j=0;j<trilist.size();j++)`
310	`{`	308	`{`
311	`Vec3i face = mesh->geometry.face(trilist[j]);`	309	`Vec3i face = mesh->geometry.face(trilist[j]);`
312	`ISectTri new_tri;`	310	`ISectTri new_tri;`
313	`new_tri.point0 = transform.mul_3D_point(mesh->geometry.vertex(face[0]));`	311	`new_tri.point0 = transform.mul_3D_point(mesh->geometry.vertex(face[0]));`
314	`new_tri.point1 = transform.mul_3D_point(mesh->geometry.vertex(face[1]));`	312	`new_tri.point1 = transform.mul_3D_point(mesh->geometry.vertex(face[1]));`
315	`new_tri.point2 = transform.mul_3D_point(mesh->geometry.vertex(face[2]));`	313	`new_tri.point2 = transform.mul_3D_point(mesh->geometry.vertex(face[2]));`
316	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`	314	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`
317	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`	315	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`
318	`new_tri.mesh_id = 0;`	316	`new_tri.mesh_id = 0;`
319	`new_tri.tri_id = trilist[j];`	317	`new_tri.tri_id = trilist[j];`
320	`isecttris.push_back(new_tri);`	318	`isecttris.push_back(new_tri);`
321	`TriAccel ta;`	319	`TriAccel ta;`
322	`create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);`	320	`create_tri_accel(new_tri.point0, new_tri.point1, new_tri.point2, ta);`
323	`ta.mesh_id = 0;`	321	`ta.mesh_id = 0;`
324	`ta.tri_id = trilist[j];`	322	`ta.tri_id = trilist[j];`
325	`triaccel.push_back(ta);`	323	`triaccel.push_back(ta);`
326	`}`	324	`}`
327		325
328	`max_objects = _max_objects;`	326	`max_objects = _max_objects;`
329	`max_level = _max_level;`	327	`max_level = _max_level;`
330	`init();`	328	`init();`
331	`}`	329	`}`
332		330
333	`void BSPTree::build()`	331	`void BSPTree::build()`
334	`{`	332	`{`
335	`if (!b_is_build)`	333	`if (!b_is_build)`
336	`{`	334	`{`
337	`vector<ISectTri*> objects;`	335	`vector<ISectTri*> objects;`
338	`vector<TriAccel*> tri_objects;`	336	`vector<TriAccel*> tri_objects;`
339	`for(unsigned int i=0;i<isecttris.size();i++)`	337	`for(unsigned int i=0;i<isecttris.size();i++)`
340	`{`	338	`{`
341	`ISectTri& tri = isecttris[i];`	339	`ISectTri& tri = isecttris[i];`
342	`TriAccel& tri_accel = triaccel[i];`	340	`TriAccel& tri_accel = triaccel[i];`
343	`objects.push_back(&tri);`	341	`objects.push_back(&tri);`
344	`tri_objects.push_back(&tri_accel);`	342	`tri_objects.push_back(&tri_accel);`
345	`}`	343	`}`
346	`subdivide_node(*root, bbox, 0, objects, tri_objects);`	344	`subdivide_node(*root, bbox, 0, objects, tri_objects);`
347	`b_is_build = true;`	345	`b_is_build = true;`
348	`}`	346	`}`
349	`isecttris.clear();`	347	`isecttris.clear();`
350	`all_objects.clear();`	348	`all_objects.clear();`
351	`make_fast_tree(root);`	349	`make_fast_tree(root);`
352	`}`	350	`}`
353		351
354	`bool BSPTree::is_build()`	352	`bool BSPTree::is_build()`
355	`{`	353	`{`
356	`return b_is_build;`	354	`return b_is_build;`
357	`}`	355	`}`
358		356
359	`void BSPTree::print(BSPNode *node, int depth)`	357	`void BSPTree::print(BSPNode *node, int depth)`
360	`{`	358	`{`
361	`if (node==0)`	359	`if (node==0)`
362	`return;`	360	`return;`
363	`for(int i=0;i<depth;i++)`	361	`for(int i=0;i<depth;i++)`
364	`cout << " ";`	362	`cout << " ";`
365	`// cout << "axis:" << node->axis_leaf << ", count:" << node->objects.size() << ", plane:" << node->plane << ", " << endl;`	363	`// cout << "axis:" << node->axis_leaf << ", count:" << node->objects.size() << ", plane:" << node->plane << ", " << endl;`
366	`print(node->left, depth+1);`	364	`print(node->left, depth+1);`
367	`print(node->right, depth+1);`	365	`print(node->right, depth+1);`
368	`}`	366	`}`
369		367
370	`int BSPTree::size(BSPNode *node)`	368	`int BSPTree::size(BSPNode *node)`
371	`{`	369	`{`
372	`if (node==0)`	370	`if (node==0)`
373	`return 0;`	371	`return 0;`
374	`int s = sizeof(BSPNode);`	372	`int s = sizeof(BSPNode);`
375	`s+= node->count * sizeof(ISectTri);`	373	`s+= node->count * sizeof(ISectTri);`
376	`s+=size(node->left);`	374	`s+=size(node->left);`
377	`s+=size(node->right);`	375	`s+=size(node->right);`
378	`return s;`	376	`return s;`
379	`}`	377	`}`
380		378
381	`int BSPTree::size()`	379	`int BSPTree::size()`
382	`{`	380	`{`
383	`return size(root);`	381	`return size(root);`
384	`}`	382	`}`
385		383
386	`/__declspec(align(16))/ static const unsigned int modulo[] = {0,1,2,0,1};`	384	`/__declspec(align(16))/ static const unsigned int modulo[] = {0,1,2,0,1};`
387		385
388	`inline bool intersect2(Ray &ray, const TriAccel &acc, double t_max)`	386	`inline bool intersect2(Ray &ray, const TriAccel &acc, double t_max)`
389	`{`	387	`{`
390	`// cout << "Acc: " << (int)&acc << " ";`	388	`// cout << "Acc: " << (int)&acc << " ";`
391	`//inline bool Intersect(TriAccel &acc,Ray &ray)`	389	`//inline bool Intersect(TriAccel &acc,Ray &ray)`
392	`#define ku modulo[acc.k+1]`	390	`#define ku modulo[acc.k+1]`
393	`#define kv modulo[acc.k+2]`	391	`#define kv modulo[acc.k+2]`
394	`// dont prefetch here, assume data has already been prefetched`	392	`// dont prefetch here, assume data has already been prefetched`
395	`// start high-latency division as early as possible`	393	`// start high-latency division as early as possible`
396	`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]);`	394	`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]);`
397	`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;`	395	`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;`
398	`// check for valid distance.`	396	`// check for valid distance.`
399	`if (!(t_max > f && f > f_eps)\|\|ray.dist<f) return false;`	397	`if (!(t_max > f && f > f_eps)\|\|ray.dist<f) return false;`
400	`// compute hitpoint positions on uv plane`	398	`// compute hitpoint positions on uv plane`
401	`const double hu = (ray.origin[ku] + f * ray.direction[ku]);`	399	`const double hu = (ray.origin[ku] + f * ray.direction[ku]);`
402	`const double hv = (ray.origin[kv] + f * ray.direction[kv]);`	400	`const double hv = (ray.origin[kv] + f * ray.direction[kv]);`
403	`// check first barycentric coordinate`	401	`// check first barycentric coordinate`
404	`const double lambda = (hu * (double)acc.b_nu + hv * (double)acc.b_nv + (double)acc.b_d);`	402	`const double lambda = (hu * (double)acc.b_nu + hv * (double)acc.b_nv + (double)acc.b_d);`
405	`if (lambda < 0.0) return false;`	403	`if (lambda < 0.0) return false;`
406	`// check second barycentric coordinate`	404	`// check second barycentric coordinate`
407	`const double mue = (hu * (double)acc.c_nu + hv * (double)acc.c_nv + (double)acc.c_d);`	405	`const double mue = (hu * (double)acc.c_nu + hv * (double)acc.c_nv + (double)acc.c_d);`
408	`if (mue < 0.0) return false;`	406	`if (mue < 0.0) return false;`
409	`// check third barycentric coordinate`	407	`// check third barycentric coordinate`
410	`if (lambda+mue > 1.0) return false;`	408	`if (lambda+mue > 1.0) return false;`
411	`// have a valid hitpoint here. store it.`	409	`// have a valid hitpoint here. store it.`
412	`ray.dist = f;`	410	`ray.dist = f;`
413	`ray.u = lambda;`	411	`ray.u = lambda;`
414	`ray.v = mue;`	412	`ray.v = mue;`
415	`ray.hit_object = (TriMesh*)acc.mesh_id;`	413	`ray.hit_object = (TriMesh*)acc.mesh_id;`
416	`ray.hit_face_id = acc.tri_id;`	414	`ray.hit_face_id = acc.tri_id;`
417	`ray.has_hit=true;`	415	`ray.has_hit=true;`
418	`return true;`	416	`return true;`
419	`}`	417	`}`
420		418
421	`bool BSPTree::intersect_node(Ray &ray, const BSPNode &node, double t_min, double t_max) const`	419	`bool BSPTree::intersect_node(Ray &ray, const BSPNode &node, double t_min, double t_max) const`
422	`{`	420	`{`
423	`//cout << node.plane << endl;`	421	`//cout << node.plane << endl;`
424	`node_calls++;`	422	`node_calls++;`
425	`static bool found;`	423	`static bool found;`
426	`static int i;`	424	`static int i;`
427		425
428	`if (node.axis_leaf==4)`	426	`if (node.axis_leaf==4)`
429	`{`	427	`{`
430	`found = false;`	428	`found = false;`
431	`for(i=0;i<node.count;i++)`	429	`for(i=0;i<node.count;i++)`
432	`{`	430	`{`
433	`// const TriAccel* tri2 = all_triaccel[node.id+i];`	431	`// const TriAccel* tri2 = all_triaccel[node.id+i];`
434	`const ISectTri* tri = all_objects[node.id+i];`	432	`const ISectTri* tri = all_objects[node.id+i];`
435	`// if (intersect2(ray, *tri2, t_max))`	433	`// if (intersect2(ray, *tri2, t_max))`
436	`// found=true;`	434	`// found=true;`
437	`if (intersect(ray, *tri, t_max))`	435	`if (intersect(ray, *tri, t_max))`
438	`found=true;`	436	`found=true;`
439	`}`	437	`}`
440	`if (found)`	438	`if (found)`
441	`return true;`	439	`return true;`
442	`else`	440	`else`
443	`return false;`	441	`return false;`
444	`}`	442	`}`
445	`else`	443	`else`
446	`{`	444	`{`
447	`BSPNode *near_node;`	445	`BSPNode *near_node;`
448	`BSPNode *far_node;`	446	`BSPNode *far_node;`
449	`if (ray.direction[node.axis_leaf]>=0)`	447	`if (ray.direction[node.axis_leaf]>=0)`
450	`{`	448	`{`
451	`near_node = node.left;`	449	`near_node = node.left;`
452	`far_node = node.right;`	450	`far_node = node.right;`
453	`}`	451	`}`
454	`else`	452	`else`
455	`{`	453	`{`
456	`near_node = node.right;`	454	`near_node = node.right;`
457	`far_node = node.left;`	455	`far_node = node.left;`
458	`}`	456	`}`
459		457
460	`// In order to avoid instability`	458	`// In order to avoid instability`
461	`double t;`	459	`double t;`
462	`if (fabs(ray.direction[node.axis_leaf])<d_eps)`	460	`if (fabs(ray.direction[node.axis_leaf])<d_eps)`
463	`t = (node.plane - ray.origin[node.axis_leaf])/d_eps;// intersect node plane;`	461	`t = (node.plane - ray.origin[node.axis_leaf])/d_eps;// intersect node plane;`
464	`else`	462	`else`
465	`t = (node.plane - ray.origin[node.axis_leaf])/ray.direction[node.axis_leaf];// intersect node plane;`	463	`t = (node.plane - ray.origin[node.axis_leaf])/ray.direction[node.axis_leaf];// intersect node plane;`
466		464
467	`if (t>t_max)`	465	`if (t>t_max)`
468	`return intersect_node(ray, *near_node, t_min, t_max);`	466	`return intersect_node(ray, *near_node, t_min, t_max);`
469		467
470	`else if (t<t_min)`	468	`else if (t<t_min)`
471	`return intersect_node(ray, *far_node, t_min, t_max);`	469	`return intersect_node(ray, *far_node, t_min, t_max);`
472	`else`	470	`else`
473	`{`	471	`{`
474	`if (intersect_node(ray, *near_node, t_min, t))`	472	`if (intersect_node(ray, *near_node, t_min, t))`
475	`return true;`	473	`return true;`
476	`else`	474	`else`
477	`return intersect_node(ray, *far_node, t, t_max);`	475	`return intersect_node(ray, *far_node, t, t_max);`
478	`}`	476	`}`
479	`}`	477	`}`
480	`}`	478	`}`
481		479
482	`bool BSPTree::intersect(Ray &ray) const`	480	`bool BSPTree::intersect(Ray &ray) const`
483	`{`	481	`{`
484	`double t_min, t_max;`	482	`double t_min, t_max;`
485	`bbox.intersect_min_max(ray, t_min, t_max);`	483	`bbox.intersect_min_max(ray, t_min, t_max);`
486	`if (t_min>t_max)`	484	`if (t_min>t_max)`
487	`return false;`	485	`return false;`
488	`if (!intersect_node(ray, *root, t_min, t_max))`	486	`if (!intersect_node(ray, *root, t_min, t_max))`
489	`return false;`	487	`return false;`
490	`// cout << "____" << endl;`	488	`// cout << "____" << endl;`
491	`// ray.reset();`	489	`// ray.reset();`
492	`// cout << "Here " << endl;`	490	`// cout << "Here " << endl;`
493	`// intersect_fast_node(ray, &fast_tree[0], t_min, t_max);`	491	`// intersect_fast_node(ray, &fast_tree[0], t_min, t_max);`
494	`// if (!ray.has_hit)`	492	`// if (!ray.has_hit)`
495	`// return false;`	493	`// return false;`
496	`else`	494	`else`
497	`{`	495	`{`
498	`// Calculate the normal at the intersection`	496	`// Calculate the normal at the intersection`
499	`ray.hit_object = trimesh[(int)ray.hit_object];`	497	`ray.hit_object = trimesh[(int)ray.hit_object];`
500	`Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);`	498	`Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);`

Subversion Repositories gelsvn

(root)/branches/SPD/src/Geometry/BSPTree.cpp – Rev 208 → 290