{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Geodesic Distance Computation on 3-D Meshes", "===========================================", "", "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) 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}$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This tour explores the computation of the geodesic distance on a mesh using", "an iterative algorithm to solve the Eikonal equation."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Triangulated Mesh", "-----------------", "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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.lighting = 1", "plot_mesh(vertex0, faces0, options)", "shading('faceted')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform sub-division to obtained a finer mesh."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.lighting = 1", "plot_mesh(vertex, faces, options)", "shading('faceted')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Useful shortcuts."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["dotp = lambda u, v: sum(u.*v, 1)", "R = lambda u: reshape(u, [1 1 length(u)])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Shortcut for the explicit inverse of a series of 2D matrices, and the multiplication with a vector."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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, : )]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Control Formulation of the Eikonal Equation", "-------------------------------------------", "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.", "", "Geodesic distance Computation Using Jacobi Iterations", "-----------------------------------------------------", "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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = ones(n, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["I = [88 2602 883 23]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Initialization of the distance map."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["U = zeros(n, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Computation of the set of indexes for $i,j,k$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["i = [faces(1, : ) faces(2, : ) faces(3, : )]", "j = [faces(2, : ) faces(3, : ) faces(1, : )]", "k = [faces(3, : ) faces(1, : ) faces(2, : )]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["|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$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["x  = vertex(: , i)", "x1 = vertex(: , j) - x", "x2 = vertex(: , k) - x"]}, {"cell_type": "markdown", "metadata": {}, "source": ["|ui| stores $U_{j}$, |uk| stores $U_k$, |w| stores $w=W_i$.", "We denote $u=(U_j,U_k) \\in \\RR^2$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["uj = U(j)", "uk = U(k)", "u = [R(uj); R(uk)]", "w = R(W(i))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["|C| stores the inner product matrix $C_{s,t} = \\dotp{x_s}{x_t}$,", "and we denote $S = C^{-1}$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["C = [R(dotp(x1, x1)) R(dotp(x1, x2)); R(dotp(x2, x1)) R(dotp(x2, x2))]", "S = Inv(C)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a = sum(sum(S))", "b = dotp(sum(S, 2), u)", "c = dotp(Mult(S, u), u) - w.^2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the reduced discriminant $\\de = b^2 - a c$, which is always", "positive."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["delta = max(b.^2 - a.*c, 0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The solution is then", "$$ d = \\frac{ b + \\sqrt{\\delta} }{ a }. $$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["d = (b + sqrt(delta))./ a"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["alpha = Mult(S, u - repmat(d, 2, 1))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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 }. $$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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)$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["U1 = accumarray(i', d, [n 1], @min);  U1(U1 = =0) = Inf"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Enforce the boundary conditions $\\forall i \\in I, \\, U_i = 0$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["U1(I) = 0"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Update the solution."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["U = U1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__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"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo1()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the fixed point energy."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["h = plot(err)", "set(h, 'LineWidth', 2)", "axis tight"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For the display of the distance."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["mycolor = lambda U, k: mod(U/ max(U), 1/ k)", "mycolor = lambda U, k: cos(2*pi*k*U/ max(U))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the distance map."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.face_vertex_color = mycolor(U, 5)", "plot_mesh(vertex, faces, options)", "colormap jet(256)", "shading interp"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Geodesic Path Extraction", "------------------------", "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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["k = 30", "pend = round(rand(k, 1)*n) + 1", "options.method = 'continuous'", "paths = compute_geodesic_mesh(U, vertex, faces, pend, options)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot_fast_marching_mesh(vertex, faces, mycolor(U, 5), paths, options)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Geodesic Distance Convergence", "-----------------------------", "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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 3000", "vertex = [2*rand(2, n)-1; zeros(1, n)]", "vertex(: , 1) = 0", "faces = compute_delaunay(vertex)", "options.name = ''", "I = 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot_mesh(vertex, faces)", "shading faceted"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Compute and display the geodesic distance."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo2()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__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)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo3()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Spacially Varying Metrics", "-------------------------", "One can define a non-constant metric on the mesh.", "", "", "Triangulation of a square."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 8000", "vertex = [rand(2, n); zeros(1, n)]", "vertex(: , 1) = [0 .3 0]'", "faces = compute_delaunay(vertex)", "I = 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate a binary metric in order to produce a reflexion effect."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = ones(n, 1)", "W(vertex(1, : ) <.5) = 1/ 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the metric."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.face_vertex_color = W", "", "plot_mesh(vertex, faces, options)", "colormap jet(256)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 4__", "", "Compute the geodesic distance for a metric $W_i$ that is not constant", "over the mesh."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo4()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}