Dijkstra and Fast Marching¶
Propagate a distance front through a weighted domain, first with a graph update and then with an eikonal update. Extract paths from the resulting distance map to see how local travel costs determine global geodesics.
Run this tour¶
Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.
# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys
working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
checkout = working / "numerical-tours-support"
if not checkout.exists():
subprocess.run(
[
"git",
"clone",
"--depth",
"1",
"--branch",
"master",
"https://github.com/gpeyre/numerical-tours.git",
str(checkout),
],
check=True,
)
python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
importlib.util.find_spec(name) is None
for name in [
"numpy",
"scipy",
"matplotlib",
"skimage",
"sklearn",
"pywt",
"ipywidgets",
"cvxpy",
"skfmm",
"autograd",
"progressbar",
"celer",
]
):
subprocess.run(
[sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
plt.rcParams.update(
{
"figure.figsize": (8, 4),
"figure.dpi": 100,
"axes.spines.top": False,
"axes.spines.right": False,
"font.size": 11,
"image.cmap": "gray",
}
)
%matplotlib inline
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tours details the implementations of Dijkstra and Fast Marching algorithms in 2-D.
The implementation are performed in Matlab, and are hence quite slow.
from nt_toolbox.general import np, plt
from nt_toolbox.signal import bilinear_interpolate, grad, imageplot, np, plt
from nt_toolbox.graph import np, perform_dijstra_fm, plt
import numpy as np
import matplotlib.pyplot as plt
from numpy import (
abs,
add,
arange,
array,
cos,
empty,
exp,
extract,
full,
gradient,
imag,
inf,
linspace,
meshgrid,
mod,
ones,
pi,
real,
sqrt,
stack,
sum,
transpose,
zeros,
)
from matplotlib.pyplot import (
axis,
clf,
close,
gray,
jet,
matplotlib,
np,
plot,
set_cmap,
step,
subplot,
time,
)
from nt_toolbox.general import np, plt
from nt_toolbox.signal import bilinear_interpolate, grad, imageplot, np, plt
from nt_toolbox.graph import np, perform_dijstra_fm, plt
def exo1(x0, W):
"""
Implement the Dijkstra algorithm by iterating these step while the
stack |I| is non empty.
Display from time to time the front that propagates.
"""
n = W.shape[0]
pstart = transpose(array([x0]))
[D, Dsvg, Ssvg] = perform_dijstra_fm(W, pstart, inf, "dijstr", "sym", n * 6)
plt.clf()
for i in arange(0, 4):
plt.subplot(2, 2, i + 1)
d = Dsvg[:, :, i]
d[d == inf] = 0
imageplot(d)
set_cmap("jet")
return D
def exo2(x0, W):
"""
Implement the FM algorithm by iterating these step while the
stack |I| is non empty.
Display from time to time the front that propagates.
"""
n = W.shape[0]
pstart = transpose(array([x0]))
[D, Dsvg, Ssvg] = perform_dijstra_fm(W, pstart, inf, "fm", "sym", n * 6)
clf
for i in arange(0, 4):
subplot(2, 2, i + 1)
d = Dsvg[:, :, i]
d[d == np.inf] = 0
imageplot(d)
set_cmap("jet")
return D
def exo3(x0, W):
"""
Compute the distance map to these starting point using the FM algorithm.
"""
n = W.shape[0]
pstart = transpose(array([x0]))
[D, Dsvg, Ssvg] = perform_dijstra_fm(W, pstart, inf, "fm", "sym", n * 6)
# display
k = 8
displ = lambda D: cos(2 * pi * k * D / max(D.flatten()))
clf
imageplot(displ(D))
set_cmap("jet")
return D
def exo4(tau, x0, x1, G):
"""
Perform the full geodesic path extraction by iterating the gradient
descent. You must be very careful when the path become close to
$x_0$, because the distance function is not differentiable at this
point. You must stop the iteration when the path is close to $x_0$.
"""
n = G.shape[0]
Geval = lambda G, x: (
bilinear_interpolate(G[:, :, 0], imag(x), real(x))
+ 1j * bilinear_interpolate(G[:, :, 1], imag(x), real(x))
)
niter = 1.5 * n / tau
# init gamma
gamma = [x1]
xtgt = x0[0] + 1j * x0[1]
for i in arange(0, niter):
g = Geval(G, gamma[-1])
gamma.append(gamma[-1] - tau * g)
if abs(gamma[-1] - xtgt) < 1:
break
gamma.append(xtgt)
return gamma
import numpy as np
import matplotlib.pyplot as plt
from numpy import (
abs,
add,
arange,
array,
cos,
empty,
exp,
extract,
full,
gradient,
imag,
inf,
linspace,
meshgrid,
mod,
ones,
pi,
real,
sqrt,
stack,
sum,
transpose,
zeros,
)
from matplotlib.pyplot import (
axis,
clf,
close,
gray,
jet,
matplotlib,
np,
plot,
set_cmap,
step,
subplot,
time,
)
%matplotlib inline
Installation¶
Navigating on the Grid¶
We use a cartesian grid of size $n \times n$, and defines operators to navigate in the grid.
We use a singe index $i \in \{1,\ldots,n^2\}$ to index a position on the 2-D grid.
Size of the grid.
n = 40
The four displacement vector to go to the four neightbors.
neigh = array([[1, -1, 0, 0], [0, 0, 1, -1]])
For simplicity of implementation, we use periodic boundary conditions.
boundary = lambda x: mod(x, n)
For a given grid index |k|, and a given neighboring index k in ${1,2,3,4}$, |Neigh(k,i)| gives the corresponding grid neighboring index.
ind2sub1 = lambda k: [np.asarray(k, dtype=int) // n, np.asarray(k, dtype=int) % n]
sub2ind1 = lambda u: int(u[0] * n + u[1])
Neigh = lambda k, i: sub2ind1(boundary(ind2sub1(k) + neigh[:, i]))
Check that these functions are indeed bijections.
print(ind2sub1(sub2ind1([13, 27])))
print(sub2ind1(ind2sub1(134)))
Dikstra Algorithm¶
The Dijkstra algorithm compute the geodesic distance on a graph. We use here a graph whose nodes are the pixels, and whose edges defines the usual 4-connectity relationship.
In the following, we use the notation $i \sim j$ to indicate that an index $j$ is a neighbor of $i$ on the graph defined by the discrete grid.
The metric $W(x)$. We use here a constant metric.
W = ones((n, n))
Set $\Ss = \{x_0\}$ of initial points.
x0 = [n / 2, n / 2]
Initialize the stack of available indexes.
I = [sub2ind1(x0)]
Initialize the distance to $+\infty$, excepted for the boundary conditions.
D = ones((n, n)) + inf
u = ind2sub1(I)
D[u[0], u[1]] = 0
Initialize the state to 0 (unexplored), excepted for the boundary point $\Ss$ (indexed by |I|) to $1$ (front).
S = zeros((n, n))
S[u[0], u[1]] = 1
Define a callbabk to use a 1-D indexing on a 2-D array.
extract = lambda x, I: x[I]
extract1d = lambda x, I: extract(x.flatten(), I)
The first step of each iteration of the method is to pop the from stack the element $i$ with smallest current distance $D_i$.
j = np.argsort(extract1d(D, I))
if np.ndim(j) == 0:
j = [j] # make sure that j is a list a not a singleton
j = j[0]
i = I[j]
a = I.pop(j)
We update its state $S$ to be dead (-1).
u = ind2sub1(i)
S[u[0], u[1]] = -1
Retrieve the list of the neighbors that are not dead and add to I those that are not yet in it.
J = []
for k in np.arange(0, 4):
j = Neigh(i, k)
if extract1d(S, j) != -1:
# add to the list of point to update
J.append(j)
if extract1d(S, j) == 0:
# add to the front
u = ind2sub1(j)
S[u[0], u[1]] = 1
I.append(j)
imageplot(S)
Update neighbor values. For each neightbo $j$ of $i$, perform the update, assuming the length of the edge between $j$ and $k$ is $W_j$. $$ D_j \leftarrow \umin{k \sim j} D_k + W_j. $$
DNeigh = lambda D, k: extract1d(D, Neigh(j, k))
for j in J:
dx = min(DNeigh(D, 0), DNeigh(D, 1))
dy = min(DNeigh(D, 2), DNeigh(D, 3))
u = ind2sub1(j)
w = extract1d(W, j)
D[u[0], u[1]] = min(dx + w, dy + w)
Worked example 1
Implement the Dijkstra algorithm by iterating these step while the stack |I| is non empty. Display from time to time the front that propagates.
D = exo1(x0, W)
Display the geodesic distance map using a cosine modulation to make the level set appears more clearly.
displ = lambda D: cos(2 * pi * 5 * D / max(D.flatten()))
imageplot(displ(D))
set_cmap("jet")
Fast Marching¶
The Dijstra algorithm suffers from a strong metrization problem, and it actually computes the $\ell^1$ distance on the grid.
The Fast Marching algorithm replace the graph update by a local resolution of the Eikonal equation. This reduces significantly the grid bias, and can be shown to converge to the underlying geodesic distance when the grid step size tends to zero.
Over a continuous domain, the distance map $D(x)$ to a set of seed points $ \Ss $ is the unique solution in the viscosity sense $$ \forall x \notin \Ss, \quad \norm{\nabla D(x)} = W(x) \qandq \forall y \in \Ss, \quad D(y) = 0. $$
The equation is then discretized on a grid of $n \times n$ pixel, and a solution $ (D_{k,\ell})_{k,\ell=1}^n \in \RR^{n \times n} $ is found by using an upwind finite difference approximation, that is faithful to the viscosity solution $$ \forall (k,\ell) \notin \tilde \Ss, \quad \norm{ (\nabla D)_{k,\ell} } = W_{k,\ell}$ \qandq \forall (k,\ell) \notin \tilde \Ss, \quad D_{k,\ell}=0, $$ where $\tilde \Ss$ is the set of discrete starting points (defined here by |x0|).
To be consisten with the viscosity solution, one needs to use a non-linear upwind gradient derivative. This corresponds to computing the norm of the gradient as $$ \norm{ (\nabla D)_{k,\ell} }^2 = \max( D_{k+1,\ell}-D_{k,\ell}, D_{k-1,\ell}-D_{k,\ell}, 0 )^2 + \max( D_{k,\ell+1}-D_{k,\ell}, D_{k,\ell-1}-D_{k,\ell}, 0 )^2. $$
A each step of the FM propagation, one update $ D_{k,\ell} \leftarrow d $ by solving the eikonal equation with respect to $D_{k,\ell}$ alone. This is equivalent to solving the quadratic equation $$ (d-d_x)^2 + (d-d_y)^2 = w^2 \qwhereq w=W_{k,\ell}. $$ and where $$ d_x = \min(D_{k+1,\ell},D_{k-1,\ell}) \qandq d_y = \min(D_{k,\ell+1},D_{k,\ell-1}). $$
The update is thus defined as $$ d = \choice{ \frac{d_x+d_y+ \sqrt{\De}}{2} \quad\text{when}\quad \De \geq 0, \\ \min(d_x,d_y)+w \quad \text{otherwise.}$ }$ \qwhereq \De = 2 w^2 - (d_x-d_y)^2. $$
Note that in the case where $\De<0$, one has to use the Dijkstra update.
Once the Dijstra algorithm is implemented, the implementation of the Fast Marching is trivial. It just corresponds to replacing the graph udpate
D[u[0], u[1]] = min(dx + w, dy + w)
by the eikonal update.
Delta = 2 * w**2 - (dx - dy) ** 2
if Delta >= 0:
D[u[0], u[1]] = (dx + dy + sqrt(Delta)) / 2
else:
D[u[0], u[1]] = min(dx + w, dy + w)
Worked example 2
Implement the Fast Marching algorithm. Display from time to time the front that propagates.
D = exo2(x0, W)
Display the geodesic distance map using a cosine modulation to make the level set appears more clearly.
imageplot(displ(D))
set_cmap("jet")
Computation of Geodesic Paths¶
We use a more complicated, non-constant metric, with a bump in the middle.
n = 100
x = linspace(-1, 1, n)
[Y, X] = meshgrid(x, x)
sigma = 0.2
W = 1 + 8 * exp(-(X**2 + Y**2) / (2 * sigma**2))
Display it.
imageplot(W)
Starting points.
x0 = [round(0.1 * n), round(0.1 * n)]
Worked example 3:
Compute the distance map to these starting point using the FM algorithm. Important: use symetric boundary conditions.
D = exo3(x0, W)
Once the geodesic distance map to $\Ss$ has been computed, the geodesic curve between any point $x_1$ and $\Ss$ extracted through gradient descent $$ \ga'(t) = - \eta_t \nabla D(\ga(t)) \qandq \ga(0)=x_1 $$ where $\eta_t>0$ controls the parameterization speed of the resulting curve.
To obtain unit speed parameterization, one can use $\eta_t = \norm{\nabla D(\ga(t))}^{-1}$ (one need to be careful when $\ga$ approaches $\Ss$ since $D$ is not smooth on $\Ss$).
Compute the gradient $G_0(x) = \nabla D(x) \in \RR^2$ of the distance map. Use centered differences.
G0 = grad(D)
Normalize the gradient to obtained $G(x) = G_0(x)/\norm{G_0(x)}$, in order to have unit speed geodesic curve (parameterized by arc length).
d = sqrt(sum(G0**2, axis=2))
U = zeros((n, n, 2))
U[:, :, 0] = d
U[:, :, 1] = d
G = G0 / U
The geodesic is then numerically computed using a discretized gradient descent, which defines a discret curve $ (\ga_k)_k $ using $$ \ga_{k+1} = \ga_k - \tau G(\ga_k) $$ where $\ga_k \in \RR^2$ is an approximation of $\ga(t)$ at time $t=k\tau$, and the step size $\tau>0$ should be small enough.
Step size $\tau$ for the gradient descent.
tau = 0.8
Initialize the path with the ending point.
x1 = round(0.9 * n) + 1j * round(0.88 * n)
gamma = [x1]
Define a shortcut to interpolate $G$ at a 2-D points. Warning: the |interp2| switches the role of the axis ...
Geval = lambda G, x: (
bilinear_interpolate(G[:, :, 0], imag(x), real(x))
+ 1j * bilinear_interpolate(G[:, :, 1], imag(x), real(x))
)
Compute the gradient at the last point in the path, using interpolation.
g = Geval(G, gamma[-1])
Perform the descent and add the new point to the path.
gamma.append(gamma[-1] - tau * g)
Worked example 4
Perform the full geodesic path extraction by iterating the gradient descent. You must be very careful when the path become close to $x_0$, because the distance function is not differentiable at this point. You must stop the iteration when the path is close to $x_0$.
gamma = exo4(tau, x0, x1, G)
Display the geodesic curve.
clf
imageplot(W)
set_cmap("gray")
h = plot(imag(gamma), real(gamma), ".b", linewidth=2)
h = plot(x0[1], x0[0], ".r", markersize=20)
h = plot(imag(x1), real(x1), ".g", markersize=20)
References and further reading¶
James A. Sethian. A Fast Marching Level Set Method for Monotonically Advancing Fronts. 1996, PNAS 93(4), 1591–1595. A causal numerical solver for the eikonal equation.
Edsger W. Dijkstra. A Note on Two Problems in Connexion with Graphs. 1959, Numerische Mathematik 1, 269–271. Shortest paths with nonnegative edge weights.
Stanley Osher and James A. Sethian. Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton–Jacobi Formulations. 1988, Journal of Computational Physics 79(1), 12–49. Implicit front propagation using level-set functions.
Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. A Global Geometric Framework for Nonlinear Dimensionality Reduction. 2000, Science 290, 2319–2323. Geodesic distances and multidimensional scaling for manifold learning.
Gabriel Peyré. Advanced Signal, Image and Surface Processing. 2010, course notes. A mathematical companion to the Numerical Tours.