tutorial/tutorials.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Libigl Tutorials\n",
    "[![PyPI version](https://badge.fury.io/py/libigl.svg)](https://pypi.org/project/libigl/)\n",
    "\n",
    "[![buildwheels](https://github.com/libigl/libigl-python-bindings/actions/workflows/wheels.yml/badge.svg)](https://github.com/libigl/libigl-python-bindings/actions/workflows/wheels.yml?query=branch%3Amain)\n",
    "\n",
    "![](images/libigl-logo.jpg)\n",
    "\n",
    "Libigl is an open source C++ library for geometry processing research and development. The python bindings combine the rapid prototyping familiar to Matlab to Python programmers with the performance and versatility of C++. The tutorial is a self-contained, hands-on introduction to libigl in Python.  Via interactive, step-by-step examples, we demonstrate how to accomplish common geometry processing tasks such as computation of differential quantities and operators, real-time deformation, parametrization, numerical optimization and remeshing. Each section of the lecture notes contains a simple Python example application."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 0\n",
    "\n",
    "We introduce libigl with a series of self-contained examples. The purpose of\n",
    "each example is to showcase a feature of libigl while applying to a practical\n",
    "problem in geometry processing. In this chapter, we will present the basic\n",
    "concepts of libigl.\n",
    "\n",
    "### Libigl design principles\n",
    "\n",
    "Before getting into the examples, we summarize the two main design principles in\n",
    "libigl:\n",
    "\n",
    "1. **No complex data types.** We mostly use `numpy` or `scipy` matrices and vectors. This greatly\n",
    "  favors code reusability and interoperability and forces the function authors to expose all the\n",
    "  parameters used by the algorithm.\n",
    "\n",
    "2. **Function encapsulation.** Every function is contained in a unique Python function.\n",
    "\n",
    "\n",
    "### Downloading Libigl\n",
    "Libigl can be downloaded from [PyPI](https://pypi.org/project/libigl/):\n",
    "```\n",
    "python -m pip install libigl \n",
    "```\n",
    "\n",
    "\n",
    "All of libigl functionality depends only on `numpy` and `scipy`. For the visualization in this tutorial we use [meshplot](https://github.com/skoch9/meshplot) which can be easily installed from Conda:\n",
    "```\n",
    "python -m pip install https://github.com/skoch9/meshplot/archive/0.4.0.tar.gz \n",
    "```\n",
    "\n",
    "\n",
    "To start using libigl (with the plots) you just need to import it together with the `numpy`, `scipy`, and `meshplot`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import igl\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "from meshplot import plot, subplot, interact\n",
    "\n",
    "import os\n",
    "root_folder = os.getcwd()\n",
    "#root_folder = os.path.join(os.getcwd(), \"tutorial\")\n",
    "data_folder = os.path.join(root_folder,\"/data\")\n",
    "# pip install gitpython\n",
    "from git import Repo\n",
    "if not os.path.isdir(data_folder):\n",
    "    Repo.clone_from(\"https://github.com/libigl/libigl-tutorial-data.git\", data_folder)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh representation\n",
    "\n",
    "Libigl uses `numpy` to encode vectors and matrices and `scipy` for sparse matrices.\n",
    "\n",
    "A triangular mesh is encoded as a pair of matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v: np.array\n",
    "f: np.array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`v` is a #N by 3 matrix which stores the coordinates of the vertices. Each\n",
    "row stores the coordinate of a vertex, with its x, y and z coordinates in the first,\n",
    "second and third column, respectively. The matrix `f` stores the triangle\n",
    "connectivity: each line of `f` denotes a triangle whose 3 vertices are\n",
    "represented as indices pointing to rows of `f`.\n",
    "\n",
    "![A simple mesh made of 2 triangles and 4 vertices.](images/VF.png )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.array([\n",
    "    [0., 0, 0],\n",
    "    [1, 0, 0],\n",
    "    [1, 1, 1],\n",
    "    [2, 1, 0]\n",
    "])\n",
    "\n",
    "F = np.array([\n",
    "    [0, 1, 2],\n",
    "    [1, 3, 2]\n",
    "])\n",
    "\n",
    "plot(V, F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the order of the vertex indices in `f` determines the orientation of\n",
    "the triangles and it should thus be consistent for the entire surface.\n",
    "This simple representation has many advantages:\n",
    "\n",
    "1. It is memory efficient and cache friendly\n",
    "2. The use of indices instead of pointers greatly simplifies debugging\n",
    "3. The data can be trivially copied and serialized\n",
    "\n",
    "Libigl provides input and output functions to read and write many common mesh formats.\n",
    "The IO functions are igl.read_\\* and igl.write_\\*.\n",
    "\n",
    "Reading a mesh from a file requires a single libigl function call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load a mesh in OFF format\n",
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"bunny.off\"))\n",
    "\n",
    "## Print the vertices and faces matrices \n",
    "print(\"Vertices: \", len(v))\n",
    "print(\"Faces: \", len(f))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function reads the mesh bumpy.off and returns the `v` and `f` matrices.\n",
    "Similarly, a mesh can be written to an OBJ file using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Save the mesh in OBJ format\n",
    "ret = igl.write_triangle_mesh(os.path.join(root_folder, \"data\", \"bunny_out.obj\"), v, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 1: Discrete Geometric Quantities and Operators\n",
    "This chapter illustrates a few discrete quantities that libigl can compute on a mesh and the libigl functions that construct popular discrete differential geometry operators. It also provides an introduction to basic drawing and coloring routines of our viewer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian curvature\n",
    "\n",
    "Gaussian curvature on a continuous surface is defined as the product of the\n",
    "principal curvatures:\n",
    "\n",
    " $k_G = k_1 k_2.$\n",
    "\n",
    "As an _intrinsic_ measure, it depends on the metric and\n",
    "not the surface's embedding.\n",
    "\n",
    "Intuitively, Gaussian curvature tells how locally spherical or _elliptic_ the\n",
    "surface is ( $k_G>0$ ), how locally saddle-shaped or _hyperbolic_ the surface\n",
    "is ( $k_G<0$ ), or how locally cylindrical or _parabolic_ ( $k_G=0$ ) the\n",
    "surface is.\n",
    "\n",
    "In the discrete setting, one definition for a \"discrete Gaussian curvature\"\n",
    "on a triangle mesh is via a vertex's _angular deficit_:\n",
    "\n",
    " $k_G(v_i) = 2π - \\sum\\limits_{j\\in N(i)}θ_{ij},$\n",
    "\n",
    "where $N(i)$ are the triangles incident on vertex $i$ and $θ_{ij}$ is the angle\n",
    "at vertex $i$ in triangle $j$ <cite data-cite=\"meyer2003\">(Meyer, 2003)</cite>.\n",
    "\n",
    "Just like the continuous analog, our discrete Gaussian curvature reveals\n",
    "elliptic, hyperbolic and parabolic vertices on the domain.\n",
    "\n",
    "Let's compute Gaussian curvature and visualize it in pseudocolor. First, calculate the curvature with libigl and then plot it in pseudocolors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"bumpy.off\"))\n",
    "k = igl.gaussian_curvature(v, f)\n",
    "plot(v, f, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, compute the massmatrix and divide the gaussian curvature values by area to get the integral average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)\n",
    "minv = sp.sparse.diags(1 / m.diagonal())\n",
    "kn = minv.dot(k)\n",
    "plot(v, f, kn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Curvature directions\n",
    "The two principal curvatures $(k_1,k_2)$ at a point on a surface measure how\n",
    "much the surface bends in different directions. The directions of maximum and\n",
    "minimum (signed) bending are called principal directions and are always\n",
    "orthogonal.\n",
    "\n",
    "Mean curvature is defined as the average of principal curvatures:\n",
    "\n",
    " $H = \\frac{1}{2}(k_1 + k_2).$\n",
    "\n",
    "One way to extract mean curvature is by examining the Laplace-Beltrami operator\n",
    "applied to the surface positions. The result is a so-called mean-curvature\n",
    "normal:\n",
    "\n",
    "  $-\\Delta \\mathbf{x} = H \\mathbf{n}.$\n",
    "\n",
    "It is easy to compute this on a discrete triangle mesh in libigl using the\n",
    "cotangent Laplace-Beltrami operator <cite data-cite=\"meyer2003\">(Meyer, 2003)</cite>. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = igl.cotmatrix(v, f)\n",
    "m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)\n",
    "\n",
    "minv = sp.sparse.diags(1 / m.diagonal())\n",
    "\n",
    "hn = -minv.dot(l.dot(v))\n",
    "h = np.linalg.norm(hn, axis=1)\n",
    "plot(v, f, h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combined with the angle defect definition of discrete Gaussian curvature, one\n",
    "can define principal curvatures and use least squares fitting to find\n",
    "directions  <cite data-cite=\"meyer2003\">(Meyer, 2003)</cite>.\n",
    "\n",
    "Alternatively, a robust method for determining principal curvatures is via\n",
    "quadric fitting  <cite data-cite=\"panozzo2010\">(Panozzo, 2010)</cite>. In the neighborhood around every vertex, a\n",
    "best-fit quadric is found and principal curvature values and directions are\n",
    "analytically computed on this quadric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1, v2, k1, k2 = igl.principal_curvature(v, f)\n",
    "h2 = 0.5 * (k1 + k2)\n",
    "p = plot(v, f, h2, shading={\"wireframe\": False}, return_plot=True)\n",
    "\n",
    "avg = igl.avg_edge_length(v, f) / 2.0\n",
    "p.add_lines(v + v1 * avg, v - v1 * avg, shading={\"line_color\": \"red\"})\n",
    "p.add_lines(v + v2 * avg, v - v2 * avg, shading={\"line_color\": \"green\"});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient\n",
    "Scalar functions on a surface can be discretized as a piecewise linear function\n",
    "with values defined at each mesh vertex:\n",
    "\n",
    " $f(\\mathbf{x}) \\approx \\sum\\limits_{i=1}^n \\phi_i(\\mathbf{x})\\, f_i,$\n",
    "\n",
    "where $\\phi_i$ is a piecewise linear hat function defined by the mesh so that\n",
    "for each triangle $\\phi_i$ is _the_ linear function which is one only at\n",
    "vertex $i$ and zero at the other corners.\n",
    "\n",
    "![Hat function $\\phi_i$ is one at vertex $i$, zero at all other vertices, and linear on incident triangles.](images/hat-function.jpg)\n",
    "\n",
    "Thus gradients of such piecewise linear functions are simply sums of gradients\n",
    "of the hat functions:\n",
    "\n",
    " $\\nabla f(\\mathbf{x}) \\approx\n",
    " \\nabla \\sum\\limits_{i=1}^n \\phi_i(\\mathbf{x})\\, f_i =\n",
    " \\sum\\limits_{i=1}^n \\nabla \\phi_i(\\mathbf{x})\\, f_i.$\n",
    "\n",
    "This reveals that the gradient is a linear function of the vector of $f_i$\n",
    "values. Because the $\\phi_i$ are linear in each triangle, their gradients are\n",
    "_constant_ in each triangle. Thus our discrete gradient operator can be written\n",
    "as a matrix multiplication taking vertex values to triangle values:\n",
    "\n",
    " $\\nabla f \\approx \\mathbf{G}\\,\\mathbf{f},$\n",
    "\n",
    "where $\\mathbf{f}$ is $n\\times 1$ and $\\mathbf{G}$ is an $md\\times n$ sparse\n",
    "matrix. This matrix $\\mathbf{G}$ can be derived geometrically <cite data-cite=\"jacobson2013\">(Jacobson, 2013)</cite>.\n",
    "\n",
    "Libigl's `grad` function computes $\\mathbf{G}$ for\n",
    "triangle and tetrahedral meshes. \n",
    "Let's see how this works. First load a mesh and a corresponding surface function.\n",
    "Next, compute the gradient operator g (#F*3 x #V) on the triangle mesh, apply it to the surface function and extract the magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"cheburashka.off\"))\n",
    "u = igl.read_dmat(os.path.join(root_folder, \"data\", \"cheburashka-scalar.dmat\"))\n",
    "\n",
    "g = igl.grad(v, f)\n",
    "gu = g.dot(u).reshape(f.shape, order=\"F\")\n",
    "\n",
    "gu_mag = np.linalg.norm(gu, axis=1)\n",
    "p = plot(v, f, u, shading={\"wireframe\":False}, return_plot=True)\n",
    "\n",
    "max_size = igl.avg_edge_length(v, f) / np.mean(gu_mag)\n",
    "bc = igl.barycenter(v, f)\n",
    "bcn = bc + max_size * gu\n",
    "p.add_lines(bc, bcn, shading={\"line_color\": \"black\"});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laplacian\n",
    "\n",
    "The discrete Laplacian is an essential geometry processing tool. Many\n",
    "interpretations and flavors of the Laplace and Laplace-Beltrami operator exist.\n",
    "\n",
    "In open Euclidean space, the _Laplace_ operator is the usual divergence of\n",
    "gradient (or equivalently the Laplacian of a function is the trace of its\n",
    "Hessian):\n",
    "\n",
    " $\\Delta f =\n",
    " \\frac{\\partial^2 f}{\\partial x^2} +\n",
    " \\frac{\\partial^2 f}{\\partial y^2} +\n",
    " \\frac{\\partial^2 f}{\\partial z^2}.$\n",
    "\n",
    "The _Laplace-Beltrami_ operator generalizes this to surfaces.\n",
    "\n",
    "When considering piecewise-linear functions on a triangle mesh, a discrete\n",
    "Laplacian may be derived in a variety of ways. The most popular in geometry\n",
    "processing is the so-called \"cotangent Laplacian\" $\\mathbf{L}$, arising\n",
    "simultaneously from FEM, DEC and applying divergence theorem to vertex\n",
    "one-rings. As a linear operator taking vertex values to vertex values, the\n",
    "Laplacian $\\mathbf{L}$ is a $n\\times n$ matrix with elements:\n",
    "\n",
    "$L_{ij} = \\begin{cases}j \\in N(i) &\\cot \\alpha_{ij} + \\cot \\beta_{ij},\\\\\n",
    "j \\notin N(i) & 0,\\\\\n",
    "i = j & -\\sum\\limits_{k\\neq i} L_{ik},\n",
    "\\end{cases}$\n",
    "\n",
    "where $N(i)$ are the vertices adjacent to (neighboring) vertex $i$, and\n",
    "$\\alpha_{ij},\\beta_{ij}$ are the angles opposite to edge ${ij}$.\n",
    "\n",
    "Libigl implements discrete \"cotangent Laplacians\" for triangles meshes and\n",
    "tetrahedral meshes, building both with fast geometric rules rather than \"by the\n",
    "book\" FEM construction which involves many (small) matrix inversions <cite data-cite=\"sharf_2007\">(Sharf, 2007)</cite>.\n",
    "\n",
    "First, load a triangle mesh and then calculate the Laplace-Beltrami operator, visualize the normals as pseudocolors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse.linalg import spsolve\n",
    "\n",
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"cow.off\"))\n",
    "l = igl.cotmatrix(v, f)\n",
    "\n",
    "n = igl.per_vertex_normals(v, f)*0.5+0.5\n",
    "c = np.linalg.norm(n, axis=1)\n",
    "p = plot(v, f, c, shading={\"wireframe\": False}, return_plot=True)\n",
    "\n",
    "vs = [v]\n",
    "cs = [c]\n",
    "for i in range(10):\n",
    "    m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_BARYCENTRIC)\n",
    "    s = (m - 0.001 * l)\n",
    "    b = m.dot(v)\n",
    "    v = spsolve(s, m.dot(v))\n",
    "    n = igl.per_vertex_normals(v, f)*0.5+0.5\n",
    "    c = np.linalg.norm(n, axis=1)\n",
    "    vs.append(v)\n",
    "    cs.append(c)\n",
    "\n",
    "@interact(level=(0, 9))\n",
    "def mcf(level=0):\n",
    "    p.update_object(vertices=vs[level], colors=cs[level])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The operator applied to mesh vertex positions amounts to smoothing by _flowing_\n",
    "the surface along the mean curvature normal direction. Note that this is equivalent to minimizing surface area. The following example computes conformalized mean curvature flow using the cotangent Laplacian <cite data-cite=\"kazhdan_2012\">(Kazhdan, 2012)</cite> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mass matrix\n",
    "The mass matrix $\\mathbf{M}$ is another $n \\times n$ matrix which takes vertex\n",
    "values to vertex values. From an FEM point of view, it is a discretization of\n",
    "the inner-product: it accounts for the area around each vertex. Consequently,\n",
    "$\\mathbf{M}$ is often a diagonal matrix, such that $M_{ii}$ is the barycentric\n",
    "or voronoi area around vertex $i$ in the mesh <cite data-cite=\"meyer_2003\">(Meyer, 2003)</cite>. The inverse of this matrix is also very useful as it transforms integrated quantities into point-wise quantities, e.g.:\n",
    "\n",
    " $\\Delta f \\approx \\mathbf{M}^{-1} \\mathbf{L} \\mathbf{f}.$\n",
    "\n",
    "In general, when encountering squared quantities integrated over the surface,\n",
    "the mass matrix will be used as the discretization of the inner product when\n",
    "sampling function values at vertices:\n",
    "\n",
    " $\\int_S x\\, y\\ dA \\approx \\mathbf{x}^T\\mathbf{M}\\,\\mathbf{y}.$\n",
    "\n",
    "An alternative mass matrix $\\mathbf{T}$ is a $md \\times md$ matrix which takes\n",
    "triangle vector values to triangle vector values. This matrix represents an\n",
    "inner-product accounting for the area associated with each triangle (i.e. the\n",
    "triangles true area).\n",
    "\n",
    "### Alternative construction of Laplacian\n",
    "\n",
    "An alternative construction of the discrete cotangent Laplacian is by\n",
    "\"squaring\" the discrete gradient operator. This may be derived by applying\n",
    "Green's identity (ignoring boundary conditions for the moment):\n",
    "\n",
    "  $\\int_S \\|\\nabla f\\|^2 dA = \\int_S f \\Delta f dA$\n",
    "\n",
    "Or in matrix form which is immediately translatable to code:\n",
    "\n",
    "  $\\mathbf{f}^T \\mathbf{G}^T \\mathbf{T} \\mathbf{G} \\mathbf{f} =\n",
    "  \\mathbf{f}^T \\mathbf{M} \\mathbf{M}^{-1} \\mathbf{L} \\mathbf{f} =\n",
    "  \\mathbf{f}^T \\mathbf{L} \\mathbf{f}.$\n",
    "\n",
    "So we have that $\\mathbf{L} = \\mathbf{G}^T \\mathbf{T} \\mathbf{G}$. This also\n",
    "hints that we may consider $\\mathbf{G}^T$ as a discrete _divergence_ operator,\n",
    "since the Laplacian is the divergence of the gradient. Naturally, $\\mathbf{G}^T$ is\n",
    "a $n \\times md$ sparse matrix which takes vector values stored at triangle faces\n",
    "to scalar divergence values at vertices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"cow.off\"))\n",
    "l = igl.cotmatrix(v, f)\n",
    "g = igl.grad(v, f)\n",
    "\n",
    "d_area = igl.doublearea(v, f)\n",
    "t = sp.sparse.diags(np.hstack([d_area, d_area, d_area]) * 0.5)\n",
    "\n",
    "k = -g.T.dot(t).dot(g)\n",
    "print(\"|k-l|: %s\"%sp.sparse.linalg.norm(k-l))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exact Discrete Geodesic Distances\n",
    "\n",
    "The discrete geodesic distance between two points is the length of the shortest\n",
    "path between then restricted to the surface. For triangle meshes, such a path is\n",
    "made of a set of segments which can be either edges of the mesh or crossing a\n",
    "triangle.\n",
    "\n",
    "Libigl includes a wrapper for the exact geodesic algorithm <cite data-cite=\"mitchell_1987\">(Mitchell, 1987)</cite>\n",
    "developed by Danil Kirsanov (https://code.google.com/archive/p/geodesic/),\n",
    "exposing it through an Eigen-based API. The function \n",
    "```python\n",
    "d = igl.exact_geodesic(v, f, vs, fs, vt, ft)\n",
    "```\n",
    "computes the closest geodesic distances of each vertex in vt or face in ft, from\n",
    "the source vertices vs or faces fs of the input mesh v, f. The output is written\n",
    "in the vector d, which lists first the distances for the vertices in vt, and\n",
    "then for the faces in ft. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"armadillo.obj\"))\n",
    "\n",
    "## Select a vertex from which the distances should be calculated\n",
    "vs = np.array([0])\n",
    "##All vertices are the targets\n",
    "vt = np.arange(v.shape[0])\n",
    "\n",
    "d = igl.exact_geodesic(v, f, vs, vt)#, fs, ft)\n",
    "\n",
    "strip_size = 0.1\n",
    "##The function should be 1 on each integer coordinate\n",
    "c = np.abs(np.sin((d / strip_size * np.pi)))\n",
    "plot(v, f, c, shading={\"wireframe\": False})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 2: Matrices and Linear Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laplace equation\n",
    "A common linear system in geometry processing is the Laplace equation:\n",
    "\n",
    " $∆z = 0$\n",
    "\n",
    "subject to some boundary conditions, for example Dirichlet boundary conditions\n",
    "(fixed value):\n",
    "\n",
    " $\\left.z\\right|_{\\partial{S}} = z_{bc}$\n",
    "\n",
    "In the discrete setting, the linear system can be written as:\n",
    "\n",
    " $\\mathbf{L} \\mathbf{z} = \\mathbf{0}$\n",
    "\n",
    "where $\\mathbf{L}$ is the $n \\times n$ discrete Laplacian and $\\mathbf{z}$ is a\n",
    "vector of per-vertex values. Most of $\\mathbf{z}$ correspond to interior\n",
    "vertices and are unknown, but some of $\\mathbf{z}$ represent values at boundary\n",
    "vertices. Their values are known so we may move their corresponding terms to\n",
    "the right-hand side.\n",
    "\n",
    "Conceptually, this is very easy if we have sorted $\\mathbf{z}$ so that interior\n",
    "vertices come first and then boundary vertices:\n",
    "\n",
    "$$\n",
    " \\left(\\begin{array}{cc}\n",
    " \\mathbf{L}_{in,in} & \\mathbf{L}_{in,b}\\\\\n",
    " \\mathbf{L}_{b,in} & \\mathbf{L}_{b,b}\\end{array}\\right)\n",
    " \\left(\\begin{array}{c}\n",
    " \\mathbf{z}_{in}\\\\\n",
    " \\mathbf{z}_{b}\\end{array}\\right) =\n",
    " \\left(\\begin{array}{c}\n",
    " \\mathbf{0}_{in}\\\\\n",
    " \\mathbf{z}_{bc}\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "The bottom block of equations is no longer meaningful so we'll only consider\n",
    "the top block:\n",
    "\n",
    "$$\n",
    " \\left(\\begin{array}{cc}\n",
    " \\mathbf{L}_{in,in} & \\mathbf{L}_{in,b}\\end{array}\\right)\n",
    " \\left(\\begin{array}{c}\n",
    " \\mathbf{z}_{in}\\\\\n",
    " \\mathbf{z}_{b}\\end{array}\\right) =\n",
    " \\mathbf{0}_{in}\n",
    "$$\n",
    "\n",
    "We can move the known values to the right-hand side:\n",
    "\n",
    "$$\n",
    " \\mathbf{L}_{in,in}\n",
    " \\mathbf{z}_{in} = -\n",
    " \\mathbf{L}_{in,b}\n",
    " \\mathbf{z}_{b}\n",
    "$$\n",
    "\n",
    "Finally we can solve this equation for the unknown values at interior vertices\n",
    "$\\mathbf{z}_{in}$.\n",
    "\n",
    "However, our vertices will often not be sorted in this way. One option would be to sort `V`,\n",
    "then proceed as above and then _unsort_ the solution `Z` to match `V`. However,\n",
    "this solution is not very general.\n",
    "\n",
    "With array slicing no explicit sort is needed. Instead we can _slice-out_\n",
    "submatrix blocks ($\\mathbf{L}_{in,in}$, $\\mathbf{L}_{in,b}$, etc.) and follow\n",
    "the linear algebra above directly. Then we can slice the solution _into_ the\n",
    "rows of `Z` corresponding to the interior vertices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse.linalg import spsolve\n",
    "\n",
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"camelhead.off\"))\n",
    "\n",
    "## Find boundary vertices\n",
    "e = igl.boundary_facets(f)\n",
    "v_b = np.unique(e)\n",
    "\n",
    "## List of all vertex indices\n",
    "v_all = np.arange(v.shape[0])\n",
    "\n",
    "## List of interior indices\n",
    "v_in = np.setdiff1d(v_all, v_b)\n",
    "\n",
    "## Construct and slice up Laplacian\n",
    "l = igl.cotmatrix(v, f)\n",
    "l_ii = l[v_in, :]\n",
    "l_ii = l_ii[:, v_in]\n",
    "\n",
    "l_ib = l[v_in, :]\n",
    "l_ib = l_ib[:, v_b]\n",
    "\n",
    "## Dirichlet boundary conditions from z-coordinate\n",
    "z = v[:, 2]\n",
    "bc = z[v_b]\n",
    "\n",
    "## Solve PDE\n",
    "z_in = spsolve(-l_ii, l_ib.dot(bc))\n",
    "\n",
    "plot(v, f, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quadratic energy minimization\n",
    "\n",
    "The same Laplace equation may be equivalently derived by minimizing Dirichlet\n",
    "energy subject to the same boundary conditions:\n",
    "\n",
    " $\\mathop{\\text{minimize }}_z \\frac{1}{2}\\int\\limits_S \\|\\nabla z\\|^2 dA$\n",
    "\n",
    "On our discrete mesh, recall that this becomes\n",
    "\n",
    " $\\mathop{\\text{minimize }}_\\mathbf{z}  \\frac{1}{2}\\mathbf{z}^T \\mathbf{G}^T \\mathbf{D}\n",
    " \\mathbf{G} \\mathbf{z} \\rightarrow \\mathop{\\text{minimize }}_\\mathbf{z} \\mathbf{z}^T \\mathbf{L} \\mathbf{z}$\n",
    "\n",
    "The general problem of minimizing some energy over a mesh subject to fixed\n",
    "value boundary conditions is so wide spread that libigl has a dedicated api for\n",
    "solving such systems.\n",
    "\n",
    "Let us consider a general quadratic minimization problem subject to different\n",
    "common constraints:\n",
    "\n",
    "$$\n",
    " \\mathop{\\text{minimize }}_\\mathbf{z}  \\frac{1}{2}\\mathbf{z}^T \\mathbf{Q} \\mathbf{z} +\n",
    " \\mathbf{z}^T \\mathbf{B} + \\text{constant},\n",
    "$$\n",
    "\n",
    " subject to\n",
    "\n",
    "$$\n",
    " \\mathbf{z}_b = \\mathbf{z}_{bc} \\text{ and } \\mathbf{A}_{eq} \\mathbf{z} =\n",
    " \\mathbf{B}_{eq},\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "  - $\\mathbf{Q}$ is a (usually sparse) $n \\times n$ positive semi-definite\n",
    "    matrix of quadratic coefficients (Hessian),\n",
    "  - $\\mathbf{B}$ is a $n \\times 1$ vector of linear coefficients,\n",
    "  - $\\mathbf{z}_b$ is a $|b| \\times 1$ portion of\n",
    "$\\mathbf{z}$ corresponding to boundary or _fixed_ vertices,\n",
    "  - $\\mathbf{z}_{bc}$ is a $|b| \\times 1$ vector of known values corresponding to\n",
    "    $\\mathbf{z}_b$,\n",
    "  - $\\mathbf{A}_{eq}$ is a (usually sparse) $m \\times n$ matrix of linear\n",
    "    equality constraint coefficients (one row per constraint), and\n",
    "  - $\\mathbf{B}_{eq}$ is a $m \\times 1$ vector of linear equality constraint\n",
    "    right-hand side values.\n",
    "\n",
    "This specification is overly general as we could write $\\mathbf{z}_b =\n",
    "\\mathbf{z}_{bc}$ as rows of $\\mathbf{A}_{eq} \\mathbf{z} =\n",
    "\\mathbf{B}_{eq}$, but these fixed value constraints appear so often that they\n",
    "merit a dedicated function: `min_quad_with_fixed`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear equality constraints\n",
    "We saw above that `min_quad_with_fixed` in libigl provides a compact way to\n",
    "solve general quadratic programs. Let's consider another example, this time\n",
    "with active linear equality constraints. Specifically let's solve the\n",
    "`bi-Laplace equation` or equivalently minimize the Laplace energy:\n",
    "\n",
    "$$\n",
    " \\Delta^2 z = 0 \\leftrightarrow \\mathop{\\text{minimize }}\\limits_z \\frac{1}{2}\n",
    " \\int\\limits_S (\\Delta z)^2 dA\n",
    "$$\n",
    "\n",
    "subject to fixed value constraints and a linear equality constraint:\n",
    "\n",
    " $z_{a} = 1, z_{b} = -1$ and $z_{c} = z_{d}$.\n",
    "\n",
    "Notice that we can rewrite the last constraint in the familiar form from above:\n",
    "\n",
    " $z_{c} - z_{d} = 0.$\n",
    "\n",
    "Now we can assembly `Aeq` as a $1 \\times n$ sparse matrix with a coefficient\n",
    "$1$ in the column corresponding to vertex $c$ and a $-1$ at $d$. The right-hand\n",
    "side `Beq` is simply zero.\n",
    "\n",
    "Internally, `min_quad_with_fixed` solves using the Lagrange Multiplier\n",
    "method. This method adds additional variables for each linear constraint (in\n",
    "general a $m \\times 1$ vector of variables $\\lambda$) and then solves the\n",
    "saddle problem:\n",
    "\n",
    "$$\n",
    "  \\mathop{\\text{find saddle }}_{\\mathbf{z},\\lambda}\\, \\frac{1}{2}\\mathbf{z}^T \\mathbf{Q} \\mathbf{z} +\n",
    "  \\mathbf{z}^T \\mathbf{B} + \\text{constant} + \\lambda^T\\left(\\mathbf{A}_{eq}\n",
    " \\mathbf{z} - \\mathbf{B}_{eq}\\right)\n",
    "$$\n",
    "\n",
    "This can be rewritten in a more familiar form by stacking $\\mathbf{z}$ and\n",
    "$\\lambda$ into one $(m+n) \\times 1$ vector of unknowns:\n",
    "\n",
    "$$\n",
    " \\mathop{\\text{find saddle }}_{\\mathbf{z},\\lambda}\\,\n",
    " \\frac{1}{2}\n",
    " \\left(\n",
    "  \\mathbf{z}^T\n",
    "  \\lambda^T\n",
    " \\right)\n",
    " \\left(\n",
    "  \\begin{array}{cc}\n",
    "  \\mathbf{Q}      & \\mathbf{A}_{eq}^T\\\\\n",
    "  \\mathbf{A}_{eq} & 0\n",
    "  \\end{array}\n",
    " \\right)\n",
    " \\left(\n",
    "  \\begin{array}{c}\n",
    "  \\mathbf{z}\\\\\n",
    "  \\lambda\n",
    "  \\end{array}\n",
    " \\right) +\n",
    " \\left(\n",
    "  \\mathbf{z}^T\n",
    "  \\lambda^T\n",
    " \\right)\n",
    " \\left(\n",
    "  \\begin{array}{c}\n",
    "  \\mathbf{B}\\\\\n",
    "  -\\mathbf{B}_{eq}\n",
    "  \\end{array}\n",
    "  \\right)\n",
    "  + \\text{constant}\n",
    "$$\n",
    "\n",
    "Differentiating with respect to $\\left( \\mathbf{z}^T \\lambda^T \\right)$ reveals\n",
    "a linear system and we can solve for $\\mathbf{z}$ and $\\lambda$. The only\n",
    "difference from the straight quadratic _minimization_ system, is that this\n",
    "saddle problem system will not be positive definite. Thus, we must use a\n",
    "different factorization technique (LDLT rather than LLT): libigl's\n",
    "`min_quad_with_fixed` automatically chooses the correct solver in\n",
    "the presence of linear equality constraints.\n",
    "\n",
    "The following example first solves with just fixed value constraints (left: 1 and -1 on the left hand and foot respectively), then solves with an additional linear equality constraint (right: points on right hand and foot constrained to be equal).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"cheburashka.off\"))\n",
    "\n",
    "## Two fixed points: Left hand, left foot should have values 1 and -1\n",
    "b = np.array([4331, 5957])\n",
    "bc = np.array([1., -1.])\n",
    "B = np.zeros((v.shape[0], 1))\n",
    "\n",
    "## Construct Laplacian and mass matrix\n",
    "L = igl.cotmatrix(v, f)\n",
    "M = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)\n",
    "Minv = sp.sparse.diags(1 / M.diagonal())\n",
    "\n",
    "## Bi-Laplacian\n",
    "Q = L @ (Minv @ L)\n",
    "\n",
    "## Solve with only equality constraints\n",
    "Aeq = sp.sparse.csc_matrix((0, 0))\n",
    "Beq = np.array([])\n",
    "_, z1 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)\n",
    "\n",
    "## Solve with equality and linear constraints\n",
    "Aeq = sp.sparse.csc_matrix((1, v.shape[0]))\n",
    "Aeq[0,6074] = 1\n",
    "Aeq[0, 6523] = -1\n",
    "Beq = np.array([0.])\n",
    "_, z2 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)\n",
    "\n",
    "## Normalize colors to same range\n",
    "min_z = min(np.min(z1), np.min(z2))\n",
    "max_z = max(np.max(z1), np.max(z2))\n",
    "z = [(z1 - min_z) / (max_z - min_z), (z2 - min_z) / (max_z - min_z)]\n",
    "\n",
    "## Plot the functions\n",
    "p = plot(v, f, z1, shading={\"wireframe\":False}, return_plot=True)\n",
    "\n",
    "@interact(function=[('z0', 0), ('z1', 1)])\n",
    "def sf(function):\n",
    "    p.update_object(colors=z[function])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quadratic programming\n",
    "\n",
    "We can generalize the quadratic optimization in the previous section even more\n",
    "by allowing inequality constraints. Specifically box constraints (lower and\n",
    "upper bounds):\n",
    "\n",
    " $\\mathbf{l} \\le \\mathbf{z} \\le \\mathbf{u},$\n",
    "\n",
    "where $\\mathbf{l},\\mathbf{u}$ are $n \\times 1$ vectors of lower and upper\n",
    "bounds\n",
    "and general linear inequality constraints:\n",
    "\n",
    " $\\mathbf{A}_{ieq} \\mathbf{z} \\le \\mathbf{B}_{ieq},$\n",
    "\n",
    "where $\\mathbf{A}_{ieq}$ is a $k \\times n$ matrix of linear coefficients and\n",
    "$\\mathbf{B}_{ieq}$ is a $k \\times 1$ matrix of constraint right-hand sides.\n",
    "\n",
    "Again, we are overly general as the box constraints could be written as\n",
    "rows of the linear inequality constraints, but bounds appear frequently enough\n",
    "to merit a dedicated api.\n",
    "\n",
    "Libigl implements its own active set routine for solving _quadratric programs_\n",
    "(QPs). This algorithm works by iteratively \"activating\" violated inequality\n",
    "constraints by enforcing them as equalities and \"deactivating\" constraints\n",
    "which are no longer needed.\n",
    "\n",
    "After deciding which constraints are active at each iteration, the problem\n",
    "reduces to a quadratic minimization subject to linear _equality_ constraints,\n",
    "and the method from the previous section is invoked. This is repeated until convergence.\n",
    "\n",
    "Currently the implementation is efficient for box constraints and sparse\n",
    "non-overlapping linear inequality constraints.\n",
    "\n",
    "Unlike alternative interior-point methods, the active set method benefits from\n",
    "a warm-start (initial guess for the solution vector $\\mathbf{z}$).\n",
    "\n",
    "The following example uses an active set solver to optimize discrete biharmonic kernels <cite data-cite=\"rustamov_2011\">(Rustamov, 2011)</cite> at multiple scales:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO: Check why results differ, add interactivity\n",
    "\n",
    "v, f, _ = igl.read_off(os.path.join(root_folder, \"data\", \"cheburashka.off\"))\n",
    "\n",
    "# One fixed point on belly\n",
    "b = np.array([[2556]])\n",
    "bc = np.array([[1.0]])\n",
    "\n",
    "# Construct Laplacian and mass matrix\n",
    "l = igl.cotmatrix(v, f)\n",
    "m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)\n",
    "minv = sp.sparse.diags(1 / m.diagonal())\n",
    "\n",
    "# Bi-Laplacian\n",
    "q = l @ (minv @ l)\n",
    "\n",
    "# Zero linear term\n",
    "bz = np.zeros((v.shape[0], 1))\n",
    "\n",
    "# Lower and upper bound\n",
    "lx = np.zeros((v.shape[0], 1))\n",
    "ux = np.ones((v.shape[0], 1))\n",
    "\n",
    "# Equality constraint constrains solution to sum to 1\n",
    "beq = np.array([[0.08]])\n",
    "aeq = sp.sparse.csc_matrix(m.diagonal())\n",
    "\n",
    "# Empty inequality constraints\n",
    "aieq = sp.sparse.csc_matrix((0, 0))\n",
    "bieq = np.array([])\n",
    "\n",
    "z = igl.active_set(q, bz, b, bc, aeq, beq, aieq, bieq, lx, ux, max_iter=8)\n",
    "plot(v, f, z[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigen Decomposition\n",
    "\n",
    "Libigl has rudimentary support for extracting eigen pairs of a generalized\n",
    "eigen value problem:\n",
    "\n",
    " $Ax = \\lambda B x$\n",
    "\n",
    "where $A$ is a sparse symmetric matrix and $B$ is a sparse positive definite\n",
    "matrix. Most commonly in geometry processing, we let $A=L$ the cotangent\n",
    "Laplacian and $B=M$ the per-vertex mass matrix <cite data-cite=\"vallet_2008\">(Vallet, 2008)</cite>.\n",
    "Typically applications will make use of the _low frequency_ eigen modes.\n",
    "Analogous to the Fourier decomposition, a function $f$ on a surface can be\n",
    "represented via its spectral decomposition of the eigen modes of the\n",
    "Laplace-Beltrami:\n",
    "\n",
    " $f = \\sum\\limits_{i=1}^\\infty a_i \\phi_i$\n",
    "\n",
    "where each $\\phi_i$ is an eigen function satisfying: $\\Delta \\phi_i = \\lambda_i\n",
    "\\phi_i$ and $a_i$ are scalar coefficients. For a discrete triangle mesh, a\n",
    "completely analogous decomposition exists, albeit with finite sum:\n",
    "\n",
    " $\\mathbf{f} = \\sum\\limits_{i=1}^n a_i \\phi_i$\n",
    "\n",
    "where now a column vector of values at vertices $\\mathbf{f} \\in \\mathcal{R}^n$\n",
    "specifies a piecewise linear function and $\\phi_i \\in \\mathcal{R}^n$ is an\n",
    "eigen vector satisfying:\n",
    "\n",
    "$\\mathbf{L} \\phi_i = \\lambda_i \\mathbf{M} \\phi_i$.\n",
    "\n",
    "Note that Vallet &amp; Levy <cite data-cite=\"vallet_2008\">(Vallet, 2008)</cite> propose solving a symmetrized\n",
    "_standard_ eigen problem $\\mathbf{M}^{-1/2}\\mathbf{L}\\mathbf{M}^{-1/2} \\phi_i\n",
    "= \\lambda_i \\phi_i$. Libigl implements a generalized eigen problem solver so\n",
    "this unnecessary symmetrization can be avoided.\n",
    "\n",
    "Often the sum above is _truncated_ to the first $k$ eigen vectors. If the low\n",
    "frequency modes are chosen, i.e. those corresponding to small $\\lambda_i$\n",
    "values, then this truncation effectively _regularizes_ $\\mathbf{f}$ to smooth,\n",
    "slowly changing functions over the mesh <cite data-cite=\"hildebrandt_2011\">(Hildebrandt, 2011)</cite>. Modal\n",
    "analysis and model subspaces have been used frequently in real-time deformation\n",
    "<cite data-cite=\"barbic_2012\">(Barbic, 2005)</cite>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following example, the first k eigen vectors of the discrete Laplace-Beltrami operator are computed and displayed in\n",
    "pseudocolors atop the beetle. \n",
    "Low frequency eigen vectors of the discrete Laplace-Beltrami operator vary smoothly and slowly over the model.\n",
    "At first, calculate the Laplace-Betrami operator and solve the generalized Eigen problem with scipy/arpack. \n",
    "Then, rescale the Eigen vectors and visualize them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"beetle.off\"))\n",
    "l = -igl.cotmatrix(v, f)\n",
    "m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)\n",
    "\n",
    "k = 10\n",
    "d, u = sp.sparse.linalg.eigsh(l, k, m, sigma=0, which=\"LM\")\n",
    "\n",
    "u = (u - np.min(u)) / (np.max(u) - np.min(u))\n",
    "bbd = 0.5 * np.linalg.norm(np.max(v, axis=0) - np.min(v, axis=0))\n",
    "\n",
    "p = plot(v, f, bbd * u[:, 0], shading={\"wireframe\":False, \"flat\": False}, return_plot=True)\n",
    "\n",
    "@interact(ev=[(\"EV %i\"%i, i) for i in range(k)])\n",
    "def sf(ev):\n",
    "    p.update_object(colors=u[:, ev])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 3: Shape deformation\n",
    "Modern mesh-based shape deformation methods satisfy user deformation\n",
    "constraints at handles (selected vertices or regions on the mesh) and propagate\n",
    "these handle deformations to the rest of the shape _smoothly_ and _without removing\n",
    "or distorting details_. Libigl provides implementations of a variety of\n",
    "state-of-the-art deformation techniques, ranging from quadratic mesh-based\n",
    "energy minimizers, to skinning methods, to non-linear elasticity-inspired\n",
    "techniques."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Biharmonic deformation\n",
    "The period of research between 2000 and 2010 produced a collection of\n",
    "techniques that cast the problem of handle-based shape deformation as a\n",
    "quadratic energy minimization problem or equivalently the solution to a linear\n",
    "partial differential equation.\n",
    "\n",
    "There are many flavors of these techniques, but a prototypical subset are those\n",
    "that consider solutions to the bi-Laplace equation, that is a biharmonic\n",
    "function <cite data-cite=\"botsch_2004\">(Botsch, 2004)</cite>. This fourth-order PDE provides sufficient\n",
    "flexibility in boundary conditions to ensure $C^1$ continuity at handle\n",
    "constraints in the limit under refinement <cite data-cite=\"jacobson_mixed_2010\">(Jacobson, 2010)</cite>.\n",
    "\n",
    "#### Biharmonic surfaces\n",
    "Let us first begin our discussion of biharmonic _deformation_, by considering\n",
    "biharmonic _surfaces_. We will casually define biharmonic surfaces as surface\n",
    "whose _position functions_ are biharmonic with respect to some initial\n",
    "parameterization:\n",
    "\n",
    " $\\Delta^2 \\mathbf{x}' = 0$\n",
    "\n",
    "and subject to some handle constraints, conceptualized as \"boundary\n",
    "conditions\":\n",
    "\n",
    " $\\mathbf{x}'_{b} = \\mathbf{x}_{bc}.$\n",
    "\n",
    "where $\\mathbf{x}'$ is the unknown 3D position of a point on the surface. So we\n",
    "are asking that the bi-Laplacian of each of spatial coordinate function to be\n",
    "zero.\n",
    "\n",
    "In libigl, one can solve a biharmonic problem with `harmonic`\n",
    "and setting $k=2$ (_bi_-harmonic).\n",
    "\n",
    "This produces a smooth surface that interpolates the handle constraints, but all\n",
    "original details on the surface will be _smoothed away_. Most obviously, if the\n",
    "original surface is not already biharmonic, then giving all handles the\n",
    "identity deformation (keeping them at their rest positions) will **not**\n",
    "reproduce the original surface. Rather, the result will be the biharmonic\n",
    "surface that does interpolate those handle positions.\n",
    "\n",
    "Thus, we may conclude that this is not an intuitive technique for shape\n",
    "deformation.\n",
    "\n",
    "#### Biharmonic deformation fields\n",
    "Now we know that one useful property for a deformation technique is \"rest pose\n",
    "reproduction\": applying no deformation to the handles should apply no\n",
    "deformation to the shape.\n",
    "\n",
    "To guarantee this by construction we can work with _deformation fields_ (ie.\n",
    "displacements)\n",
    "$\\mathbf{d}$ rather\n",
    "than directly with positions $\\mathbf{x}$. Then the deformed positions can be\n",
    "recovered as\n",
    "\n",
    " $\\mathbf{x}' = \\mathbf{x}+\\mathbf{d}.$\n",
    "\n",
    "A smooth deformation field $\\mathbf{d}$ which interpolates the deformation\n",
    "fields of the handle constraints will impose a smooth deformed shape\n",
    "$\\mathbf{x}'$. Naturally, we consider _biharmonic deformation fields_:\n",
    "\n",
    " $\\Delta^2 \\mathbf{d} = 0$\n",
    "\n",
    "subject to the same handle constraints, but rewritten in terms of their implied\n",
    "deformation field at the boundary (handles).\n",
    "\n",
    " $\\mathbf{d}_b = \\mathbf{x}_{bc} - \\mathbf{x}_b.$\n",
    "\n",
    "Again we can use `harmonic` with $k=2$, but this time solve for the\n",
    "deformation field and then recover the deformed positions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"decimated-max.obj\"))\n",
    "v[:,[0, 2]] = v[:,[2, 0]] # Swap X and Z axes\n",
    "u = v.copy()\n",
    "\n",
    "s = igl.read_dmat(os.path.join(root_folder, \"data\", \"decimated-max-selection.dmat\"))\n",
    "b = np.array([[t[0] for t in [(i, s[i]) for i in range(0, v.shape[0])] if t[1] >= 0]]).T\n",
    "\n",
    "## Boundary conditions directly on deformed positions\n",
    "u_bc = np.zeros((b.shape[0], v.shape[1]))\n",
    "v_bc = np.zeros((b.shape[0], v.shape[1]))\n",
    "\n",
    "for bi in range(b.shape[0]):\n",
    "    v_bc[bi] = v[b[bi]]\n",
    "\n",
    "    if s[b[bi]] == 0: # Don't move handle 0\n",
    "        u_bc[bi] = v[b[bi]]\n",
    "    elif s[b[bi]] == 1: # Move handle 1 down\n",
    "        u_bc[bi] = v[b[bi]] + np.array([[0, -50, 0]])\n",
    "    else: # Move other handles forward\n",
    "        u_bc[bi] = v[b[bi]] + np.array([[-25, 0, 0]])\n",
    "\n",
    "p = plot(v, f, s, shading={\"wireframe\": False, \"colormap\": \"tab10\"}, return_plot=True)\n",
    "\n",
    "@interact(deformation_field=True, step=(0.0, 2.0))\n",
    "def update(deformation_field, step=0.0):\n",
    "    # Determine boundary conditions\n",
    "    u_bc_anim = v_bc + step * (u_bc - v_bc)\n",
    "\n",
    "    if deformation_field:\n",
    "        d_bc = u_bc_anim - v_bc\n",
    "        d = igl.harmonic(v, f, b, d_bc, 2)\n",
    "        u = v + d\n",
    "    else:\n",
    "        u = igl.harmonic(v, f, b, u_bc_anim, 2)\n",
    "    p.update_object(vertices=u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Relationship to \"differential coordinates\" and Laplacian surface editing\n",
    "Biharmonic functions (whether positions or displacements) are solutions to the\n",
    "bi-Laplace equation, but also minimizers of the \"Laplacian energy\". For\n",
    "example, for displacements $\\mathbf{d}$, the energy reads\n",
    "\n",
    " $\\int\\limits_S \\|\\Delta \\mathbf{d}\\|^2 dA,$\n",
    "\n",
    "where we define $\\Delta \\mathbf{d}$ to simply apply the Laplacian\n",
    "coordinate-wise.\n",
    "\n",
    "By linearity of the Laplace(-Beltrami) operator we can reexpress this energy in\n",
    "terms of the original positions $\\mathbf{x}$ and the unknown positions\n",
    "$\\mathbf{x}' = \\mathbf{x} - \\mathbf{d}$:\n",
    "\n",
    " $\\int\\limits_S \\|\\Delta (\\mathbf{x}' - \\mathbf{x})\\|^2 dA = \\int\\limits_S\n",
    " \\|\\Delta \\mathbf{x}' - \\Delta \\mathbf{x})\\|^2 dA.$\n",
    "\n",
    "In the early work of Sorkine et al., the quantities $\\Delta \\mathbf{x}'$ and\n",
    "$\\Delta \\mathbf{x}$ were dubbed \"differential coordinates\"  <cite data-cite=\"sorkine_2002\">(Sorkine, 2004)</cite>.\n",
    "Their deformations (without linearized rotations) is thus equivalent to\n",
    "biharmonic deformation fields."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polyharmonic deformation\n",
    "We can generalize biharmonic deformation by considering different powers of\n",
    "the Laplacian, resulting in a series of PDEs of the form:\n",
    "\n",
    " $\\Delta^k \\mathbf{d} = 0.$\n",
    "\n",
    "with $k\\in{1,2,3,\\dots}$. The choice of $k$ determines the level of continuity\n",
    "at the handles. In particular, $k=1$ implies $C^0$ at the boundary, $k=2$\n",
    "implies $C^1$, $k=3$ implies $C^2$ and in general $k$ implies $C^{k-1}$.\n",
    "\n",
    "The following example deforms a flat domain (left) into a bump as a solution to various $k$-harmonic PDEs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"bump-domain.obj\"))\n",
    "u = v.copy()\n",
    "\n",
    "# Find boundary vertices outside annulus\n",
    "vrn = np.linalg.norm(v, axis = 1)\n",
    "is_outer = [vrn[i] - 1.00 > -1e-15 for i in range(v.shape[0])]\n",
    "is_inner = [vrn[i] - 0.15 < 1e-15 for i in range(v.shape[0])]\n",
    "in_b = [is_outer[i] or is_inner[i] for i in range(len(is_outer))]\n",
    "\n",
    "b = np.array([i for i in range(v.shape[0]) if (in_b[i])]).T\n",
    "bc = np.zeros(b.size)\n",
    "\n",
    "for bi in range(b.size):\n",
    "    bc[bi] = 0.0 if is_outer[b[bi]] else 1.0\n",
    "\n",
    "c = np.array(is_outer)\n",
    "\n",
    "p = plot(u, f, c, shading={\"wire_width\": 0.01, \"colormap\": \"tab10\"}, return_plot=True)                \n",
    "\n",
    "@interact(z_max=(0.0, 1.0), k=(1, 4))\n",
    "def update(z_max, k):\n",
    "    z = igl.harmonic(v, f, b, bc, int(k))\n",
    "    u[:, 2] = z_max * z\n",
    "    p.update_object(vertices=u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### As-rigid-as-possible\n",
    "\n",
    "Skinning and other linear methods for deformation are inherently limited.\n",
    "Difficulties arise especially when large rotations are imposed by the handle constraints.\n",
    "\n",
    "In the context of energy-minimization approaches, the problem stems from\n",
    "comparing positions (our displacements) in the coordinate frame of the\n",
    "undeformed shape. These quadratic energies are at best invariant to global\n",
    "rotations of the entire shape, but not smoothly varying local rotations. Thus\n",
    "linear techniques will not produce non-trivial bending and twisting.\n",
    "\n",
    "Furthermore, when considering solid shapes (e.g. discretized with tetrahedral\n",
    "meshes) linear methods struggle to maintain local volume, and they often suffer from\n",
    "shrinking and bulging artifacts.\n",
    "\n",
    "Non-linear deformation techniques present a solution to these problems.\n",
    "They work by comparing the deformation of a mesh\n",
    "vertex to its rest position _rotated_ to a new coordinate frame which best\n",
    "matches the deformation. The non-linearity stems from the mutual dependence of\n",
    "the deformation and the best-fit rotation. These techniques are often labeled\n",
    "\"as-rigid-as-possible\" as they penalize the sum of all local deformations'\n",
    "deviations from rotations.\n",
    "\n",
    "To arrive at such an energy, let's consider a simple per-triangle energy:\n",
    "\n",
    " $E_\\text{linear}(\\mathbf{X}') = \\sum\\limits_{t \\in T} a_t \\sum\\limits_{\\{i,j\\}\n",
    " \\in t} w_{ij} \\left\\|\n",
    " \\left(\\mathbf{x}'_i - \\mathbf{x}'_j\\right) -\n",
    " \\left(\\mathbf{x}_i - \\mathbf{x}_j\\right)\\right\\|^2$\n",
    "\n",
    "where $\\mathbf{X}'$ are the mesh's unknown deformed vertex positions, $t$ is a\n",
    "triangle in a list of triangles $T$, $a_t$ is the area of triangle $t$ and\n",
    "$\\{i,j\\}$ is an edge in triangle $t$. Thus, this energy measures the norm of\n",
    "change between an edge vector in the original mesh $\\left(\\mathbf{x}_i -\n",
    "\\mathbf{x}_j\\right)$ and the unknown mesh $\\left(\\mathbf{x}'_i -\n",
    "\\mathbf{x}'_j\\right)$.\n",
    "\n",
    "This energy is **not** rotation invariant. If we rotate the mesh by 90 degrees\n",
    "the change in edge vectors not aligned with the axis of rotation will be large,\n",
    "despite the overall deformation being perfectly rigid.\n",
    "\n",
    "So, the \"as-rigid-as-possible\" solution is to append auxiliary variables\n",
    "$\\mathbf{R}_t$\n",
    "for each triangle $t$ which are constrained to be rotations. Then the energy is\n",
    "rewritten, this time comparing deformed edge vectors to their rotated rest\n",
    "counterparts:\n",
    "\n",
    "\n",
    " $E_\\text{arap}(\\mathbf{X}',\\{\\mathbf{R}_1,\\dots,\\mathbf{R}_{|T|}\\}) = \\sum\\limits_{t \\in T} a_t \\sum\\limits_{\\{i,j\\}\n",
    " \\in t} w_{ij} \\left\\|\n",
    " \\left(\\mathbf{x}'_i - \\mathbf{x}'_j\\right)-\n",
    " \\mathbf{R}_t\\left(\\mathbf{x}_i - \\mathbf{x}_j\\right)\\right\\|^2.$\n",
    "\n",
    "The separation into the primary vertex position variables $\\mathbf{X}'$ and the\n",
    "rotations $\\{\\mathbf{R}_1,\\dots,\\mathbf{R}_{|T|}\\}$ lead to strategy for\n",
    "optimization, too. If the rotations $\\{\\mathbf{R}_1,\\dots,\\mathbf{R}_{|T|}\\}$\n",
    "are held fixed then the energy is quadratic in the remaining variables\n",
    "$\\mathbf{X}'$ and can be optimized by solving a (sparse) global linear system.\n",
    "Alternatively, if $\\mathbf{X}'$ are held fixed then each rotation is the\n",
    "solution to a localized _Procrustes_ problem (found via $3 \\times 3$ SVD or\n",
    "polar decompostion). These two steps---local and global---each weakly decrease\n",
    "the energy, thus we may safely iterate them until convergence.\n",
    "\n",
    "The different flavors of \"as-rigid-as-possible\" depend on the dimension and\n",
    "codimension of the domain and the edge-sets $T$. The proposed surface\n",
    "manipulation technique by Sorkine and Alexa  <cite data-cite=\"sorkine_2007\">(Sorkine, 2007)</cite>, considers $T$ to\n",
    "be the set of sets of edges emanating from each vertex (spokes). Later, Chao et\n",
    "al.  derived the relationship between \"as-rigid-as-possible\" mesh energies and\n",
    "co-rotational elasticity considering 0-codimension elements as edge-sets:\n",
    "triangles in 2D and tetrahedra in 3D  <cite data-cite=\"chao_2010\">(Chao, 2010)</cite>. They also showed how\n",
    "Sorkine and Alexa's edge-sets are not a discretization of a continuous energy,\n",
    "proposing instead edge-sets for surfaces containing all edges of elements\n",
    "incident on a vertex (spokes and rims). They show that this amounts to\n",
    "measuring bending, albeit in a discretization-dependent way.\n",
    "\n",
    "Libigl, supports these common flavors. Selecting one is a matter of setting the energy type before the precompuation phase.\n",
    "\n",
    "```python\n",
    "#arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES;\n",
    "#arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES_AND_RIMS;\n",
    "#arap_data.energy = igl::ARAP_ENERGY_TYPE_ELEMENTS;\n",
    "arap = igl.ARAP(v, f, dimension, b)\n",
    "```\n",
    "Just like `igl.min_quad_with_fixed_*`, this precomputation phase only depends on the mesh, fixed vertex indices `b` and the energy parameters. To solve with certain constraints on the positions of vertices in `b`, we may call:\n",
    "\n",
    "```python\n",
    "vn = arap.solve(bc, v)\n",
    "```\n",
    "\n",
    "which uses `v` as an initial guess and then computes the solution into it.\n",
    "\n",
    "Libigl's implementation of as-rigid-as-possible deformation takes advantage of the highly optimized singular value decomposition code from McAdams et al. <cite data-cite=\"mcadams_2011\">(McAdams, 2011)</cite> which leverages SSE intrinsics.\n",
    "\n",
    "The following example deforms a surface as if it were made of an elastic material. The concept of local rigidity will be revisited shortly in the context of surface parameterization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"decimated-knight.off\"))\n",
    "s = igl.read_dmat(os.path.join(root_folder, \"data\", \"decimated-knight-selection.dmat\"))\n",
    "\n",
    "# Vertices in selection\n",
    "b = np.array([[t[0] for t in [(i, s[i]) for i in range(0, v.shape[0])] \n",
    "      if t[1] >= 0]]).T\n",
    "\n",
    "# Centroid\n",
    "mid = 0.5 * (np.max(v, axis=0) + np.min(v, axis=0))\n",
    "\n",
    "# Precomputation\n",
    "arap = igl.ARAP(v, f, 3, b)\n",
    "\n",
    "# Set color based on selection\n",
    "c = np.ones_like(f) * np.array([1.0, 228/255, 58/255])\n",
    "for fi in range(0, f.shape[0]):\n",
    "    if s[f[fi, 0]] >= 0 and s[f[fi, 1]] >= 0 and s[f[fi, 2]] >= 0:\n",
    "        c[fi] = np.array([80/255, 64/255, 1.0])\n",
    "\n",
    "# Plot the mesh with pseudocolors\n",
    "p = plot(v, f, c, return_plot=True)\n",
    "\n",
    "@interact(t=(0.0, 10.0))\n",
    "def update(t=1.0):\n",
    "    bc = np.zeros((b.size, v.shape[1]))\n",
    "    for i in range(0, b.size):\n",
    "        bc[i] = v[b[i]]\n",
    "        if s[b[i]] == 0:\n",
    "            r = mid[0] * 0.25\n",
    "            bc[i, 0] += r * np.sin(0.5 * t * 2 * np.pi)\n",
    "            bc[i, 1] = bc[i, 1] - r + r * np.cos(np.pi + 0.5 * t * 2 * np.pi)\n",
    "        elif s[b[i]] == 1:\n",
    "            r = mid[1] * 0.15\n",
    "            bc[i, 1] = bc[i, 1] + r + r * np.cos(np.pi + 0.15 * t * 2 * np.pi)\n",
    "            bc[i, 2] -= r * np.sin(0.15 * t * 2 * np.pi)\n",
    "        elif s[b[i]] == 2:\n",
    "            r = mid[1] * 0.15\n",
    "            bc[i, 2] = bc[i, 2] + r + r * np.cos(np.pi + 0.35 * t * 2 * np.pi)\n",
    "            bc[i, 0] += r * np.sin(0.35 * t * 2 * np.pi)\n",
    "\n",
    "    vn = arap.solve(bc, v)\n",
    "    p.update_object(vertices=vn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 4: Parametrization\n",
    "\n",
    "In computer graphics, we denote as surface parametrization a map from the\n",
    "surface to \\\\(\\mathbf{R}^2\\\\). It is usually encoded by a new set of 2D\n",
    "coordinates for each vertex of the mesh (and possibly also by a new set of\n",
    "faces in one to one correspondence with the faces of the original surface).\n",
    "Note that\n",
    "this definition is the *inverse* of the classical differential geometry\n",
    "definition.\n",
    "\n",
    "A parametrization has many applications, ranging from texture mapping to\n",
    "surface remeshing. Many algorithms have been proposed, and they can be broadly\n",
    "divided in four families:\n",
    "\n",
    "1. **Single patch, fixed boundary**: these algorithm can parametrize a\n",
    "disk-like part of the surface given fixed 2D positions for its boundary. These\n",
    "algorithms are efficient and simple, but they usually produce high-distortion maps due to the fixed boundary.\n",
    "\n",
    "2. **Single patch, free boundary:** these algorithms let the boundary\n",
    "deform freely, greatly reducing the map distortion. Care should be taken to\n",
    "prevent the border to self-intersect.\n",
    "\n",
    "3. **Global parametrization**: these algorithms work on meshes with arbitrary\n",
    "genus. They initially cut the mesh in multiple patches that can be separately parametrized. The generated maps are discontinuous on the cuts (often referred as *seams*).\n",
    "\n",
    "4. **Global seamless parametrization**: these are global parametrization algorithm that hides the seams, making the parametrization \"continuous\", under specific assumptions that we will discuss later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Harmonic parametrization\n",
    "\n",
    "Harmonic parametrization <cite data-cite=\"eck_2005\">(Eck, 2005)</cite> is a single patch, fixed boundary parametrization\n",
    "algorithm that computes the 2D coordinates of the flattened mesh as two\n",
    "harmonic functions.\n",
    "\n",
    "The algorithm is divided in 3 steps:\n",
    "\n",
    "1. Detection of the boundary vertices\n",
    "2. Map the boundary vertices to a circle\n",
    "3. Compute two harmonic functions (one for u and one for the v coordinate). The harmonic functions use the fixed vertices on the circle as boundary constraints.\n",
    "\n",
    "The algorithm is coded with libigl in the following example. `bnd` contains the indices of the boundary vertices, bnd_uv their position on the UV plane, and \"1\" denotes that we want to compute an harmonic function (2 will be for biharmonic, 3 for triharmonic, etc.). Note that each of the three\n",
    "functions is designed to be reusable in other parametrization algorithms.\n",
    "The UV coordinates are then used to apply a procedural checkerboard texture to the mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f  = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"camelhead.off\"))\n",
    "## Find the open boundary\n",
    "bnd = igl.boundary_loop(f)\n",
    "\n",
    "## Map the boundary to a circle, preserving edge proportions\n",
    "bnd_uv = igl.map_vertices_to_circle(v, bnd)\n",
    "\n",
    "## Harmonic parametrization for the internal vertices\n",
    "uv = igl.harmonic(v, f, bnd, bnd_uv, 1)\n",
    "v_p = np.hstack([uv, np.zeros((uv.shape[0],1))])\n",
    "\n",
    "p = plot(v, f, uv=uv, shading={\"wireframe\": False, \"flat\": False}, return_plot=True)\n",
    "\n",
    "@interact(mode=['3D','2D'])\n",
    "def switch(mode):\n",
    "    if mode == \"3D\":\n",
    "        plot(v, f, uv=uv, shading={\"wireframe\": False, \"flat\": False}, plot=p)\n",
    "    if mode == \"2D\":\n",
    "        plot(v_p, f, uv=uv, shading={\"wireframe\": True, \"flat\": False}, plot=p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least squares conformal maps\n",
    "\n",
    "Least squares conformal maps parametrization <cite data-cite=\"levy_2002\">(Levy, 2002)</cite> minimizes the\n",
    "conformal (angular) distortion of the parametrization. Differently from\n",
    "harmonic parametrization, it does not need to have a fixed boundary.\n",
    "\n",
    "LSCM minimizes the following energy:\n",
    "\n",
    "\\\\[ E_{LSCM}(\\mathbf{u},\\mathbf{v}) = \\int_X \\frac{1}{2}| \\nabla \\mathbf{u}^{\\perp} - \\nabla \\mathbf{v} |^2 dA \\\\]\n",
    "\n",
    "which can be rewritten in matrix form as <cite data-cite=\"mullen_2008\">(Mullen, 2008)</cite>:\n",
    "\n",
    "\\\\[ E_{LSCM}(\\mathbf{u},\\mathbf{v}) = \\frac{1}{2} [\\mathbf{u},\\mathbf{v}]^t (L_c - 2A) [\\mathbf{u},\\mathbf{v}] \\\\]\n",
    "\n",
    "where $L_c$ is the cotangent Laplacian matrix and $A$ is a matrix such that\n",
    "$[\\mathbf{u},\\mathbf{v}]^t A  [\\mathbf{u},\\mathbf{v}]$ is equal to the [vector\n",
    "area](http://en.wikipedia.org/wiki/Vector_area) of the mesh.\n",
    "\n",
    "Using libigl, this matrix energy can be written in a few lines of code. The\n",
    "cotangent matrix can be computed using `igl.cotmatrix`. Note that we want to apply the Laplacian matrix to the u and v coordinates at the same time, thus we need to extend it taking the left Kronecker product with a 2x2 identity matrix. The area matrix is computed with `igl.vector_area_matrix`.\n",
    "\n",
    "The final energy matrix is $L_{flat} - 2A$. Note that in this\n",
    "case we do not need to fix the boundary. To remove the null space of the energy and make the minimum unique, it is sufficient to fix two arbitrary\n",
    "vertices to two arbitrary positions. The full source code is provided in the following LSCM parametrization example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"camelhead.off\"))\n",
    "\n",
    "# Fix two points on the boundary\n",
    "b = np.array([2, 1])\n",
    "\n",
    "bnd = igl.boundary_loop(f)\n",
    "b[0] = bnd[0]\n",
    "b[1] = bnd[int(bnd.size / 2)]\n",
    "\n",
    "bc = np.array([[0.0, 0.0], [1.0, 0.0]])\n",
    "\n",
    "# LSCM parametrization\n",
    "_, uv = igl.lscm(v, f, b, bc)\n",
    "\n",
    "p = plot(v, f, uv=uv, shading={\"wireframe\": False, \"flat\": False}, return_plot=True)\n",
    "\n",
    "@interact(mode=['3D','2D'])\n",
    "def switch(mode):\n",
    "    if mode == \"3D\":\n",
    "        plot(v, f, uv=uv, shading={\"wireframe\": False, \"flat\": False}, plot=p)\n",
    "    if mode == \"2D\":\n",
    "        plot(uv, f, uv=uv, shading={\"wireframe\": True, \"flat\": False}, plot=p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### As-rigid-as-possible parametrization\n",
    "\n",
    "As-rigid-as-possible parametrization <cite data-cite=\"liu\">(Liu, 2008)</cite> is a powerful single-patch, non-linear algorithm to compute a parametrization that strives to preserve\n",
    "distances (and thus angles). The idea is very similar to ARAP surface\n",
    "deformation: each triangle is mapped to the plane trying to preserve its\n",
    "original shape, up to a rigid rotation.\n",
    "\n",
    "The algorithm can be implemented reusing the functions discussed in the\n",
    "deformation chapter: `igl.ARAP` and `arap.solve`. The only\n",
    "difference is that the optimization has to be done in 2D instead of 3D and that\n",
    "we need to compute a starting point. While for 3D deformation the optimization\n",
    "is bootstrapped with the original mesh, this is not the case for ARAP\n",
    "parametrization since the starting point must be a 2D mesh. \n",
    "\n",
    "In the following example, we initialize the optimization with harmonic\n",
    "parametrization. Similarly to LSCM, the boundary is free to deform to minimize\n",
    "the distortion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f  = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"camelhead.off\"))\n",
    "\n",
    "## Find the open boundary\n",
    "bnd = igl.boundary_loop(f)\n",
    "\n",
    "## Map the boundary to a circle, preserving edge proportions\n",
    "bnd_uv = igl.map_vertices_to_circle(v, bnd)\n",
    "\n",
    "## Harmonic parametrization for the internal vertices\n",
    "uv = igl.harmonic(v, f, bnd, bnd_uv, 1)\n",
    "\n",
    "arap = igl.ARAP(v, f, 2, np.zeros(0))\n",
    "uva = arap.solve(np.zeros((0, 0)), uv)\n",
    "\n",
    "p = plot(v, f, uv=uva, shading={\"wireframe\": False, \"flat\": False}, return_plot=True)\n",
    "\n",
    "@interact(mode=['3D','2D'])\n",
    "def switch(mode):\n",
    "    if mode == \"3D\":\n",
    "        plot(v, f, uv=uva, shading={\"wireframe\": False, \"flat\": False}, plot=p)\n",
    "    if mode == \"2D\":\n",
    "        plot(uva, f, uv=uva, shading={\"wireframe\": True, \"flat\": False}, plot=p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### N-rotationally symmetric tangent fields\n",
    "\n",
    "The design of tangent fields is a basic tool used to design guidance fields for uniform quadrilateral and hexahedral remeshing. Libigl contains an implementation of all the state-of-the-art algorithms to design N-RoSy fields and their generalizations.\n",
    "\n",
    "In libigl, tangent unit-length vector fields are piece-wise constant on the faces of a triangle mesh, and they are described by one or more vectors per-face. The function\n",
    "\n",
    "```python\n",
    "output_field, output_singularities = igl.nrosy(v, f, b, bc, b_soft, b_soft_weight, bc_soft, n, 0.5)\n",
    "```\n",
    "\n",
    "creates a smooth unit-length vector field (n=1) starting from a sparse set of constrained faces, whose indices are listed in b and their constrained value is specified in bc. The functions supports soft_constraints (b_soft, b_soft_weight, bc_soft), and returns the interpolated field for each face of the triangle mesh (output_field), plus the singularities of the field (output_singularities).\n",
    "\n",
    "The singularities are vertices where the field vanishes (highlighted in red in the figure above). `igl.nrosy` can also generate N-RoSy fields <cite data-cite=\"levy_2008\">(Levy, 2008)</cite>, which are a generalization of vector fields where in every face the vector is defined up to a constant rotation of $2\\pi / N$. As can be observed in the following figure, the singularities of the fields generated with different N are of different types and they appear in different positions.\n",
    "\n",
    "We demonstrate how to call and plot N-RoSy fields in the following example, where the degree of the field can be changed. `igl.nrosy` implements the algorithm proposed in <cite data-cite=\"bommes_2009\">(Bommes, 2009)</cite>. N-RoSy fields can also be interpolated with many other algorithms, see the library [libdirectional](https://github.com/avaxman/libdirectional) for a reference implementation of the most popular ones. For a complete categorization of fields used in various applications see <cite data-cite=\"vaxman_2016\">(Vaxman, 2016)</cite>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Converts a representative vector per face in the full set of vectors that describe an N-RoSy field\n",
    "# def representative_to_nrosy(v, f, r, n):\n",
    "#     b1, b2, b3 = igl.local_basis(v, f)\n",
    "#     ym = np.zeros((f.shape[0] * n, 3))\n",
    "\n",
    "#     for i in range(f.shape[0]):\n",
    "#         x = r[i] * b1[i].T\n",
    "#         y = r[i] * b2[i].T\n",
    "#         angle = np.arctan2(y[0], x[0])\n",
    "\n",
    "#         for j in range(0, n):\n",
    "#             anglej = angle + 2 * np.pi * j / float(n)\n",
    "#             xj = np.cos(anglej)\n",
    "#             yj = np.sin(anglej)\n",
    "#             ym[i * n + j] = xj * b1[i] + yj * b2[i]\n",
    "#     return ym\n",
    "\n",
    "# # Load a mesh in OFF format and plot it\n",
    "# v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"bumpy.off\"))\n",
    "#     #print(points_n)\n",
    "    \n",
    "# # Constrained faces id\n",
    "# b = np.array([[0]])\n",
    "\n",
    "# # Highlight in red the constrained faces\n",
    "# c = np.ones((f.shape[0], 3))\n",
    "# for i in range(b.size):\n",
    "#     c[b[i]] = np.array([1, 0, 0])\n",
    "\n",
    "# p = plot(v, f, c)\n",
    "# pn_id = pp_id = e_id = None\n",
    "\n",
    "\n",
    "# # Constrained faces representative vector\n",
    "# bcv = np.array([[1.0, 1.0, 1.0]])\n",
    "# avg = igl.avg_edge_length(v, f)\n",
    "# bc = igl.barycenter(v, f)\n",
    "\n",
    "# # Plots the mesh with an N-RoSy field and its singularities on top\n",
    "# # The constrained faces (b) are colored in red.\n",
    "# @interact(n=(1, 10))\n",
    "# def plot_mesh_nrosy(n=1):\n",
    "#     global pn_id, pp_id, e_id, bcv, avg, b, bc\n",
    "\n",
    "#     r, s = igl.nrosy(v, f, b, bcv, np.array([[]], dtype=np.int64), np.array([[]]), np.array([[]]), n, 0.5)\n",
    "    \n",
    "#     # Expand the representative vectors in the full vector set and plot them as lines\n",
    "#     y = representative_to_nrosy(v, f, r, n)\n",
    "#     be = np.zeros((bc.shape[0] * n, 3))\n",
    "#     for i in range(bc.shape[0]):\n",
    "#         for j in range(n):\n",
    "#             be[i * n + j] = bc[i]\n",
    "\n",
    "#     if e_id:\n",
    "#         p.remove_object(e_id)\n",
    "#     e_id = p.add_lines(be, be + y * (avg / 2))\n",
    "\n",
    "#     # Plot the singularities as colored dots (red for negative, blue for positive)\n",
    "#     points_n = []\n",
    "#     points_p = []\n",
    "#     for i in range(0, s.size):\n",
    "#         if s[i] < -0.001:\n",
    "#             points_n.append(v[i])\n",
    "#         elif s[i] > 0.001:\n",
    "#             points_p.append(v[i])\n",
    "    \n",
    "#     if pp_id and pn_id:\n",
    "#         p.remove_object(pn_id)\n",
    "#         p.remove_object(pp_id)\n",
    "#     if len(points_n) > 0:\n",
    "#         pn_id = p.add_points(np.array(points_n), c=\"red\", shading={\"point_size\": 2.0})\n",
    "#     if len(points_p) > 0:\n",
    "#         pp_id = p.add_points(np.array(points_p), c=\"blue\", shading={\"point_size\": 2.0})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Planarization\n",
    "\n",
    "A quad mesh can be transformed in a planar quad mesh with Shape-Up <cite data-cite=\"bouaziz_2012\">(Bouaziz, 2012)</cite>, a local/global approach that uses the global step to enforce surface continuity and the local step to enforce planarity.\n",
    "\n",
    "The following example planarizes a quad mesh until it satisfies a user-given planarity threshold. The colors represent the planarity of the quads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load a quad mesh generated by a conjugate field\n",
    "vqc, fqc, _ = igl.read_off(os.path.join(root_folder, \"data\", \"inspired_mesh_quads_Conjugate.off\"))\n",
    "\n",
    "# Convert it to a triangle mesh\n",
    "fqc_tri = np.zeros((fqc.shape[0] * 2, 3), dtype=\"int64\")\n",
    "fqc_tri[:fqc.shape[0]] = fqc[:, :3]\n",
    "fqc_tri[fqc.shape[0]:, 0] = fqc[:, 2]\n",
    "fqc_tri[fqc.shape[0]:, 1] = fqc[:, 3]\n",
    "fqc_tri[fqc.shape[0]:, 2] = fqc[:, 0]\n",
    "\n",
    "# Planarize it\n",
    "vqc_p = igl.planarize_quad_mesh(vqc, fqc, 100, 0.005)\n",
    "\n",
    "# Calculate a color to each quad that corresponds to its planarity\n",
    "planarity = igl.quad_planarity(vqc, fqc)\n",
    "planarity_p = igl.quad_planarity(vqc_p, fqc)\n",
    "\n",
    "c = np.concatenate([planarity, planarity])\n",
    "c_p = np.concatenate([planarity_p, planarity_p])\n",
    "\n",
    "p = plot(vqc, fqc_tri, c, shading={\"normalize\":[min(np.min(c), np.min(c_p)), max(np.max(c), np.max(c_p))]})\n",
    "\n",
    "@interact(mode=['Curved','Planar'])\n",
    "def switch(mode):\n",
    "    if mode == \"Curved\":\n",
    "        p.update_object(colors=c)\n",
    "    if mode == \"Planar\":\n",
    "        p.update_object(vertices=vqc_p, colors=c_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 5: External libraries\n",
    "\n",
    "An additional positive side effect of using matrices as basic types is that it\n",
    "is easy to exchange data between libigl and other software and libraries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baking ambient occlusion\n",
    "\n",
    "[Ambient occlusion](http://en.wikipedia.org/wiki/Ambient_occlusion) is a\n",
    "rendering technique used to calculate the exposure of each point in a surface\n",
    "to ambient lighting. It is usually encoded as a scalar (normalized between 0\n",
    "and 1) associated with the vertice of a mesh.\n",
    "\n",
    "Formally, ambient occlusion is defined as:\n",
    "\n",
    "\\\\[ A_p = \\frac{1}{\\pi} \\int_\\omega V_{p,\\omega}(n \\cdot \\omega) d\\omega \\\\]\n",
    "\n",
    "where $V_{p,\\omega}$ is the visibility function at  p, defined to be zero if p\n",
    "is occluded in the direction $\\omega$ and one otherwise, and $d\\omega$ is the\n",
    "infinitesimal solid angle step of the integration variable $\\omega$.\n",
    "\n",
    "The integral is usually approximated by casting rays in random directions\n",
    "around each vertex. This approximation can be computed using the function:\n",
    "\n",
    "```\n",
    "ao = igl.ambient_occlusion(v, f, v_samples, n_samples, 500)\n",
    "```\n",
    "\n",
    "that given a scene described in `v` and `f`, computes the ambient occlusion of\n",
    "the points in `v_samples` whose associated normals are `n_samples`. The\n",
    "number of casted rays can be controlled (usually at least 300-500 rays are\n",
    "required to get a smooth result) and the result is returned in `ao`, as a\n",
    "single scalar for each sample.\n",
    "\n",
    "Ambient occlusion can be used to darken the surface colors, as shown in the following example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"fertility.off\"))\n",
    "\n",
    "n = igl.per_vertex_normals(v, f)\n",
    "\n",
    "# Compute ambient occlusion factor using embree\n",
    "ao = igl.ambient_occlusion(v, f, v, n, 50)\n",
    "ao = 1.0 - ao\n",
    "\n",
    "plot(v, f, ao, shading={\"colormap\": \"gist_gray\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 6: Miscellaneous\n",
    "\n",
    "Libigl contains a _wide_ variety of geometry processing tools and functions for\n",
    "dealing with meshes and the linear algebra related to them: far too many to\n",
    "discuss in this introductory tutorial. We've pulled out a couple of the\n",
    "interesting functions in this chapter to highlight."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh Statistics\n",
    "\n",
    "Libigl contains various mesh statistics, including face angles, face areas and\n",
    "the detection of singular vertices, which are vertices with more or less than 6\n",
    "neighbours in triangulations or 4 in quadrangulations.\n",
    "\n",
    "The example computes these quantities and\n",
    "does a basic statistic analysis that allows to estimate the isometry and\n",
    "regularity of a mesh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"horse_quad.obj\"))\n",
    "\n",
    "## Count the number of irregular vertices, the border is ignored\n",
    "irregular = igl.is_irregular_vertex(v, f) \n",
    "v_count = v.shape[0]\n",
    "irregular_v_count = np.sum(irregular)\n",
    "irregular_ratio = irregular_v_count / v_count\n",
    "\n",
    "print(\"Irregular vertices: \\n%d/%d (%.2f%%)\\n\"%(irregular_v_count, v_count, irregular_ratio * 100))\n",
    "\n",
    "## Compute areas, min, max and standard deviation\n",
    "area = igl.doublearea(v, f) / 2.0\n",
    "\n",
    "area_avg = np.mean(area)\n",
    "area_min = np.min(area) / area_avg\n",
    "area_max = np.max(area) / area_avg\n",
    "area_ns = (area - area_avg) / area_avg\n",
    "area_sigma = np.sqrt(np.mean(np.square(area_ns)))\n",
    "\n",
    "print(\"Areas (Min/Max)/Avg_Area Sigma: \\n%.2f/%.2f (%.2f)\\n\"%(area_min, area_max, area_sigma))\n",
    "\n",
    "## Compute per face angles, min, max and standard deviation\n",
    "angles = igl.internal_angles(v, f)\n",
    "angles = 360.0 * (angles / (2 * np.pi))\n",
    "\n",
    "angle_avg = np.mean(angles)\n",
    "angle_min = np.min(angles)\n",
    "angle_max = np.max(angles)\n",
    "angle_ns = angles - angle_avg\n",
    "angle_sigma = np.sqrt(np.mean(np.square(angle_ns)))\n",
    "\n",
    "print(\"Angles in degrees (Min/Max) Sigma: \\n%.2f/%.2f (%.2f)\\n\"%(angle_min, angle_max, angle_sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first row contains the number and percentage of irregular vertices, which\n",
    "is particularly important for quadrilateral meshes when they are used to define\n",
    "subdivision surfaces: every singular point will result in a point of the\n",
    "surface that is only C^1.\n",
    "\n",
    "The second row reports the area of the minimal element, maximal element and the\n",
    "standard deviation.  These numbers are normalized by the mean area, so in the\n",
    "example above 5.33 max area means that the biggest face is 5 times larger than\n",
    "the average face. An ideal isotropic mesh would have both min and max area\n",
    "close to 1.\n",
    "\n",
    "The third row measures the face angles, which should be close to 60 degrees (90\n",
    "for quads) in a perfectly regular triangulation. For FEM purposes, the closer\n",
    "the angles are to 60 degrees the more stable will the optimization be. In this\n",
    "case, it is clear that the mesh is of bad quality and it will probably result\n",
    "in artifacts if used for solving PDEs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subdivision surfaces\n",
    "\n",
    "Given a coarse mesh (aka cage) with vertices `V` and faces `F`, one can createa\n",
    "higher-resolution mesh with more vertices and faces by _subdividing_ every\n",
    "face. That is, each coarse triangle in the input is replaced by many smaller\n",
    "triangles. Libigl has three different methods for subdividing a triangle mesh.\n",
    "\n",
    "An \"in plane\" subdivision method will not change the point set or carrier\n",
    "surface of the mesh. New vertices are added on the planes of existing triangles\n",
    "and vertices surviving from the original mesh are not moved.\n",
    "\n",
    "By adding new faces, a subdivision algorithm changes the _combinatorics_ of the\n",
    "mesh. The change in combinatorics and the formula for positioning the\n",
    "high-resolution vertices is called the \"subdivision rule\".\n",
    "\n",
    "For example, in the _in plane_ subdivision method of `igl.upsample`, vertices\n",
    "are added at the midpoint of every edge: $v_{ab} = \\frac{1}{2}(v_a + v_b)$ and\n",
    "each triangle $(i_a,i_b,i_c)$ is replaced with four triangles:\n",
    "$(i_a,i_{ab},i_{ca})$, $(i_b,i_{bc},i_{ab})$, $(i_{ab},i_{bc},i_{ca})$, and\n",
    "$(i_{bc},i_{c},i_{ca})$. This process may be applied recursively, resulting in\n",
    "a finer and finer mesh.\n",
    "\n",
    "The subdivision method of `igl.loop` is not in plane. The vertices of the\n",
    "refined mesh are moved to weight combinations of their neighbors: the mesh is\n",
    "smoothed as it is refined <cite data-cite=\"loop_1987\">(Loop, 1987)</cite>. This and other _smooth subdivision_\n",
    "methods can be understood as generalizations of spline curves to surfaces. In\n",
    "particular the Loop subdivision method will converge to a $C^1$ surface as we\n",
    "consider the limit of recursive applications of subdivision. Away from\n",
    "\"irregular\" or \"extraordinary\" vertices (vertices of the original cage with\n",
    "valence not equal to 6), the surface is $C^2$. The combinatorics (connectivity\n",
    "and number of faces) of `igl.loop` and `igl.upsample` are identical: the only\n",
    "difference is that the vertices have been smoothed in `igl.loop`.\n",
    "\n",
    "Finally, libigl also implements a form of _in plane_ \"false barycentric\n",
    "subdivision\" in `igl.false_barycentric_subdivision`. This method simply adds\n",
    "the barycenter of every triangle as a new vertex $v_{abc}$ and replaces each\n",
    "triangle with three triangles $(i_a,i_b,i_{abc})$, $(i_b,i_c,i_{abc})$, and\n",
    "$(i_c,i_a,i_{abc})$. In contrast to `igl.upsample`, this method will create\n",
    "triangles with smaller and smaller internal angles and new vertices will sample\n",
    "the carrier surfaces with extreme bias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ov, of = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"decimated-knight.off\"))\n",
    "uv, uf = igl.upsample(ov, of)\n",
    "lv, lf = igl.loop(ov, of)\n",
    "\n",
    "p = plot(ov, of, shading={\"wireframe\": True})\n",
    "\n",
    "@interact(mode=['Coarse','Upsample', 'Loop'])\n",
    "def switch(mode):\n",
    "    if mode == \"Coarse\":\n",
    "        plot(ov, of, shading={\"wireframe\": True}, plot=p)\n",
    "    if mode == \"Upsample\":\n",
    "        plot(uv, uf, shading={\"wireframe\": True}, plot=p)\n",
    "        #p.update_object(vertices=uv, faces=uf)\n",
    "    if mode == \"Loop\":\n",
    "        plot(lv, lf, shading={\"wireframe\": True}, plot=p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data smoothing\n",
    "\n",
    "A noisy function $f$ defined on a surface $\\Omega$ can be smoothed using an energy minimization that balances a smoothing term $E_S$ with a quadratic fitting term:\n",
    "\n",
    "$u = \\operatorname{argmin}_u \\alpha E_S(u) + (1-\\alpha)\\int_\\Omega ||u-f||^2 dx$\n",
    "\n",
    "The parameter $\\alpha$ determines how aggressively the function is smoothed.\n",
    "\n",
    "A classical choice for the smoothness energy is the Laplacian energy of the function with zero Neumann boundary conditions, which is a form of the biharmonic energy. It is constructed using the cotangent Laplacian `L` and\n",
    "the mass matrix `M`: `QL = L'*(M\\L)`. Because of the implicit zero Neumann boundary conditions however, the function behavior is significantly warped at the boundary if $f$ does not have zero normal gradient at the boundary.\n",
    "\n",
    "In <cite data-cite=\"stein_2017\">(Stein, 2017)</cite> it is suggested to use the Biharmonic energy with natural\n",
    "Hessian boundary conditions instead, which corresponds to the hessian energy with the matrix `QH = H'*(M2\\H)`, where `H` is a finite element Hessian and `M2` is a stacked mass matrix. The matrices `H` and `QH` are implemented in\n",
    "libigl as `igl.hessian` and `igl.hessian_energy` respectively. \n",
    "\n",
    "In the following example the differences between the Laplacian energy with zero Neumann boundary conditions and the Hessian energy can be clearly seen: whereas the zero Neumann boundary condition in the third image bias the isolines\n",
    "of the function to be perpendicular to the boundary, the Hessian energy gives an unbiased result.\n",
    "\n",
    "The following example shows a function on the beetle mesh, the function with added noise, the result of smoothing with the Laplacian energy and zero Neumann boundary conditions, and the result of smoothing with the Hessian energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh(os.path.join(root_folder, \"data\", \"beetle.off\"))\n",
    "e = igl.edges(f)\n",
    "\n",
    "# Constructing an exact function to smooth\n",
    "z_exact = v[0:, 2] + 0.5 * v[0:, 1] + v[0:, 1] * v[0:, 1] + v[0:, 2] * v[0:, 2] * v[0:, 2]\n",
    "    \n",
    "# Make the exact function noisy\n",
    "s = 0.2 * (np.max(z_exact) - np.min(z_exact))\n",
    "np.random.seed(5)\n",
    "z_noisy = z_exact + s * np.random.rand(*z_exact.shape)\n",
    "\n",
    "# Constructing the squared Laplacian and squared Hessian energy\n",
    "l = igl.cotmatrix(v, f)\n",
    "m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_BARYCENTRIC)\n",
    "\n",
    "m_inv_l = sp.sparse.linalg.spsolve(m, l)\n",
    "ql = l.T @ m_inv_l\n",
    "qh = igl.hessian_energy(v, f)\n",
    "\n",
    "# Solve to find Laplacian-smoothed and Hessian-smoothed solutions\n",
    "al = 8e-4;\n",
    "zl = sp.sparse.linalg.spsolve(al * ql + (1 - al) * m, al * m.dot(z_noisy))\n",
    "ah = 5e-6;\n",
    "zh = sp.sparse.linalg.spsolve(ah * qh + (1 - ah) * m, ah * m.dot(z_noisy))\n",
    "\n",
    "# Calculate isolines\n",
    "ilx_v, ilx_e = igl.isolines(v, f, z_exact, 30)\n",
    "iln_v, iln_e = igl.isolines(v, f, z_noisy, 30)\n",
    "ill_v, ill_e = igl.isolines(v, f, zl, 30)\n",
    "ilh_v, ilh_e = igl.isolines(v, f, zh, 30)\n",
    "\n",
    "p = plot(v, f, z_exact, return_plot=True)\n",
    "e_id = p.add_edges(ilx_v, ilx_e)\n",
    "\n",
    "@interact(mode=['Original', 'Noisy', 'Biharmonic smoothing (0-Neumann)', 'Biharmonic smoothing (Natural Hessian)'])\n",
    "def switch(mode):\n",
    "    global e_id\n",
    "    p.remove_object(e_id)\n",
    "    if mode == \"Original\":\n",
    "        p.update_object(colors=z_exact)\n",
    "        e_id = p.add_edges(ilx_v, ilx_e)\n",
    "    if mode == \"Noisy\":\n",
    "        p.update_object(colors=z_noisy)\n",
    "        e_id = p.add_edges(iln_v, iln_e)\n",
    "    if mode == \"Biharmonic smoothing (0-Neumann)\":\n",
    "        p.update_object(colors=zl)\n",
    "        e_id = p.add_edges(ill_v, ill_e)\n",
    "    if mode == \"Biharmonic smoothing (Natural Hessian)\":\n",
    "        p.update_object(colors=zh)\n",
    "        e_id = p.add_edges(ilh_v, ilh_e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marching Tetrahedra\n",
    "\n",
    "Often 3D data is captured as scalar field defined over space $f(\\mathbf{x}) : \\mathcal{R}^3 \\rightarrow \\mathcal{R}$. Lurking within this field, _iso-surfaces_ of the scalar field are often salient geometric objects. The\n",
    "iso-surface at value $v$ is composed of all points $\\mathbf{x}$ in $\\mathcal{R}^3$ such that $f(\\mathbf{x}) = v$. A core problem in geometry processing is to extract an iso-surface as a triangle mesh for further mesh-based processing or visualization. This is referred to as iso-contouring.\n",
    "\n",
    "\"Marching Tetrahedra\" <cite data-cite=\"treece_1999\">(Treece, 1999)</cite> is a [famous method](https://en.wikipedia.org/wiki/Marching_tetrahedra) for iso-contouring tri-linear functions $f$ on a 3D simplicial complex (aka a tet mesh). The core idea of this method is to contour the iso-surface passing through each cell  (if it does at all) with a predefined topology (aka connectivity) chosen from a look up tabledepending on the function values at each vertex of the cell. The method\n",
    "iterates (\"marches\") over all cells (\"tetrahedra\") in the complex and stitches together the final mesh.\n",
    "\n",
    "In libigl, `igl.marching_tets` constructs a triangle mesh `(v,f)` approximating the iso-level set for the value `isovalue` from an input scalar field `s` sampled at the vertices of a tet mesh locations `(tv, tt)`:\n",
    "\n",
    "```python\n",
    "v, f = igl.marching_tets(tv, tt, s, isovalue)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tv = np.load(os.path.join(root_folder, \"data\", \"marching_cube_tv.npy\"))\n",
    "tt = np.load(os.path.join(root_folder, \"data\", \"marching_cube_tt.npy\"))\n",
    "s = np.linalg.norm(tv, axis=1)\n",
    "\n",
    "svs = []\n",
    "sfs = []\n",
    "for i in np.linspace(0.05, 0.75, 15):\n",
    "    sv, sf, _, _ = igl.marching_tets(tv, tt, s, i)\n",
    "    svs.append(sv)\n",
    "    sfs.append(sf)\n",
    "\n",
    "p = plot(sv, sf, return_plot=True)\n",
    "oid = 0\n",
    "\n",
    "@interact(t=(0, 14))\n",
    "def update(t=0):\n",
    "    global oid\n",
    "    p.remove_object(oid)\n",
    "    oid = p.add_mesh(svs[t], sfs[t])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "<!-- Chapter 2 -->\n",
    "\n",
    "[^jacobson_thesis_2013]: Alec Jacobson, [_Algorithms and Interfaces for Real-Time Deformation of 2D and 3D Shapes_](https://www.google.com/search?q=Algorithms+and+Interfaces+for+Real-Time+Deformation+of+2D+and+3D+Shapes), 2013.\n",
    "[^kazhdan_2012]: Michael Kazhdan, Jake Solomon, Mirela Ben-Chen, [Can Mean-Curvature Flow Be Made Non-Singular](https://www.google.com/search?q=Can+Mean-Curvature+Flow+Be+Made+Non-Singular), 2012.\n",
    "[^meyer_2003]: Mark Meyer, Mathieu Desbrun, Peter Schröder and Alan H.  Barr, [Discrete Differential-Geometry Operators for Triangulated 2-Manifolds](https://www.google.com/search?q=Discrete+Differential-Geometry+Operators+for+Triangulated+2-Manifolds), 2003.\n",
    "[^mitchell_1987]: Joseph S. B. Mitchell, David M. Mount, Christos H. Papadimitriou. [The Discrete Geodesic Problem](https://www.google.com/search?q=The+Discrete+Geodesic+Problem), 1987\n",
    "[^panozzo_2010]: Daniele Panozzo, Enrico Puppo, Luigi Rocca, [Efficient Multi-scale Curvature and Crease Estimation](https://www.google.com/search?q=Efficient+Multi-scale+Curvature+and+Crease+Estimation), 2010.\n",
    "[^sharf_2007]: Andrei Sharf, Thomas Lewiner, Gil Shklarski, Sivan Toledo, and Daniel Cohen-Or. [Interactive topology-aware surface reconstruction](https://www.google.com/search?q=Interactive+topology-aware+surface+reconstruction), 2007.\n",
    "\n",
    "<!-- Chapter 3 -->\n",
    "\n",
    "[^barbic_2005]: Jernej Barbic and Doug James. [Real-Time Subspace Integration for St.Venant-Kirchhoff Deformable Models](https://www.google.com/search?q=Real-Time+Subspace+Integration+for+St.Venant-Kirchhoff+Deformable+Models), 2005.\n",
    "[^hildebrandt_2011]: Klaus Hildebrandt, Christian Schulz, Christoph von Tycowicz, and Konrad Polthier. [Interactive Surface Modeling using Modal Analysis](https://www.google.com/search?q=Interactive+Surface+Modeling+using+Modal+Analysis), 2011.\n",
    "[^rustamov_2011]: Raid M. Rustamov, [Multiscale Biharmonic Kernels](https://www.google.com/search?q=Multiscale+Biharmonic+Kernels), 2011.\n",
    "[^vallet_2008]: Bruno Vallet and Bruno Lévy. [Spectral Geometry Processing with Manifold Harmonics](https://www.google.com/search?q=Spectral+Geometry+Processing+with+Manifold+Harmonics), 2008.\n",
    "\n",
    "<!-- Chapter 4 -->\n",
    "\n",
    "[^botsch_2004]: Matrio Botsch and Leif Kobbelt. [An Intuitive Framework for Real-Time Freeform Modeling](https://www.google.com/search?q=An+Intuitive+Framework+for+Real-Time+Freeform+Modeling), 2004.\n",
    "[^chao_2010]: Isaac Chao, Ulrich Pinkall, Patrick Sanan, Peter Schröder. [A Simple Geometric Model for Elastic Deformations](https://www.google.com/search?q=A+Simple+Geometric+Model+for+Elastic+Deformations), 2010.\n",
    "[^jacobson_2011]: Alec Jacobson, Ilya Baran, Jovan Popović, and Olga Sorkine. [Bounded Biharmonic Weights for Real-Time Deformation](https://www.google.com/search?q=Bounded+biharmonic+weights+for+real-time+deformation), 2011.\n",
    "[^jacobson_2012]: Alec Jacobson, Ilya Baran, Ladislav Kavan, Jovan Popović, and Olga Sorkine. [Fast Automatic Skinning Transformations](https://www.google.com/search?q=Fast+Automatic+Skinning+Transformations), 2012.\n",
    "[^jacobson_mixed_2010]: Alec Jacobson, Elif Tosun, Olga Sorkine, and Denis Zorin. [Mixed Finite Elements for Variational Surface Modeling](https://www.google.com/search?q=Mixed+Finite+Elements+for+Variational+Surface+Modeling), 2010.\n",
    "[^jacobson_skinning_course_2014]: Alec Jacobson, Zhigang Deng, Ladislav Kavan, J.P. Lewis. [_Skinning: Real-Time Shape Deformation_](https://www.google.com/search?q=Skinning+Real-Time+Shape+Deformation), 2014.\n",
    "[^kavan_2008]: Ladislav Kavan, Steven Collins, Jiri Zara, and Carol O'Sullivan. [Geometric Skinning with Approximate Dual Quaternion Blending](https://www.google.com/search?q=Geometric+Skinning+with+Approximate+Dual+Quaternion+Blending), 2008.\n",
    "[^mcadams_2011]: Alexa McAdams, Andrew Selle, Rasmus Tamstorf, Joseph Teran, Eftychios Sifakis. [Computing the Singular Value Decomposition of 3x3 matrices with minimal branching and elementary floating point operations](https://www.google.com/search?q=Computing+the+Singular+Value+Decomposition+of+3x3+matrices+with+minimal+branching+and+elementary+floating+point+operations), 2011.\n",
    "[^sorkine_2004]: Olga Sorkine, Yaron Lipman, Daniel Cohen-Or, Marc Alexa, Christian Rössl and Hans-Peter Seidel. [Laplacian Surface Editing](https://www.google.com/search?q=Laplacian+Surface+Editing), 2004.\n",
    "[^sorkine_2007]: Olga Sorkine and Marc Alexa. [As-rigid-as-possible Surface Modeling](https://www.google.com/search?q=As-rigid-as-possible+Surface+Modeling), 2007.\n",
    "[^wang_bc_2015]: Yu Wang, Alec Jacobson, Jernej Barbic, Ladislav Kavan. [Linear Subspace Design for Real-Time Shape Deformation](https://www.google.com/search?q=Linear+Subspace+Design+for+Real-Time+Shape+Deformation), 2015\n",
    "\n",
    "<!-- Chapter 5 -->\n",
    "\n",
    "[^bommes_2009]: David Bommes, Henrik Zimmer, Leif Kobbelt. [Mixed-integer quadrangulation](http://www-sop.inria.fr/members/David.Bommes/publications/miq.pdf), 2009.\n",
    "[^bouaziz_2012]: Sofien Bouaziz, Mario Deuss, Yuliy Schwartzburg, Thibaut Weise, Mark Pauly [Shape-Up: Shaping Discrete Geometry with Projections](http://lgg.epfl.ch/publications/2012/shapeup.pdf), 2012\n",
    "[^eck_2005]: Matthias Eck, Tony DeRose, Tom Duchamp, Hugues Hoppe, Michael Lounsbery, Werner Stuetzle.  [Multiresolution Analysis of Arbitrary Meshes](http://research.microsoft.com/en-us/um/people/hoppe/mra.pdf), 2005.\n",
    "[^levy_2002]: Bruno Lévy, Sylvain Petitjean, Nicolas Ray, Jérome Maillot. [Least Squares Conformal Maps, for Automatic Texture Atlas Generation](http://www.cs.jhu.edu/~misha/Fall09/Levy02.pdf), 2002.\n",
    "[^levy_2008]: Nicolas Ray, Bruno Vallet, Wan Chiu Li, Bruno Lévy. [N-Symmetry Direction Field Design](http://alice.loria.fr/publications/papers/2008/DGF/NSDFD-TOG.pdf), 2008.\n",
    "[^liu_2008]: Ligang Liu, Lei Zhang, Yin Xu, Craig Gotsman, Steven J. Gortler. [A Local/Global Approach to Mesh Parameterization](http://cs.harvard.edu/~sjg/papers/arap.pdf), 2008.\n",
    "[^mullen_2008]: Patrick Mullen, Yiying Tong, Pierre Alliez, Mathieu Desbrun. [Spectral Conformal Parameterization](http://www.geometry.caltech.edu/pubs/MTAD08.pdf), 2008.\n",
    "[^panozzo_2014]: Daniele Panozzo, Enrico Puppo, Marco Tarini, Olga Sorkine-Hornung.  [Frame Fields: Anisotropic and Non-Orthogonal Cross Fields](http://cs.nyu.edu/~panozzo/papers/frame-fields-2014.pdf), 2014.\n",
    "[^vaxman_2016]: Amir Vaxman, Marcel Campen, Olga Diamanti, Daniele Panozzo, David Bommes, Klaus Hildebrandt, Mirela Ben-Chen. [Directional Field Synthesis, Design, and Processing](https://www.google.com/search?q=Directional+Field+Synthesis+Design+and+Processing), 2016\n",
    "\n",
    "<!-- Chapter 6 -->\n",
    "\n",
    "[^schuller_2013]: Christian Schüller, Ladislav Kavan, Daniele Panozzo, Olga Sorkine-Hornung.  [Locally Injective Mappings](http://igl.ethz.ch/projects/LIM/), 2013.\n",
    "[^zhou_2016]: Qingnan Zhou, Eitan Grinspun, Denis Zorin. [Mesh Arrangements for Solid Geometry](https://www.google.com/search?q=Mesh+Arrangements+for+Solid+Geometry), 2016\n",
    "\n",
    "<!-- Chapter 7 -->\n",
    "\n",
    "[^baerentzen_2005]: J Andreas Baerentzen and Henrik Aanaes. [Signed distance computation using the angle weighted pseudonormal](https://www.google.com/search?q=Signed+distance+computation+using+the+angle+weighted+pseudonormal), 2005.\n",
    "[^bouaziz_2012]: Sofien Bouaziz, Mario Deuss, Yuliy Schwartzburg, Thibaut Weise, Mark Pauly [Shape-Up: Shaping Discrete Geometry with Projections](http://lgg.epfl.ch/publications/2012/shapeup.pdf), 2012\n",
    "[^garg_2016]: Akash Garg, Alec Jacobson, Eitan Grinspun. [Computational Design of Reconfigurables](https://www.google.com/search?q=Computational+Design+of+Reconfigurables), 2016\n",
    "[^hoppe_1996]: Hugues Hoppe. [Progressive Meshes](https://www.google.com/search?q=Progressive+meshes), 1996\n",
    "[^jacobson_2013]: Alec Jacobson, Ladislav Kavan, and Olga Sorkine. [Robust Inside-Outside Segmentation using Generalized Winding Numbers](https://www.google.com/search?q=Robust+Inside-Outside+Segmentation+using+Generalized+Winding+Numbers), 2013.\n",
    "[^loop_1987]: Charles Loop. [Smooth Subdivision Surfaces Based on Triangles](https://www.google.com/search?q=smooth+subdivision+surfaces+based+on+triangles), 1987.\n",
    "[^lorensen_1987]: W.E. Lorensen and Harvey E. Cline. [Marching cubes: A high resolution 3d surface construction algorithm](https://www.google.com/search?q=Marching+cubes:+A+high+resolution+3d+surface+construction+algorithm), 1987.\n",
    "[^rabinovich_2016]: Michael Rabinovich, Roi Poranne, Daniele Panozzo, Olga Sorkine-Hornung. [Scalable Locally Injective Mappings](http://cs.nyu.edu/~panozzo/papers/SLIM-2016.pdf), 2016.\n",
    "[^schroeder_1994]: William J. Schroeder, William E. Lorensen, and Steve Linthicum. [Implicit Modeling of Swept Surfaces and Volumes](https://www.google.com/search?q=implicit+modeling+of+swept+surfaces+and+volumes), 1994.\n",
    "[^takayama14]: Kenshi Takayama, Alec Jacobson, Ladislav Kavan, Olga Sorkine-Hornung. [A Simple Method for Correcting Facet Orientations in Polygon Meshes Based on Ray Casting](https://www.google.com/search?q=A+Simple+Method+for+Correcting+Facet+Orientations+in+Polygon+Meshes+Based+on+Ray+Casting), 2014.\n",
    "[^treece_1999]: G.M. Treece, R.W. Prager, and A.H.Gee [Regularised marching tetrahedra: improved iso-surface extraction](https://www.sciencedirect.com/science/article/pii/S009784939900076X), 1999.\n",
    "[^crane_2013]: Keenan Crane, Clarisse Weischedel, and Max Wardetzky. [Geodesics in Heat: A New Approach to Computing Distance Based on Heat Flow](https://www.google.com/search?q=geodesics+in+heat+a+new+approach+to+computing+distance+based+on+heat+flow), 2013.\n",
    "[^bobenko_2005]: Alexander I. Bobenko and Boris A. Springborn. [A discrete Laplace-Beltrami operator for simplicial surfaces](https://www.google.com/search?q=a+discrete+laplace-beltrami+operator+for+simplicial+surfaces), 2005.\n",
    "[^jiang_2017]: Zhongshi Jiang, Scott Schaefer, Daniele Panozzo. [SCAF: Simplicial Complex Augmentation Framework for Bijective Maps](https://doi.org/10.1145/3130800.3130895), 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}