{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ShapeDNA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ShapeDNA is an n-dimensional intrinsic shape descriptor (see Reuter et al., CAD Journal, 2006). It can be used to compare two geometric objects independent of their pose or posture as the ShapeDNA is not affected by (near)-isometric deformations. This tutorial shows how you compute, normalize and re-weight Laplace-Beltrami spectra to obtain the ShapeDNA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "from lapy import TetMesh, TriaMesh, shapedna"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we load some data: a tria mesh representing the boundary of a cube and a tetrahedral mesh representing the full cube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--> VTK format         ... \n",
      " --> DONE ( V: 2402 , T: 4800 )\n",
      "\n",
      "--> VTK format         ... \n",
      " --> DONE ( V: 9261 , T: 48000 )\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# load data\n",
    "tria = TriaMesh.read_vtk(\"../data/cubeTria.vtk\")\n",
    "tet = TetMesh.read_vtk(\"../data/cubeTetra.vtk\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the first three eigenvalues and eigenvectors of the triangle mesh..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TriaMesh with regular Laplace-Beltrami\n",
      "Solver: spsolve (LU decomposition) ...\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-4.0165149e-05,  4.1696410e+00,  4.1704664e+00], dtype=float32)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute eigenvalues and eigenvectors for tria mesh\n",
    "ev = shapedna.compute_shapedna(tria, k=3)\n",
    "ev['Eigenvectors']\n",
    "ev['Eigenvalues']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we perform a normalization of the eigenvalues using the method \"geometry\" which is equal to surface area normalization for 2d meshes. The resulting eigenvalues are the same as when computing them on the same shape with unit surface area (=1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.4099089e-04,  2.5017845e+01,  2.5022799e+01], dtype=float32)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# volume / surface / geometry normalization of tria eigenvalues\n",
    "shapedna.normalize_ev(tria, ev[\"Eigenvalues\"], method=\"geometry\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For surfaces, eigenvalues increase linearly with their ordering. In order to reduce the influence of higher (and probably more noise affected) eigenvalues it is common practice to perform a linear re-weighting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-4.01651487e-05,  2.08482051e+00,  1.39015547e+00])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# linear reweighting of tria eigenvalues\n",
    "shapedna.reweight_ev(ev[\"Eigenvalues\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normalized and re-weighted eigenvalues are called the ShapeDNA. We can now compute the distance between two shapes by comparing their ShapeDNA. The default is the Euclidean distance between two ShapeDNA vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute distance for tria eigenvalues (trivial case)\n",
    "shapedna.compute_distance(ev[\"Eigenvalues\"], ev[\"Eigenvalues\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, that usually more eigenvalues are used (in the order of 15 to 50) for shape comparison. Also you can do other analyses, e.g. find clusters in this shape space or project it via PCA for visualization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now repeat the above steps for a tetrahedral mesh, again computing the first three eigenvalues and -vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TetMesh with regular Laplace\n",
      "Solver: spsolve (LU decomposition) ...\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([8.4440224e-05, 9.8897915e+00, 9.8898811e+00], dtype=float32)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute eigenvalues and eigenvectors for tet mesh\n",
    "evTet = shapedna.compute_shapedna(tet, k=3)\n",
    "evTet['Eigenvectors']\n",
    "evTet['Eigenvalues']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For 3d meshes the \"geometry\" normalization defaults to unit volume normalization. Since the cube is already unit volume, nothing happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 4800 triangles on boundary.\n",
      "Searched mesh after 79 flooding iterations (0.012834310531616211 sec).\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([8.4440224e-05, 9.8897915e+00, 9.8898811e+00], dtype=float32)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# volume / surface / geometry normalization of tet eigenvalues\n",
    "shapedna.normalize_ev(tet, evTet[\"Eigenvalues\"], method=\"geometry\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we perform linear re-weighting. This is only meaningful for small eigenvalues as the asymtotic trend of eigenvalues of 3d solids is not linear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([8.44402239e-05, 4.94489574e+00, 3.29662704e+00])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# linear reweighting of tet eigenvalues\n",
    "shapedna.reweight_ev(evTet[\"Eigenvalues\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the ShapeDNA of the 3D solid cube, we can compare it to other ShapeDNA (or to itself, which of course yields zero)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute distance for tria eigenvalues (trivial case)\n",
    "shapedna.compute_distance(evTet[\"Eigenvalues\"], evTet[\"Eigenvalues\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}