matlab/meshproc_7_geodesic_poisson/index.html


<!DOCTYPE html
  PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN">
<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd">
   <head>
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
   
      <!--
This HTML is auto-generated from an M-file.
To make changes, update the M-file and republish this document.
      --><script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML"></script><p style="font-size:0px">
         \[
         \newcommand{\NN}{\mathbb{N}}
         \newcommand{\CC}{\mathbb{C}}
         \newcommand{\GG}{\mathbb{G}}
         \newcommand{\LL}{\mathbb{L}}
         \newcommand{\PP}{\mathbb{P}}
         \newcommand{\QQ}{\mathbb{Q}}
         \newcommand{\RR}{\mathbb{R}}
         \newcommand{\VV}{\mathbb{V}}
         \newcommand{\ZZ}{\mathbb{Z}}
         \newcommand{\FF}{\mathbb{F}}
         \newcommand{\KK}{\mathbb{K}}
         \newcommand{\UU}{\mathbb{U}}
         \newcommand{\EE}{\mathbb{E}}
         
         \newcommand{\Aa}{\mathcal{A}}
         \newcommand{\Bb}{\mathcal{B}}
         \newcommand{\Cc}{\mathcal{C}}
         \newcommand{\Dd}{\mathcal{D}}
         \newcommand{\Ee}{\mathcal{E}}
         \newcommand{\Ff}{\mathcal{F}}
         \newcommand{\Gg}{\mathcal{G}}
         \newcommand{\Hh}{\mathcal{H}}
         \newcommand{\Ii}{\mathcal{I}}
         \newcommand{\Jj}{\mathcal{J}}
         \newcommand{\Kk}{\mathcal{K}}
         \newcommand{\Ll}{\mathcal{L}}
         \newcommand{\Mm}{\mathcal{M}}
         \newcommand{\Nn}{\mathcal{N}}
         \newcommand{\Oo}{\mathcal{O}}
         \newcommand{\Pp}{\mathcal{P}}
         \newcommand{\Qq}{\mathcal{Q}}
         \newcommand{\Rr}{\mathcal{R}}
         \newcommand{\Ss}{\mathcal{S}}
         \newcommand{\Tt}{\mathcal{T}}
         \newcommand{\Uu}{\mathcal{U}}
         \newcommand{\Vv}{\mathcal{V}}
         \newcommand{\Ww}{\mathcal{W}}
         \newcommand{\Xx}{\mathcal{X}}
         \newcommand{\Yy}{\mathcal{Y}}
         \newcommand{\Zz}{\mathcal{Z}}
         
         \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}
         \renewcommand{\epsilon}{\varepsilon}
         
         \newcommand{\Calpha}{\mathrm{C}^\al}
         \newcommand{\Cbeta}{\mathrm{C}^\be}
         \newcommand{\Cal}{\text{C}^\al}
         \newcommand{\Cdeux}{\text{C}^{2}}
         \newcommand{\Cun}{\text{C}^{1}}
         \newcommand{\Calt}[1]{\text{C}^{#1}}
         
         \newcommand{\lun}{\ell^1}
         \newcommand{\ldeux}{\ell^2}
         \newcommand{\linf}{\ell^\infty}
         \newcommand{\ldeuxj}{{\ldeux_j}}
         \newcommand{\Lun}{\text{\upshape L}^1}
         \newcommand{\Ldeux}{\text{\upshape L}^2}
         \newcommand{\Lp}{\text{\upshape L}^p}
         \newcommand{\Lq}{\text{\upshape L}^q}
         \newcommand{\Linf}{\text{\upshape L}^\infty}
         \newcommand{\lzero}{\ell^0}
         \newcommand{\lp}{\ell^p}
         
         
         \renewcommand{\d}{\ins{d}}
         
         \newcommand{\Grad}{\text{Grad}}
         \newcommand{\grad}{\text{grad}}
         \renewcommand{\div}{\text{div}}
         \newcommand{\diag}{\text{diag}}
         
         \newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }
         \newcommand{\pdd}[2]{ \frac{ \partial^2 #1}{\partial #2^2} }
         
         \newcommand{\dotp}[2]{\langle #1,\,#2\rangle}
         \newcommand{\norm}[1]{|\!| #1 |\!|}
         \newcommand{\normi}[1]{\norm{#1}_{\infty}}
         \newcommand{\normu}[1]{\norm{#1}_{1}}
         \newcommand{\normz}[1]{\norm{#1}_{0}}
         \newcommand{\abs}[1]{\vert #1 \vert}
         
         
         \newcommand{\argmin}{\text{argmin}}
         \newcommand{\argmax}{\text{argmax}}
         \newcommand{\uargmin}[1]{\underset{#1}{\argmin}\;}
         \newcommand{\uargmax}[1]{\underset{#1}{\argmax}\;}
         \newcommand{\umin}[1]{\underset{#1}{\min}\;}
         \newcommand{\umax}[1]{\underset{#1}{\max}\;}
         
         \newcommand{\pa}[1]{\left( #1 \right)}
         \newcommand{\choice}[1]{ \left\{  \begin{array}{l} #1 \end{array} \right. }
         
         \newcommand{\enscond}[2]{ \left\{ #1 \;:\; #2 \right\} }
         
         \newcommand{\qandq}{ \quad \text{and} \quad }
         \newcommand{\qqandqq}{ \qquad \text{and} \qquad }
         \newcommand{\qifq}{ \quad \text{if} \quad }
         \newcommand{\qqifqq}{ \qquad \text{if} \qquad }
         \newcommand{\qwhereq}{ \quad \text{where} \quad }
         \newcommand{\qqwhereqq}{ \qquad \text{where} \qquad }
         \newcommand{\qwithq}{ \quad \text{with} \quad }
         \newcommand{\qqwithqq}{ \qquad \text{with} \qquad }
         \newcommand{\qforq}{ \quad \text{for} \quad }
         \newcommand{\qqforqq}{ \qquad \text{for} \qquad }
         \newcommand{\qqsinceqq}{ \qquad \text{since} \qquad }
         \newcommand{\qsinceq}{ \quad \text{since} \quad }
         \newcommand{\qarrq}{\quad\Longrightarrow\quad}
         \newcommand{\qqarrqq}{\quad\Longrightarrow\quad}
         \newcommand{\qiffq}{\quad\Longleftrightarrow\quad}
         \newcommand{\qqiffqq}{\qquad\Longleftrightarrow\qquad}
         \newcommand{\qsubjq}{ \quad \text{subject to} \quad }
         \newcommand{\qqsubjqq}{ \qquad \text{subject to} \qquad }
         \]
         
      </p>
      <title>Geodesic Distance with Poisson Equation</title>
      <NOSCRIPT>
         <DIV STYLE="color:#CC0000; text-align:center"><B>Warning: <A HREF="http://www.math.union.edu/locate/jsMath">jsMath</A> 
               	requires JavaScript to process the mathematics on this page.<BR> 
               	If your browser supports JavaScript, be sure it is enabled.</B></DIV>
         <HR>
      </NOSCRIPT>
      <meta name="generator" content="MATLAB 8.4">
      <meta name="date" content="2014-11-27">
      <meta name="m-file" content="index">
      <LINK REL="stylesheet" HREF="../style.css" TYPE="text/css">
   </head>
   <body>
      <div class="content">
         <h1>Geodesic Distance with Poisson Equation</h1>
         <introduction>
            <p>This tour explores the method detailed in the paper <a href="#biblio">[CraneWeischedelWardetzky13]</a> to compute geodesic distance by simply solving a Poisson equation.
            </p>
         </introduction>
         <h2>Contents</h2>
         <div>
            <ul>
               <li><a href="#1">Installing toolboxes and setting up the path.</a></li>
               <li><a href="#8">Gradient, Divergence and Laplacian on Surfaces</a></li>
               <li><a href="#30">Heat Diffusion and Time Stepping</a></li>
               <li><a href="#37">Geodesic in Heat Method</a></li>
               <li><a href="#48">Bibliography</a></li>
            </ul>
         </div>
         <h2>Installing toolboxes and setting up the path.<a name="1"></a></h2>
         <p>You need to download the following files: <a href="../toolbox_signal.zip">signal toolbox</a>, <a href="../toolbox_general.zip">general toolbox</a> and <a href="../toolbox_graph.zip">graph toolbox</a>.
         </p>
         <p>You need to unzip these toolboxes in your working directory, so that you have <tt>toolbox_signal</tt>, <tt>toolbox_general</tt> and <tt>toolbox_graph</tt> in your directory.
         </p>
         <p><b>For Scilab user:</b> you must replace the Matlab comment '%' by its Scilab counterpart '//'.
         </p>
         <p><b>Recommandation:</b> You should create a text file named for instance <tt>numericaltour.sce</tt> (in Scilab) or <tt>numericaltour.m</tt> (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run <tt>exec('numericaltour.sce');</tt> (in Scilab) or <tt>numericaltour;</tt> (in Matlab) to run the commands.
         </p>
         <p>Execute this line only if you are using Matlab.</p><pre class="codeinput">getd = @(p)path(p,path); <span class="comment">% scilab users must *not* execute this</span>
</pre><p>Then you can add the toolboxes to the path.</p><pre class="codeinput">getd(<span class="string">'toolbox_signal/'</span>);
getd(<span class="string">'toolbox_general/'</span>);
getd(<span class="string">'toolbox_graph/'</span>);
</pre><pre class="codeoutput">Warning: Name is nonexistent or not a directory: toolbox_signal 
Warning: Name is nonexistent or not a directory: toolbox_general 
Warning: Name is nonexistent or not a directory: toolbox_graph 
</pre><h2>Gradient, Divergence and Laplacian on Surfaces<a name="8"></a></h2>
         <p>The topology of a triangulation is defined via a set of indexes \(\Vv = \{1,\ldots,n\}\) that indexes the \(n\) vertices,
            a set of edges \(\Ee \subset \Vv \times \Vv\) and a set of \(m\) faces \(\Ff \subset \Vv  \times \Vv \times \Vv\).
         </p>
         <p>We load a mesh. The set of faces \(\Ff\) is stored in a matrix \(F \in \{1,\ldots,n\}^{3 \times m}\). The positions \(x_i
            \in \RR^3\), for \(i \in V\), of the \(n\) vertices are stored in a matrix \(X_0 = (x_{0,i})_{i=1}^n \in \RR^{3 \times n}\).
         </p><pre class="codeinput">clear <span class="string">options</span>;
name = <span class="string">'elephant-50kv'</span>;
options.name = name; <span class="comment">% useful for displaying</span>
[X,F] = read_mesh(name);
</pre><p>Number \(n\) of vertices and number \(m\) of faces.</p><pre class="codeinput">n = size(X,2);
m = size(F,2);
</pre><p>Display the mesh in 3-D.</p><pre class="codeinput">options.lighting = 1;
clf;
plot_mesh(X,F,options);
axis(<span class="string">'tight'</span>);
</pre><img vspace="5" hspace="5" src="index_01.png"> <p>In this tour, we use piecewise linear finite element to compute the gradient operators, which in turns allows us to compute
            the divergence (its transposed operator) and the Laplacian (as the composition of the divergence and the gradient).
         </p>
         <p>The gradient operator \(\nabla\) can be understood as a collection of 3 sparse matrices \((\nabla_s)_{s=1,2,3}\) of size \((m,n)\)
            that computes each coordinate of \(\nabla u=(\nabla_s u)_{s=1,2,3}\) through the formula, for each face \(f\), \[ (\nabla
            u)_f = \frac{1}{2A_f} \sum_{i \in f} u_i (N_f \wedge e_i) \] where \(A_f\) is the area of face \(f\), \(N_f\) is the normal
            to the face, \(e_i\) is the edge opposite to vertex \(i\), and \(\wedge\) is the cross product in \(\RR^3\).
         </p>
         <p>Callback to get the coordinates of all the vertex of index \(i=1,2,3\) in all faces.</p><pre class="codeinput">XF = @(i)X(:,F(i,:));
</pre><p>Compute un-normalized normal through the formula \(e_1 \wedge e_2 \) where \(e_i\) are the edges.</p><pre class="codeinput">Na = cross( XF(2)-XF(1), XF(3)-XF(1) );
</pre><p>Compute the area of each face as half the norm of the cross product.</p><pre class="codeinput">amplitude = @(X)sqrt( sum( X.^2 ) );
A = amplitude(Na)/2;
</pre><p>Compute the set of unit-norm normals to each face.</p><pre class="codeinput">normalize = @(X)X ./ repmat(amplitude(X), [3 1]);
N = normalize(Na);
</pre><p>Populate the sparse entries of the matrices for the operator implementing \( \sum_{i \in f} u_i (N_f \wedge e_i) \).</p><pre class="codeinput">I = []; J = []; V = []; <span class="comment">% indexes to build the sparse matrices</span>
<span class="keyword">for</span> i=1:3
    <span class="comment">% opposite edge e_i indexes</span>
    s = mod(i,3)+1;
    t = mod(i+1,3)+1;
    <span class="comment">% vector N_f^e_i</span>
    wi = cross(XF(t)-XF(s),N);
    <span class="comment">% update the index listing</span>
    I = [I, 1:m];
    J = [J, F(i,:)];
    V = [V, wi];
<span class="keyword">end</span>
</pre><p>Sparse matrix with entries \(1/(2A_f)\).</p><pre class="codeinput">dA = spdiags(1./(2*A(:)),0,m,m);
</pre><p>Compute gradient.</p><pre class="codeinput">GradMat = {};
<span class="keyword">for</span> k=1:3
    GradMat{k} = dA*sparse(I,J,V(k,:),m,n);
<span class="keyword">end</span>
</pre><p>\(\nabla\) gradient operator.</p><pre class="codeinput">Grad = @(u)[GradMat{1}*u, GradMat{2}*u, GradMat{3}*u]';
</pre><p>Compute divergence matrices as transposed of grad for the face area inner product.</p><pre class="codeinput">dAf = spdiags(2*A(:),0,m,m);
DivMat = {GradMat{1}'*dAf, GradMat{2}'*dAf, GradMat{3}'*dAf};
</pre><p>Div operator.</p><pre class="codeinput">Div = @(q)DivMat{1}*q(1,:)' + DivMat{2}*q(2,:)' + DivMat{3}*q(3,:)';
</pre><p>Laplacian operator as the composition of grad and div.</p><pre class="codeinput">Delta = DivMat{1}*GradMat{1} + DivMat{2}*GradMat{2} + DivMat{3}*GradMat{3};
</pre><p>Cotan of an angle between two vectors.</p><pre class="codeinput">cota = @(a,b)cot( acos( dot(normalize(a),normalize(b)) ) );
</pre><p>Compute cotan weights Laplacian.</p><pre class="codeinput">I = []; J = []; V = []; <span class="comment">% indexes to build the sparse matrices</span>
Ia = []; Va = []; <span class="comment">% area of vertices</span>
<span class="keyword">for</span> i=1:3
    <span class="comment">% opposite edge e_i indexes</span>
    s = mod(i,3)+1;
    t = mod(i+1,3)+1;
    <span class="comment">% adjacent edge</span>
    ctheta = cota(XF(s)-XF(i), XF(t)-XF(i));
    <span class="comment">% ctheta = max(ctheta, 1e-2); % avoid degeneracy</span>
    <span class="comment">% update the index listing</span>
    I = [I, F(s,:), F(t,:)];
    J = [J, F(t,:), F(s,:)];
    V = [V, ctheta, ctheta];
    <span class="comment">% update the diagonal with area of face around vertices</span>
    Ia = [Ia, F(i,:)];
    Va = [Va, A];
<span class="keyword">end</span>
<span class="comment">% Aread diagonal matrix</span>
Ac = sparse(Ia,Ia,Va,n,n);
<span class="comment">% Cotan weights</span>
Wc = sparse(I,J,V,n,n);
<span class="comment">% Laplacian with cotan weights.</span>
DeltaCot = spdiags(full(sum(Wc))', 0, n,n) - Wc;
</pre><p>Check that the Laplacian with cotan weights is actually equal to the composition of divergence and gradient.</p><pre class="codeinput">fprintf(<span class="string">'Should be 0: %e\n'</span>, norm(Delta-DeltaCot, <span class="string">'fro'</span>)/norm(Delta, <span class="string">'fro'</span>));
</pre><pre class="codeoutput">Should be 0: 1.604879e-10
</pre><p>Display a function \(f\) on the mesh.</p><pre class="codeinput">f = X(2,:);
options.face_vertex_color = f(:);
clf; plot_mesh(X,F,options);
axis(<span class="string">'tight'</span>);
colormap <span class="string">parula(256)</span>;
</pre><img vspace="5" hspace="5" src="index_02.png"> <p>Display its Laplacian.</p><pre class="codeinput">g = Delta*f(:);
g = clamp(g, -3*std(g), 3*std(g));
options.face_vertex_color = rescale(g);
clf; plot_mesh(X,F,options);
axis(<span class="string">'tight'</span>);
colormap <span class="string">parula(256)</span>;
</pre><img vspace="5" hspace="5" src="index_03.png"> <h2>Heat Diffusion and Time Stepping<a name="30"></a></h2>
         <p>The method developped in <a href="#biblio">[CraneWeischedelWardetzky13]</a> relies on the fact that the level set of the geodesic distance function to a starting point \(i\) agrees with the level set
            of the solution of the heat diffusion when the time of diffusion tends to zero. This fundamental result is proved in <a href="#biblio">[Varadhan67]</a>.
         </p>
         <p>In fact, the same result holds true when replacing the heat diffusion solution by a single Euler implicit step in time, with
            time step \(t\). This means one should consider the solution \(u\) to the equation \[ (\text{Id}+t \Delta) u = \delta_i \]
            where \(\delta_i\) is the Dirac vector at vertex index \(i\).
         </p>
         <p>Select index \(i\).</p><pre class="codeinput">i = 21000;
</pre><p>Set time \(t\).</p><pre class="codeinput">t = 1000;
</pre><p>Solve the linear system.</p><pre class="codeinput">delta = zeros(n,1);
delta(i) = 1;
u = (Ac+t*Delta)\delta;
</pre><p>Display this solution.</p><pre class="codeinput">options.face_vertex_color = u;
clf; plot_mesh(X,F,options);
axis(<span class="string">'tight'</span>);
colormap <span class="string">parula(256)</span>;
</pre><img vspace="5" hspace="5" src="index_04.png"> <p><i>Exercice 1:</i> (<a href="../missing-exo/">check the solution</a>) Solve the heat diffusion equation \[ \frac{\partial g}{\partial t} = -\Delta g \] by performing several implicit time stepping.
         </p><pre class="codeinput">exo1;
</pre><img vspace="5" hspace="5" src="index_05.png"> <h2>Geodesic in Heat Method<a name="37"></a></h2>
         <p>The main point of the method <a href="#biblio">[CraneWeischedelWardetzky13]</a> is to retrieve an approximation of the distance function \(\phi\) from the level set of implicit heat diffusion step \(u\).
         </p>
         <p>This is achieved by using the fact that \(\norm{\nabla \phi}=1\), i.e. one should have \[ \nabla \phi \approx -\frac{\nabla
            u}{\norm{\nabla u}}. \] Solving this equation in the least square sense leads to the resolution of a Poisson equation.
         </p>
         <p>Compute the solution \(u\) with explicit time stepping.</p><pre class="codeinput">t = .1;
u = (Ac+t*DeltaCot)\delta;
</pre><p>Compute the gradient field.</p><pre class="codeinput">g = Grad(u);
</pre><p>Normalize it to obtain \[ h = -\frac{\nabla u}{\norm{\nabla u}}. \]</p><pre class="codeinput">h = -normalize(g);
</pre><p>Integrate it back by solving \[ \Delta \phi = \text{div}(h).  \]</p><pre class="codeinput">phi = Delta \ Div(h);
</pre><p>Display.</p><pre class="codeinput">options.face_vertex_color = phi;
clf; plot_mesh(X,F,options);
axis(<span class="string">'tight'</span>);
colormap <span class="string">parula(256)</span>;
</pre><img vspace="5" hspace="5" src="index_06.png"> <p>Functions for displaying the level sets of the distance.</p><pre class="codeinput">p = 30;
DispFunc = @(phi)cos(2*pi*p*phi);
</pre><p>Same, but using a different colormap.</p><pre class="codeinput">options.face_vertex_color = DispFunc(phi);
clf; plot_mesh(X,F,options);
axis(<span class="string">'tight'</span>);
colormap <span class="string">parula(256)</span>;
</pre><img vspace="5" hspace="5" src="index_07.png"> <p><i>Exercice 2:</i> (<a href="../missing-exo/">check the solution</a>) Display the result obtained for several time \(t\).
         </p><pre class="codeinput">exo2;
</pre><img vspace="5" hspace="5" src="index_08.png"> <p><i>Exercice 3:</i> (<a href="../missing-exo/">check the solution</a>) Compute distances from an increasing number of starting points that are computed using a farthest point sampling.
         </p><pre class="codeinput">exo3;
</pre><img vspace="5" hspace="5" src="index_09.png"> <h2>Bibliography<a name="48"></a></h2>
         <p><a name="biblio"></a></p>
         <div>
            <ul>
               <li>[CraneWeischedelWardetzky13] K. Crane, C. Weischedel, M. Wardetzky, <i>Geodesics in heat: A new approach to computing distance based on heat flow</i>, ACM Transactions on Graphics , vol. 32, no. 5, pp. 152:1-152:11, 2013.
               </li>
               <li>[Varadhan67] S.R.S. Varadhan, <i>On the behavior of the fundamental solution of the heat equation with variable coefficients</i> Communications on Pure and Applied Mathematics 20, 2, 431-455, 1967
               </li>
            </ul>
         </div>
         <p class="footer"><br>
            Copyright  (c) 2010 Gabriel Peyre<br></p>
      </div>
      <!--
##### SOURCE BEGIN #####
%% Geodesic Distance with Poisson Equation
% This tour explores the method detailed in the paper <#biblio [CraneWeischedelWardetzky13]> to compute geodesic distance by 
% simply solving a Poisson equation.

%% Installing toolboxes and setting up the path.

%%
% You need to download the following files: 
% <../toolbox_signal.zip signal toolbox>, 
% <../toolbox_general.zip general toolbox> and 
% <../toolbox_graph.zip graph toolbox>.

%%
% You need to unzip these toolboxes in your working directory, so
% that you have 
% |toolbox_signal|, 
% |toolbox_general| and 
% |toolbox_graph|
% in your directory.

%%
% *For Scilab user:* you must replace the Matlab comment '%' by its Scilab
% counterpart '//'.

%%
% *Recommandation:* You should create a text file named for instance |numericaltour.sce| (in Scilab) or |numericaltour.m| (in Matlab) to write all the
% Scilab/Matlab command you want to execute. Then, simply run |exec('numericaltour.sce');| (in Scilab) or |numericaltour;| (in Matlab) to run the commands. 

%%
% Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

%%
% Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');
getd('toolbox_graph/');

%% Gradient, Divergence and Laplacian on Surfaces
% The topology of a triangulation is defined via a set of indexes \(\Vv = \{1,\ldots,n\}\)
% that indexes the \(n\) vertices, a set of edges \(\Ee \subset \Vv \times \Vv\)
% and a set of \(m\) faces \(\Ff \subset \Vv  \times \Vv \times \Vv\).

%%
% We load a mesh. The set of faces \(\Ff\) is stored in a matrix \(F \in
% \{1,\ldots,n\}^{3 \times m}\).
% The positions \(x_i \in \RR^3\), for \(i \in V\), of the \(n\) vertices
% are stored in a matrix \(X_0 = (x_{0,i})_{i=1}^n \in \RR^{3 \times n}\).

clear options;
name = 'elephant-50kv';
options.name = name; % useful for displaying
[X,F] = read_mesh(name);

%% 
% Number \(n\) of vertices and number \(m\) of faces.

n = size(X,2);
m = size(F,2);

%%
% Display the mesh in 3-D.

options.lighting = 1;
clf;
plot_mesh(X,F,options); 
axis('tight');

%%
% In this tour, we use piecewise linear finite element to compute the gradient
% operators, which in turns allows us to compute the divergence (its
% transposed operator) and the Laplacian 
% (as the composition of the divergence and the gradient).

%%
% The gradient operator \(\nabla\) can be understood as a collection of 3 sparse 
% matrices \((\nabla_s)_{s=1,2,3}\) of size \((m,n)\)
% that computes each coordinate of \(\nabla u=(\nabla_s u)_{s=1,2,3}\)
% through the formula, for each face \(f\), 
% \[ (\nabla u)_f = \frac{1}{2A_f} \sum_{i \in f} u_i (N_f \wedge e_i) \]
% where \(A_f\) is the area of face \(f\), \(N_f\) is the normal to the
% face, \(e_i\) is the edge opposite to vertex \(i\), and
% \(\wedge\) is the cross product in \(\RR^3\).

%%
% Callback to get the coordinates of all the vertex of index \(i=1,2,3\) in
% all faces.

XF = @(i)X(:,F(i,:));

%%
% Compute un-normalized normal through the formula 
% \(e_1 \wedge e_2 \) where \(e_i\) are the edges.

Na = cross( XF(2)-XF(1), XF(3)-XF(1) );

%%
% Compute the area of each face as half the norm of the cross product.

amplitude = @(X)sqrt( sum( X.^2 ) );
A = amplitude(Na)/2;

%%
% Compute the set of unit-norm normals to each face.

normalize = @(X)X ./ repmat(amplitude(X), [3 1]);
N = normalize(Na);

%%
% Populate the sparse entries of the matrices for the 
% operator implementing \( \sum_{i \in f} u_i (N_f \wedge e_i) \).

I = []; J = []; V = []; % indexes to build the sparse matrices
for i=1:3
    % opposite edge e_i indexes
    s = mod(i,3)+1;
    t = mod(i+1,3)+1;
    % vector N_f^e_i
    wi = cross(XF(t)-XF(s),N);    
    % update the index listing
    I = [I, 1:m];
    J = [J, F(i,:)]; 
    V = [V, wi];
end

%% 
% Sparse matrix with entries \(1/(2A_f)\).

dA = spdiags(1./(2*A(:)),0,m,m); 

%%
% Compute gradient.

GradMat = {};
for k=1:3
    GradMat{k} = dA*sparse(I,J,V(k,:),m,n);
end

%%
% \(\nabla\) gradient operator.

Grad = @(u)[GradMat{1}*u, GradMat{2}*u, GradMat{3}*u]';

%%
% Compute divergence matrices as transposed of grad for the face area inner product.

dAf = spdiags(2*A(:),0,m,m); 
DivMat = {GradMat{1}'*dAf, GradMat{2}'*dAf, GradMat{3}'*dAf};

%%
% Div operator. 

Div = @(q)DivMat{1}*q(1,:)' + DivMat{2}*q(2,:)' + DivMat{3}*q(3,:)';

%%
% Laplacian operator as the composition of grad and div. 

Delta = DivMat{1}*GradMat{1} + DivMat{2}*GradMat{2} + DivMat{3}*GradMat{3};

%%
% Cotan of an angle between two vectors.

cota = @(a,b)cot( acos( dot(normalize(a),normalize(b)) ) );

%%
% Compute cotan weights Laplacian. 

I = []; J = []; V = []; % indexes to build the sparse matrices
Ia = []; Va = []; % area of vertices
for i=1:3
    % opposite edge e_i indexes
    s = mod(i,3)+1;
    t = mod(i+1,3)+1;
    % adjacent edge
    ctheta = cota(XF(s)-XF(i), XF(t)-XF(i));
    % ctheta = max(ctheta, 1e-2); % avoid degeneracy
    % update the index listing
    I = [I, F(s,:), F(t,:)];
    J = [J, F(t,:), F(s,:)]; 
    V = [V, ctheta, ctheta];
    % update the diagonal with area of face around vertices
    Ia = [Ia, F(i,:)];
    Va = [Va, A]; 
end
% Aread diagonal matrix
Ac = sparse(Ia,Ia,Va,n,n);
% Cotan weights
Wc = sparse(I,J,V,n,n);
% Laplacian with cotan weights.
DeltaCot = spdiags(full(sum(Wc))', 0, n,n) - Wc;

%%
% Check that the Laplacian with cotan weights is actually equal to 
% the composition of divergence and gradient.

fprintf('Should be 0: %e\n', norm(Delta-DeltaCot, 'fro')/norm(Delta, 'fro'));

%%
% Display a function \(f\) on the mesh.

f = X(2,:);
options.face_vertex_color = f(:);
clf; plot_mesh(X,F,options); 
axis('tight');
colormap parula(256);

%%
% Display its Laplacian.

g = Delta*f(:);
g = clamp(g, -3*std(g), 3*std(g));
options.face_vertex_color = rescale(g);
clf; plot_mesh(X,F,options); 
axis('tight');
colormap parula(256);

%% Heat Diffusion and Time Stepping
% The method developped in <#biblio [CraneWeischedelWardetzky13]> relies on
% the fact that the level set of the geodesic distance function to a
% starting point \(i\) agrees with the level set of the solution of
% the heat diffusion when the time of diffusion tends to zero. This
% fundamental result is proved in <#biblio [Varadhan67]>.

%%
% In fact, the
% same result holds true when replacing the heat diffusion solution by a
% single Euler implicit step in time, with time step \(t\). This means one
% should consider the solution \(u\) to the equation 
% \[ (\text{Id}+t \Delta) u = \delta_i \]
% where \(\delta_i\) is the Dirac vector at vertex index \(i\).

%%
% Select index \(i\).

i = 21000; 

%%
% Set time \(t\).

t = 1000;

%%
% Solve the linear system. 

delta = zeros(n,1); 
delta(i) = 1;
u = (Ac+t*Delta)\delta;

%%
% Display this solution.

options.face_vertex_color = u;
clf; plot_mesh(X,F,options); 
axis('tight');
colormap parula(256);

%%
% _Exercice 1:_ (<../missing-exo/ check the solution>)
% Solve the heat diffusion equation 
% \[ \frac{\partial g}{\partial t} = -\Delta g \]
% by performing several implicit time stepping.

exo1;

%% Geodesic in Heat Method
% The main point of the method <#biblio [CraneWeischedelWardetzky13]> is to
% retrieve an approximation of the distance function \(\phi\) from the level set of
% implicit heat diffusion step \(u\). 

%%
% This is achieved by using the fact
% that \(\norm{\nabla \phi}=1\), i.e. one should have 
% \[ \nabla \phi \approx -\frac{\nabla u}{\norm{\nabla u}}. \]
% Solving this equation in the least square sense leads to the resolution
% of a Poisson equation. 

%%
% Compute the solution \(u\) with explicit time stepping.

t = .1;
u = (Ac+t*DeltaCot)\delta;

%%
% Compute the gradient field.

g = Grad(u); 

%%
% Normalize it to obtain 
% \[ h = -\frac{\nabla u}{\norm{\nabla u}}. \]

h = -normalize(g);

%%
% Integrate it back by solving
% \[ \Delta \phi = \text{div}(h).  \]

phi = Delta \ Div(h);

%%
% Display.

options.face_vertex_color = phi;
clf; plot_mesh(X,F,options); 
axis('tight');
colormap parula(256);

%% 
% Functions for displaying the level sets of the distance.

p = 30;
DispFunc = @(phi)cos(2*pi*p*phi);

%%
% Same, but using a different colormap.

options.face_vertex_color = DispFunc(phi);
clf; plot_mesh(X,F,options); 
axis('tight');
colormap parula(256);

%%
% _Exercice 2:_ (<../missing-exo/ check the solution>)
% Display the result obtained for several time \(t\).

exo2;

%%
% _Exercice 3:_ (<../missing-exo/ check the solution>)
% Compute distances from an increasing number of starting points
% that are computed using a farthest point sampling.

exo3;

%% Bibliography
% <html><a name="biblio"></a></html>

%%
% * [CraneWeischedelWardetzky13] K. Crane, C. Weischedel, M. Wardetzky, _Geodesics in heat: A new approach to computing distance based on heat flow_, ACM Transactions on Graphics , vol. 32, no. 5, pp. 152:1-152:11, 2013.
% * [Varadhan67] S.R.S. Varadhan, _On the behavior of the fundamental solution of the heat equation with variable coefficients_ Communications on Pure and Applied Mathematics 20, 2, 431-455, 1967

	
##### SOURCE END #####
-->
   </body>
</html>