example_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1 - Imports and defining functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from pyFM.mesh import TriMesh\n",
    "from pyFM.functional import FunctionalMapping\n",
    "\n",
    "import meshplot as mp\n",
    "\n",
    "def plot_mesh(myMesh,cmap=None):\n",
    "    mp.plot(myMesh.vertlist, myMesh.facelist,c=cmap)\n",
    "    \n",
    "def double_plot(myMesh1,myMesh2,cmap1=None,cmap2=None):\n",
    "    d = mp.subplot(myMesh1.vertlist, myMesh1.facelist, c=cmap1, s=[2, 2, 0])\n",
    "    mp.subplot(myMesh2.vertlist, myMesh2.facelist, c=cmap2, s=[2, 2, 1], data=d)\n",
    "\n",
    "def visu(vertices):\n",
    "    min_coord,max_coord = np.min(vertices,axis=0,keepdims=True),np.max(vertices,axis=0,keepdims=True)\n",
    "    cmap = (vertices-min_coord)/(max_coord-min_coord)\n",
    "    return cmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2- Loading and processing a mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic Mesh methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A TriMesh class can be created from a path (to a .off or a .obj file) or simply an array of vertices and an optional array of faces.\n",
    "\n",
    "The mesh can be centered, area-normalized, rotated or translated when loading.\n",
    "\n",
    "\n",
    "Vertices and faces are stored in the 'vertlist' and 'facelist' attributes. One can also use 'mesh.vertices' and 'mesh.faces' to access them. While these notations can feel non-intuitive they result in clearer functions as it avoids expressions of the form ```mesh.vertices - vertices```.\n",
    "\n",
    "A TriMesh class possess multiple attributes like edges, per-face area, per-vertex area, per-face normals, per-vertex normals, ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh1 = TriMesh('data/cat-00.off', area_normalize=True, center=False)\n",
    "mesh2 = TriMesh(mesh1.vertlist, mesh1.facelist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Attributes are computed on the fly and cached\n",
    "edges = mesh1.edges\n",
    "\n",
    "area = mesh1.area\n",
    "\n",
    "face_areas = mesh1.face_areas\n",
    "vertex_areas = mesh1.vertex_areas\n",
    "face_normals = mesh1.normals\n",
    "\n",
    "# AREA WEIGHTED VERTEX NORMALS\n",
    "vertex_normals_a = mesh1.vertex_normals\n",
    "\n",
    "# UNIFORM WEIGHTED VERTEX NORMALS\n",
    "mesh1.set_vertex_normal_weighting('uniform')\n",
    "vertex_normals_u = mesh1.vertex_normals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geodesics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We propose three versions to compute geodesics :\n",
    "- Heat method - based on [potpourri3d](https://github.com/nmwsharp/potpourri3d) using robust laplacian (recommended)\n",
    "- Heat method - pure python implementation from pyFM (not robust but control on the whole code)\n",
    "- Dijkstra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Geodesic distance from a given index\n",
    "# Set robust to False to obtain result from the Python implementation\n",
    "dists = mesh1.geod_from(1000, robust=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d837b8a9816540bb965fab8486718de5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/7207 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "S1_geod = mesh1.get_geodesic(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laplacian and functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectrum of the LBO can be computed easily.\n",
    "\n",
    "Eigenvalues and eigenvectors are stored in the ```mesh.eigenvalues``` and ```mesh.eigenvectors``` attributes.\n",
    "\n",
    "Gradient and divergence can be computed using the associated methods. Using the ```mesh.project``` and ```mesh.unproject``` functions allows to switch between seeing a function in the LBO basis or on the complete shape.\n",
    "\n",
    "The squared $L^2$ norm and $H^1_0$ norm can be computed via the ```mesh.l2_sqnorm``` and ```mesh.h1_sqnorm``` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing 100 eigenvectors\n",
      "\tDone in 1.08 s\n"
     ]
    }
   ],
   "source": [
    "# By default does not use the intrinsic delaunay Laplacian\n",
    "mesh1.process(k=100, intrinsic=False, verbose=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4efe32749d2d4b52ba2d593af6e6b7e5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.001617…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the third eigenfunction\n",
    "plot_mesh(mesh1, mesh1.eigenvectors[:,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 - Computing the functional map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Loading data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mesh 1 : 7207 vertices, 14410 faces\n",
      "Mesh 2 : 5000 vertices,  9996 faces\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9c9620e3996348aca2cdc95f54c64015",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Output(), Output()))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c3ce767353ed4e6e9a23221795b56edd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Output(), Output()))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mesh1 = TriMesh('data/cat-00.off')\n",
    "mesh2 = TriMesh('data/lion-00.off')\n",
    "print(f'Mesh 1 : {mesh1.n_vertices:4d} vertices, {mesh1.n_faces:5d} faces\\n'\n",
    "      f'Mesh 2 : {mesh2.n_vertices:4d} vertices, {mesh2.n_faces:5d} faces')\n",
    "\n",
    "double_plot(mesh1,mesh2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Computing descriptors**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Computing Laplacian spectrum\n",
      "Computing 200 eigenvectors\n",
      "\tDone in 2.27 s\n",
      "Computing 200 eigenvectors\n",
      "\tDone in 1.54 s\n",
      "\n",
      "Computing descriptors\n",
      "\tNormalizing descriptors\n",
      "\n",
      "\t120 out of 600 possible descriptors kept\n"
     ]
    }
   ],
   "source": [
    "process_params = {\n",
    "    'n_ev': (35,35),  # Number of eigenvalues on source and Target\n",
    "    'landmarks': np.loadtxt('data/landmarks.txt',dtype=int)[:5],  # loading 5 landmarks\n",
    "    'subsample_step': 5,  # In order not to use too many descriptors\n",
    "    'descr_type': 'WKS',  # WKS or HKS\n",
    "}\n",
    "\n",
    "model = FunctionalMapping(mesh1,mesh2)\n",
    "model.preprocess(**process_params,verbose=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Fitting the model**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\RR}{\\mathbb{R}}$\n",
    "$\\newcommand{\\Ss}{\\mathcal{S}}$\n",
    "$\\newcommand{\\uargmin}[1]{\\underset{#1}{\\text{argmin}}\\;}$\n",
    "$\\newcommand{\\uargmax}[1]{\\underset{#1}{\\text{argmax}}\\;}$\n",
    "$\\def\\*#1{\\mathbf{#1}}$\n",
    "\n",
    "In pyFM, we always consider functional maps $\\*C:\\Ss_1\\to\\Ss_2$ and pointwise maps $T:\\Ss_2\\to\\Ss_1$ going in opposite directions, with $\\*C$ always going from shape 1 to shape 2 !\n",
    "\n",
    "Optimization problem is\n",
    "\\begin{equation}\n",
    "\\uargmin{\\*C\\in\\RR^{k_2\\times k_1}} w_{descr}\\|\\*C\\*A - \\*B\\|^2 + w_{lap}\\|\\*C\\Delta_1 - \\Delta_2\\*C\\|^2 + w_{\\text{d- comm}}\\sum_i \\|\\*C\\Gamma_1^i - \\Gamma_2^i\\*C\\|^2 + w_{\\text{orient}}\\sum_i \\|\\*C\\Lambda_1^i - \\Lambda_2^i\\*C\\|^2\n",
    "\\end{equation}\n",
    "\n",
    "with $\\Gamma_1^i$ and $\\Gamma_2^i$ [multipliative operators](http://www.lix.polytechnique.fr/~maks/papers/fundescEG17.pdf) associated to the $i$-th descriptors, $\\Lambda_1^i$ and $\\Lambda_2^i$ [orientation preserving operators](https://arxiv.org/abs/1806.04455) associated to the $i$-th descriptors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing commutativity operators\n",
      "\tScaling LBO commutativity weight by 3.9e-09\n",
      "\n",
      "Optimization :\n",
      "\t35 Ev on source - 35 Ev on Target\n",
      "\tUsing 120 Descriptors\n",
      "\tHyperparameters :\n",
      "\t\tDescriptors preservation :1.0e+00\n",
      "\t\tDescriptors commutativity :1.0e-01\n",
      "\t\tLaplacian commutativity :1.0e-02\n",
      "\t\tOrientation preservation :0.0e+00\n",
      "\n",
      "\tTask : CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH, funcall : 47, nit : 44, warnflag : 0\n",
      "\tDone in 0.55 seconds\n"
     ]
    }
   ],
   "source": [
    "fit_params = {\n",
    "    'w_descr': 1e0,\n",
    "    'w_lap': 1e-2,\n",
    "    'w_dcomm': 1e-1,\n",
    "    'w_orient': 0\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "model.fit(**fit_params, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Visualizing the associated point to point map**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "p2p_21 = model.get_p2p(n_jobs=1)\n",
    "cmap1 = visu(mesh1.vertlist); cmap2 = cmap1[p2p_21]\n",
    "double_plot(mesh1,mesh2,cmap1,cmap2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4 - Refining the Functional Map\n",
    "```model.FM``` returns the current state of functional map. One can change which one is returned by using ```model.change_FM_type(FM_type)```, as one can see below. \n",
    "\n",
    "**ICP**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.icp_refine(verbose=True)\n",
    "p2p_21_icp = model.get_p2p()\n",
    "cmap1 = visu(mesh1.vertlist); cmap2 = cmap1[p2p_21_icp]\n",
    "double_plot(mesh1,mesh2,cmap1,cmap2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Zoomout**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "model.change_FM_type('classic') # We refine the first computed map, not the icp-refined one\n",
    "model.zoomout_refine(nit=15, step = 1, verbose=True)\n",
    "print(model.FM.shape)\n",
    "p2p_21_zo = model.get_p2p()\n",
    "cmap1 = visu(mesh1.vertlist); cmap2 = cmap1[p2p_21_zo]\n",
    "double_plot(mesh1,mesh2,cmap1,cmap2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluating Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyFM.eval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute geodesic distance matrix on the cat mesh\n",
    "A_geod = mesh1.get_geodesic(verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load an approximate ground truth map\n",
    "gt_p2p = np.loadtxt('data/lion2cat',dtype=int)\n",
    "\n",
    "acc_base = pyFM.eval.accuracy(p2p_21, gt_p2p, A_geod, sqrt_area=mesh1.sqrtarea)\n",
    "\n",
    "acc_icp = pyFM.eval.accuracy(p2p_21_icp, gt_p2p, A_geod, sqrt_area=np.sqrt(mesh1.area))\n",
    "\n",
    "acc_zo = pyFM.eval.accuracy(p2p_21_zo, gt_p2p, A_geod, sqrt_area=np.sqrt(mesh1.area))\n",
    "\n",
    "print(f'Accuracy results\\n'\n",
    "      f'\\tBasic FM : {1e3*acc_base:.2f}\\n'\n",
    "      f'\\tICP refined : {1e3*acc_icp:.2f}\\n'\n",
    "      f'\\tZoomOut refined : {1e3*acc_zo:.2f}\\n')"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}