{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Direct Bundle Registration\n\nThis example explains how you can register two bundles from two different\nsubjects directly in the space of streamlines [Garyfallidis15]_, [Garyfallidis14]_.\n\nTo show the concept we will use two pre-saved cingulum bundles. The algorithm\nused here is called Streamline-based Linear Registration (SLR) [Garyfallidis15]_.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.viz import window, actor\nfrom time import sleep\nfrom dipy.data import two_cingulum_bundles\n\ncb_subj1, cb_subj2 = two_cingulum_bundles()\n\nfrom dipy.align.streamlinear import StreamlineLinearRegistration\nfrom dipy.tracking.streamline import set_number_of_points"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "An important step before running the registration is to resample the\nstreamlines so that they both have the same number of points per streamline.\nHere we will use 20 points. This step is not optional. Inputting streamlines\nwith a different number of points will break the theoretical advantages of using\nthe SLR as explained in [Garyfallidis15]_.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cb_subj1 = set_number_of_points(cb_subj1, 20)\ncb_subj2 = set_number_of_points(cb_subj2, 20)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's say now that we want to move the ``cb_subj2`` (moving) so that it can be\naligned with ``cb_subj1`` (static). Here is how this is done.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "srr = StreamlineLinearRegistration()\n\nsrm = srr.optimize(static=cb_subj1, moving=cb_subj2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After the optimization is finished we can apply the transformation to\n``cb_subj2``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cb_subj2_aligned = srm.transform(cb_subj2)\n\n\ndef show_both_bundles(bundles, colors=None, show=True, fname=None):\n\n    scene = window.Scene()\n    scene.SetBackground(1., 1, 1)\n    for (i, bundle) in enumerate(bundles):\n        color = colors[i]\n        lines_actor = actor.streamtube(bundle, color, linewidth=0.3)\n        lines_actor.RotateX(-90)\n        lines_actor.RotateZ(90)\n        scene.add(lines_actor)\n    if show:\n        window.show(scene)\n    if fname is not None:\n        sleep(1)\n        window.record(scene, n_frames=1, out_path=fname, size=(900, 900))\n\n\nshow_both_bundles([cb_subj1, cb_subj2],\n                  colors=[window.colors.orange, window.colors.red],\n                  show=False,\n                  fname='before_registration.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: before_registration.png\n   :align: center\n\n   Before bundle registration.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "show_both_bundles([cb_subj1, cb_subj2_aligned],\n                  colors=[window.colors.orange, window.colors.red],\n                  show=False,\n                  fname='after_registration.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: after_registration.png\n   :align: center\n\n   After bundle registration.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As you can see the two cingulum bundles are well aligned although they contain\nmany streamlines of different lengths and shapes.\n\nStreamline-based Linear Registration (SLR) is a method which given two sets of\nstreamlines (fixed and moving) and a streamline-based cost function, will\nminimize the cost function and transform the moving set of streamlines (target)\nto the fixed (reference), so that they maximally overlap under the condition\nthat the transform stays linear.\n\nWe denote a single streamline with s and a set of streamlines with S.\nA streamline s is an ordered sequence of line segments connecting 3D vector\npoints $\\mathbf{x}_{k} \\in \\mathbb{R}^{3}$ with $k \\in[1, K]$ where K is the\ntotal number of points of streamline s. Given two bundles(two sets of\nstreamlines), we denote $S_{a}=\\left\\{s_{1}^{a}, \\ldots, S_{A}^{a}\\right\\}$ and\n$S_{b}=\\left\\{s_{1}^{b}, \\ldots, s_{B}^{b}\\right\\}$, where A and B are the\ntotal numbers of streamlines in each set respectively. We want to minimize a\ncost function so that we can align the two sets together. For this purpose, we\nintroduce a new cost function, the Bundle-based Minimum Distance (BMD), which\nis defined as:\n\n\\begin{align}\\operatorname{BMD}\\left(S_{a}, S_{b}\\right)=\\frac{1}{4}\\left(\\frac{1}{A} \\\\\n        \\sum_{i=1}^{A} \\min _{j} D(i, j)+\\frac{1}{B} \\sum_{j=1}^{B} \\\\\n        \\min _{i} D(i, j)\\right)^{2}\\end{align}\n\n\nwhere D is the rectangular matrix given by all pairwise Minimum average\nDirect-Flip (MDF) streamline distances (Garyfallidis et al., 2012). Therefore,\nevery element of matrix D is equal to\n$D_{i j}=M D F\\left(s^{a}{ }_{i}, s^{b}{ }_{j}\\right)$.\n\nNotice, how in Eq. (1), the most similar streamlines from one streamline set\nto the other are weighted more by averaging the minimum values of the rows\nand columns of matrix D. This makes our method robust to fanning streamlines\nnear endpoints of bundles and spurious streamlines if any in the bundle. The\nMDF is a symmetric distance between two individual streamlines. It was\nprimarily used for clustering (Garyfallidis et al., 2010; Visser et al., 2011)\nand tractography simplification (see Garyfallidis et al., 2012). This distance\ncan be applied only when both streamlines have the same number of points.\nTherefore we assume from now on that an initial interpolation of streamlines\nhas been applied, so that all streamlines have the same number of points K,\nand all segments of each streamline have equal length. The length of each\nsegment is equal to the length of the streamline divided by the number of\nsegments $(K-1)$. This is achieved by a simple linear interpolation with the\nstarting and ending points of the streamlines intact. When K is small, the\ninterpolation provides a rough representation of the streamline, but as K\nbecomes larger and larger the shape of the interpolated streamline becomes\nidentical with the shape of the initial streamline. Under this assumption, the\nMDF for two streamlines $S_{a}$ and $S_{b}$ is defined as:\n\n\n\\begin{align}\\operatorname{MDF}\\left(s_{i}^{a}, s_{j}^{b}\\right)=\\min \\\\\n        \\left(d_{\\text {direct }}\\left(s_{i}^{a}, s_{j}^{b}\\right), \\\\\n        d_{\\text {flipped }}\\left(s_{i}^{a}, s_{j}^{b}\\right)\\right)\\end{align}\n\n\nwhere $d_{\\text {direct }}$ is the direct distance which is defined as:\n\n\\begin{align}d_{\\text {direct }}\\left(s_{i}^{a}, s_{j}^{b}\\right)=\\frac{1}{K} \\\\\n        \\sum_{k=1}^{K}\\left\\|\\mathbf{x}_{k}^{a}-\\mathbf{x}_{k}^{b}\\right\\|_{2}\\end{align}\n\nwhere $x_{k}^{a}$ is the k-th point of streamline $S_{i}^{a}$ and $x_{k}^{b}$\nis the k-th point of streamline $S_{j}^{b}$. $d_{\\text {flipped }}$ is the one\nof the streamlines flipped and it is defined as:\n\n\\begin{align}d_{\\text {flipped }}\\left(s_{i}^{a}, s_{j}^{b}\\right)=\\frac{1}{K} \\\\\n        \\sum_{k=1}^{K}\\left\\|\\mathbf{x}_{k}^{a}-\\mathbf{x}_{K-k+1}^{b}\\\\\n        \\right\\|_{2}\\end{align}\n\nand K is the total number of points in $x^{a}$ and $x^{b}$.\nThe MDF has two very useful properties. First, it takes into consideration that\nstreamlines have no preferred orientation. Second, it is a mathematically sound\nmetric distance in the space of streamlines as proved in Garyfallidis et al.\n(2012). This means that the MDF is nonnegative, 0 only when both streamlines\nare identical, symmetric and it satisfies the triangle inequality. Now that we\nhave defined our cost function in Eq. (1) we can formulate the following\noptimization problem. Given a fixed bundle S and a moving bundle M we would\nlike to find the vector of parameters t which transforms M to S using a linear\ntransformation T so that BMD is minimum:\n\n\\begin{align}\\operatorname{SLR}(S, M)=\\\\underset{\\mathbf{t}}{\\operatorname{argmin}} \\\\\n        \\operatorname{BMD}(S, T(M, \\mathbf{t}))\\end{align}\n\n\nHere, $\\mathbf{t}$ is a vector in $\\mathbb{R}^{n}$ holding the parameters of\nthe linear transform where n = 12 for affine or n = 6 for rigid registration.\nFrom this vector we can then compose the transformation matrix which is applied\nto all the points of bundle M.\n\n\n## References\n\n.. [Garyfallidis15] Garyfallidis et al., \"Robust and efficient linear\n                    registration of white-matter fascicles in the space\n                    of streamlines\", Neuroimage, 117:124-140, 2015.\n.. [Garyfallidis14] Garyfallidis et al., \"Direct native-space fiber bundle\n                    alignment for group comparisons\", ISMRM, 2014.\n\n\n\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.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}