$\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} \renewcommand{\phi}{\varphi} \renewcommand{\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.
%load_ext pymatbridge
Starting MATLAB on ZMQ socket ipc:///tmp/pymatbridge Send 'exit' command to kill the server .MATLAB started and connected!
%%matlab
warning off
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/fastmarching_4bis_geodesic_mesh')
You need to download numerical_tours and install the IPython notebook to run the code.
You must also install the python-matlab-bridge.
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.
%%matlab
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.
%%matlab
clf;
options.lighting = 1;
plot_mesh(vertex0,faces0,options);
shading('faceted');
Perform sub-division to obtained a finer mesh.
%%matlab
nsub = 2;
% number of subdivision steps
options.sub_type = 'loop';
options.verb = 0;
% [vertex,faces] = perform_mesh_subdivision(vertex0,faces0,nsub,options);
vertex = vertex0;
faces = faces0;
n = size(vertex,2);
Display.
%%matlab
clf;
options.lighting = 1;
plot_mesh(vertex,faces,options);
shading('faceted');
Useful shortcuts.
%%matlab
dotp = @(u,v)sum(u.*v,1);
R = @(u)reshape(u, [1 1 length(u)]);
Shortcut for the explicit inverse of a series of 2D matrices, and the multiplication with a vector.
%%matlab
Inv1 = @(M,d)[M(2,2,:)./d -M(1,2,:)./d; -M(2,1,:)./d M(1,1,:)./d];
Inv = @(M)Inv1(M, M(1,1,:).*M(2,2,:) - M(1,2,:).*M(2,1,:));
Mult = @(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.
%%matlab
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.
%%matlab
I = [88 2602 883 23];
Initialization of the distance map.
%%matlab
U = zeros(n,1);
Computation of the set of indexes for $i,j,k$.
%%matlab
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$.
%%matlab
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$.
%%matlab
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}$.
%%matlab
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$.
%%matlab
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.
%%matlab
delta = max( b.^2 - a.*c, 0);
The solution is then $$ d = \frac{ b + \sqrt{\delta} }{ a }. $$
%%matlab
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$.
%%matlab
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 }. $$
%%matlab
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)$.
%%matlab
U1 = accumarray(i', d, [n 1], @min); U1(U1==0) = Inf;
Enforce the boundary conditions $\forall i \in I, \, U_i = 0$.
%%matlab
U1(I) = 0;
Update the solution.
%%matlab
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
%%matlab
exo1()
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) <ipython-input-24-9872be35b9e1> in <module>() ----> 1 get_ipython().run_cell_magic('matlab', '', 'exo1()') /Users/gpeyre/anaconda/envs/py3k/lib/python3.3/site-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell) 2160 magic_arg_s = self.var_expand(line, stack_depth) 2161 with self.builtin_trap: -> 2162 result = fn(magic_arg_s, cell) 2163 return result 2164 /Users/gpeyre/anaconda/envs/py3k/lib/python3.3/site-packages/pymatbridge/matlab_magic.py in matlab(self, line, cell, local_ns) /Users/gpeyre/anaconda/envs/py3k/lib/python3.3/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k) 191 # but it's overkill for just that one bit of state. 192 def magic_deco(arg): --> 193 call = lambda f, *a, **k: f(*a, **k) 194 195 if callable(arg): /Users/gpeyre/anaconda/envs/py3k/lib/python3.3/site-packages/pymatbridge/matlab_magic.py in matlab(self, line, cell, local_ns) 215 e_s += "\n-----------------------" 216 e_s += "\nAre you sure Matlab is started?" --> 217 raise RuntimeError(e_s) 218 219 RuntimeError: There was an error running the code: exo1() ----------------------- Are you sure Matlab is started?
%%matlab
%% Insert your code here.
Display the fixed point energy.
%%matlab
clf;
h = plot(err);
set(h, 'LineWidth', 2);
axis tight;
For the display of the distance.
%%matlab
mycolor = @(U,k)mod(U/max(U), 1/k);
mycolor = @(U,k)cos(2*pi*k*U/max(U));
Display the distance map.
%%matlab
clf;
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.
%%matlab
k = 30;
pend = round(rand(k,1)*n)+1;
options.method = 'continuous';
paths = compute_geodesic_mesh(U, vertex, faces, pend, options);
Display.
%%matlab
clf;
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.
%%matlab
n = 3000;
vertex = [2*rand(2,n)-1; zeros(1,n)];
vertex(:,1) = 0;
faces = compute_delaunay(vertex);
options.name = '';
I = 1;
Display.
%%matlab
clf
plot_mesh(vertex,faces);
shading faceted;
Exercise 2
Compute and display the geodesic distance.
%%matlab
exo2()
%%matlab
%% 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).
%%matlab
exo3()
%%matlab
%% Insert your code here.
One can define a non-constant metric on the mesh.
Triangulation of a square.
%%matlab
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.
%%matlab
W = ones(n,1);
W(vertex(1,:)<.5) = 1/2;
Display the metric.
%%matlab
options.face_vertex_color = W;
clf;
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.
%%matlab
exo4()
%%matlab
%% Insert your code here.