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


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

Subversion Repositories gelsvn

(root)/trunk/src/Geometry/BSPTree.cpp – Rev 541 → 595