fix ray intersection bugs by avoiding dividing zero

include/cubvh/bounding_box.cuh

@@ -149,15 +149,15 @@ struct BoundingBox {
     }
 
     __host__ __device__ Eigen::Vector2f ray_intersect(Eigen::Ref<const Eigen::Vector3f> pos, Eigen::Ref<const Eigen::Vector3f> dir) const {
-        float tmin = (min.x() - pos.x()) / dir.x();
-        float tmax = (max.x() - pos.x()) / dir.x();
+        float tmin = safe_divide(min.x() - pos.x(), dir.x());
+        float tmax = safe_divide(max.x() - pos.x(), dir.x());
 
         if (tmin > tmax) {
             host_device_swap(tmin, tmax);
         }
 
-        float tymin = (min.y() - pos.y()) / dir.y();
-        float tymax = (max.y() - pos.y()) / dir.y();
+        float tymin = safe_divide(min.y() - pos.y() / dir.y());
+        float tymax = safe_divide(max.y() - pos.y() / dir.y());
 
         if (tymin > tymax) {
             host_device_swap(tymin, tymax);
@@ -175,8 +175,8 @@ struct BoundingBox {
             tmax = tymax;
         }
 
-        float tzmin = (min.z() - pos.z()) / dir.z();
-        float tzmax = (max.z() - pos.z()) / dir.z();
+        float tzmin = safe_divide(min.z() - pos.z(), dir.z());
+        float tzmax = safe_divide(max.z() - pos.z(), dir.z());
 
         if (tzmin > tzmax) {
             host_device_swap(tzmin, tzmax);

include/cubvh/common.h

@@ -95,4 +95,14 @@ __device__ __host__ Eigen::Vector3f fibonacci_dir(uint32_t i, const Eigen::Vecto
 	return cylindrical_to_dir(Eigen::Vector2f{fractf((i+epsilon) / (N_DIRS-1+2*epsilon) + offset.x()), fractf(i / GOLDEN_RATIO + offset.y())});
 }
 
+__host__ __device__ float safe_divide(float numerator, float denominator, float epsilon = 1e-6f) {
+	if (fabs(denominator) < epsilon) {
+		if (denominator <= 0)
+			return -(numerator / epsilon); 
+		else
+			return numerator / epsilon; 
+	}
+	return numerator / denominator;
+}
+
 }
\ No newline at end of file

include/cubvh/triangle.cuh

@@ -30,7 +30,7 @@ struct Triangle {
         Eigen::Vector3f rov0 = ro - a;
         n = v1v0.cross( v2v0 );
         Eigen::Vector3f q = rov0.cross( rd );
-        float d = 1.0f/rd.dot( n );
+        float d = safe_divide(1.0f, rd.dot(n));
         float u = d*-q.dot( v2v0 );
         float v = d* q.dot( v1v0 );
         float t = d*-n.dot( rov0 );