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<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) 
  {
//  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 > 0.001)||ray.dist<f) return false;
    // compute hitpoint positions on uv plane
    const double hu = (ray.origin[ku] + f * ray.direction[ku]);
    const double hv = (ray.origin[kv] + f * ray.direction[kv]);
    // check first barycentric coordinate
    const double lambda = (hu * (double)acc.b_nu + hv * (double)acc.b_nv + (double)acc.b_d);
    if (lambda < 0.0) return false;
    // check second barycentric coordinate
    const double mue = (hu * (double)acc.c_nu + hv * (double)acc.c_nv + (double)acc.c_d);
    if (mue < 0.0) return false;
    // check third barycentric coordinate
    if (lambda+mue > 1.0) return false;
    // have a valid hitpoint here. store it.
    ray.dist = f;
    ray.u = lambda;
    ray.v = mue;
    ray.hit_object = (TriMesh*)acc.mesh_id;
    ray.hit_face_id = acc.tri_id;
    ray.has_hit=true;
    return true;
  }
 
  bool BSPTree::intersect_node(Ray &ray, const BSPNode &node, double t_min, double t_max) const 
  {
//    cout << node.plane << endl;
    node_calls++;    
 
 
 
    if (node.axis_leaf==4) 
    {
      bool found = false; 
      for(int 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;
 
//    printf("%.2f %.2f %.2f\n", ray.dist, t_min, t_max);
    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 
    {
      //printf("HIT\n");
 
      // Calculate the normal at the intersection
      ray.id = reinterpret_cast<int>(ray.hit_object);
      ray.hit_object = trimesh[ray.id];
      
      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 = 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;
/*
      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
    {
      //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))
#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
    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 < 0.0)
      return false;
    if(fabs(t) < f_eps)
      return false;
 
 
 
 
 
    if(t_max < t)
      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 428	Rev 429
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		122
123	`int new_axis=-1;`	123	`int new_axis=-1;`
124	`double min_cost=CGLA::BIG;`	124	`double min_cost=CGLA::BIG;`
125	`int new_pos = -1;`	125	`int new_pos = -1;`
126		126
127	`for(i=0;i<3;i++)`	127	`for(i=0;i<3;i++)`
128	`{`	128	`{`
129	`for(int k=1;k<TESTS;k++)`	129	`for(int k=1;k<TESTS;k++)`
130	`{`	130	`{`
131	`BBox left_bbox = bbox;`	131	`BBox left_bbox = bbox;`
132	`BBox right_bbox = bbox;`	132	`BBox right_bbox = bbox;`
133		133
134	`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];`
135	`node.plane = center;`	135	`node.plane = center;`
136		136
137	`left_bbox.max_corner[i] = center;`	137	`left_bbox.max_corner[i] = center;`
138	`right_bbox.min_corner[i] = center;`	138	`right_bbox.min_corner[i] = center;`
139		139
140	`// Try putting the triangles in the left and right boxes`	140	`// Try putting the triangles in the left and right boxes`
141	`int left_count = 0;`	141	`int left_count = 0;`
142	`int right_count = 0;`	142	`int right_count = 0;`
143	`for(unsigned int j=0;j<objects.size();j++)`	143	`for(unsigned int j=0;j<objects.size();j++)`
144	`{`	144	`{`
145	`ISectTri* tri = objects[j];`	145	`ISectTri* tri = objects[j];`
146	`left_count += left_bbox.intersect_triangle(*tri);`	146	`left_count += left_bbox.intersect_triangle(*tri);`
147	`right_count += right_bbox.intersect_triangle(*tri);`	147	`right_count += right_bbox.intersect_triangle(*tri);`
148	`}`	148	`}`
149		149
150	`//double len = bbox.max_corner[i] - bbox.min_corner[i];`	150	`//double len = bbox.max_corner[i] - bbox.min_corner[i];`
151	`double cost = left_countleft_bbox.area() + right_countright_bbox.area(); // - len*len;`	151	`double cost = left_countleft_bbox.area() + right_countright_bbox.area(); // - len*len;`
152	`if(cost < min_cost)`	152	`if(cost < min_cost)`
153	`{`	153	`{`
154	`min_cost = cost;`	154	`min_cost = cost;`
155	`new_axis = i;`	155	`new_axis = i;`
156	`new_pos = k;`	156	`new_pos = k;`
157	`right_zero = (right_count==0);`	157	`right_zero = (right_count==0);`
158	`left_zero = (left_count==0);`	158	`left_zero = (left_count==0);`
159	`}`	159	`}`
160	`}`	160	`}`
161	`}`	161	`}`
162	`node.axis_leaf = new_axis;`	162	`node.axis_leaf = new_axis;`
163	`left_node->axis_leaf = static_cast<unsigned char>(-1);`	163	`left_node->axis_leaf = static_cast<unsigned char>(-1);`
164	`right_node->axis_leaf = static_cast<unsigned char>(-1);`	164	`right_node->axis_leaf = static_cast<unsigned char>(-1);`
165		165
166	`// Now chose the right splitting plane`	166	`// Now chose the right splitting plane`
167	`BBox left_bbox = bbox;`	167	`BBox left_bbox = bbox;`
168	`BBox right_bbox = bbox;`	168	`BBox right_bbox = bbox;`
169		169
170	`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];`
171	`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];`
172	`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;`
173		173
174	`if (left_zero)`	174	`if (left_zero)`
175	`{`	175	`{`
176	`// Find min position of all triangle vertices and place the center there`	176	`// Find min position of all triangle vertices and place the center there`
177	`center = bbox.max_corner[node.axis_leaf];`	177	`center = bbox.max_corner[node.axis_leaf];`
178	`for(unsigned int j=0;j<objects.size();j++)`	178	`for(unsigned int j=0;j<objects.size();j++)`
179	`{`	179	`{`
180	`ISectTri* tri = objects[j];`	180	`ISectTri* tri = objects[j];`
181	`if (tri->point0[node.axis_leaf]<center)`	181	`if (tri->point0[node.axis_leaf]<center)`
182	`center=tri->point0[node.axis_leaf];`	182	`center=tri->point0[node.axis_leaf];`
183	`if (tri->point1[node.axis_leaf]<center)`	183	`if (tri->point1[node.axis_leaf]<center)`
184	`center=tri->point1[node.axis_leaf];`	184	`center=tri->point1[node.axis_leaf];`
185	`if (tri->point2[node.axis_leaf]<center)`	185	`if (tri->point2[node.axis_leaf]<center)`
186	`center=tri->point2[node.axis_leaf];`	186	`center=tri->point2[node.axis_leaf];`
187	`}`	187	`}`
188	`center -= diff;`	188	`center -= diff;`
189	`}`	189	`}`
190	`if (right_zero)`	190	`if (right_zero)`
191	`{`	191	`{`
192	`// Find max position of all triangle vertices and place the center there`	192	`// Find max position of all triangle vertices and place the center there`
193	`center = bbox.min_corner[node.axis_leaf];`	193	`center = bbox.min_corner[node.axis_leaf];`
194	`for(unsigned int j=0;j<objects.size();j++)`	194	`for(unsigned int j=0;j<objects.size();j++)`
195	`{`	195	`{`
196	`ISectTri* tri = objects[j];`	196	`ISectTri* tri = objects[j];`
197	`if (tri->point0[node.axis_leaf]>center)`	197	`if (tri->point0[node.axis_leaf]>center)`
198	`center=tri->point0[node.axis_leaf];`	198	`center=tri->point0[node.axis_leaf];`
199	`if (tri->point1[node.axis_leaf]>center)`	199	`if (tri->point1[node.axis_leaf]>center)`
200	`center=tri->point1[node.axis_leaf];`	200	`center=tri->point1[node.axis_leaf];`
201	`if (tri->point2[node.axis_leaf]>center)`	201	`if (tri->point2[node.axis_leaf]>center)`
202	`center=tri->point2[node.axis_leaf];`	202	`center=tri->point2[node.axis_leaf];`
203	`}`	203	`}`
204	`center += diff;`	204	`center += diff;`
205	`}`	205	`}`
206		206
207	`node.plane = center;`	207	`node.plane = center;`
208	`left_bbox.max_corner[node.axis_leaf] = center;`	208	`left_bbox.max_corner[node.axis_leaf] = center;`
209	`right_bbox.min_corner[node.axis_leaf] = center;`	209	`right_bbox.min_corner[node.axis_leaf] = center;`
210		210
211	`// Now put the triangles in the right and left node`	211	`// Now put the triangles in the right and left node`
212	`for(i=0;i<objects.size();i++)`	212	`for(i=0;i<objects.size();i++)`
213	`{`	213	`{`
214	`ISectTri* tri = objects[i];`	214	`ISectTri* tri = objects[i];`
215	`TriAccel *tri_accel = tri_objects[i];`	215	`TriAccel *tri_accel = tri_objects[i];`
216	`if (left_bbox.intersect_triangle(*tri))`	216	`if (left_bbox.intersect_triangle(*tri))`
217	`{`	217	`{`
218	`left_objects.push_back(tri);`	218	`left_objects.push_back(tri);`
219	`tri_left_objects.push_back(tri_accel);`	219	`tri_left_objects.push_back(tri_accel);`
220	`}`	220	`}`
221	`if (right_bbox.intersect_triangle(*tri))`	221	`if (right_bbox.intersect_triangle(*tri))`
222	`{`	222	`{`
223	`right_objects.push_back(tri);`	223	`right_objects.push_back(tri);`
224	`tri_right_objects.push_back(tri_accel);`	224	`tri_right_objects.push_back(tri_accel);`
225	`}`	225	`}`
226	`}`	226	`}`
227	`//if (left_zero\|\|right_zero)`	227	`//if (left_zero\|\|right_zero)`
228	`// cout << left_objects.size() << "," << right_objects.size() << "," << level << endl;`	228	`// cout << left_objects.size() << "," << right_objects.size() << "," << level << endl;`
229		229
230	`objects.clear();`	230	`objects.clear();`
231	`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);`
232	`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);`
233	`}`	233	`}`
234	`}`	234	`}`
235		235
236	`void BSPTree::init()`	236	`void BSPTree::init()`
237	`{`	237	`{`
238	`root = new BSPNode();`	238	`root = new BSPNode();`
239	`bbox.compute_bbox(isecttris);`	239	`bbox.compute_bbox(isecttris);`
240	`bbox.min_corner-=Vec3f(1.0);`	240	`bbox.min_corner-=Vec3f(1.0);`
241	`bbox.max_corner+=Vec3f(1.0);`	241	`bbox.max_corner+=Vec3f(1.0);`
242	`}`	242	`}`
243		243
244	`void BSPTree::init(vector<const TriMesh*>& _trimesh,`	244	`void BSPTree::init(vector<const TriMesh*>& _trimesh,`
245	`vector<Mat4x4f>& _transforms,`	245	`vector<Mat4x4f>& _transforms,`
246	`int _max_objects, int _max_level)`	246	`int _max_objects, int _max_level)`
247	`{`	247	`{`
248	`trimesh = _trimesh;`	248	`trimesh = _trimesh;`
249	`transforms = _transforms;`	249	`transforms = _transforms;`
250	`for(unsigned int i=0;i<trimesh.size();i++)`	250	`for(unsigned int i=0;i<trimesh.size();i++)`
251	`{`	251	`{`
252	`const TriMesh *mesh = trimesh[i];`	252	`const TriMesh *mesh = trimesh[i];`
253	`// Loop through all triangles and add them to intersection structure`	253	`// Loop through all triangles and add them to intersection structure`
254	`for(int j=0;j<mesh->geometry.no_faces();j++)`	254	`for(int j=0;j<mesh->geometry.no_faces();j++)`
255	`{`	255	`{`
256	`Vec3i face = mesh->geometry.face(j);`	256	`Vec3i face = mesh->geometry.face(j);`
257	`ISectTri new_tri;`	257	`ISectTri new_tri;`
258	`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]));`
259	`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]));`
260	`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]));`
261	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`	261	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`
262	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`	262	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`
263	`new_tri.mesh_id = i;`	263	`new_tri.mesh_id = i;`
264	`new_tri.tri_id = j;`	264	`new_tri.tri_id = j;`
265	`isecttris.push_back(new_tri);`	265	`isecttris.push_back(new_tri);`
266	`TriAccel ta;`	266	`TriAccel ta;`
267	`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);`
268	`ta.mesh_id = i;`	268	`ta.mesh_id = i;`
269	`ta.tri_id = j;`	269	`ta.tri_id = j;`
270	`triaccel.push_back(ta);`	270	`triaccel.push_back(ta);`
271	`}`	271	`}`
272	`}`	272	`}`
273		273
274	`max_objects = _max_objects;`	274	`max_objects = _max_objects;`
275	`max_level = _max_level;`	275	`max_level = _max_level;`
276	`init();`	276	`init();`
277	`}`	277	`}`
278		278
279	`void BSPTree::init(const TriMesh* mesh, Mat4x4f transform,`	279	`void BSPTree::init(const TriMesh* mesh, Mat4x4f transform,`
280	`vector<int> &trilist,`	280	`vector<int> &trilist,`
281	`int _max_objects, int _max_level)`	281	`int _max_objects, int _max_level)`
282	`{`	282	`{`
283	`trimesh.push_back(mesh);`	283	`trimesh.push_back(mesh);`
284	`transforms.push_back(transform);`	284	`transforms.push_back(transform);`
285	`// Loop through all triangles and add them to intersection structure`	285	`// Loop through all triangles and add them to intersection structure`
286	`for(unsigned int j=0;j<trilist.size();j++)`	286	`for(unsigned int j=0;j<trilist.size();j++)`
287	`{`	287	`{`
288	`Vec3i face = mesh->geometry.face(trilist[j]);`	288	`Vec3i face = mesh->geometry.face(trilist[j]);`
289	`ISectTri new_tri;`	289	`ISectTri new_tri;`
290	`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]));`
291	`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]));`
292	`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]));`
293	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`	293	`new_tri.edge0 = new_tri.point1 - new_tri.point0;`
294	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`	294	`new_tri.edge1 = new_tri.point2 - new_tri.point0;`
295	`new_tri.mesh_id = 0;`	295	`new_tri.mesh_id = 0;`
296	`new_tri.tri_id = trilist[j];`	296	`new_tri.tri_id = trilist[j];`
297	`isecttris.push_back(new_tri);`	297	`isecttris.push_back(new_tri);`
298	`TriAccel ta;`	298	`TriAccel ta;`
299	`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);`
300	`ta.mesh_id = 0;`	300	`ta.mesh_id = 0;`
301	`ta.tri_id = trilist[j];`	301	`ta.tri_id = trilist[j];`
302	`triaccel.push_back(ta);`	302	`triaccel.push_back(ta);`
303	`}`	303	`}`
304		304
305	`max_objects = _max_objects;`	305	`max_objects = _max_objects;`
306	`max_level = _max_level;`	306	`max_level = _max_level;`
307	`init();`	307	`init();`
308	`}`	308	`}`
309		309
310	`void BSPTree::build()`	310	`void BSPTree::build()`
311	`{`	311	`{`
312	`if (!b_is_build)`	312	`if (!b_is_build)`
313	`{`	313	`{`
314	`vector<ISectTri*> objects;`	314	`vector<ISectTri*> objects;`
315	`vector<TriAccel*> tri_objects;`	315	`vector<TriAccel*> tri_objects;`
316	`for(unsigned int i=0;i<isecttris.size();i++)`	316	`for(unsigned int i=0;i<isecttris.size();i++)`
317	`{`	317	`{`
318	`ISectTri* tri = &isecttris[i];`	318	`ISectTri* tri = &isecttris[i];`
319	`TriAccel* tri_accel = &triaccel[i];`	319	`TriAccel* tri_accel = &triaccel[i];`
320	`objects.push_back(tri);`	320	`objects.push_back(tri);`
321	`tri_objects.push_back(tri_accel);`	321	`tri_objects.push_back(tri_accel);`
322	`}`	322	`}`
323	`subdivide_node(*root, bbox, 0, objects, tri_objects);`	323	`subdivide_node(*root, bbox, 0, objects, tri_objects);`
324	`b_is_build = true;`	324	`b_is_build = true;`
325	`}`	325	`}`
326	`make_fast_tree(root);`	326	`make_fast_tree(root);`
327	`}`	327	`}`
328		328
329	`bool BSPTree::is_build()`	329	`bool BSPTree::is_build()`
330	`{`	330	`{`
331	`return b_is_build;`	331	`return b_is_build;`
332	`}`	332	`}`
333		333
334	`void BSPTree::print(BSPNode *node, int depth)`	334	`void BSPTree::print(BSPNode *node, int depth)`
335	`{`	335	`{`
336	`if (node==0)`	336	`if (node==0)`
337	`return;`	337	`return;`
338	`for(int i=0;i<depth;i++)`	338	`for(int i=0;i<depth;i++)`
339	`cout << " ";`	339	`cout << " ";`
340	`// cout << "axis:" << node->axis_leaf << ", count:" << node->objects.size() << ", plane:" << node->plane << ", " << endl;`	340	`// cout << "axis:" << node->axis_leaf << ", count:" << node->objects.size() << ", plane:" << node->plane << ", " << endl;`
341	`print(node->left, depth+1);`	341	`print(node->left, depth+1);`
342	`print(node->right, depth+1);`	342	`print(node->right, depth+1);`
343	`}`	343	`}`
344		344
345	`int BSPTree::size(BSPNode *node)`	345	`int BSPTree::size(BSPNode *node)`
346	`{`	346	`{`
347	`if (node==0)`	347	`if (node==0)`
348	`return 0;`	348	`return 0;`
349	`int s = sizeof(BSPNode);`	349	`int s = sizeof(BSPNode);`
350	`s+= node->count * sizeof(ISectTri);`	350	`s+= node->count * sizeof(ISectTri);`
351	`s+=size(node->left);`	351	`s+=size(node->left);`
352	`s+=size(node->right);`	352	`s+=size(node->right);`
353	`return s;`	353	`return s;`
354	`}`	354	`}`
355		355
356	`int BSPTree::size()`	356	`int BSPTree::size()`
357	`{`	357	`{`
358	`return size(root);`	358	`return size(root);`
359	`}`	359	`}`
360		360
361	`/__declspec(align(16))/ static const unsigned int modulo[] = {0,1,2,0,1};`	361	`/__declspec(align(16))/ static const unsigned int modulo[] = {0,1,2,0,1};`
362		362
363	`inline bool intersect2(Ray &ray, const TriAccel &acc, double t_max)`	363	`inline bool intersect2(Ray &ray, const TriAccel &acc, double t_max)`
364	`{`	364	`{`
365	`// cout << "Acc: " << (int)&acc << " ";`	365	`// cout << "Acc: " << (int)&acc << " ";`
366	`//inline bool Intersect(TriAccel &acc,Ray &ray)`	366	`//inline bool Intersect(TriAccel &acc,Ray &ray)`
367	`#define ku modulo[acc.k+1]`	367	`#define ku modulo[acc.k+1]`
368	`#define kv modulo[acc.k+2]`	368	`#define kv modulo[acc.k+2]`
369	`// dont prefetch here, assume data has already been prefetched`	369	`// dont prefetch here, assume data has already been prefetched`
370	`// start high-latency division as early as possible`	370	`// start high-latency division as early as possible`
371	`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 nd = 1.0/((double)ray.direction[acc.k] + (double)acc.n_u * (double)ray.direction[ku] + (double)acc.n_v * (double)ray.direction[kv]);`
372	`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	`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;`
373	`// check for valid distance.`	373	`// check for valid distance.`
374	`if (!(t_max > f && f > 0.001)\|\|ray.dist<f) return false;`	374	`if (!(t_max > f && f > 0.001)\|\|ray.dist<f) return false;`
375	`// compute hitpoint positions on uv plane`	375	`// compute hitpoint positions on uv plane`
376	`const double hu = (ray.origin[ku] + f * ray.direction[ku]);`	376	`const double hu = (ray.origin[ku] + f * ray.direction[ku]);`
377	`const double hv = (ray.origin[kv] + f * ray.direction[kv]);`	377	`const double hv = (ray.origin[kv] + f * ray.direction[kv]);`
378	`// check first barycentric coordinate`	378	`// check first barycentric coordinate`
379	`const double lambda = (hu * (double)acc.b_nu + hv * (double)acc.b_nv + (double)acc.b_d);`	379	`const double lambda = (hu * (double)acc.b_nu + hv * (double)acc.b_nv + (double)acc.b_d);`
380	`if (lambda < 0.0) return false;`	380	`if (lambda < 0.0) return false;`
381	`// check second barycentric coordinate`	381	`// check second barycentric coordinate`
382	`const double mue = (hu * (double)acc.c_nu + hv * (double)acc.c_nv + (double)acc.c_d);`	382	`const double mue = (hu * (double)acc.c_nu + hv * (double)acc.c_nv + (double)acc.c_d);`
383	`if (mue < 0.0) return false;`	383	`if (mue < 0.0) return false;`
384	`// check third barycentric coordinate`	384	`// check third barycentric coordinate`
385	`if (lambda+mue > 1.0) return false;`	385	`if (lambda+mue > 1.0) return false;`
386	`// have a valid hitpoint here. store it.`	386	`// have a valid hitpoint here. store it.`
387	`ray.dist = f;`	387	`ray.dist = f;`
388	`ray.u = lambda;`	388	`ray.u = lambda;`
389	`ray.v = mue;`	389	`ray.v = mue;`
390	`ray.hit_object = (TriMesh*)acc.mesh_id;`	390	`ray.hit_object = (TriMesh*)acc.mesh_id;`
391	`ray.hit_face_id = acc.tri_id;`	391	`ray.hit_face_id = acc.tri_id;`
392	`ray.has_hit=true;`	392	`ray.has_hit=true;`
393	`return true;`	393	`return true;`
394	`}`	394	`}`
395		395
396	`bool BSPTree::intersect_node(Ray &ray, const BSPNode &node, double t_min, double t_max) const`	396	`bool BSPTree::intersect_node(Ray &ray, const BSPNode &node, double t_min, double t_max) const`
397	`{`	397	`{`
398	`// cout << node.plane << endl;`	398	`// cout << node.plane << endl;`
399	`node_calls++;`	399	`node_calls++;`
400	`static bool found;`	-
401	`static int i;`	-
402		-
403	`if (node.axis_leaf==4)`	400	`if (node.axis_leaf==4)`
404	`{`	401	`{`
405	`found = false;`	402	`bool found = false;`
406	`for(i=0;i<node.count;i++)`	403	`for(int i=0; i < node.count; ++i)`
407	`{`	404	`{`
408	`// const TriAccel* tri2 = all_triaccel[node.id+i];`	405	`// const TriAccel* tri2 = all_triaccel[node.id+i];`
409	`const ISectTri* tri = all_objects[node.id+i];`	406	`const ISectTri* tri = all_objects[node.id+i];`
410	`// if (intersect2(ray, *tri2, t_max))`	407	`// if (intersect2(ray, *tri2, t_max))`
411	`// found=true;`	408	`// found=true;`
412	`if (intersect(ray, *tri, t_max))`	409	`if (intersect(ray, *tri, t_max))`
413	`found=true;`	410	`found=true;`
414	`}`	411	`}`
415	`if (found)`	412	`if (found)`
416	`return true;`	413	`return true;`
417	`else`	414	`else`
418	`return false;`	415	`return false;`
419	`}`	416	`}`
420	`else`	417	`else`
421	`{`	418	`{`
422	`BSPNode *near_node;`	419	`BSPNode *near_node;`
423	`BSPNode *far_node;`	420	`BSPNode *far_node;`
424	`if (ray.direction[node.axis_leaf]>=0)`	421	`if (ray.direction[node.axis_leaf]>=0)`
425	`{`	422	`{`
426	`near_node = node.left;`	423	`near_node = node.left;`
427	`far_node = node.right;`	424	`far_node = node.right;`
428	`}`	425	`}`
429	`else`	426	`else`
430	`{`	427	`{`
431	`near_node = node.right;`	428	`near_node = node.right;`
432	`far_node = node.left;`	429	`far_node = node.left;`
433	`}`	430	`}`
434		431
435	`// In order to avoid instability`	432	`// In order to avoid instability`
436	`double t;`	433	`double t;`
437	`if (fabs(ray.direction[node.axis_leaf])<d_eps)`	434	`if (fabs(ray.direction[node.axis_leaf])<d_eps)`
438	`t = (node.plane - ray.origin[node.axis_leaf])/d_eps;// intersect node plane;`	435	`t = (node.plane - ray.origin[node.axis_leaf])/d_eps;// intersect node plane;`
439	`else`	436	`else`
440	`t = (node.plane - ray.origin[node.axis_leaf])/ray.direction[node.axis_leaf];// intersect node plane;`	437	`t = (node.plane - ray.origin[node.axis_leaf])/ray.direction[node.axis_leaf];// intersect node plane;`
441		438
442	`if (t>t_max)`	439	`if (t>t_max)`
443	`return intersect_node(ray, *near_node, t_min, t_max);`	440	`return intersect_node(ray, *near_node, t_min, t_max);`
444	`else if (t<t_min)`	441	`else if (t<t_min)`
445	`return intersect_node(ray, *far_node, t_min, t_max);`	442	`return intersect_node(ray, *far_node, t_min, t_max);`
446	`else`	443	`else`
447	`{`	444	`{`
448	`if (intersect_node(ray, *near_node, t_min, t))`	445	`if (intersect_node(ray, *near_node, t_min, t))`
449	`return true;`	446	`return true;`
450	`else`	447	`else`
451	`return intersect_node(ray, *far_node, t, t_max);`	448	`return intersect_node(ray, *far_node, t, t_max);`
452	`}`	449	`}`
453	`}`	450	`}`
454	`}`	451	`}`
455		452
456	`bool BSPTree::intersect(Ray &ray) const`	453	`bool BSPTree::intersect(Ray &ray) const`
457	`{`	454	`{`
458	`double t_min, t_max;`	455	`double t_min, t_max;`
459	`bbox.intersect_min_max(ray, t_min, t_max);`	456	`bbox.intersect_min_max(ray, t_min, t_max);`
460	`if (t_min>t_max)`	457	`if (t_min>t_max)`
461	`return false;`	458	`return false;`
462		459
463	`// printf("%.2f %.2f %.2f\n", ray.dist, t_min, t_max);`	460	`// printf("%.2f %.2f %.2f\n", ray.dist, t_min, t_max);`
464	`if (!intersect_node(ray, *root, t_min, t_max))`	461	`if (!intersect_node(ray, *root, t_min, t_max))`
465	`return false;`	462	`return false;`
466	`// cout << "____" << endl;`	463	`// cout << "____" << endl;`
467	`// ray.reset();`	464	`// ray.reset();`
468	`// cout << "Here " << endl;`	465	`// cout << "Here " << endl;`
469	`// intersect_fast_node(ray, &fast_tree[0], t_min, t_max);`	466	`// intersect_fast_node(ray, &fast_tree[0], t_min, t_max);`
470	`// if (!ray.has_hit)`	467	`// if (!ray.has_hit)`
471	`//return false;`	468	`//return false;`
472	`else`	469	`else`
473	`{`	470	`{`
474	`//printf("HIT\n");`	471	`//printf("HIT\n");`
475		472
476	`// Calculate the normal at the intersection`	473	`// Calculate the normal at the intersection`
477	`ray.id = reinterpret_cast<int>(ray.hit_object);`	474	`ray.id = reinterpret_cast<int>(ray.hit_object);`
478	`ray.hit_object = trimesh[ray.id];`	475	`ray.hit_object = trimesh[ray.id];`
479		476
480	`Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);`	477	`Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);`
481	`Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);`	478	`Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);`
482	`Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);`	479	`Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);`
483	`Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);`	480	`Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);`
484	`ray.hit_normal = transforms[ray.id].mul_3D_vector(`	481	`ray.hit_normal = transforms[ray.id].mul_3D_vector(`
485	`normalize(normal0*(1 - ray.u - ray.v)`	482	`normalize(normal0*(1 - ray.u - ray.v)`
486	`+normal1*ray.u`	483	`+normal1*ray.u`
487	`+normal2*ray.v));`	484	`+normal2*ray.v));`
488		485
489	`ray.hit_pos = ray.origin + ray.direction*ray.dist;`	486	`ray.hit_pos = ray.origin + ray.direction*ray.dist;`
490	`/*`	487	`/*`
491	`Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);`	488	`Vec3i face = ray.hit_object->normals.face(ray.hit_face_id);`
492	`Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);`	489	`Vec3f normal0 = ray.hit_object->normals.vertex(face[0]);`
493	`Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);`	490	`Vec3f normal1 = ray.hit_object->normals.vertex(face[1]);`
494	`Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);`	491	`Vec3f normal2 = ray.hit_object->normals.vertex(face[2]);`
495	`ray.hit_normal = normalize(normal0*(1 - ray.u - ray.v)`	492	`ray.hit_normal = normalize(normal0*(1 - ray.u - ray.v)`
496	`+normal1*ray.u`	493	`+normal1*ray.u`
497	`+normal2*ray.v);`	494	`+normal2*ray.v);`
498	`ray.hit_pos = ray.origin + ray.direction*ray.dist;`	495	`ray.hit_pos = ray.origin + ray.direction*ray.dist;`
499	`*/`	496	`*/`
500	`return true;`	497	`return true;`

Subversion Repositories gelsvn

(root)/trunk/src/Geometry/BSPTree.cpp – Rev 428 → 429