correspondence.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correspondence\n",
    "\n",
    "Solution to find matching triangles between two meshes\n",
    "from the paper:<br>\n",
    "- <i>Robert W. Sumner, Jovan Popovic.</i> Deformation Transfer for Triangle Meshes. ACM Transactions on Graphics. 23, 3. August 2004.\n",
    "    http://graphics.csail.mit.edu/~sumner/research/deftransfer/\n",
    "\n",
    "\n",
    "\n",
    "<div style=\"display:table;\">\n",
    "<div style=\"display: table-row;\">\n",
    "<div style=\"display: table-cell; min-width:100px;\"><b>Task:</b></div>\n",
    "<div style=\"display: table-cell;\">Find a triangle mapping of a source mesh to a target mesh.</div>\n",
    "</div>\n",
    "<div style=\"display: table-row;\">\n",
    "<div style=\"display: table-cell;min-width:100px;\"><b>Idea:</b></div>\n",
    "<div style=\"display: table-cell;\">Repeated closest point transformations controlled by user-defined markers</div>\n",
    "</div>\n",
    "</div>\n",
    "<br/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorithm:\n",
    "1. Transform the source mesh to match user-defined markers while keeping the identity.\n",
    "2. Search for closest vertices in the target.\n",
    "3. Transform the source mesh again but additionally with the closest points.\n",
    "4. Repeat from step 2.\n",
    "\n",
    "\n",
    "The transformations in each step start from the initial source mesh.\n",
    "It is updated iteratively by the transformation\n",
    "that is created by the least squares solution of cost functions.\n",
    "Instead for searching for a transformation, the costs are solved directly for the transformed mesh.\n",
    "The user-defined markers forces the equation system\n",
    "to solve for the target shape.\n",
    "The smoothness cost $E_S$ and identity cost $E_I$ try to preserve the source mesh.\n",
    "By enforcing the markers these cost functions find a \"best\" fit for the source\n",
    "into the target mesh.\n",
    "The closest point cost $E_C$ is the third cost function which\n",
    "weighting is increase iteratively and so it slowly adapts the target mesh.\n",
    "The closest points are re-computed after each iteration.\n",
    "In the first iteration there are no closest points used (since there are none yet)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Triangle transformation\n",
    "\n",
    "A triangle (or face) without position can be described in the form of span vectors written in columns\n",
    "that are produced by the three vertices $\\vec{a}, \\vec{b}, \\vec{c}$ and a fourth square-rooted normal $\\vec{d}$.\n",
    "The vertex $\\vec{d}$ allows lateral movement of the face $V_i$ when a transformation is applied.\n",
    "\n",
    "$$V_i = \\left[\\vec{b}-\\vec{a}, \\vec{c}-\\vec{a}, \\vec{d}-\\vec{a}\\right]$$\n",
    "\n",
    "The fourth span vector $\\vec{d}$ is the square-root normalized normal of the face, given by:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\vec{d} = \\vec{a} + \\frac{\\left(\\vec{b}-\\vec{a}\\right)\\times\\left(\\vec{c}-\\vec{a}\\right)}\n",
    "    {\\sqrt{\\left|\\left(\\vec{b}-\\vec{a}\\right)\\times\\left(\\vec{c}-\\vec{a}\\right)\\right|}}\n",
    "\\end{equation}\n",
    "\n",
    "A face (or triangle) transformation $T_i$ for a face-index $i$ is then given by:\n",
    "\n",
    "$$\\widetilde{V}_i = T_i\\cdot V_i \\longrightarrow T_i=\\widetilde{V}_i V_i^{-1}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformation Definition\n",
    "\n",
    "The transformation to the target mesh is controlled by a cost function\n",
    "that can be solved directly for the vertices of the transformed mesh including the fourth face normal.\n",
    "It searches for the least square solution for the vertices $x = (\\vec{v}_1, ..., \\vec{v}_n)$, $x \\in \\mathbb{R}^{3 \\times n}$ which contains\n",
    "the vertices of the mesh as well as the artifical fourth face normals. </br>\n",
    "The solution can be expressed as the sum of the cost functions weighted by $w_S$, $w_I$, and $w_C$.\n",
    "\n",
    "\\begin{equation}\n",
    "    \\widetilde{x} = \\underset{\\widetilde{V}}{\\min}\\left( w_S E_S + w_I  E_I + w_C E_C \\right)\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\n",
    "Name | Equation\n",
    ":---- | :-----\n",
    "Smoothness | $$E_S(v_1,...,v_n)=\\sum_{T_i}\\sum_{T_j\\in\\text{adj}(T_i)}\\left\\|T_i - T_j\\right\\|^2_F$$\n",
    "Identity | $$ E_I(v_1,...,v_n) = \\sum_{T_i}\\left\\|T_i-\\mathbf{I}\\right\\|^2_F $$\n",
    "Closest Valid Point | $$E_C(v_1,...,v_n; c_1,...,c_n)=\\sum_{i}\\left\\|v_i-c_i\\right\\|^2$$\n",
    "\n",
    "Each cost function can be expanded and expressed as a sum of individual terms\n",
    "and merged back into a single Frobenius norm (like an accordion).</br>\n",
    "For example the identity cost can be converted as followed:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frobenius Norm\n",
    "#### Merging equations\n",
    "\n",
    "The equations can be merged into a single Frobenius Norm which can be solved by common solvers.\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_{k}\\left\\|Z_k\\right\\|^2_F &= \\sum_{k}z_{k1}^2+z_{k2}^2+...+z_{kl}^2\n",
    "\\\\\n",
    " &= (z_{11}^2+z_{12}^2+...+z_{1l}^2)+ (z_{21}^2+z_{22}^2+...+z_{2l}^2) + ... + (z_{k1}^2+z_{k2}^2+...+z_{kl}^2)\n",
    " \\\\\n",
    " &= \\left\\|\\sum_{k}Z_k\\right\\|^2_F\n",
    "\\end{align}\n",
    "\n",
    "with $l$ element in $Z$.\n",
    "\n",
    "#### Transposed\n",
    "\n",
    "The frobenius norm and its minimization does not change if the equation is transposed.\n",
    "\n",
    "\\begin{equation}\n",
    "\\left\\|Z\\right\\|^2_F = z_{1}^2+z_{2}^2+...+z_{l}^2 = \\left\\|Z_k^T\\right\\|^2_F\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformation Structure\n",
    "\n",
    "A equation system of squared Frobenius norms can be reordered into another structure.\n",
    "This can be used to structure a matrix equation system that can be solved efficiently.\n",
    "\n",
    "\n",
    "\n",
    "### Matrix Structure\n",
    "The triangle transformation $T_i$ can be expressed as a equation on all vertices $x \\in \\mathbb{R}^{3\\times n}$\n",
    "and $\\hat{V}^{-1} \\in \\mathbb{R}^{n\\times 3}$:\n",
    "\n",
    "\\begin{equation}\n",
    "T_i = \\widetilde{V} V^{-1} \\quad \\xrightarrow[]{\\text{Expand to } x} \\quad \\hat{T}_i = x \\hat{V}^{-1}\n",
    "\\end{equation}\n",
    "\n",
    "<!--\n",
    "% \\quad \\xrightarrow[]{Transpose} \\quad \\hat{T}_i^T=(\\hat{V}^{-1})^Tx^T\n",
    "-->\n",
    "\n",
    "This applied to any $T_i$, it follows with $C_i \\in \\mathbb{R}^{3\\times 3}$ as constant:\n",
    "\n",
    "\\begin{align}\n",
    "\\left\\|T_i-C_i\\right\\|^2_F\n",
    "&= \\left\\|\\widetilde{V}_i V_i^{-1}-C_i\\right\\|^2_F\n",
    "& \\xrightarrow[]{\\text{Expand to } x} \\quad& \\left\\|x \\hat{V}_i^{-1}-C_i\\right\\|^2_F\n",
    "\\end{align}\n",
    "\n",
    "Transposing and merging the equations for $m$ triangle transformations $T_i$\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_{i=1}^{m}\\left\\|x\\hat{V}^{-1}-C_i\\right\\|^2_F\n",
    "&= \\sum_{i=1}^{m}\\left\\| \\left( x\\hat{V}^{-1}-C_i \\right)^T \\right\\|^2_F\n",
    "\\\\&=\\sum_{i=1}^{m}\\left\\| \\left(\\hat{V}_i^{-1}\\right)^T x^T - C_i^T  \\right\\|^2_F\n",
    "\\\\&= \\left\\| \\begin{pmatrix}\n",
    "\\left(\\hat{V}_1^{-1}\\right)^T \\\\\n",
    "\\left(\\hat{V}_2^{-1}\\right)^T \\\\\n",
    "... \\\\\n",
    "\\left(\\hat{V}_m^{-1}\\right)^T\n",
    "\\end{pmatrix}\n",
    "x^T - \\begin{pmatrix}\n",
    "C_1^T \\\\\n",
    "C_2^T \\\\\n",
    "... \\\\\n",
    "C_m^T\n",
    "\\end{pmatrix} \\right\\|^2_F\n",
    "&&= \\left\\|Ax^T-b\\right\\|^2_F\n",
    "\\end{align}\n",
    "\n",
    "### Cost Matrices\n",
    "\n",
    "Each cost builds a separate matrix in which each row is minimized with a least square method.\n",
    "\n",
    "#### Identity cost matrix\n",
    "\n",
    "\\begin{equation}\n",
    "E_I(v_1,...,v_n) = \\sum_m^{\\left|T_i\\right|}\\left\\|T_i-\\mathbf{I}\\right\\|^2_F\n",
    "= \\left\\| \\begin{pmatrix}\n",
    "\\left(\\hat{V}_1^{-1}\\right)^T \\\\\n",
    "\\left(\\hat{V}_2^{-1}\\right)^T \\\\\n",
    "... \\\\\n",
    "\\left(\\hat{V}_m^{-1}\\right)^T\n",
    "\\end{pmatrix}\n",
    "x^T - \\begin{pmatrix}\n",
    "\\mathbf{I} \\\\\n",
    "\\mathbf{I} \\\\\n",
    "... \\\\\n",
    "\\mathbf{I}\n",
    "\\end{pmatrix} \\right\\|^2_F\n",
    "= \\left\\|A_Ix^T-b_I\\right\\|^2_F\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "where $A_I \\in \\mathbb{R}^{3m\\times n}$  and $b_i \\in \\mathbb{R}^{3m\\times 3}$\n",
    "\n",
    "#### Smoothness cost matrix\n",
    "\n",
    "\\begin{equation}\n",
    "E_S(v_1,...,v_n)=\\sum_m^{\\left|T_i\\right|}\\sum_{T_j\\in\\text{adj}(T_i)}\\left\\|T_i - T_j\\right\\|^2_F\n",
    "= \\left\\| \\begin{pmatrix}\n",
    "\\left(\\hat{V}_1^{-1} - \\hat{V}_{j_1}^{-1}\\right)^T \\\\\n",
    "\\left(\\hat{V}_1^{-1} - \\hat{V}_{j_2}^{-1}\\right)^T \\\\\n",
    "... \\\\\n",
    "\\left(\\hat{V}_1^{-1} - \\hat{V}_{j_{\\left|\\text{adj}(T_i)\\right|}}^{-1}\\right)^T \\\\\n",
    "\\left(\\hat{V}_2^{-1} - \\hat{V}_{j_1}^{-1}\\right)^T \\\\\n",
    "... \\\\\n",
    "\\end{pmatrix}\n",
    "x^T - 0\n",
    "\\right\\|^2_F\n",
    "= \\left\\|A_Sx^T-b_S\\right\\|^2_F\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "where $A_S \\in \\mathbb{R}^{3q\\times n}$, $b_i \\in \\mathbb{R}^{3q\\times 3}$, and\n",
    "$q = \\sum_m^{\\left|T_i\\right|}\\left|\\text{adj}(T_i)\\right|$.\n",
    "\n",
    "\n",
    "### Closest point cost matrix\n",
    "\n",
    "The closest point is the 2-Norm and can be handled the same way like the Frobenius norm.\n",
    "\n",
    "\\begin{equation}\n",
    "E_C(v_1,...,v_n; c_1,...,c_n)=\\sum_{i}\\left\\|v_i-c_i\\right\\|^2\n",
    "= \\left\\| \\mathbf{I}\n",
    "x^T - \\begin{pmatrix}\n",
    "c_1^T \\\\\n",
    "c_2^T \\\\\n",
    "... \\\\\n",
    "c_n^T\n",
    "\\end{pmatrix}\\right\\|^2_F\n",
    "= \\left\\|A_Cx^T-b_C\\right\\|^2_F\n",
    "\\end{equation}\n",
    "\n",
    "where $A_C \\in \\mathbb{R}^{n\\times n}$  is the identity, and $b_c \\in \\mathbb{R}^{n\\times 3}$ are the transposed\n",
    "closes points.\n",
    "\n",
    "\n",
    "### Final equation for least-squares solver\n",
    "Each cost is weighted, which needs to be reflected into the equation system.\n",
    "We can do this by weighting each row.\n",
    "\n",
    "\\begin{equation}\n",
    "E = w_S E_S + w_I  E_I + w_C E_C\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "\\widetilde{x} = \\underset{x}{\\text{arg}\\min}\\left\\|\n",
    "\\begin{pmatrix}\n",
    "w_I A_I \\\\\n",
    "w_S A_S \\\\\n",
    "w_C A_C \\\\\n",
    "\\end{pmatrix}x^T -\n",
    "\\begin{pmatrix}\n",
    "w_I b_I \\\\\n",
    "w_S b_S \\\\\n",
    "w_C b_C \\\\\n",
    "\\end{pmatrix}\n",
    "\\right\\|^2_F\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Before the least-squares solution is searched, the columns for all marker vertices\n",
    "are resolved and subtracted from the right side $b_I, b_S, b_C$.\n",
    "Then the marker positions on the left side $A_I, A_S, A_C$ are set to 0.\n",
    "This forces the equation to find a solution with the marker vertices fixed.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transformation Vertex Indicies\n",
    "In the expansion of V to $\\hat{V}$ for the big $x$ vertices list each transformation is indexed as followed:\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{T} &= \\widetilde{V}V^{-1}\n",
    "    = \\left[\\vec{b}-\\vec{a}, \\vec{c}-\\vec{a}, \\vec{d}-\\vec{a}\\right] V^{-1}  \\\\\n",
    "    \\\\\n",
    "    &= \\begin{pmatrix}\n",
    "        b_{1} - a_{1} & c_{1} - a_{1} & d_{1} - a_{1} \\\\\n",
    "        b_{2} - a_{2} & c_{2} - a_{2} & d_{2} - a_{2} \\\\\n",
    "        b_{3} - a_{3} & c_{3} - a_{3} & d_{3} - a_{3}\n",
    "    \\end{pmatrix}\n",
    "    \\begin{pmatrix}\n",
    "        \\overline{v}_{1,1} & \\overline{v}_{1,2} & \\overline{v}_{1,3} \\\\\n",
    "        \\overline{v}_{2,1} & \\overline{v}_{2,2} & \\overline{v}_{2,3} \\\\\n",
    "        \\overline{v}_{3,1} & \\overline{v}_{3,2} & \\overline{v}_{3,3}\n",
    "    \\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "The resulting matrix $\\mathbb{R}^{3\\times 3}$ has the following equation for each entry:\n",
    "\n",
    "\\begin{align}\n",
    "    w_{i,j}\n",
    "    &= (b_{i} - a_{i}) \\overline{v}_{1,j}\n",
    "    + (c_{i} - a_{i}) \\overline{v}_{2,j}\n",
    "    + (d_{i} - a_{i}) \\overline{v}_{3,j}\n",
    "    \\\\\n",
    "    &= b_{i}\\overline{v}_{1,j}  - a_{i}\\overline{v}_{1,j}\n",
    "    + c_{i}\\overline{v}_{2,j} - a_{i}\\overline{v}_{2,j}\n",
    "    + d_{i}\\overline{v}_{3,j} - a_{i}\\overline{v}_{3,j}\n",
    "    \\\\\n",
    "    &=\n",
    "    - (\\overline{v}_{1,j} + \\overline{v}_{2,j}+ \\overline{v}_{3,j})a_{i}\n",
    "    +\\overline{v}_{1,j}b_{i}\n",
    "    + \\overline{v}_{2,j}c_{i}\n",
    "    + \\overline{v}_{3,j}d_{i}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implemenation\n",
    "\n",
    "Numpy + Scipy + Plotly (for plotting)\n",
    "\n",
    "### Optimization\n",
    "\n",
    "Since the three XYZ components of the vertices are independent, the minimization problem can be solved independetly\n",
    "by component. One optimization is to create a diagonal block matrix and solve it single-threaded;\n",
    "or use multiple threads to solve the least-square problem for each XYZ component."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}