{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Icosahedron matching\n\nCredit: C Ambroise\n\nA simple example on how to match two icosahedrons of the same order.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nimport warnings\nimport numpy as np\nfrom scipy.spatial import transform\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import axes3d\nfrom surfify.plotting import plot_trisurf\nfrom surfify.utils import icosahedron, ico2ico"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first build the reference icosahedron.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "order = 3\nvertices_norm, triangles_norm = icosahedron(order, standard_ico=True)\nprint(vertices_norm.shape, triangles_norm.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then we fetch freesurfer's icosahedron of the same order.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vertices, triangles = icosahedron(order)\nprint(vertices.shape, triangles.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We try to find the optimal rotation between the two icosahedrons using\nthe scipy module.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rotation, rmse = transform.Rotation.align_vectors(vertices_norm, vertices)\nprint(rmse)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Okay, does not seem to be working, because the rmse is supposed to be very\nclose or equal to zero\n\nWe print the vertices to try to find the issue here: it seems that the order\nof the vertices is not the same in the two matrices. That is why the previous\nalgorithm did not work properly, since it can only match to the corresponding\nrow in the other matrix.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(vertices_norm)\nprint(vertices)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we plot the sufaces together to show that they have the same structure\nbut their vertices are not at the same places.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1, 1, subplot_kw={\n        \"projection\": \"3d\", \"aspect\": \"auto\"}, figsize=(10, 10))\nplot_trisurf(vertices_norm, triangles_norm, fig=fig, ax=ax, alpha=0.3,\n             edgecolors=\"blue\")\nplot_trisurf(vertices, triangles, fig=fig, ax=ax, alpha=0.3,\n             edgecolors=\"green\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To compute the rotation between the two structures, we do not need all the\nvertices. So we consider a small subset of 4 points that have correspondances\nin both icosahedrons (for instance they have the same absolute\nvalues, only sign differs, for each dimension).\nThe 4 firsts work for the reference icosahedron.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vertices_of_interest_norm = vertices_norm[:4]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we search for 4 similar vertices in the FreeSurfer icosahedron.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for i in range(len(vertices)):\n    coords_of_interest = vertices[i]\n    idx_of_interest = (np.abs(vertices) == np.abs(coords_of_interest)).all(1)\n    if idx_of_interest.sum() == 4:\n        vertices_of_interest = vertices[idx_of_interest]\n        fs_row_idx = i\n        break\nprint(fs_row_idx)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we need to find a rotation between these two set of points.\nMany can be possible, depending on the ordering of the points. To do this,\nwe compute the optimal rotation matrix between our reference points and\nthe others variously permuted, until we find a rotation that works.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import itertools\npermutations = itertools.permutations(range(4))\nn_permutations = np.math.factorial(4)\nit = 0\nbest_rmse = rmse\nbest_rotation = rotation\nwhile rmse > 0 and it < n_permutations:\n    it += 1\n    order = np.array(next(permutations))\n    with warnings.catch_warnings():\n        warnings.simplefilter(\"ignore\", category=UserWarning)\n        rotation, rmse = transform.Rotation.align_vectors(\n                vertices_of_interest_norm, vertices_of_interest[order])\n    if rmse < best_rmse:\n        best_rmse = rmse\n        best_rotation = rotation\nprint(\"Number of permutations tested {}/{}\".format(it, np.math.factorial(4)))\nprint(best_rotation.as_matrix())\nprint(best_rmse)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we found a rotation that works, we can simply apply it to the icosahedron\nso it matches the reference one.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1, 1, subplot_kw={\n        \"projection\": \"3d\", \"aspect\": \"auto\"}, figsize=(10, 10))\nplot_trisurf(vertices_norm, triangles_norm, fig=fig, ax=ax, alpha=0.3,\n             edgecolors=\"blue\")\nplot_trisurf(best_rotation.apply(vertices), triangles, fig=fig, ax=ax,\n             alpha=0.3, edgecolors=\"green\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To easily solve the issue outlined in this example, you can find a\nfunction in the `surfify.utils` module. `ico2ico` allows you to find a\nproper rotation between a reference icosahedron and another one.\nWe plot only half of the triangles of each icosahedron so it clearly appears\nthat they are the same.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rotation = ico2ico(vertices, vertices_norm)\nfig, ax = plt.subplots(1, 1, subplot_kw={\n        \"projection\": \"3d\", \"aspect\": \"auto\"}, figsize=(10, 10))\nplot_trisurf(vertices_norm, triangles_norm[::2], fig=fig, ax=ax, alpha=0.3,\n             edgecolors=\"blue\")\nplot_trisurf(rotation.apply(vertices), triangles[::2], fig=fig, ax=ax,\n             alpha=0.3, edgecolors=\"green\")\nplt.show()"
      ]
    }
  ],
  "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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}