ray_trace ray doesn't start at given position

[yportier]

Hi, are there any restrictions to the velocity.values array ?
I am perhaps trying to abuse the library with highly heterogeneous velocities and sometimes, the ray that is returned by ray_trace simply doesn't start at the given start_idx.

[malcolmw]

Hi, @yportier,

The algorithm should be unconditionally stable, so arbitrarily complex velocity fields are fine. The ray tracer uses the Euler method to integrate the gradient of the traveltime field in reverse direction. That means it starts from the end point and traces its way back to the source. Integration terminates when the ray takes an uphill step, meaning it has overshot the source location. There is, however, no guarantee that the tracer will terminate exactly at the source location. Error is accumulated along the integration path and so the tracer misses by some amount. The amount that it misses by can be reduced by decreasing the interval between nodes on the computational grid.

Take a look at https://github.com/malcolmw/pykonal/blob/master/jupyter/figure_6.ipynb for more on this.

- Malcolm

[yportier]

Hello Malcolm, Thanks for your feedback. I was not able to figure out my problem. Here is a colab notebook that shows my issue.

[malcolmw]

@yportier, it looks to me like the ray tracer is hitting the corner of one your "walls", taking an uphill step as a result, and terminating integration. What happens if you double the sampling density of your velocity model?

[yportier]

I tried setting intervals=(2, 2), here's the output:

[malcolmw]

What I mean is actually doubling the number of samples in your velocity model. So if you velocity model has node intervals (1, 1, 1) and values

[[1, 2, 3],
[2, 3, 4],
[3, 4, 5]]

it would become

[[1, 1, 2, 2, 3, 3],
[1, 1, 2, 2, 3, 3],
[2, 2, 3, 3, 4, 4],
[2, 2, 3, 3, 4, 4],
[3, 3, 4, 4, 5, 5],
[3, 3, 4, 4, 5, 5]]

with node intervals (1/2, 1/2, 1/2).

[yportier]

ok, sorry, I tried "doubling" the array and using intervals (0.5, 0.5, 0.5) as suggested, it only get a tad closer to the wall:

[malcolmw]

A couple of things:

1. Be careful when switching between node indices and physical coordinates.
2. I set your the velocity inside your walls to be 0.
3. I used a smoothing filter to smooth your velocity model.

See the code below for a working example of what I think you are trying to do.

import scipy.ndimage

speedmap[np.nonzero(speedmap == 1)] = 0
speedmap = scipy.ndimage.gaussian_filter(speedmap, 1/2)

# Computation

ORIGIN = [20, 80]
DESTIN = [80, 45]

src_idx = ORIGIN

speedmap_ = np.repeat(np.repeat(speedmap, 2, axis=0), 2, axis=1)

solver = pykonal.EikonalSolver(coord_sys='cartesian')
solver.velocity.min_coords = 0, 0, 0
solver.velocity.node_intervals = (0.5, 0.5, 0.5)
solver.velocity.npts = *speedmap_.shape, 1
solver.velocity.values = np.expand_dims(speedmap_, axis=-1)
src_idx = *src_idx, 0
solver.traveltime.values[src_idx] = 0
solver.unknown[src_idx] = False
solver.trial.push(*src_idx)
solver.solve()

# Display

fig, axes = plt.subplots(ncols=2, figsize=(16, 8))

ax = axes[0]

ax.pcolormesh(
    solver.vv.nodes[:, :, 0, 0], 
    solver.vv.nodes[:, :, 0, 1], 
    solver.vv.values[:, :, 0]
)
ax.scatter(*solver.vv.nodes[ORIGIN[0], ORIGIN[1], 0, :2])
ax.scatter(DESTIN[0], DESTIN[1])
ax.plot(path[:, 0], path[:, 1], color="red")

ax = axes[1]
ax.pcolormesh(
    solver.tt.nodes[:, :, 0, 0], 
    solver.tt.nodes[:, :, 0, 1], 
    solver.tt.values[:, :, 0]
)
ax.scatter(*solver.tt.nodes[ORIGIN[0], ORIGIN[1], 0, :2])
ax.scatter(DESTIN[0], DESTIN[1])
ax.plot(path[:, 0], path[:, 1], color="red")

[yportier]

Thanks!

I've been able to reproduce your output (it's only missing the definition of path, but I assumed it did not change). However, the origin is still not correct. It is not at 20, 80 but around 10, 40, so I am assuming it needs to be transformed according to the "densification" of the input. By doing that, it seems to be working for me (see below - I acknowledge your warning about confusing indices and coordinates, I'm just doing a dirty display to check the consistency of my path...).

Also, note that I get the following warning at execution:

/usr/local/lib/python3.7/dist-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/local/lib/python3.7/dist-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/local/lib/python3.7/dist-packages/numpy/lib/function_base.py:1096: RuntimeWarning: invalid value encountered in subtract out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n