fix nan error in raytracing when using honor_grid

keurfonluu · Oct 20, 2021 · 4d1018b · 4d1018b
1 parent 03de066
commit 4d1018b
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 9 deletions.
diff --git a/fteikpy/_fteik/_ray2d.py b/fteikpy/_fteik/_ray2d.py
@@ -50,19 +50,21 @@ def _ray2d(z, x, zgrad, xgrad, zend, xend, zsrc, xsrc, stepsize, max_step, honor
         if honor_grid:
             fac = shrink(pcur, delta, lower, upper)
             pcur -= fac * delta
+            pcur[0] = min(max(pcur[0], z[0]), z[-1])
+            pcur[1] = min(max(pcur[1], x[0]), x[-1])
 
             if fac < 1.0:
+                # Handle precision issues due to fac
+                pcur[0] = numpy.round(pcur[0], 8)
+                pcur[1] = numpy.round(pcur[1], 8)
+
                 i = numpy.searchsorted(z, pcur[0], side="right") - 1
                 j = numpy.searchsorted(x, pcur[1], side="right") - 1
                 lower[0] = z[max(i - 1, 0)] if pcur[0] == z[i] else z[i]
                 lower[1] = x[max(j - 1, 0)] if pcur[1] == x[j] else x[j]
                 upper[0] = z[i + 1]
                 upper[1] = x[j + 1]
 
-                # Handle precision issues due to fac
-                pcur[0] = numpy.round(pcur[0], 8)
-                pcur[1] = numpy.round(pcur[1], 8)
-
                 if (pcur != ray[count - 1]).any():
                     ray[count] = pcur.copy()
                     count += 1
@@ -75,6 +77,9 @@ def _ray2d(z, x, zgrad, xgrad, zend, xend, zsrc, xsrc, stepsize, max_step, honor
 
         else:
             pcur -= delta
+            pcur[0] = min(max(pcur[0], z[0]), z[-1])
+            pcur[1] = min(max(pcur[1], x[0]), x[-1])
+
             ray[count] = pcur.copy()
             count += 1
 

diff --git a/fteikpy/_fteik/_ray3d.py b/fteikpy/_fteik/_ray3d.py
@@ -74,8 +74,16 @@ def _ray3d(
         if honor_grid:
             fac = shrink(pcur, delta, lower, upper)
             pcur -= fac * delta
+            pcur[0] = min(max(pcur[0], z[0]), z[-1])
+            pcur[1] = min(max(pcur[1], x[0]), x[-1])
+            pcur[2] = min(max(pcur[2], y[0]), y[-1])
 
             if fac < 1.0:
+                # Handle precision issues due to fac
+                pcur[0] = numpy.round(pcur[0], 8)
+                pcur[1] = numpy.round(pcur[1], 8)
+                pcur[2] = numpy.round(pcur[2], 8)
+
                 i = numpy.searchsorted(z, pcur[0], side="right") - 1
                 j = numpy.searchsorted(x, pcur[1], side="right") - 1
                 k = numpy.searchsorted(y, pcur[2], side="right") - 1
@@ -86,11 +94,6 @@ def _ray3d(
                 upper[1] = x[j + 1]
                 upper[2] = y[k + 1]
 
-                # Handle precision issues due to fac
-                pcur[0] = numpy.round(pcur[0], 8)
-                pcur[1] = numpy.round(pcur[1], 8)
-                pcur[2] = numpy.round(pcur[2], 8)
-
                 if (pcur != ray[count - 1]).any():
                     ray[count] = pcur.copy()
                     count += 1
@@ -103,6 +106,10 @@ def _ray3d(
 
         else:
             pcur -= delta
+            pcur[0] = min(max(pcur[0], z[0]), z[-1])
+            pcur[1] = min(max(pcur[1], x[0]), x[-1])
+            pcur[2] = min(max(pcur[2], y[0]), y[-1])
+
             ray[count] = pcur.copy()
             count += 1