Is there an example demonstrating how to compute shortest-traveltime paths?

[r03ert0]

Maybe code to reproduce the figure in the README.md? thank you!

[malcolmw]

Hi, @r03ert0,

Try out the code below to generate a figure similar to the attached.

Out of curiosity, what is your intended use case? Your profile suggests you are working on something different than the typical seismological applications I originally designed the tool for.

import matplotlib.pyplot as plt
import numpy as np
import pykonal
import scipy.ndimage

n_srcs = 12

velocity = pykonal.fields.ScalarField3D(coord_sys="cartesian")
velocity.min_coords = 0, 0, 0
velocity.node_intervals = 0.12524461839530332, 0.12598425196850394, 0.1
velocity.npts = 512, 128, 1

for iy in range(velocity.npts[1]):
    depth = velocity.nodes[0, iy, 0, 1]
    velocity.values[:, iy] = 3.5 + 0.1 * depth

velocity.values += 20 * np.random.randn(*velocity.npts)
velocity.values = scipy.ndimage.gaussian_filter(velocity.values, 8)

velocity.values[116:204, 32:96] += 2
velocity.values[308:396, 32:96] -= 2
velocity.values = scipy.ndimage.gaussian_filter(velocity.values, 2)

velocity.values = np.clip(velocity.values, 1, 8)

plt.close("all")
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

qmesh = ax.pcolormesh(
    velocity.nodes[:, :, 0, 0],
    velocity.nodes[:, :, 0, 1],
    velocity.values[:, :, 0],
    cmap=plt.get_cmap("nipy_spectral_r"),
    vmin=1,
    vmax=8
)
cbar = fig.colorbar(qmesh, ax=ax, orientation="horizontal")
cbar.set_label("Velocity (km/s)")
cbar.ax.xaxis.set_minor_locator(plt.MultipleLocator(0.25))

nodes = velocity.nodes

for isrc in range(n_srcs):
    solver = pykonal.solver.PointSourceSolver(coord_sys="cartesian")
    solver.velocity.min_coords = velocity.min_coords
    solver.velocity.npts = velocity.npts
    solver.velocity.node_intervals = velocity.node_intervals
    solver.velocity.values = velocity.values
    solver.src_loc = np.random.rand(3) * [64, 16, 0]
    solver.solve()

    for rx in range(0, 65, 4):
        ax.scatter(
            rx, 
            0,
            marker=((0, 0), (-2, 4), (2, 4), (0, 0)),
            edgecolor="k",
            facecolor="w",
            linewidth=0.5,
            s=128,
            clip_on=False,
            zorder=3
        )
        ray = solver.trace_ray(np.array([rx, 0., 0.]))
        ax.plot(ray[:, 0], ray[:, 1], color="k", linewidth=0.5)

    ax.scatter(
        *solver.src_loc[:2],
        marker="*",
        edgecolor="k",
        facecolor="w",
        s=256,
        zorder=3
    )

ax.invert_yaxis()
ax.set_xlabel("Horizontal offset (km)")
ax.set_ylabel("Depth (km)")
ax.set_aspect(1)
ax.set_xticks(range(0, 65, 8))
ax.set_yticks(range(0, 17, 4))

[r03ert0]

hey! Thank you very much, that's awesome! 🚀 I was looking for a nice replacement for an old fast marching algorithm that I have been using for some time. I looked 1st into scikit-fmm, but yours is nicer : ) (in particular because I wanted to have those geodesic paths). I have two ideas where I would like to try fast marching. In the first one, I'm trying to match the vertices in two surfaces. One is the outer surface of a brain, and the 2nd is an internal surface computed using a sign distance function. I found a nice library for optimal transport, but sometimes the vertices follow incorrect paths. So, the first idea is to try to use fast marching to create a "geodesic" cost matrix for optimal transport. Should work! The 2nd idea is more fancy. I've been trying to come up with a model for the development of brain connections. Right now I have a quite nice model for the development of brain folding. Now, I would like to make fibres grow as the brain folds. Many fibres start to grow and travel across the brain at the time its surface folds. Fibres are constrained by the geometry, but also by the density and compression of the brain tissue. I think I can model all that as a cost function inside the brain model. Then, I would like to use your library to make fibres grow. I don't know what you see in your nice figure with the heat map, the stars and the paths: I see a cortex, with different densities, and neurones growing their axons 🤓 .

[malcolmw]

@r03ert0, Sounds very interesting! And I am glad that the code is useful to communities beyond my own!

The figure above is intended to represent the path that seismic waves travel between earthquakes (stars) and seismometers on the surface (triangles).

Let me know if you run into any trouble applying the code for your purposes and I will try my best to help. And keep me posted if you obtain any significant results; I will be very interested to see them.

Closing this issue as answered, but feel free to open another if you have other questions.

[r03ert0]

thank you @malcolmw ! I'll keep you posted :)