Important: Please read the installation page for details about how to install the toolboxes. $\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{\umin}[1]{\underset{#1}{\min}\;}$ $\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{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This tour explores the computation of the geodesic distance on a mesh using an iterative algorithm to solve the Eikonal equation.
from __future__ import division
import nt_toolbox as nt
from nt_solutions import fastmarching_4bis_geodesic_mesh as solutions
%matplotlib inline
%load_ext autoreload
%autoreload 2
We consider a 3-D surface discretized using a triangular mesh with $N$ vertices. We denote $ \{x_i\}_{i=1}^{N} $ the set of vertices, where $x_i \in \RR^3$. We denote $\Ff \subset \{1, \ldots, N\}^3$ the set of faces.
First load a mesh.
name = 'beetle'; % 1k
name = 'camel'; % 600 v
name = 'cow'; % 700v
name = 'venus'; % 700 v
name = 'mannequin'; % 400 v
name = 'nefertiti'; % 300 v
[vertex0, faces0] = read_mesh(name)
clear options; options.name = name
Display.
options.lighting = 1
plot_mesh(vertex0, faces0, options)
shading('faceted')
Perform sub-division to obtained a finer mesh.
nsub = 2; % number of subdivision steps
options.sub_type = 'loop'
options.verb = 0
[vertex, faces] = perform_mesh_subdivision(vertex0, faces0, nsub, options)
n = size(vertex, 2)
Display.
options.lighting = 1
plot_mesh(vertex, faces, options)
shading('faceted')
Useful shortcuts.
dotp = lambda u, v: sum(u.*v, 1)
R = lambda u: reshape(u, [1 1 length(u)])
Shortcut for the explicit inverse of a series of 2D matrices, and the multiplication with a vector.
Inv1 = lambda M, d: [M(2, 2, : )./ d -M(1, 2, : )./ d; -M(2, 1, : )./ d M(1, 1, : )./ d]
Inv = lambda M: Inv1(M, M(1, 1, : ).*M(2, 2, : ) - M(1, 2, : ).*M(2, 1, : ))
Mult = lambda M, u: [M(1, 1, : ).*u(1, 1, : ) + M(1, 2, : ).*u(2, 1, : ); M(2, 1, : ).*u(1, 1, : ) + M(2, 2, : ).*u(2, 1, : )]
Our goal is to compute the geodesic distance to a set of starting points.
We consider an isotropic metric on the mesh $W_i > 0$.
Given a set $I$ of starting point, we want to compute the discrete geodesic distance map $ U \in \RR^N $ to the starting points $ \{x_i\}_{i \in I} $.
Let's consider a geodesic distance $d(x,y)$ on a manifold. A general mathematical result is that the geodesic distance $$ d(x,A) = \umin{a \in A} d(x,a) $$ to a set $A$ is the unique solution to the following equation $$ \forall \, x \notin A, \quad d(x,A) = \umin{ y \in B(x) } d(y,A) + d(x,y) $$ where $ B(x) $ is a small disk around $x$ that does not contain $A$.
In a discrete setting, one can use this property with $B(x_i)$ being the one ring, made of the edges connecting vertices that are around $x_i$. Using a piecewise linear interpolation of $U$ on the mesh, one obtain that $U$ is a solution to a fixed point equation $$ U = \Ga(U) $$ where the update operator $ \Ga : \RR^N \rightarrow \RR^N $ is defined as $$ \Ga(U)_i = \umin{ (i,j,k) \in \Ff } d_{i,j,k} $$ where $$ d_{i,j,k} = \umin{0 \leq t \leq 1} t U_{j} + (1-t) U_k + W_i \norm{ t x_j + (1-t) x_k - x }. $$ The quantities $ d_{i,j,k} $ can in fact be computed by solving a second order polynomial, as we describe in the following section.
The mapping $\Ga$ is monotone $$ U \leq V \qarrq \Ga(U) \leq \Ga(V). $$ If one uses an initialization $U^{(0)}$ such that $\Ga(U^{(0)}) \geq U^{(0)} $ (for instance setting $U^{(0)}=0$), then the iterates $$ U^{(\ell+1)} = \Ga(U^{(\ell)}) $$ are increasing. One can prove that they are bounded, and hence they converge to the (unique) fixed point solution of $\Ga(U)=U$.
Define the metric $W_i > 0$ on the mesh. We start with the uniform (constant) metric.
W = ones(n, 1)
Indexes $I \subset \{1,\ldots,N\}$ of starting points. Warning: this list works for the mesh |nefertitit| refined twice. For other meshes, use the function |select3dtool| to retrive the indexes of a few vertices.
I = [88 2602 883 23]
Initialization of the distance map.
U = zeros(n, 1)
Computation of the set of indexes for $i,j,k$.
i = [faces(1, : ) faces(2, : ) faces(3, : )]
j = [faces(2, : ) faces(3, : ) faces(1, : )]
k = [faces(3, : ) faces(1, : ) faces(2, : )]
|x| stores the position of all the $x_i$ vertices (point where the value of $d=d_{i,j,k}$ is computed), while |x1| and |x2| store respectively $x_j-x_i$ and $x_k-x_i$.
x = vertex(: , i)
x1 = vertex(: , j) - x
x2 = vertex(: , k) - x
|ui| stores $U_{j}$, |uk| stores $U_k$, |w| stores $w=W_i$. We denote $u=(U_j,U_k) \in \RR^2$.
uj = U(j)
uk = U(k)
u = [R(uj); R(uk)]
w = R(W(i))
|C| stores the inner product matrix $C_{s,t} = \dotp{x_s}{x_t}$, and we denote $S = C^{-1}$.
C = [R(dotp(x1, x1)) R(dotp(x1, x2)); R(dotp(x2, x1)) R(dotp(x2, x2))]
S = Inv(C)
If we denote $g = \al_1 x_1 + \al_2 x_2$ the (unknown) update direction, the geodesic distance (which is affine in the triangle), is written as $ U(z) = \dotp{g}{z-x} + d $ where $d \in \RR$ is unknown.
We need to compute both $g$ and $d$. Since $U(x_j)=u_1$ and $U(x_k)=u_2$, one obtains $$ \al = S ( u - d \mathbb{I} ), $$ where $\mathbb{I} = (1,1) \in \RR^2$. Making use of the Eikonal equation $ \norm{\nabla U(x)} = \norm{g} = w $, one obtains that $d$ is the (maximum) solution of the following second order equation $$ a d^2 - 2 b d + c = 0 \qwhereq \choice{ a = \dotp{S \mathbb{I}}{\mathbb{I}} \\ b = \dotp{S \mathbb{I}}{u} \\ c = \dotp{S u}{u} - w^2 } $$
Compute the values of $a, b, c$.
a = sum(sum(S))
b = dotp(sum(S, 2), u)
c = dotp(Mult(S, u), u) - w.^2
Compute the reduced discriminant $\de = b^2 - a c$, which is always positive.
delta = max(b.^2 - a.*c, 0)
The solution is then $$ d = \frac{ b + \sqrt{\delta} }{ a }. $$
d = (b + sqrt(delta))./ a
To ensure correctness of the scheme, the update should come from within the triangle (which correspond to the condition $$t \in [0,1]$$ in the original definition of $d_{i,j,k}$ using the control formulation). This corresponds to having $ \al_1<=0 $ and $\al_2 \leq 0$ where $ \al = S ( u - d \mathbb{I} ) $.
Compute $\al \in \RR^2$.
alpha = Mult(S, u - repmat(d, 2, 1))
For the index $i \in J$ where the update comes from outside the edge $ [x_j,x_k] $, we use an update along the edge $$ d = \min\pa{ U_j + \norm{x-x_j}W_i, U_k + \norm{x-x_k}W_i }. $$
J = find(alpha(1, 1, : ) >0 | alpha(2, 1, : ) >0)
d1 = sqrt(dotp(x1, x1)); d1 = d1(: ).*w(: ) + uj(: )
d2 = sqrt(dotp(x2, x2)); d2 = d2(: ).*w(: ) + uk(: )
d = d(: )
d(J) = min(d1(J), d2(J))
Now that |d| stores the value of |d_{i,j,k}| for a bunch of faces assign in |U1(i)| the mininimum between the previous |U1(i)| and all the entries in |d| that corresponds to a face $(i,j,k)$.
U1 = accumarray(i', d, [n 1], @min); U1(U1 = =0) = Inf
Enforce the boundary conditions $\forall i \in I, \, U_i = 0$.
U1(I) = 0
Update the solution.
U = U1
Exercise 1
Perform the full iterative algorithm until convergence. Store in |err(l)| the fixed point error $ \norm{ U^{(\ell+1)} - U^{(\ell)} } $. Note: you might want to put outside of the loop all the quantities that do not depend on $u$, e.g. $S, a, $ etc. isplay
solutions.exo1()
## Insert your code here.
Display the fixed point energy.
h = plot(err)
set(h, 'LineWidth', 2)
axis tight
For the display of the distance.
mycolor = lambda U, k: mod(U/ max(U), 1/ k)
mycolor = lambda U, k: cos(2*pi*k*U/ max(U))
Display the distance map.
options.face_vertex_color = mycolor(U, 5)
plot_mesh(vertex, faces, options)
colormap jet(256)
shading interp
A geodesic minimal path from any point on the mesh to the starting point can be computed using a gradient descent. Since the computed geodesic distance is affine on each triangle, it makes sense to define a discretized path that is a line segment on each face.
Extract minimal paths from all over the mesh, starting from the boundary.
k = 30
pend = round(rand(k, 1)*n) + 1
options.method = 'continuous'
paths = compute_geodesic_mesh(U, vertex, faces, pend, options)
Display.
plot_fast_marching_mesh(vertex, faces, mycolor(U, 5), paths, options)
The discrete geodesic distance converges to the continuous distance when the mesh is refined.
Compute a Delaunay triangulation of random point on a square with the first point in the center.
n = 3000
vertex = [2*rand(2, n)-1; zeros(1, n)]
vertex(: , 1) = 0
faces = compute_delaunay(vertex)
options.name = ''
I = 1
Display.
plot_mesh(vertex, faces)
shading faceted
Exercise 2
Compute and display the geodesic distance.
solutions.exo2()
## Insert your code here.
Exercise 3
Display the convergence of the computed geodesic distance to the the true geodesic distance (which is the Euclidean distance $ \norm{x_i} $) as $n$ increases. Note: the triangulation with increasing number of points should be refining (i.e. a finer triangulation should contains all the other ones).
solutions.exo3()
## Insert your code here.
One can define a non-constant metric on the mesh.
Triangulation of a square.
n = 8000
vertex = [rand(2, n); zeros(1, n)]
vertex(: , 1) = [0 .3 0]'
faces = compute_delaunay(vertex)
I = 1
Generate a binary metric in order to produce a reflexion effect.
W = ones(n, 1)
W(vertex(1, : ) <.5) = 1/ 2
Display the metric.
options.face_vertex_color = W
plot_mesh(vertex, faces, options)
colormap jet(256)
Exercise 4
Compute the geodesic distance for a metric $W_i$ that is not constant over the mesh.
solutions.exo4()
## Insert your code here.