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

Subversion Repositories gelsvn

(root)/trunk/src/Geometry/BSPTree.cpp – Rev 311 → 313