In this notebook, we demonstrate anisotropic fast marching with asymmetric quadratic metrics, in two and three dimensions. The metrics considered in the notebook are a generalization of Riemannian metrics, featuring an additional non-symmetric term. They are also a special case of non-symmetric Finslerian metrics. An asymmetric quadratic metric measures vectors according to the formula: $$ F_x(v) := \sqrt{ \|v\|^2_{M(x)} + \max\{0,<\omega(x), v>\}^2} $$ where $M$ is a field of symmetric positive definite tensors, and $\omega$ is a vector field. As evidenced by the above formula, the role of the vector field $\omega$ is to further penalize motion in its direction.
The HFM software computes the distance associated to a given asymmetric quadratic metric, and the corresponding minimal paths, by solving a variant of the eikonal PDE. Namely for all $x$ within a domain $\Omega$ $$ \sqrt{\|\nabla u(x)\|^2_{D(x)} + \max\{0,<\eta(x),\nabla u(x)>\}^2} = 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.
Technical note The two and three dimensional implementations use different discretization schemes: semi-Lagrangian, vs Eulerian. The latter one in addition involves a relaxation parameter. As a result two dimensional implementation can reach high levels of accuracy (including second/third order), whereas the three dimensional one is best used in contexts where speed and qualitative behavior are most important.
References The experiments presented in this notebook, or related variants, are presented in the following publications.
Mirebeau, J.-M. (2014). Efficient fast marching with Finsler metrics. Numerische Mathematik, 126(3), 515–557. link
Duits, R., Meesters, S. P., Mirebeau, J.-M., & Portegies, J. M. (2018). Optimal paths for variants of the 2D and 3D Reeds-Shepp car with applications in image analysis. Journal of Mathematical Imaging and Vision, 1–33. http://doi.org/ https://doi.org/10.1007/s10851-018-0795-z
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('AsymmetricQuadratic','FMM'))
from agd import Eikonal
from agd.Metrics import AsymQuad,Riemann
from agd.Plotting import savefig; #savefig.dirName = 'Figures/AsymmetricQuadratic'
from agd import FiniteDifferences as fd
from agd import AutomaticDifferentiation as ad
import numpy as np; xp = np
import matplotlib.pyplot as plt
Uncomment the following line to use the GPU accelerated eikonal solver.
#xp,Eikonal,plt = map(ad.cupy_friendly,(xp,Eikonal,plt))
In order to illustrate the concept of asymmetric quadratic metric, we compute the distance map with respect to a metric independent of the position $x$ within the domain. It is defined by $$ F(v)^2 = \|v\|_D^2 + \max\{0,<\omega,v>\}^2, $$ where $D$ is the identity tensor, and $\omega = (1,1)$.
hfmIn = Eikonal.dictIn({
'model':'AsymmetricQuadratic2',
'exportValues':1,
'seed':[0,0],
})
# Define the domain
n=201
hfmIn.SetRect(sides=[[-1,1],[-1,1]],dimx=n)
X,Y = hfmIn.Grid()
hfmIn.SetUniformTips((6,6))
hfmIn['metric']= AsymQuad(xp.eye(2),[1,1]) # Inputs are : D, w
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.014662 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
Since the metric is constant, the distance map from the origin is simply $u(x) = F(x)$.
The additional cost of motion in the direction $\omega = (1,1)$ is clearly visible. Each level line of $u$ is built of two half ellipses, defined by the tensors $D$ and $D+\omega \omega^T$.
plt.figure(figsize=[5,4]); plt.title('Distance map, asymmetric quadratic norm'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values']);
plt.colorbar();
Since the metric is constant, minimal geodesics are straight lines toward the origin.
fig = plt.figure(figsize=[4,4]); plt.title('Minimal geodesics, for a constant metric'); plt.axis('equal');
for geo in hfmOut['geodesics']: plt.plot(*geo)
Non-symmetric metrics are well fit to extract structures possessing a preferred orientation. For instance, the contours of two dimensional objects are naturally oriented, by rotating the outward normal clockwise. See the notebook VI - Rander metrics for a discussion on image sub-domain segmentation.
In the following experiment, we illustrate a different use of non-symmetric metrics, to avoid the shortcut's problem in tubular structure segmentation. A similar effect can be achieved with the use of non-holonomic metrics penalizing curvature, such as Euler-Mumford elasticae, see IV - Curvature penalized planar paths. Curvature penalization models can be simpler to use, but are much more computationally expensive.
def gamma(t):
return ad.array((xp.cos(t),0.5*xp.sin(2*t)))
def gamma_tgt(t):
tX,tY = -xp.sin(t),xp.cos(2*t)
tN = xp.sqrt(tX**2+tY**2)
return ad.array((tX/tN,tY/tN))
def gamma_normal(t):
tX,tY = gamma_tgt(t)
return ad.array((tY,-tX))
T=xp.linspace(0,2*np.pi,100)
plt.axis('equal')
plt.plot(*gamma(T));
The construction of our example involves running a first fast marching.
isoIn = Eikonal.dictIn({
'model':'Isotropic2',
'exportValues':1,
'seeds':gamma(T).T, # Second T is transposition
'speed':1,
'seedValueVariation':gamma_tgt(T) # Interpolates data defined at the seeds
})
# Define the domain
n=201
isoIn.SetRect(sides=[[-1.1,1.1],[-0.6,0.6]],dimx=n)
X,Y = isoIn.Grid()
isoOut = isoIn.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.003738 s.
We computed the distance to the curve, and an extension of its tangents.
plt.contourf(X,Y,ad.remove_ad(isoOut['values'])); plt.axis('equal');
var = isoOut['values'].gradient()
plt.quiver(X[::5,::5],Y[::5,::5],var[0,::5,::5],var[1,::5,::5]); plt.axis('equal');
distGamma = ad.remove_ad(isoOut['values'])
baseSpeed = np.exp(-(distGamma/0.1)**2)+0.1
baseCost = 1/baseSpeed
tgtGammaX,tgtGammaY = var
tgtGammaN = np.sqrt(tgtGammaX**2+tgtGammaY**2)
tgtGammaX,tgtGammaY = tgtGammaX/tgtGammaN,tgtGammaY/tgtGammaN
plt.contourf(X,Y,baseCost); plt.axis('equal');
hfmIn = Eikonal.dictIn({
'model':'AsymmetricQuadratic2',
'exportValues':1,
'seed':[0.75,-0.5],
'tip':[-0.75,-0.5]
})
# Define the domain
n=201
hfmIn.SetRect(sides=[[-1.1,1.1],[-0.6,0.6]],dimx=n)
X,Y = hfmIn.Grid()
tubularMetric = AsymQuad(baseCost*fd.as_field(xp.eye(2),baseCost.shape), # D=baseCost*Id
[5*tgtGammaX,5*tgtGammaY]) # w
hfmIn['metric'] = tubularMetric
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.013311 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
The front propagates along the tubular structure. This propagation is asymmetric, since the direction proportional to the tangent is strongly penalized.
plt.contourf(X,Y,hfmOut['values']); plt.axis('equal');
plt.contourf(X,Y,ad.remove_ad(isoOut['values'])); plt.axis('equal')
plt.plot(*hfmOut['geodesic'],color='red')
[<matplotlib.lines.Line2D at 0x12f2ed610>]
Exchanging the role of the seeds and tips yields the other half of the tubular structure, thanks to the metric asymmetry.
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.012945 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.contourf(X,Y,ad.remove_ad(isoOut['values'])); plt.axis('equal')
plt.plot(*hfmOut['geodesic'],color='red');
In contrast, a symmetric metric cannot recover such a structure, due to the shortcuts problem.
hfmIn['metric'] = AsymQuad(hfmIn['metric'].m,[0,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.008679 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.contourf(X,Y,ad.remove_ad(isoOut['values'])); plt.axis('equal')
plt.plot(*hfmOut['geodesic'],color='red');
# Save for future use
hfmIn['metric']=tubularMetric
tubularIn = hfmIn
In this section, we present a highly tunable class of metrics, strongly asymmetric and anisotropic, referred to as asymmetric Rander metrics. They generalize both Rander metrics and asymmetric quadratic metrics, and take the following form $$ F_x(\dot x) := \sqrt{ \|\dot x\|^2_{M(x)} + \max\{0,<u(x), \dot x>\}^2+\max\{0,<v(x), \dot x>\}^2}+<w(x), \dot x>, $$ where $M$ is a field of positive definite matrices, and $u,v,w$ are vector fields, with $w$ small enough.
Clearly, asymmetric Rander metrics reduce to:
Using several non-zero vector fields among $u,v,w$ one can achieve anisotropies that are not realizable with the simpler metric classes, and whose uses cases are yet to be determined.
Limitations. Fast marching with Asymmetric Rander metrics is currently only available in two dimensions, and is not GPU accelerated.
if xp is not np: raise ad.DeliberateNotebookError("Sorry, asymmetric Rander metrics are not supported by the GPU Eikonal solver")
from agd.Metrics.asym_rander import AsymRander
hfmIn = Eikonal.dictIn({
'model':'AsymRander2',
'seed':(0,0),
'exportValues':1,
'factoringRadius':-1,
'exportGeodesicFlow':1,
})
hfmIn.SetRect([[-1,1],[-1,1]],dimx=101)
X = hfmIn.Grid()
Since we use a large factoring radius, constant metrics are reproduced up to numerical precision.
def check(constant_metric):
# Run
hfmIn['metric'] = constant_metric
hfmOut = hfmIn.Run()
# Display
plt.axis('equal')
plt.contourf(*X,hfmOut['values'])
plt.scatter(0,0,color='red') # Seed position
# Check :
assert np.allclose(hfmOut['values'],constant_metric.norm(X))
assert np.all(np.isfinite(hfmOut['flow']))
assert np.sum((constant_metric.norm(hfmOut['flow']) - 1)>1e-8)<12
Choosing $u$ and $v$ as large vectors, while $w=0$, yields a metric that almost forbids motion except in the cone defined as $$ <u,\dot x> \leq 0 \text{ and } <v,\dot x> \leq 0 $$
Constructor of AsymRander
. Input arguments are $m,u,v,w$.
check( AsymRander( np.eye(2), [4.,0], [0.,4.], None ) )
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field order defaults to 1 Field seedRadius defaults to 2 Field factoringPointChoice defaults to Seed Field exportFactoring defaults to 0 Fast marching solver completed in 0.008158 s.
Here we choose $u=0$, while $v$ and $w$ are non-zero.
check( AsymRander( np.eye(2), None, [0.,4.], [0.5,0.] ) )
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field order defaults to 1 Field seedRadius defaults to 2 Field factoringPointChoice defaults to Seed Field exportFactoring defaults to 0 Fast marching solver completed in 0.009008 s.
Obviously, all three of $u,v,w$ may be zero, producing even stranger anisotropy shapes.
Alternatively, as discussed above, we may reproduce an asymmetric quadratic metric by choosing $v=w=0$.
asymQuadOut = tubularIn.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.01285 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
tubularIn2 = tubularIn.copy()
tubularIn2.update({
'model':'AsymRander2',
'metric':AsymRander(tubularMetric.m,tubularMetric.w,None,None),
})
asymRanderOut = tubularIn2.Run()
Field verbosity defaults to 1 Field cosAngleMin defaults to 0.5 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.015457 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45
Most values (95%) are reproduced up to machine precision.
assert np.sum(np.abs(asymRanderOut['values']-asymQuadOut['values'])>1e-8)/tubularIn.size < 0.05
The difference between the non-reproduced values is rather small, yet non-zero.
valueDiff = np.abs(asymRanderOut['values']-asymQuadOut['values'])
assert np.max(valueDiff) < 0.06
assert np.mean(valueDiff) < 3e-4
The reason some values are not exactly reproduced is because the numerical scheme involves adaptive and anisotropic stencils, which depend discontinuously on the metric. Due to numerical discretization errors, the stencil may differ at one place, and this creates a discrepancy between the numerical solutions, which is propagated in a subdomain. This behavior may be avoided by using fixed stencils, but in that case the scheme becomes more numerically expensive and less accurate.
plt.title("Places where the cost functions differ")
plt.contourf(valueDiff>1e-8);