Replace ray triangle intersection with MOLLER TRUMBORE

Author: tigert1998

ray-tracing-cuda/ray.cu

@@ -10,6 +10,6 @@ __device__ __host__ Ray::Ray(vec3 position, vec3 direction)
   direction_ = normalize(direction);
 }
 
-__device__ __host__ vec3 Ray::position() const { return position_; }
+__device__ __host__ const vec3 Ray::position() const { return position_; }
 
-__device__ __host__ vec3 Ray::direction() const { return direction_; }
+__device__ __host__ const vec3 Ray::direction() const { return direction_; }

ray-tracing-cuda/ray.cuh

@@ -6,8 +6,8 @@ class Ray {
  public:
   __device__ __host__ Ray();
   __device__ __host__ Ray(glm::vec3 position, glm::vec3 direction);
-  __device__ __host__ glm::vec3 position() const;
-  __device__ __host__ glm::vec3 direction() const;
+  __device__ __host__ const glm::vec3 position() const;
+  __device__ __host__ const glm::vec3 direction() const;
 
  private:
   glm::vec3 position_, direction_;

ray-tracing-cuda/ray_tracing.cu

@@ -6,20 +6,16 @@
 #include "utils.cuh"
 
 using glm::vec3;
-using std::pair;
 
 constexpr int TRACE_DEPTH_LIMIT = 10;
 
 __device__ glm::vec3 Trace(HitableList *world, Ray ray, curandState *states) {
-  struct Layer {
-    vec3 emitted;
-    vec3 attenuation;
-  };
-
   Layer storage[TRACE_DEPTH_LIMIT];
 
   int storage_size = 0;
 
+  bool should_debug = false;
+
   vec3 result;
   for (int depth = 0;; depth++) {
     HitRecord record;
@@ -39,11 +35,18 @@ __device__ glm::vec3 Trace(HitableList *world, Ray ray, curandState *states) {
       result = emitted;
       break;
     }
+    storage[storage_size].t = record.t;
+    storage[storage_size].pos = ray.position();
+    storage[storage_size].target = hit_point;
     storage[storage_size].emitted = emitted;
     storage[storage_size++].attenuation = attenuation;
     ray = reflection;
   }
 
+  if (should_debug) {
+    DebugTracePath(storage, storage_size);
+  }
+
   for (int i = storage_size - 1; i >= 0; i--) {
     result = storage[i].emitted + storage[i].attenuation * result;
   }

ray-tracing-cuda/ray_tracing.cuh

@@ -6,6 +6,14 @@
 #include "hitable_list.cuh"
 #include "ray.cuh"
 
+struct Layer {
+  glm::vec3 pos;
+  double t;
+  glm::vec3 target;
+  glm::vec3 emitted;
+  glm::vec3 attenuation;
+};
+
 __device__ glm::vec3 Trace(HitableList *world, Ray ray);
 
 __global__ void RayTracing(HitableList *world, Camera *camera, int height,

ray-tracing-cuda/utils.cu

@@ -1,15 +1,20 @@
 #ifndef STB_IMAGE_WRITE_IMPLEMENTATION
 #define STB_IMAGE_WRITE_IMPLEMENTATION
 #endif
+#ifndef GLM_ENABLE_EXPERIMENTAL
+#define GLM_ENABLE_EXPERIMENTAL
+#endif
 
 #include <glog/logging.h>
 #include <stb_image_write.h>
 
 #include <fstream>
+#include <glm/gtx/intersect.hpp>
 #include <string>
 
 #include "hitable.cuh"
 #include "ray.cuh"
+#include "ray_tracing.cuh"
 #include "utils.cuh"
 
 using glm::cross;
@@ -22,26 +27,38 @@ __global__ void CudaRandomInit(uint64_t seed, curandState *state) {
   curand_init(seed, idx, 0, &state[idx]);
 }
 
-__device__ bool TriangleHit(vec3 p[3], const Ray &ray, double t_from,
+__device__ bool TriangleHit(const vec3 p[3], const Ray &ray, double t_from,
                             double t_to, double *out_t, vec3 *out_normal) {
-  vec3 n = normalize(cross(p[2] - p[0], p[1] - p[2]));
-  double den = dot(ray.direction(), n);
-  double num = dot(p[0] - ray.position(), n);
-  double t = num / den;
-  if (std::isnan(t) || std::isinf(t)) return false;
-  if (t < t_from || t > t_to) return false;
-  vec3 hit_point = ray.position() + (float)t * ray.direction();
-  {
-    double dot_values[3];
-    for (int i = 0; i < 3; i++) {
-      vec3 a = p[i], b = p[(i + 1) % 3];
-      dot_values[i] = dot(cross(a - hit_point, b - hit_point), n);
-    }
-    for (int i = 0; i < 3; i++)
-      if (dot_values[i] * dot_values[(i + 1) % 3] < 0) return false;
+  //  MOLLER_TRUMBORE
+  double eps = 1e-7;
+  vec3 v0v1 = p[1] - p[0];
+  vec3 v0v2 = p[2] - p[0];
+  vec3 pvec = cross(ray.direction(), v0v2);
+  float det = dot(v0v1, pvec);
+
+  // ray and triangle are parallel if det is close to 0
+  if (fabs(det) < eps) return false;
+
+  float invDet = 1 / det;
+
+  vec3 tvec = ray.position() - p[0];
+  float u = dot(tvec, pvec) * invDet;
+  if (u < 0 || u > 1) return false;
+
+  vec3 qvec = cross(tvec, v0v1);
+  float v = dot(ray.direction(), qvec) * invDet;
+  if (v < 0 || u + v > 1) return false;
+
+  float t = dot(v0v2, qvec) * invDet;
+
+  if (!(t_from <= t && t <= t_to)) {
+    return false;
   }
+
   *out_t = t;
-  *out_normal = dot(n, ray.direction()) < 0 ? n : -n;
+  vec3 n = normalize(cross(v0v1, v0v2));
+  *out_normal = dot(ray.direction(), n) < 0 ? n : -n;
+
   return true;
 }
 
@@ -55,3 +72,16 @@ void WriteImage(std::vector<vec3> &pixels, int height, int width,
   LOG(INFO) << "Writing image to " << path << "...";
   stbi_write_jpg(path.c_str(), width, height, 3, data.data(), 100);
 }
+
+__device__ void DebugTracePath(Layer *layers, int n) {
+  for (int i = 0; i < n; i++) {
+    printf(
+        "[%d] (%.3f, %.3f, %.3f) -> (%.3f, %.3f, %.3f) with t: %.3f, "
+        "attenuation: (%.3f, %.3f, %.3f), emitted: (%.3f, %.3f, %.3f)\n",
+        i, layers[i].pos[0], layers[i].pos[1], layers[i].pos[2],
+        layers[i].target[0], layers[i].target[1], layers[i].target[2],
+        layers[i].t, layers[i].attenuation[0], layers[i].attenuation[1],
+        layers[i].attenuation[2], layers[i].emitted[0], layers[i].emitted[1],
+        layers[i].emitted[2]);
+  }
+}

ray-tracing-cuda/utils.cuh

@@ -23,7 +23,7 @@ __global__ void CudaRandomInit(uint64_t seed, curandState *state);
 void WriteImage(std::vector<glm::vec3> &pixels, int height, int width,
                 const std::string &path);
 
-__device__ bool TriangleHit(glm::vec3 p[3], const Ray &ray, double t_from,
+__device__ bool TriangleHit(const glm::vec3 p[3], const Ray &ray, double t_from,
                             double t_to, double *out_t, glm::vec3 *out_normal);
 
 template <typename T>
@@ -48,3 +48,6 @@ __device__ void QuickSort(
   if (j > low) QuickSort(array, low, j, comp);
   if (i < high) QuickSort(array, i, high, comp);
 }
+
+struct Layer;
+__device__ void DebugTracePath(Layer *layers, int n);