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

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/Geometry/BSPTree.cpp – Rev 601 → 630