In this notebook, we demonstrate anisotropic fast marching with Rander metrics, in two dimensions.
Rander metrics are a generalization of Riemannian metrics, featuring an additional linear term. They are also a special case of non-symmetric Finslerian metrics. A Rander metric measures vectors according to the formula: $$ F_x(v) := \|v\|_{M(x)} + <\omega(x), v> $$ where $M$ is a field of symmetric positive definite tensors, and $\omega$ is a vector field. Applications of Rander metrics include:
The HFM software computes the distance associated to a given Rander metric, and the corresponding minimal paths, by solving a variant of the eikonal PDE. Namely for all $x$ within a domain $\Omega$ $$ \|\nabla u(x)\|_{D(x)} + <\eta(x),\nabla u(x)> = 1, % \|\nabla u(x) - \omega(x)\|_{M(x)^{-1}} = 1. $$ where $(D,\eta)$ is the dual metric. Some algebraic formulas allow to express the dual metric in terms of $(M,\omega)$, the primal metric, see the first two references below.
References The experiments presented in this notebook, or close variants, are presented in the following publications.
Mirebeau, J.-M. (2014). Efficient fast marching with Finsler metrics. Numerische Mathematik, 126(3), 515–557. link
Mirebeau, J.-M. (2017, April 12). Riemannian fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms. HAL (Preprint). link
Mirebeau, J.-M., Cohen, L. D., Chen, D., & Mirebeau, J.-M. (2016). Finsler Geodesics Evolution Model for Region based Active Contours. In E. R. H. Richard C Wilson & W. A. P. Smith (Eds.), (pp. 22.1–22.12). Presented at the Proceedings of the British Machine Vision Conference (BMVC), BMVA Press. http://doi.org/10.5244/C.30.22
Summary of volume Fast Marching Methods, this series of notebooks.
Main summary of the Adaptive Grid Discretizations book of notebooks, including the other volumes.
This Python® notebook is intended as documentation and testing for the HamiltonFastMarching (HFM) library, which also has interfaces to the Matlab® and Mathematica® languages. More information on the HFM library in the manuscript:
Copyright Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
import sys; sys.path.insert(0,"..") # Allow import of agd from parent directory (useless if conda package installed)
#from Miscellaneous import TocTools; print(TocTools.displayTOC('Rander','FMM'))
from agd import Eikonal
from agd.Plotting import savefig, quiver; #savefig.dirName = 'Figures/Rander'
from agd.Metrics import Rander
import numpy as np; xp=np
%matplotlib inline
import matplotlib.pyplot as plt
Uncomment the following line use the GPU eikonal solver. (Comment it for the CPU eikonal solver.)
#import agd.AutomaticDifferentiation as ad; xp,quiver,plt,Eikonal = map(ad.cupy_friendly,(xp,quiver,plt,Eikonal))
We compute the travel time and shortest path for a vehicle with unit velocity in all directions, but suject to a drift. In order words, the vehicle trajectory $\gamma$ obeys at all times $$ \| \gamma'(t) - V(\gamma(t)) \|_{D(\gamma(t))^{-1}} = 1, $$ where $\eta$ is the given drift, e.g. ocean current. The positive definite tensor field $D$ is chosen here equal to the identity, but could be modified to account for similar problems posed on manifolds.
Zermelo's navigation is locally controllable iff for all $x$ within the domain $\Omega$ one has $$ \|V(x)\|_{D(x)^{-1}} < 1. $$ This condition is a pre-requisite for our eikonal solver to function.
hfmIn = Eikonal.dictIn({
'model':'Rander2',
'exportValues':1,
'seed':[0.,0.],
})
# Define the domain
n=201
hfmIn.SetRect(sides=[[-0.5,0.5],[-0.5,0.5]],dimx=n)
hfmIn.SetUniformTips((6,6))
X = hfmIn.Grid()
R = np.linalg.norm(X,axis=0)
driftMult = 0.9*np.sin(4*np.pi*X[0])*np.sin(4.*np.pi*X[1])
drift = (driftMult/R) * X
The drift, illustrated in the next figure, increases or decreases the vehicle motion depending on its position and orientation. In particular, we see that close to the image center, the drift field
fig = plt.figure(figsize=[4,4]); plt.axis('equal');
plt.title('Drift vector field');
quiver(*X,*drift,subsampling=(5,5));
savefig(fig,"DriftVectorField.png")
The Rander metric associated to Zermelo's problem is not completely trivial. The dual metric of Zermelo's navigation problem is $(D,-V)$, whereas the primal metric has a more complex expression, for which we refer to the papers cited in the introduction. In any case, the relevant formulas are implemented in the Rander class.
hfmIn['metric'] = Rander.from_Zermelo(xp.eye(2),drift) # Riemannian metric, and drift
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.025412 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
fig = plt.figure(figsize=[4,4]); plt.axis('equal');
plt.title('Distance map, Zermelo\'s navigation problem');
plt.contourf(*X,hfmOut['values']);
savefig(fig,"ZermeloDistance.png")
fig = plt.figure(figsize=[4,4]); plt.axis('equal');
plt.title('Minimal geodesics, Zermelo\'s navigation problem');
for geo in hfmOut['geodesics']: plt.plot(*geo)
savefig(fig,"ZermeloPaths.png")
Variational image segmentation methods define image sub-domains by minimizing an energy, comprised of region terms and a boundary terms. The most well known model, the Mumford-Shah model, is also among the most complex and will not be addressed here (unkown number of regions, whose energy contribution is determined by a PDE, cracks, etc ...).
In the case of binary image segmentation, one is minimizing $$ E(U) = F(U) + G(\partial U) $$ among all sub-domains $U \subset \Omega$, the second region being $\Omega \setminus U$.
In the Chan-Vese model, the region term is linear and defined by a function $f$, while the boundary term is a measure of length, possibly defined by a Riemannian metric $M$. $$ \begin{aligned} F(U) &= \int_U f, & G(\partial U) &= \mathrm{len}_M(\partial U) \end{aligned} $$
Assume that $\partial U$ is parametrized counter-clockwise by a path $\gamma : [0,1] \to \Omega$, and that $f = \mathrm{div} V$ in the neighborhood of $\gamma$ (more on the computation of $V$ below). Then the energy rewrites, using Stokes theorem for the region term, $$ E(U) = \int_0^1 \| \gamma'(t) \|_{M(\gamma(t))} + < \gamma'(t), V(t)^\perp> \mathrm{d}t. $$ In other words, $U$ is optimal if and only iff its boundary is a minimal path w.r.t. the Rander metric of parameters $(M,V^\perp)$.
Application examples. Due to lack of time, only a synthetic problem is presented in this notebook. See the reference in the introduction for application examples.
An iterative or non-iterative method ? In ideal conditions, the above characterization allows to extract the optimal sub-domain $U$ in a single run of the fast marching algorithm, with the suitable Rander metric. In practice, an iterative approach is nevertheless required for the following reasons:
Note on numerical experiments We content ourselves with a purely synthetic instance, and refer to the reference cited in the introduction for application inspired examples.
We check the soundness of the approach on a completely analytic instance. Contour length is measured using the Euclidean metric $M = \mathrm{Id}$. On the other hand, the region term is defined by the integral of $$ \begin{aligned} f(x) = 4(\|x\|^2-1) e^{-\|x\|^2}. \end{aligned} $$ This function is chosen so that $f=\mathrm{div} V$ with $$ \begin{aligned} V(x) = - 2 x e^{-\|x\|^2}. \end{aligned} $$
In view of the problem symmetries, the energy is maximized on some disk centered at the origin. The energy of a disk of radius $r$ is $$ R(D_r) = 2 \pi \int_0^r f(r) r \mathrm{d}r + 2 \pi r. $$ This energy is extremal when $1+r f(r)=0$, which occurs for $r=0.64..$ (minimum), and $r=0.30..$ (local maximum).
# Roots of 1+r f(r)
r=xp.linspace(0,1)
def f(r):
return 4*(r**2-1)*np.exp(-r**2)
plt.plot(r,1+r*f(r),
r,0.*r);
hfmIn = Eikonal.dictIn({
'model':'Rander2',
'exportValues':1,
# Ideal case : the seed and tip are
# already on the optimal region boundary
'seed':[0.64,0.],
'tip': [-0.64,0.],
})
# Define the domain
n=201
hfmIn.SetRect(sides=[[-1,1],[-1,1]],dimx=n)
X = hfmIn.Grid()
R = np.linalg.norm(X,axis=0)
v = -2*X*np.exp(-R**2)
hfmIn['metric'] = Rander(xp.eye(2),[v[1],-v[0]])
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.032072 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
plt.figure(figsize=[4,4]); plt.title('Distance map, Isoperimetric problem'); plt.axis('equal');
plt.contourf(*X,hfmOut['values']);
# As expected, the extracted boundary is an arc of circle.
fig = plt.figure(figsize=[4,4]); plt.title('Minimal geodesic, Isoperimetric problem'); plt.axis('equal');
plt.plot(*hfmOut['geodesic']);
Extracting the second half of the boundary
By exchanging the roles of the seeds and tips, we extract the second half of the boundary.
hfmIn['seed'],hfmIn['tip'] = hfmIn['tip'],hfmIn['seed']
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field refineStencilAtWallBoundary defaults to 0 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.034914 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
plt.figure(figsize=[4,4]); plt.title('Distance map, Isoperimetric problem'); plt.axis('equal');
plt.contourf(*X,hfmOut['values']);
# As expected, the extracted boundary is an arc of circle.
fig = plt.figure(figsize=[4,4]); plt.title('Minimal geodesics, Isoperimetric problem'); plt.axis('equal');
plt.plot(*hfmOut['geodesic']);