ray: Nudge ray direction of axis when intersecting

The formula used to calculate the inverse of the direction vector
doesn't handle the direction vector aligning with an axis. Depending on
the SIMD (or not SIMD) implementation used, a axis aligned vector would
either remain the same, or e.g. end up with NaN components messing up
any future calculations.

Fixing the math to handle this is non-trivial, so for now work around
this by nudging the direction vector slightly off axis so that it has a
better hand of hitting the right box even when the direction is axis
aligned.

Closes: #214

src/graphene-ray.c

@@ -467,6 +467,22 @@ graphene_ray_intersects_sphere (const graphene_ray_t    *r,
   return graphene_ray_intersect_sphere (r, s, NULL) != GRAPHENE_RAY_INTERSECTION_KIND_NONE;
 }
 
+static inline float
+nudge_off_axis (float v)
+{
+  if (graphene_approx_val (v, 0.f))
+    {
+      if (v < 0.f)
+        return -2 * FLT_EPSILON;
+      else
+        return 2 * FLT_EPSILON;
+    }
+  else
+    {
+      return v;
+    }
+}
+
 /**
  * graphene_ray_intersect_box:
  * @r: a #graphene_ray_t
@@ -485,10 +501,18 @@ graphene_ray_intersect_box (const graphene_ray_t *r,
                             const graphene_box_t *b,
                             float                *t_out)
 {
+  graphene_vec3_t safe_direction;
   graphene_vec3_t inv_dir;
+  float d[3];
+
+  graphene_vec3_to_float (&r->direction, d);
+  graphene_vec3_init (&safe_direction,
+                      nudge_off_axis (d[0]),
+                      nudge_off_axis (d[1]),
+                      nudge_off_axis (d[2]));
 
   /* FIXME: Needs a graphene_vec3_reciprocal() */
-  inv_dir.value = graphene_simd4f_reciprocal (r->direction.value);
+  inv_dir.value = graphene_simd4f_reciprocal (safe_direction.value);
 
   graphene_vec3_t inv_min;
   graphene_vec3_subtract (&(b->min), &r->origin, &inv_min);

tests/ray.c

@@ -166,6 +166,72 @@ ray_intersect_triangle (void)
                  NULL);
 }
 
+static void
+ray_intersects_box (void)
+{
+  graphene_point3d_t min;
+  graphene_point3d_t max;
+  graphene_point3d_t origin;
+  graphene_vec3_t direction;
+  graphene_box_t box;
+  graphene_ray_t ray;
+
+  /* Off center box */
+
+  graphene_point3d_init (&min, 41.843132f, 27.356903f, -50.368336f);
+  graphene_point3d_init (&max, 51.698078f, 29.080172f, -50.368336f);
+  graphene_box_init (&box, &min, &max);
+
+  /* Ray from (0, 0, 0) along an axis *NOT* hitting the above box
+   */
+
+  graphene_point3d_init (&origin, 0, 0, 0);
+  graphene_vec3_init (&direction, 0, 0.495176f, -0.868793f);
+  graphene_ray_init (&ray, &origin, &direction);
+
+  mutest_expect ("intersection kind should be NONE",
+                 mutest_int_value (graphene_ray_intersects_box (&ray, &box)),
+                 mutest_to_be_false,
+                 NULL);
+
+  /* Nudged variant of the above ray */
+
+  graphene_vec3_init (&direction, 0 + 0.0001f, 0.495176f, -0.868793f);
+  graphene_ray_init (&ray, &origin, &direction);
+
+  mutest_expect ("intersection kind should still be NONE",
+                 mutest_int_value (graphene_ray_intersects_box (&ray, &box)),
+                 mutest_to_be_false,
+                 NULL);
+
+  /* Centered box */
+
+  graphene_point3d_init (&min, -5.654480f, 27.356903f, -50.368336f);
+  graphene_point3d_init (&max, 5.654475f, 29.080172f, -50.368336f);
+  graphene_box_init (&box, &min, &max);
+
+  /* Ray from (0, 0, 0) along the axis hitting the above box */
+
+  graphene_point3d_init (&origin, 0, 0, 0);
+  graphene_vec3_init (&direction, 0, 0.495176f, -0.868793f);
+  graphene_ray_init (&ray, &origin, &direction);
+
+  mutest_expect ("intersection kind should be ENTER",
+                 mutest_int_value (graphene_ray_intersects_box (&ray, &box)),
+                 mutest_to_be_true,
+                 NULL);
+
+  /* Nudged variant of the above ray */
+
+  graphene_vec3_init (&direction, 2 * FLT_EPSILON, 0.495176f, -0.868793f);
+  graphene_ray_init (&ray, &origin, &direction);
+
+  mutest_expect ("intersection kind should still be ENTER",
+                 mutest_int_value (graphene_ray_intersects_box (&ray, &box)),
+                 mutest_to_be_true,
+                 NULL);
+}
+
 static void
 ray_suite (void)
 {
@@ -175,6 +241,7 @@ ray_suite (void)
   mutest_it ("can compute the closest point to a point on the ray", ray_closest_point_to_point);
   mutest_it ("can be transformed", ray_matrix_transform);
   mutest_it ("can intersect triangles", ray_intersect_triangle);
+  mutest_it ("can intersect on axis", ray_intersects_box);
 }
 
 MUTEST_MAIN (