This notebook shows how the HFM library can solve isotropic eikonal equations and extract the related minimal paths. The assumption of isotropy makes this task rather simple and classical, and more complex models are addressed in the subsequent notebooks. We present some two and three dimensional test cases, featuring obstacles, position and eventually time-dependent velocities, various stopping criteria, ...
The functionality presented in this notebook, in the context of isotropic fast marching, is also applicable to the more complex anisotropic models implemented in the HFM software.
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('Isotropic','FMM'))
from agd import Eikonal
from agd.Plotting import savefig, quiver; # savefig.dirName = 'Figures/Isotropic'
# We will also need some standard python libraries.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
# Required for 3D plotting
from mpl_toolkits.mplot3d import Axes3D # Plots geodesics
useMayavi = False
if useMayavi:
from mayavi import mlab # Plots implicit surfaces
Uncomment the following line to use the GPU eikonal solver. (Comment it for the CPU eikonal solver.)
#from agd import AutomaticDifferentiation as ad; plt,quiver,Eikonal = map(ad.cupy_friendly,(plt,quiver,Eikonal))
In this document, we solve the standard isotropic eikonal equation, on domain denoted $\Omega$. A cost function $c : \Omega \to ]0,\infty[$ and some boundary data $\sigma : \partial \Omega \to ]-\infty,\infty]$ are given. This partial differential equation (PDE) reads \begin{align} \forall x \in \Omega, \|\nabla u(x)\| &= c(x), & \forall x \in \partial \Omega, u(x) &= \sigma(x). \end{align}
This PDE admits a unique viscosity solution to this PDE, defined as follows. The solution value $u(x)$ is minimal length (plus initial value penalization) of a path $\gamma : [0,1] \to \overline \Omega$ from the domain boundary $\partial \Omega$ to $x$. \begin{equation*} u(x) = \inf_{\gamma(1)=x} \sigma(\gamma(0)) + \int_0^1 c(\gamma(t)) \|\gamma'(t)\| \, \mathrm d t. \end{equation*} Equivalently, the solution $u$ can be regarded as the arrival time at $x \in \Omega$ of a front originating from $\partial \Omega$ at time $\sigma : \partial \Omega \to ]-\infty,\infty]$, and propagating at speed $s = 1/c$ where $c : \Omega \to ]0,\infty[$.
hfmIn = Eikonal.dictIn({ # Create a dictionary-like structure to store the input data
'model':'Isotropic2', # Isotropic two-dimensional eikonal equation
'order':2., # Use a second order scheme, so as to increase accuracy
'cost':1. # Unit cost for now. A position dependent cost function is also considered below.
#'speed':1. # One may equivalently provide a speed function. In that case cost = 1/speed.
})
The above 'model', defined as 'Isotropic2', also determines the nature of the enclosing domain, that we denote $\Omega_0$. It is a two dimensional box, with outflow boundary conditions on $\partial \Omega_0$. (As opposed to e.g. periodic.)
We will use the rectangle $[-1,1] \times [0,1]$, discretized on a $2n \times n$ grid where $n=100$.
hfmIn.SetRect(sides=[[-1,1],[0,1]],gridScale=0.01)
For the eikonal equation to be well posed, we introduce at set $S$ of 'seeds' in the domain, to which we attach initial values $(\sigma_s)_{s\in S}$. The PDE domain $\Omega\subset \Omega_0$ its boundary $\partial \Omega$, and the boundary conditions $\sigma : \partial \Omega \to ]-\infty,\infty]$ are thus as follows.
\begin{align} \Omega &= \Omega_0 \setminus S, & \partial \Omega &= \partial \Omega_0 \cup S, & \sigma(x) = \begin{cases} \sigma_x & \text{ if } x \in S\\ +\infty & \text{ if } x \in \partial \Omega_0. \end{cases} \end{align}hfmIn.update({
# Introduce two seeds, at positions (-0.5,0.3) and (0.5,0.8)
'seeds':[[-0.5,0.3],[0.5,0.8]],
# Boundary conditions imposed at the seeds.
'seedValues':[0.,0.5],
# 'seedValues' defaults [0.,0.] if unspecified.
})
The HFM software rounds the seed positions to the closest point of the discretization grid. (If high accuracy is needed, please use source factorization methods, which not require the point to be on the grid. See notebook on high accuracy.)
The followinf method allows to get the multi-index of one or several points in the discretization grid.
index,rounding_error = hfmIn.IndexFromPoint(hfmIn['seeds'])
print(f"Seed positions rounded to {index}, with rounding errors {rounding_error}")
Seed positions rounded to [[ 50 29] [150 80]], with rounding errors [[-0.5 0.5] [-0.5 -0.5]]
Conversely, the position associated to a given index, with integer or floating point entries, can be recovered.
print(f"Rounded positions {hfmIn.PointFromIndex(index)}, original positions {hfmIn.PointFromIndex(index+rounding_error)}")
Rounded positions [[-0.495 0.295] [ 0.505 0.805]], original positions [[-0.5 0.3] [ 0.5 0.8]]
One last step, before running the Fast-Marching solver, is to indicate which outputs are desired. Here we request:
hfmIn['exportValues']=1. # Ask for the PDE solution
hfmIn['exportGeodesicFlow']=1 # Ask for the geodesic flow
We also request three backtracked geodesics, from tips that we specify. In contrast with the seeds, defined above, the geodesic tips are not rounded to a point of the discretization grid.
The HFM software implements two geodesic backtracking techniques, which usually produce similar results, and are referred to as 'ODE' (integration of the geodesicFlow using the Euler midpoint scheme), and 'Discrete' (based on automatic differentiation). For details see the paper.
hfmIn['tips'] = [[0.,0.6],[-0.9,0.5],[0.8,0.8]] # Ask for the geodesics from these three points
#hfmInput['geodesicSolver']='ODE' # Choose the backtracking method. Defaults to 'Discrete' if unspecified. (CPU only)
For demonstration purposes, we also add a dummy, useless input key
hfmIn['dummyKey']='dummyValue'
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.003668 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45 ***** Warning ! ***** Unused fields from user: dummyKey ********************
Unused keys are frequently mispelled fields. That is why they are shown with a warning. Defaulted inputs might also be of interest, as well as the list of outputs.
print("Output keys : ", hfmOut.keys(),"\n")
Output keys : dict_keys(['FMCPUTime', 'GeodesicCPUTime', 'MaxStencilWidth', 'StencilCPUTime', 'defaulted', 'flow', 'geodesics', 'nAccepted', 'unusedFromCompute', 'unusedFromUser', 'values', 'visitedUnset', 'retcode'])
if hfmIn.mode=='cpu':
print("Unused input keys : ", hfmOut['unusedFromUser'],"\n")
print("Defaulted input keys : ", hfmOut['defaulted'],"\n")
elif hfmIn.mode=='gpu':
print("Available detail on keys",hfmOut['keys'].keys())
print("Unused input keys : ", hfmOut['keys']['unused'],"\n")
print("Defaulted input keys : ", hfmOut['keys']['default'],"\n")
Unused input keys : dummyKey Defaulted input keys : exportActiveNeighs exportActiveOffsets geodesicSolver geodesicStep geodesicVolumeBound geodesicWeightThreshold refineStencilAtWallBoundary seedRadius showProgress verbosity
hfmIn.pop('dummyKey',None)
'dummyValue'
The level lines of the computed distance map can be displayed using the contourf function.
X,Y = hfmIn.Grid() # Create a coordinate system
plt.figure(figsize=[6,3]); plt.title('Distance function, for a uniform cost'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values']); # Display the value function.
The geodesics, three in our case, are returned as a list of points, and a list of lengths. The following line splits them into separate arrays, for more convenient manipulation. We then plot.
fig = plt.figure(figsize=[6,3]); plt.title('Distance function and minimal paths'); plt.axis('equal');
# Draw the geodesics
for geo in hfmOut['geodesics']: plt.plot(*geo)
# Overlay the contour plot of the distance map
plt.contourf(X,Y,hfmOut['values'],cmap='Greys');
savefig(fig,'Dist_Path_UnifCost.png')
The minimal geodesics associated with our path planning problem obey an Ordinary Differential Equation (ODE), defined in terms of the solution $u$ of the eikonal PDE. More precisely, the minimal geodesic reaching the point $x \in \Omega$ at the time $T = u(x)$ obeys \begin{align*} \forall t \leq T, \gamma'(t) &= V(\gamma(t)), & \gamma(T) &= x. \end{align*} In practice, this ODE is solved backwards in time, as evidenced by terminal boundary condition $\gamma(T) = x$..
The vector field $V : \Omega \to {\mathbb R}^2$, displayed in the next cell, is referred to as the geodesic flow direction, and is expressed in terms of the gradient of the arrival times function $u : \Omega \to \mathbb R$. In the case of an isotropic metric, defined by a cost function $c: \Omega \to ]0,\infty[$ \begin{equation*} V(x) = c(x)^{-2} \nabla u(x), \end{equation*}
fig = plt.figure(figsize=[6,3]); plt.title('Geodesic flow vector field'); plt.axis('equal');
flowX, flowY = hfmOut['flow']
# Quiver plot of the geodesic flow. We subsample by a factor 5, using subscript [::5,::5], for better visualisation.
quiver(X,Y,flowX,flowY,subsampling=(5,5))
savefig(fig,'GeodesicFlow.png')
By construction, the geodesic flow has unit speed, as measured by the metric of the problem. In the present case of an isotropic metric, this is easily checked formally, since we get for each point $x \in \Omega$ \begin{equation*} c(x) \|V(x)\| = c(x)^{-1} \| \nabla u(x)\| = 1. \end{equation*} Note that this identity is not satisfied at the two seed points, which strictly speaking are part of $\partial \Omega$. We check this identity numerically in the next cell.
fig = plt.figure(figsize=[6,3]); plt.title('Norm of the geodesic flow vector field'); plt.axis('equal');
norms = hfmIn['cost']*np.sqrt(flowX**2+flowY**2)
plt.contourf(X[::5,::5],Y[::5,::5],norms[::5,::5])
plt.colorbar();
savefig(fig,'GeodesicFlowNorm.png')
hfmIn['exportGeodesicFlow']=0. # Do not export the geodesic flow in further experiments.
A obstacle is described as an array of 0 (obstacle absent) or 1 (obstacle present). One-pixel-thin structures can be used to represent walls. In the case, considered in the following notebooks, where adaptive wide stencil discretization schemes are used, a test is done to ensure that the scheme offsets do not jump over the walls.
Outflow boundary conditions, $\sigma(x) = \infty$, are applied at each point $x$ of the obstacles.
# Let us construct a round shaped obstacle, and 1-pixel wide a barrier
disk = (X-0.3)**2 + (Y-0.3)**2 <= 0.2**2
barrier = np.logical_and(X==X[70,0], Y>=0.4)
obstacle = np.logical_or(disk,barrier)
hfmIn['walls']=obstacle # Add the obstacle to the input parameters.
fig = plt.figure(figsize=[6,3]); plt.title('Obstacles introduced in the domain'); plt.axis('equal');
plt.contourf(X,Y,obstacle);
savefig(fig,'Obstacles.png')
We next consider a position dependent cost function. Recall the eikonal equation \begin{align*} \forall x \in \Omega, \, \|\nabla u(x)\| &= c(x) & \forall x \in \partial \Omega, \, u(x) = \sigma(x). \end{align*}
hfmIn['cost'] = np.exp(-0.5*(X**2+Y**2))
fig = plt.figure(figsize=[6,3]); plt.title('Cost function $c(x)=\exp(-\|x\|^2/2)$'); plt.axis('equal');
plt.contourf(X,Y,hfmIn['cost']);
savefig(fig,'CostFunction.png')
hfmOut =hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.003627 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=[6,3]); plt.title('Distance and geodesics, with obstacles and non-uniform cost'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values'],cmap='Greys') # Display the value function.
for geo in hfmOut['geodesics']: plt.plot(*geo)
savefig(fig,'Dist_Paths_CostObstacles.png')
if hfmIn.mode=='gpu': raise ad.DeliberateNotebookError(
"Some of the functionality described in the rest of this notebook is absent from the GPU eikonal solver")
The HFM software by default terminates when eikonal PDE is solved on the entire domain $\Omega$, by propagating a front in a single pass, in a Dijkstra-like manner. However it can be interesting to abort the computation early, in order for instance to save CPU time. Various stopping criteria are available for that purpose.
We provide an option for stopping the front propagation as soon as every member of a collection of points is reached. In many applications, the minimal geodesics are more valuable than the distance map itself. In that case, one may stop the front propagation as soon as all the geodesic tips are reached. Computation time is reduced, and geodesic backtracking is still feasible.
hfmIn['stopWhenAllAccepted']=hfmIn['tips'] # Abort computation when all the tips are reached by the front
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.002661 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=[6,3]); plt.title('Stop propagation when all tips are accepted'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values'],cmap='Greys'); # Display the value function.
for geo in hfmOut['geodesics']: plt.plot(*geo)
savefig(fig,'StopAllTipsAccepted.png')
A second option stops the front propagation as soon as any member of a collection of points is reached. For illustration, we stop as soon as the first geodesic tip is reached. The geodesic originating from this particular point can be backtraced, but obviously the others cannot anymore.
hfmIn['stopWhenAnyAccepted']=hfmIn['tips']
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.000965 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=[6,3]); plt.title('Stop propagation when first tip is accepted'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values'],cmap='Greys'); # Display the value function.
for geo in hfmOut['geodesics']: plt.plot(*geo)
# We won't use these stopping criteria anymore.
hfmIn.pop('stopWhenAllAccepted',None);
hfmIn.pop('stopWhenAnyAccepted',None);
We here decide to terminate the front propagation once the value function, i.e. the distance from the seeds (plus possibly the boundary condition), reaches a provided threshold.
hfmIn['stopAtDistance']=0.7
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.002787 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=[6,3]); plt.title('Stop propagation at a prescribed distance'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values'],cmap='Greys'); # Display the value function.
for geo in hfmOut['geodesics']: plt.plot(*geo)
hfmIn.pop('stopAtDistance',None); # We won't use this stopping criterion anymore
The HFM software can compute several side products in addition to the (approximate) distance map $u : \Omega \to ]-\infty,\infty]$. In this section, we present two of them: the euclidean length of the minimal geodesics, and the Voronoi diagram associated with the seeds. Since these computations are performed simultaneously with the front propagation, they are also used within stopping criteria.
The Voronoi cell of a seed consists of all the points which are closer to this seed than to the others.
hfmIn['seedFlags']=[0,1] # Define voronoi region indices associated with the seeds.
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Field voronoiStoppingCriterion defaults to None Fast marching solver completed in 0.003823 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45 Field exportVoronoiFlags defaults to 1
fig = plt.figure(figsize=[6,3]); plt.title('Voronoi regions'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['voronoiFlags'],cmap='Greys'); # Display the Voronoi regions.
for geo in hfmOut['geodesics']: plt.plot(*geo)
savefig(fig,'VoronoiRegions.png')
In next experiment, we terminate the computation as soon as the two Voronoi regions $V_0,V_1$ meet. More precisely, as soon as they contain points $x_0,x_1$ that differ by a single unit at a single coordinate \begin{equation*} x_1-x_0 \in \{(0,\pm 1), (\pm 1,0)\}. \end{equation*} We also display the geodesics backtraced from $x_0$ and $x_1$. By construction, they can be concatenated into a single minimal geodesic joining the two seeds.
hfmIn['voronoiStoppingCriterion']='RegionsMeeting'
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.003027 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45 Field exportVoronoiFlags defaults to 1 Field voronoiDiagram_exportGeodesicFromMeetingPoint defaults to 1
print("The Voronoi regions meet at the neighbor points x0 =",hfmOut['voronoiDiagram_meetingPoint0'],
"and x1 =", hfmOut['voronoiDiagram_meetingPoint1'])
The Voronoi regions meet at the neighbor points x0 = [0.225 0.725] and x1 = [0.215 0.725]
fig = plt.figure(figsize=[6,3]); plt.title('Propagation stopped when Voronoi regions met'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['voronoiFlags'],cmap='Greys'); # Display the Voronoi regions.
for geo in hfmOut['geodesics_voronoiDiagram']: plt.plot(*geo)
hfmIn.pop('voronoiStoppingCriterion',None); # We won't use this stopping criterion anymore.
A simple addition to the numerical scheme allows to compute the euclidean length of the minimal geodesics from the seeds towards the points of the fast marching front. Optionally, we may to interrupt the computations once these lengths reach a prescribed threshold.
The scale of the grid must be specified in each direction. Optionally, this scale may vary over the domain.
gridScale=hfmIn['gridScale']
hfmIn['euclideanScale']=[gridScale,gridScale] # euclidean scale of the grid in the x and y directions
hfmIn['stopAtEuclideanLength']=1. # Optional. Stopping criterion based on euclidean length of geodesics.
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Field voronoiStoppingCriterion defaults to None Fast marching solver completed in 0.003707 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45 Field exportEuclideanLengths defaults to 1 Field euclideanLength_exportGeodesicFromStoppingPoint defaults to 1 Field exportVoronoiFlags defaults to 1
hfmOut.keys()
dict_keys(['FMCPUTime', 'GeodesicCPUTime', 'MaxStencilWidth', 'StencilCPUTime', 'defaulted', 'euclideanLength_stoppingPoint', 'euclideanLengths', 'geodesics', 'geodesics_euclideanLength', 'nAccepted', 'unusedFromCompute', 'values', 'visitedUnset', 'voronoiFlags', 'retcode'])
print("Euclidean length stopping criterion was triggered at :", hfmOut['euclideanLength_stoppingPoint'])
Euclidean length stopping criterion was triggered at : [0.445 0.145]
fig=plt.figure(figsize=[6,3]); plt.title('Euclidean length of geodesics, shown, used as stopping criterion'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['euclideanLengths'],cmap='Greys'); # Display the geodesic euclidean lengths.
# Display the point that triggered the stopping criterion, and the corresponding minimal geodesic
plt.plot(*hfmOut['geodesic_euclideanLength'])
savefig(fig,'EuclideanLength.png')
As a consistency check, we compare the value function with the euclidean length of geodesics, the latter being computed with an adequately chosen point dependent scale. The two quantities are equal up to machine precision, as expected, on the voronoi region of the left seed. On the voronoi region of the right seed, they differ by the seed boundary value, namely $0.5$.
hfmIn['euclideanScale'] = gridScale*np.stack((hfmIn['cost'],hfmIn['cost']),2)
hfmIn.pop('stopAtEuclideanLength',None)
hfmIn['order']=1 # Disable second order enhancement, so as to get an exact match of the two methods.
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Field voronoiStoppingCriterion defaults to None Fast marching solver completed in 0.003823 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45 Field exportEuclideanLengths defaults to 1 Field exportVoronoiFlags defaults to 1
plt.figure(figsize=[6,3]); plt.title('Distance and geodesic length are consistent'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['euclideanLengths']-hfmOut['values']);
# Restore parameters
hfmIn['order']=2.
hfmIn.pop('euclideanScale',None);
In this section, we consider a path planning problem in which the speed function depends on both position and time \begin{equation*} s : \Omega \times \mathbb R \to ]0,\infty[. \end{equation*} The corresponding eikonal equation reads \begin{align*} \forall x \in \Omega, s(x,u(x)) \| \nabla u(x) \| &= 1, & \forall x \in \partial \Omega, u(x) &= \sigma(x). \end{align*} Recall that the speed function is the reciprocal of the cost function, i.e. $s=1/c$. The numerical example below can trivially be adapted to the (equivalent) case of a time dependent cost function.
In order to define a time dependent speed function, one must provide an increasing family of times $t_0<t_1< \ldots <t_K$ is provided, and a family of speed functions at these time points $s_0,s_1,\ldots, s_K : \Omega \to ]0,\infty[$. In each time interval, the HFM software performs linear interpolation: \begin{align*} s(x,t) &= (1-\delta) s_k(x) + \delta s_{k+1}(x) & \text{ if } t &= (1-\delta)t_k+\delta t_{k+1}, \end{align*} for some $0 \leq k < K$, and some $\delta \in [0,1]$. The speed function is prolongated by a constant, w.r.t. time, before $t_0$ and after $t_K$ \begin{align*} s(x,t)&=s_0(x) \text{ if } t \leq t_0, & s(x,t)&=s_K(x) \text{ if } t \geq t_K. \end{align*}
hfmIn.pop('cost',None) # Remove the previous cost function
hfmIn['speed_times'] = [0.2,0.6,1.4] # Define the interpolation times t0,t1,t2
hfmIn['speed_0'] = np.exp(0.5*(X**2+Y**2)) # Initial speed function s0
hfmIn['speed_1'] = 0.2 # Intermediate speed function s1, defined as a small constant over the domain
hfmIn['speed_2'] = hfmIn['speed_0'] # Last speed function, chosen equal to the first one
hfmOut = hfmIn.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Field voronoiStoppingCriterion defaults to None Fast marching solver completed in 0.00402 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 8.45 Field exportVoronoiFlags defaults to 1
fig = plt.figure(figsize=[6,3]); plt.title('Case of a time dependent speed function'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values'],cmap='Greys'); plt.axis('equal'); # Display the value function.
for geo in hfmOut['geodesics']: plt.plot(*geo)
savefig(fig,'TimeDependent.png')
In this section, we demonstrate three dimensional isotropic fast marching. The mayavi plugin is used for 3D surface plotting. Like numpy.meshgrid, it has the convention of transposing the first two indices of the arrays.
Like numpy.meshgrid, mayavi has the convention of transposing the first two indices of the arrays.
hfmIn3D = Eikonal.dictIn({
'model':'Isotropic3',
'order':2,
'exportValues':1,
'seed':[0.3,0.5,0.7],
})
# Define the domain as the block [0,3] x [0,2] x [0,1]
hfmIn3D.SetRect(sides=[[0,3],[0,2],[0,1]],dimx = 150)
X3D,Y3D,Z3D = hfmIn3D.Grid()
hfmIn3D['speed'] = np.exp(-(X3D-1.5)**2-(Y3D-1)**2-(Z3D-0.5)**2)
hfmIn3D.SetUniformTips((3,3,3))
hfmOut = hfmIn3D.Run()
Field verbosity defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.336461 s. Field geodesicSolver defaults to Discrete Field geodesicStep defaults to 0.25 Field geodesicWeightThreshold defaults to 0.001 Field geodesicVolumeBound defaults to 10.985
plt.axes(projection='3d')
for geo in hfmOut['geodesics']: plt.plot(*geo)
if useMayavi:
mlab.contour3d(hfmOutput['values'], contours=[1,3,5,10])
mlab.show() # Displays in an external window.
This is what the previous command should display.
This section is devoted to diagonal metrics, a classical, minor but often useful, generalization of the isotropic considered above. The path cost takes the form \begin{equation*} L_M(\gamma) := \int_0^1 \| \gamma'(t)\|_{M(\gamma(t))} \,\mathrm d t, \end{equation*} where $M(x)$ is a diagonal matrix, dependent on the current position $x \in \Omega$. Recall that $\|v\|_M = \sqrt{<M v,v>}$. The cost can be rewritten in the form \begin{equation*} \| \dot x \|_{M(x)} = \sqrt{\sum_{1 \leq i \leq d} c_i(x)^2 \dot x_i^2}, \end{equation*} where $\dot x = (\dot x_1,\cdots, \dot x_d)$, and where the diagonal coefficients of the metric are $(c_1(x)^2,\cdots, c_d(x)^2)$.
The user inputs $d$ cost functions $c_1,\cdots,c_d : \overline \Omega \to ]0,\infty[$, where $d$ is the dimension, instead of a single one in the isotropic case. All the tools described in the isotropic case of course also apply to the diagonal models, as well as to the more sophisticated models described in the further notebooks.
# The three dimensional counterpart Diagonal3 is also implemented
hfmIn2 = Eikonal.dictIn({
'model':'Diagonal2',
'cost':np.stack(( np.exp(-0.5*(X**2+Y**2)), np.ones(X.shape)),axis=2),
})
for key in ['walls','dims','seeds','seedValues','tips','exportValues','origin','gridScale']:
hfmIn2[key]=hfmIn[key]
hfmOut = hfmIn2.Run()
Field verbosity defaults to 1 Field order defaults to 1 Field seedRadius defaults to 0 Fast marching solver completed in 0.003189 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=[6,3]); plt.title('Distance and minimal paths, for a diagonal metric'); plt.axis('equal');
plt.contourf(X,Y,hfmOut['values'],cmap='Greys') # Display the value function.
for geo in hfmOut['geodesics']: plt.plot(*geo)
savefig(fig,'Diagonal.png')