Simple BVH (#1)

* Wtf * Fix bvh and add header deps in ninja build
raywan · May 23, 2019 · 344273d · 344273d
1 parent 53cc1ed
commit 344273d
Show file tree

Hide file tree

Showing 9 changed files with 221 additions and 152 deletions.
diff --git a/TODO.txt b/TODO.txt
@@ -1,8 +1,10 @@
 Build:
-  [ ] Create windows target for ninja build script generator (configure.py)
+  [x] Create windows target for ninja build script generator (configure.py)
   [x] build.bat for Windows
   [x] Create Windows entry point i.e. win32_main.cpp
 
 
-[ ] Tile size independent multithreaded rendering
-      - Height and width are must be multiples of the tile size currently
+Multithreading:
+  [ ] Fix off by one error in tile compositing @BUG
+  [ ] Tile size independent multithreaded rendering
+        Height and width are must be multiples of the tile size currently
diff --git a/configure.py b/configure.py
@@ -32,6 +32,9 @@
 
     n.variable(key='ninja_required_version', value='1.9')
 
+    if sys.platform == 'win32':
+        n.variable(key='msvc_deps_prefix', value='Note: including file:')
+
     n.variable(key='cc', value=CC)
     n.variable(key='cflags', value=' '.join(CFLAGS))
     n.variable(key='project_name', value=PROJECT_NAME)
@@ -46,7 +49,7 @@
     ############################################################################
 
     if sys.platform == 'win32':
-        n.rule('compile', command='$cc $cflags -c $in -Fo$out')
+        n.rule('compile', command='$cc /showIncludes $cflags -c $in -Fo$out', deps='msvc')
     else:
         n.rule('compile',
                command='$cc $cflags -c $in -o $out -MMD -MF $out.d',

diff --git a/src/bvh.cpp b/src/bvh.cpp
@@ -1,8 +1,12 @@
 #include "bvh.h"
 #include "mesh.h"
 #include "world.h"
+#include "metrics.h"
+#include <rw/rw_th.h>
 
-#define BOUND_MESH 1
+#define BOUND_MESH 0
+
+#define DEBUG 0
 
 uint32_t xor_shift_u32(uint32_t *state) {
   uint32_t x = *state;
@@ -13,24 +17,38 @@ uint32_t xor_shift_u32(uint32_t *state) {
   return x;
 }
 
+BVHNode *bvh_leaf_node(int prim_idx, int axis, Rect3 bounds, int *total_nodes) {
+  BVHNode *node = (BVHNode *) malloc(sizeof(BVHNode));
+  node->split_axis = axis;
+  node->prim_idx = prim_idx;
+  node->bounds = bounds;
+  node->left = NULL;
+  node->right = NULL;
+  node->leaf = true;
+  node->id = (*total_nodes)++;
+  return node;
+}
+
 std::vector<BVHPrimitive> bvh_preprocess_world(World *w) {
   std::vector<BVHPrimitive> result;
+  int cur_idx = 0;
   // Spheres
   for (int i = 0; i < w->spheres.size(); i++) {
     BVHPrimitive prim;
-    prim.idx = i;
+    prim.idx = cur_idx++;
+    prim.p_idx = i;
     prim.type = PT_SPHERE;
     prim.bounds = sphere_get_bounds(&w->spheres[i]);
     prim.centroid = 0.5 * prim.bounds.min_p + 0.5f * prim.bounds.max_p;
-   result.push_back(prim);
+    result.push_back(prim);
   }
   // Meshes
   for (int i = 0; i < w->meshes.size(); i++) {
     // Loop through each triangle in the mesh
     int cur_v_idx = 0;
     Mesh *cur_mesh = w->meshes[i];
     BVHPrimitive prim;
-    prim.idx = i;
+    prim.mesh_idx = i;
     for (int j = 0; j < cur_mesh->f.size(); j++) {
       prim.f_idx = j;
       // Intialize the bounds with one point of the triangle
@@ -39,6 +57,7 @@ std::vector<BVHPrimitive> bvh_preprocess_world(World *w) {
         // Extend the bounding box using the rest of the points
         prim.bounds = rwm_r3_union_p(prim.bounds, cur_mesh->v[cur_mesh->v_idx[cur_v_idx++]]);
       }
+      prim.idx = cur_idx++;
       prim.type = PT_TRIANGLE;
       prim.centroid = 0.5 * prim.bounds.min_p + 0.5f * prim.bounds.max_p;
       result.push_back(prim);
@@ -48,53 +67,157 @@ std::vector<BVHPrimitive> bvh_preprocess_world(World *w) {
   return result;
 }
 
-BVHNode *bvh_recursive_build(std::vector<BVHPrimitive> *prims, int lo, int hi) {
-  // Calculate total bounds of the primitives
-  BVHNode *node = (BVHNode *) malloc(sizeof(BVHNode));
-  // NOTE(ray): I don't think i need to do this computation here
-  Rect3 total_bounds = (*prims)[lo].bounds;
-  for (int i = lo + 1; i < hi; i++) {
-    total_bounds = rwm_r3_union(total_bounds, (*prims)[i].bounds);
+void print_prims(BVHPrimitive *prims, int n) {
+#if DEBUG
+  for (int i = 0; i < n; i++) {
+    printf("idx %d\n", prims[i].idx);
+    printf("\tp_idx %d\n", prims[i].p_idx);
+    printf("\tmesh_idx %d\n", prims[i].mesh_idx);
+    printf("\tf_idx %d\n", prims[i].f_idx);
+    printf("\ttype: %d\n", prims[i].type);
+    rwm_v3_printf("\tbounds.min_p", &prims[i].bounds.min_p);
+    rwm_v3_printf("\tbounds.max_p", &prims[i].bounds.max_p);
+    rwm_v3_printf("\tcentroid", &prims[i].centroid);
+    puts("");
+  }
+#endif
+}
+
+void bvhn_print(BVHNode *node) {
+#if DEBUG
+  printf("node addr: %p\n", node);
+  printf("node.leaf: %d\n", node->leaf);
+  if (node->leaf) {
+    printf("node.prim_idx: %d\n", node->prim_idx);
+  } else {
+    printf("node.left: %p\n", node->left);
+    printf("node.right: %p\n", node->right);
   }
+  rwm_v3_printf("node.bounds.min_p", &node->bounds.min_p);
+  rwm_v3_printf("node.bounds.max_p", &node->bounds.max_p);
+  puts("");
+#endif
+}
 
-  int n_primitives = hi - lo;
+int x_axis_comp(const void *a, const void *b) {
+  BVHPrimitive *left = (BVHPrimitive *) a;
+  BVHPrimitive *right = (BVHPrimitive *) b;
+  if (left->bounds.min_p.x - right->bounds.min_p.x < 0.0) return -1;
+  return 1;
+}
+
+int y_axis_comp(const void *a, const void *b) {
+  BVHPrimitive *left = (BVHPrimitive *) a;
+  BVHPrimitive *right = (BVHPrimitive *) b;
+  if (left->bounds.min_p.y - right->bounds.min_p.y < 0.0) return -1;
+  return 1;
+}
 
-  if (n_primitives == 1) {
-    node->leaf = true;
-    node->left = NULL;
+int z_axis_comp(const void *a, const void *b) {
+  BVHPrimitive *left = (BVHPrimitive *) a;
+  BVHPrimitive *right = (BVHPrimitive *) b;
+  if (left->bounds.min_p.z - right->bounds.min_p.z < 0.0) return -1;
+  return 1;
+}
+
+
+BVHNode *bvh_recursive_build(BVHPrimitive *prims, int n, int *total_nodes) {
+  // Calculate total bounds of the primitives
+  // printf("n == %d\n", n);
+  BVHNode *node = (BVHNode *) malloc(sizeof(BVHNode));
+  node->leaf = false;
+  node->id = (*total_nodes)++;
+  static uint32_t rng_state = 4;
+  int axis = xor_shift_u32(&rng_state) % 3;
+  // printf("Sorting on axis: %d\n", axis);
+  if (axis == 0) {
+    qsort(prims, n, sizeof(BVHPrimitive), x_axis_comp);
+  } else if (axis == 1) {
+    qsort(prims, n, sizeof(BVHPrimitive), y_axis_comp);
+  } else {
+    qsort(prims, n, sizeof(BVHPrimitive), z_axis_comp);
+  }
+  #if 1
+  if (n == 1) {
+    node->left = bvh_leaf_node(prims[0].idx, axis, prims[0].bounds, total_nodes);
     node->right = NULL;
-    node->bounds = total_bounds;
-    node->prim_idx = lo;
+    bvhn_print(node->left);
+  } else if (n == 2) {
+    node->left = bvh_leaf_node(prims[0].idx, axis, prims[0].bounds, total_nodes);
+    bvhn_print(node->left);
+    node->right = bvh_leaf_node(prims[1].idx, axis, prims[1].bounds, total_nodes);
+    bvhn_print(node->right);
   } else {
-    // Partition based on the centroid
-    Rect3 centroid_bounds = rwm_r3_init_p((*prims)[lo].centroid);
-    for (int i = lo + 1; i < hi; i++) {
-      centroid_bounds = rwm_r3_union_p(centroid_bounds, (*prims)[i].centroid);
-    }
-    int axis = rwm_r3_max_extent(centroid_bounds);
-    int mid = (hi - lo)/2;
-    if (centroid_bounds.max_p.e[axis] == centroid_bounds.min_p.e[axis]) {
-      node->leaf = true;
-      node->left = NULL;
-      node->right = NULL;
-      node->bounds = total_bounds;
-      node->prim_idx = lo;
-    } else {
-      node->left = bvh_recursive_build(prims, lo, mid);
-      node->right = bvh_recursive_build(prims, mid, hi);
-      node->bounds = rwm_r3_union(node->left->bounds, node->right->bounds);
-      node->split_axis = axis;
-    }
+    node->left = bvh_recursive_build(prims, n/2, total_nodes);
+    node->right = bvh_recursive_build(prims+n/2, n-n/2, total_nodes);
+  }
+  if (node->left && node->right) {
+    node->bounds = rwm_r3_union(node->left->bounds, node->right->bounds);
+  } else if (node->left) {
+    node->bounds = node->left->bounds;
+  } else if (node->right) {
+    node->bounds = node->right->bounds;
   }
+
+  bvhn_print(node);
+  #endif
   return node;
 }
 
 BVHNode *bvh_build(World *world) {
-  std::vector<BVHPrimitive> prims = bvh_preprocess_world(world);
-  BVHNode *root = bvh_recursive_build(&prims, 0, prims.size());
+  world->bvh_prims = bvh_preprocess_world(world);
+  std::vector<BVHPrimitive> bvh_prims_work_copy = world->bvh_prims;
+  printf("total prims: %llu\n", world->bvh_prims.size());
+  print_prims(world->bvh_prims.data(), world->bvh_prims.size());
+  int total_nodes = 0;
+  BVHNode *root = bvh_recursive_build(bvh_prims_work_copy.data(), world->bvh_prims.size(), &total_nodes);
+  printf("total nodes: %d\n", total_nodes);
   return root;
 }
 
-bool bvh_intersect(Ray *r, IntersectInfo *ii) {
+bool bvh_intersect(World *world, BVHNode *root, Ray *orig_ray, IntersectInfo *out_ii, Ray *out_r) {
+  if (!root) return false;
+  if (bb_intersect(&root->bounds, orig_ray, out_ii)) {
+    if (root->leaf) {
+      BVHPrimitive bvhp = world->bvh_prims[root->prim_idx];
+      switch (bvhp.type) {
+        case PT_SPHERE:
+          if (sphere_intersect(&(world->spheres[bvhp.p_idx]), out_r, out_ii)) {
+            rwth_atomic_add_i64((int64_t volatile *) &mtr_num_sphere_isect, 1);
+            out_ii->material = &world->sphere_materials[bvhp.p_idx];
+            return true;
+          }
+          break;
+        case PT_TRIANGLE: {
+          int f_idx = bvhp.f_idx;
+          Mesh *cur_mesh  = world->meshes[bvhp.mesh_idx];
+          Triangle triangle;
+          triangle.v0 = cur_mesh->v[cur_mesh->v_idx[f_idx * 3]];
+          triangle.v1 = cur_mesh->v[cur_mesh->v_idx[f_idx * 3 + 1]];
+          triangle.v2 = cur_mesh->v[cur_mesh->v_idx[f_idx * 3 + 2]];
+          if (triangle_intersect(&triangle, out_r, out_ii)) {
+            rwth_atomic_add_i64((int64_t volatile *) &mtr_num_triangle_isect, 1);
+            // Get surface properties: texture coordinates, vertex normal
+            Vec3 col = {triangle.u , triangle.v, triangle.w};
+            out_ii->color = col;
+            Vec2 uv0 = cur_mesh->uv[cur_mesh->uv_idx[f_idx * 3]];
+            Vec2 uv1 = cur_mesh->uv[cur_mesh->uv_idx[f_idx * 3 + 1]];
+            Vec2 uv2 = cur_mesh->uv[cur_mesh->uv_idx[f_idx * 3 + 2]];
+            out_ii->tex_coord = (triangle.w * uv0) + (triangle.u * uv1) + (triangle.v * uv2);
+            Vec3 n0 = cur_mesh->n[cur_mesh->n_idx[f_idx * 3]];
+            Vec3 n1 = cur_mesh->n[cur_mesh->n_idx[f_idx * 3 + 1]];
+            Vec3 n2 = cur_mesh->n[cur_mesh->n_idx[f_idx * 3 + 2]];
+            out_ii->normal = (triangle.w * n0) + (triangle.u * n1) + (triangle.v * n2);
+            out_ii->material = &world->mesh_materials[bvhp.mesh_idx];
+            return true;
+          }
+        } break;
+      }
+    } else {
+      bool left_result = bvh_intersect(world, root->left, orig_ray, out_ii, out_r);
+      bool right_result = bvh_intersect(world, root->right, orig_ray, out_ii, out_r);
+      return left_result || right_result;
+    }
+  }
   return false;
 }
diff --git a/src/bvh.h b/src/bvh.h
@@ -7,6 +7,7 @@
 struct World;
 
 struct BVHNode {
+  int id;
   Rect3 bounds;
   int prim_idx;
   int split_axis;
@@ -17,14 +18,17 @@ struct BVHNode {
 
 struct BVHPrimitive {
   int idx;
-  int f_idx; // for meshes
+  int p_idx;
   PrimitiveType type;
   Rect3 bounds;
   Point3 centroid;
+  // For meshes
+  int mesh_idx;
+  int f_idx; 
 };
 
 std::vector<BVHPrimitive> bvh_preprocess_world(World *w);
 BVHNode *bvh_build(World *world);
-bool bvh_intersect(Ray *r, IntersectInfo *ii);
+bool bvh_intersect(World *world, BVHNode *root, Ray *orig_ray, IntersectInfo *out_ii, Ray *out_r);
 
 #endif
diff --git a/src/primitive.cpp b/src/primitive.cpp
@@ -172,19 +172,33 @@ bool mesh_intersect(Mesh *mesh, Ray *r, IntersectInfo *out_ii) {
 }
 
 bool bb_intersect(Rect3 *bb, Ray *r, IntersectInfo *out_ii) {
+  #if 0
+  // TODO(ray): The optimized version doesn't work properly...
+  // Something is wrong. Come back to this later.
   float t_min = (bb->p[r->sign[0]].x - r->origin.x) * r->inv_dir.x;
   float t_max = (bb->p[1-r->sign[0]].x - r->origin.x) * r->inv_dir.x;
   float t_ymin = (bb->p[r->sign[1]].y - r->origin.y) * r->inv_dir.y;
   float t_ymax = (bb->p[1-r->sign[1]].y - r->origin.y) * r->inv_dir.y;
-  if (t_min > t_ymax || t_ymin > t_max) return false;
-  t_min = MAX(t_min, t_ymin);
-  t_max = MIN(t_max, t_ymax);
+  if (t_min > t_ymin || t_ymin > t_max) return false;
+  if (t_ymin > t_min) t_min = t_ymin;
+  if (t_ymax < t_max) t_max = t_ymax;
   float t_zmin = (bb->p[r->sign[2]].z - r->origin.z) * r->inv_dir.z;
   float t_zmax = (bb->p[1-r->sign[2]].z - r->origin.z) * r->inv_dir.z;
   if (t_min > t_zmax || t_zmin > t_max) return false;
-  t_min = MAX(t_min, t_zmin);
-  t_max = MIN(t_max, t_zmax);
-  r->at_t = t_min;
-
+  if (t_zmin > t_min) t_min = t_zmin;
+  if (t_zmax < t_max) t_max = t_zmax;
+  return (t_min < FLT_MAX) && (t_max > 0);
+#else
+  float t0 = 0;
+  float t1 = FLT_MAX;
+  for (int i = 0; i < 3; i++) {
+    float t_near = (bb->min_p.e[i] - r->origin.e[i]) * r->inv_dir.e[i];
+    float t_far = (bb->max_p.e[i] - r->origin.e[i]) * r->inv_dir.e[i];
+    if (t_near > t_far) std::swap(t_near, t_far);
+    t0 = t_near > t0 ? t_near : t0;
+    t1 = t_far < t1 ? t_far : t1;
+    if (t0 > t1) return false;
+  }
+#endif
   return true;
 }