trace_ray crashing

[viktor76525]

Trying to run trace_ray following: https://malcolmw.github.io/pykonal-docs/api/classes/fields/ScalarField3D.html#pykonal.fields.ScalarField3D.trace_ray

pykonal 0.2.3b2

Using fractional node intervals doesn't work:

import numpy
import pykonal

ny, nz, nx = 100, 100, 100
solver = pykonal.EikonalSolver(coord_sys="cartesian")
solver.velocity.min_coords = 0, 0, 0
solver.velocity.node_intervals = 0.2, 0.1, 0.2
solver.velocity.npts = ny, nz, nx
solver.velocity.values = np.ones((ny, nz, nx))
src_idx = 50, 80, 50
solver.traveltime.values[src_idx] = 0
solver.unknown[src_idx] = False
solver.trial.push(*src_idx)
solver.solve()
solver.trace_ray(np.array([20, 0, 20], dtype=float))
array([[20.,  0., 20.]])

It just returns the last point, not the trace. Works with solver.node_intervals = 1, 1, 1

Might be useful to add to the documentation that you have to run solve() before trace_ray() otherwise I get:

/usr/lib/python3.7/site-packages/numpy/lib/function_base.py:1068: RuntimeWarning: invalid value encountered in subtract
  out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
/usr/lib/python3.7/site-packages/numpy/lib/function_base.py:1089: RuntimeWarning: invalid value encountered in subtract
  out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
/usr/lib/python3.7/site-packages/numpy/lib/function_base.py:1096: RuntimeWarning: invalid value encountered in subtract
  out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n

[malcolmw]

Hi, @viktor76525,

Thanks for pointing out this source of confusion.

The problem here arises from an implicit change from grid-node indices to Cartesian coordinates. When the source location is initialized on lines

src_idx = 50, 80, 50
solver.traveltime.values[src_idx] = 0
solver.unknown[src_idx] = False
solver.trial.push(*src_idx)

it is correctly done via references to grid-node indices. The coordinates passed to solver.trace_ray(), however, must be Cartesian coordinates. The code that you have places the source at the grid node with indices (50, 80, 50), but the Cartesian coordinates of the endpoint of the raypath being traced are (20., 0., 20.), which lies just outside the grid. Hence the raypath not being traced and the endpoint being returned.

I recommend using the PointSourceSolver class. It is more accurate, easier to use, and more flexible. It alleviates the need for the user to transition between grid-node indices and continuous coordinates when setting the source location and tracing rays. Everything is done in continuous coordinates after the grid is initialized.

And thank you for the suggestion regarding the needed call to solve() before calling trace_ray(). I will make a point to add this to the documentation.

Let me know if any of that is unclear or if you have any other questions. I am going to close this issue as "Not a Bug", but feel free to continue to post questions on this thread or open a new one if you find any bugs.

Thanks again, Malcolm