Graphene can't handle axis aligned ray directions

[jadahl]

Experienced behavior

Using graphene_ray_intersect_box() using a ray with a direction that is axis aligned (one or more components is 0.0), the formula used to calculate whether the ray intersects with the box fails.

Expected behavior

The ray/box intersection should be correctly detected.

Steps to reproduce

Reproducer:

#include <graphene.h>
#include <stdio.h>

int
main (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;

  /* Box far to the top right */

  graphene_point3d_init (&min, 41.843132, 27.356903, -50.368336);
  graphene_point3d_init (&max, 51.698078, 29.080172, -50.368336);
  graphene_box_init (&box, &min, &max);

  /* Ray from (0, 0, 0) along the Y axis "upwards" *NOT* hitting the above box
   */

  graphene_point3d_init (&origin, 0, 0, 0);
  graphene_vec3_init (&direction, 0, 0.495176, -0.868793);
  graphene_ray_init (&ray, &origin, &direction);

  printf("1, intersects (exact): %d (should be 0)\n",
	 graphene_ray_intersects_box (&ray, &box));

  /* Nudged variant of the above ray */

  graphene_vec3_init (&direction, 0 + 0.0001, 0.495176, -0.868793);
  graphene_ray_init (&ray, &origin, &direction);
  printf("2, intersects (nudged): %d (should be 0)\n",
	 graphene_ray_intersects_box (&ray, &box));

  /* Box on the center top */

  graphene_point3d_init (&min, -5.654480, 27.356903, -50.368336);
  graphene_point3d_init (&max, 5.654475, 29.080172, -50.368336);
  graphene_box_init (&box, &min, &max);

  /* Ray from (0, 0, 0) along the Y axis "upwards" hitting the above box */

  graphene_point3d_init (&origin, 0, 0, 0);
  graphene_vec3_init (&direction, 0, 0.495176, -0.868793);
  graphene_ray_init (&ray, &origin, &direction);

  printf("3, intersects (exact): %d (should be 1)\n",
	 graphene_ray_intersects_box (&ray, &box));

  /* Nudged variant of the above ray */

  graphene_vec3_init (&direction, 2 * FLT_EPSILON, 0.495176, -0.868793);
  graphene_ray_init (&ray, &origin, &direction);
  printf("4, intersects (nudged): %d (should be 1)\n",
	 graphene_ray_intersects_box (&ray, &box));

  return 0;
}

With SSE2, the above program prints:

1, intersects (exact): 1 (should be 0)
2, intersects (nudged): 0 (should be 0)
3, intersects (exact): 1 (should be 1)
4, intersects (nudged): 1 (should be 1)

and with SSE2 turned off, it prints:

1, intersects (exact): 0 (should be 0)
2, intersects (nudged): 0 (should be 0)
3, intersects (exact): 0 (should be 1)
4, intersects (nudged): 1 (should be 1)

Operating system in use

Fedora 33.

SIMD implementation in use

See above.

[GeorgesStavracas]

I've looked at it a bit, and found an interesting oddity. In graphene_ray_intersect_box(), at line 492, the inv_dir.value = graphene_simd4f_reciprocal (r->direction.value) call results in a vec3 containing: -nan, 2.0194845200, -1.1510224342. Everything else derails from there.

[GeorgesStavracas]

I think the bug here is that inv_dir.x should have resulted in +∞ instead of NaN

[jadahl]

Indeed. I spent some time trying to work around that, and come up with a SSE 4.1 implementation that turned that -nan into FLT_MAX, which made it work on my CPU; but it wasn't very portable, and my aarch64 board isn't booting so I didn't look further into the neon variant.

Also, not sure requiring SSE4.1 is an option.

[jadahl]

Fwiw, here is the SSE4.1 variant:

#  define graphene_simd4f_reciprocal(v) \
  (__extension__ ({ \
    const graphene_simd4f_t __two = graphene_simd4f_init (2.0f, 2.0f, 2.0f, 2.0f); \
    const graphene_simd4f_t __inf = graphene_simd4f_init (FLT_MAX, FLT_MAX, FLT_MAX, FLT_MAX); \
    graphene_simd4f_t __mask = _mm_cmpeq_ps(_mm_set1_ps(0.0), (v)); \
    graphene_simd4f_t __s = _mm_andnot_ps (__mask, _mm_rcp_ps ((v))); \
    graphene_simd4f_t __r = _mm_blendv_ps (__s, __inf, __mask); \
    graphene_simd4f_mul (__r, graphene_simd4f_sub (__two, graphene_simd4f_mul ((v), __r))); \
  }))

[jadahl]

Also not sure if graphene_simd4f_t __mask = _mm_cmpeq_ps(_mm_set1_ps(0.0), (v)); is enough; as it does == 0.0 and not withit FLT_EPSILON * 2 which was needed in my tests to make all paths pass (SIMD & non-SIMD).

[GeorgesStavracas]

So, based on your code for Mutter, this seems to work here:

@@ -485,27 +486,46 @@ 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];
+
+#define V(v) (graphene_approx_val (d[v], 0.f) ? 2 * FLT_EPSILON : d[v])
+  graphene_vec3_to_float (&r->direction, d);
+  graphene_vec3_init (&safe_direction, V (0), V (1), V (2));
+#undef V
 
   /* FIXME: Needs a graphene_vec3_reciprocal() */
-  inv_dir.value = graphene_simd4f_reciprocal (r->direction.value);
+  inv_dir.value = graphene_simd4f_reciprocal (safe_direction.value);

[jadahl]

Here seems like a better place for this code than in mutter.

[jadahl]

Ping @ebassi .

[ebassi]

The solution from @GeorgesStavracas looks good to me. I wonder if he wants to open a PR, or if I should do that.

[jadahl]

I went and created #217 that works slightly the same, except it uses a static inline function and makes sure the direction stays on the correct side of the axis.