{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Streamline length and size reduction\n\nThis example shows how to calculate the lengths of a set of streamlines and\nalso how to compress the streamlines without considerably reducing their\nlengths or overall shape.\n\nA streamline in DIPY_ is represented as a numpy array of size\n$(N       imes 3)$ where each row of the array represents a 3D point of the\nstreamline. A set of streamlines is represented with a list of\nnumpy arrays of size $(N_i        imes 3)$ for $i=1:M$ where $M$ is the\nnumber of streamlines in the set.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom dipy.tracking.distances import approx_polygon_track\nfrom dipy.tracking.streamline import set_number_of_points\nfrom dipy.tracking.utils import length\nimport matplotlib.pyplot as plt\nfrom dipy.viz import window, actor"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's first create a simple simulation of a bundle of streamlines using\na cosine function.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def simulated_bundles(no_streamlines=50, n_pts=100):\n    t = np.linspace(-10, 10, n_pts)\n\n    bundle = []\n    for i in np.linspace(3, 5, no_streamlines):\n        pts = np.vstack((np.cos(2 * t/np.pi), np.zeros(t.shape) + i, t )).T\n        bundle.append(pts)\n\n    start = np.random.randint(10, 30, no_streamlines)\n    end = np.random.randint(60, 100, no_streamlines)\n\n    bundle = [10 * streamline[start[i]:end[i]]\n              for (i, streamline) in enumerate(bundle)]\n    bundle = [np.ascontiguousarray(streamline) for streamline in bundle]\n\n    return bundle\n\n\nbundle = simulated_bundles()\n\nprint('This bundle has %d streamlines' % len(bundle))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This bundle has 50 streamlines.\n\nUsing the ``length`` function we can retrieve the lengths of each streamline.\nBelow we show the histogram of the lengths of the streamlines.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lengths = list(length(bundle))\n\nfig_hist, ax = plt.subplots(1)\nax.hist(lengths, color='burlywood')\nax.set_xlabel('Length')\nax.set_ylabel('Count')\n# plt.show()\nplt.legend()\nplt.savefig('length_histogram.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: length_histogram.png\n   :align: center\n\n   **Histogram of lengths of the streamlines**\n\n``Length`` will return the length in the units of the coordinate system that\nstreamlines are currently. So, if the streamlines are in world coordinates then\nthe lengths will be in millimeters (mm). If the streamlines are for example in\nnative image coordinates of voxel size 2mm isotropic then you will need to\nmultiply the lengths by 2 if you want them to correspond to mm. In this example\nwe process simulated data without units, however this information is good to have\nin mind when you calculate lengths with real data.\n\nNext, let's find the number of points that each streamline has.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_pts = [len(streamline) for streamline in bundle]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Often, streamlines are represented with more points than what is actually\nnecessary for specific applications. Also, sometimes every streamline has a\ndifferent number of points, which could be a problem for some algorithms.\nThe function ``set_number_of_points`` can be used to set the number of points\nof a streamline at a specific number and at the same time enforce that all the\nsegments of the streamline will have equal length.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bundle_downsampled = set_number_of_points(bundle, 12)\nn_pts_ds = [len(s) for s in bundle_downsampled]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Alternatively, the function ``approx_polygon_track`` allows reducing the number\nof points so that there are more points in curvy regions and less points in\nless curvy regions. In contrast with ``set_number_of_points`` it does not\nenforce that segments should be of equal size.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bundle_downsampled2 = [approx_polygon_track(s, 0.25) for s in bundle]\nn_pts_ds2 = [len(streamline) for streamline in bundle_downsampled2]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Both, ``set_number_of_points`` and ``approx_polygon_track`` can be thought as\nmethods for lossy compression of streamlines.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Enables/disables interactive visualization\ninteractive = False\n\nscene = window.Scene()\nscene.SetBackground(*window.colors.white)\nbundle_actor = actor.streamtube(bundle, window.colors.red, linewidth=0.3)\n\nscene.add(bundle_actor)\n\nbundle_actor2 = actor.streamtube(bundle_downsampled, window.colors.red, linewidth=0.3)\nbundle_actor2.SetPosition(0, 40, 0)\n\nbundle_actor3 = actor.streamtube(bundle_downsampled2, window.colors.red, linewidth=0.3)\nbundle_actor3.SetPosition(0, 80, 0)\n\nscene.add(bundle_actor2)\nscene.add(bundle_actor3)\n\nscene.set_camera(position=(0, 0, 0), focal_point=(30, 0, 0))\nwindow.record(scene, out_path='simulated_cosine_bundle.png', size=(900, 900))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: simulated_cosine_bundle.png\n   :align: center\n\n   Initial bundle (down), downsampled at 12 equidistant points (middle),\n   downsampled with points that are not equidistant (up).\n\nFrom the figure above we can see that all 3 bundles look quite similar. However,\nwhen we plot the histogram of the number of points used for each streamline, it\nbecomes obvious that we have managed to reduce in a great amount the size of the\ninitial dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig_hist, ax = plt.subplots(1)\nax.hist(n_pts, color='r', histtype='step', label='initial')\nax.hist(n_pts_ds, color='g', histtype='step', label='set_number_of_points (12)')\nax.hist(n_pts_ds2, color='b', histtype='step', label='approx_polygon_track (0.25)')\nax.set_xlabel('Number of points')\nax.set_ylabel('Count')\n\n# plt.show()\nplt.legend()\nplt.savefig('n_pts_histogram.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: n_pts_histogram.png\n   :align: center\n\n   Histogram of the number of points of the streamlines.\n\nFinally, we can also show that the lengths of the streamlines haven't changed\nconsiderably after applying the two methods of downsampling.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lengths_downsampled = list(length(bundle_downsampled))\nlengths_downsampled2 = list(length(bundle_downsampled2))\n\nfig, ax = plt.subplots(1)\nax.plot(lengths, color='r', label='initial')\nax.plot(lengths_downsampled, color='g', label='set_number_of_points (12)')\nax.plot(lengths_downsampled2, color='b', label='approx_polygon_track (0.25)')\nax.set_xlabel('Streamline ID')\nax.set_ylabel('Length')\n\n# plt.show()\nplt.legend()\nplt.savefig('lengths_plots.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: lengths_plots.png\n   :align: center\n\n   Lengths of each streamline for every one of the 3 bundles.\n\n.. include:: ../links_names.inc\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
}