Fix DDA for rays aligned with grid edges and zero-direction rays

drjit/dda.py

@@ -177,8 +177,8 @@ def dda_fun(state: list, index: Array3i,
     t_max_v2 = dr.maximum(t_min_v, t_max_v)
 
     # Disable extent computation for dims where the ray direction is zero
-    t_min_v[inf_t] = -dr.inf
-    t_max_v[inf_t] =  dr.inf
+    t_max_v2[inf_t] = dr.inf
+    t_min_v2[inf_t] = -dr.inf
 
     # Reduce constraints to a single ray interval
     t_min = dr.maximum(dr.max(t_min_v2), 0)
@@ -187,6 +187,10 @@ def dda_fun(state: list, index: Array3i,
     # Only run the DDA algorithm if the interval is nonempty
     active = active & (t_max > t_min) & dr.isfinite(t_max) # type: ignore
 
+    # Deactivate rays that have zero direction along any axis 
+    # and whose origin along that axis is outside the grid bounds
+    active = active & dr.all(~inf_t | ((0 <= ray_o) & (ray_o <= grid_res_f)))
+
     # Advance the ray to the start of the interval
     ray_o = dr.fma(ray_d, t_min, ray_o)
     t_min, t_max = 0, t_max - t_min # type: ignore
@@ -203,7 +207,7 @@ def dda_fun(state: list, index: Array3i,
 
     # Step size to next interaction
     dt_v = dr.select(ray_d >= 0, dr.fma(-p0, rcp_d, rcp_d), -p0 * rcp_d)
-    dt_v[inf_t] = rcp_d
+    dt_v[inf_t] = dr.inf
 
     def body_fn(
         active: BoolT, state: StateT, dt_v: ArrayNfT, p0: ArrayNfT, pi: ArrayNiT, t_rem: Any,

tests/test_dda.py

@@ -40,6 +40,9 @@ def dda_bruteforce(
 
         for i in range(n):
             if ray_d[i] == 0:
+                if ray_o[i] < idx[i] or ray_o[i] > (idx[i] + 1):
+                    valid = False
+                    break
                 continue
 
             t_min_i = (idx[i] - ray_o[i]) / ray_d[i]
@@ -430,3 +433,61 @@ def test17_rev_ad_random(t):
     rng = m.PCG32(16)
     grad_val = rng.next_float32()*2-1
     check_grad(t, rng, vol, n_samples=1024, diff='rev', rtol=1e-4, grad_val=grad_val)
+
+@pytest.mark.parametrize("corner", range(8))
+@pytest.mark.parametrize("direction", range(6))
+def test18_edges(corner, direction):
+    """Test a ray that goes along the edges of the grid"""
+    ray_o = [corner&4 != 0, corner&2 != 0, corner&1 != 0]
+    ray_o = [2 * v - 1 for v in ray_o]
+
+    ray_d = [0, 0, 0]
+    ray_d[direction//2] = 2 * (direction&1) - 1
+
+    dda_check(
+        ray_o = ray_o,
+        ray_d = ray_d,
+        grid_res = (3, 3, 3),
+        grid_min = (-1, -1, -1),
+        grid_max = (1, 1, 1)
+    )
+
+@pytest.mark.parametrize("face", range(6))
+@pytest.mark.parametrize("direction", range(2))
+def test19_negative_zero_direction(face, direction):
+    """Test a ray with a negative zero direction"""
+
+    ray_o = [-0.0, -0.0, -0.0]
+    ray_o[face//2] = 2 * (face&1) - 1 
+
+    ray_d = [-0.0, -0.0, -0.0]
+    ray_d[face//2] = 2 * direction - 1
+
+    dda_check(
+        ray_o = ray_o,
+        ray_d = ray_d,
+        grid_res = (3, 3, 3),
+        grid_min = (-1, -1, -1),
+        grid_max = (1, 1, 1)
+    )
+
+@pytest.mark.parametrize("face", range(6))
+@pytest.mark.parametrize("direction", range(4))
+def test20_zero_direction_no_hit(face, direction):
+    """Test a ray with a zero direction that does not hit the grid"""
+    eps = 1e-6
+
+    ray_o = [0, 0, 0]
+    ray_o[face//2] = 2 * (face&1) - 1 + (2 * (face&1) - 1) * eps
+
+    ray_d = [0, 0, 0]
+    ray_d[(face//2 + 1) % 3] = 2 * (direction & 2 != 0) - 1
+    ray_d[(face//2 + 2) % 3] = 2 * (direction & 1 != 0) - 1
+
+    dda_check(
+        ray_o = ray_o,
+        ray_d = ray_d,
+        grid_res = (3, 3, 3),
+        grid_min = (-1, -1, -1),
+        grid_max = (1, 1, 1)
+    )