{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Linear fascicle evaluation (LiFE)\n\nEvaluating the results of tractography algorithms is one of the biggest\nchallenges for diffusion MRI. One proposal for evaluation of tractography\nresults is to use a forward model that predicts the signal from each of a set\nof streamlines, and then fit a linear model to these simultaneous predictions\n[Pestilli2014]_.\n\nWe will use streamlines generated using probabilistic tracking on CSA\npeaks. For brevity, we will include in this example only streamlines going\nthrough the corpus callosum connecting left to right superior frontal\ncortex. The process of tracking and finding these streamlines is fully\ndemonstrated in the `streamline_tools` example. If this example has been\nrun, we can read the streamlines from file. Otherwise, we'll run that example\nfirst, by importing it. This provides us with all of the variables that were\ncreated in that example:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from os.path import join as pjoin\n\nimport numpy as np\nimport matplotlib\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.axes_grid1 import AxesGrid\n\nimport dipy.core.optimize as opt\nfrom dipy.data import fetch_stanford_tracks\nfrom dipy.io.streamline import load_trk\nimport dipy.tracking.life as life\nfrom dipy.viz import window, actor, colormap as cmap\n\n    # We'll need to know where the corpus callosum is from these variables:\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti_data, load_nifti\n\nhardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\nlabel_fname = get_fnames('stanford_labels')\nt1_fname = get_fnames('stanford_t1')\n\ndata, affine, hardi_img = load_nifti(hardi_fname, return_img=True)\nlabels = load_nifti_data(label_fname)\nt1_data = load_nifti_data(t1_fname)\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\ncc_slice = labels == 2\n\n# Let's now fetch a set of streamlines from the Stanford HARDI dataset.\n# Those streamlines were generated during the :ref:`streamline_tools` example.\n# Read the candidates from file in voxel space:\n\nstreamlines_files = fetch_stanford_tracks()\nlr_superiorfrontal_path = pjoin(streamlines_files[1],\n                                'hardi-lr-superiorfrontal.trk')\n\ncandidate_sl_sft = load_trk(lr_superiorfrontal_path, 'same')\ncandidate_sl_sft.to_vox()\ncandidate_sl = candidate_sl_sft.streamlines"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The streamlines that are entered into the model are termed 'candidate\nstreamlines' (or a 'candidate connectome'):\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's visualize the initial candidate group of streamlines in 3D, relative to\nthe anatomical structure of this brain:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Enables/disables interactive visualization\ninteractive = False\n\ncandidate_streamlines_actor = actor.streamtube(candidate_sl,\n                                               cmap.line_colors(candidate_sl))\ncc_ROI_actor = actor.contour_from_roi(cc_slice, color=(1., 1., 0.),\n                                      opacity=0.5)\n\nvol_actor = actor.slicer(t1_data)\n\nvol_actor.display(x=40)\nvol_actor2 = vol_actor.copy()\nvol_actor2.display(z=35)\n\n# Add display objects to canvas\nscene = window.Scene()\nscene.add(candidate_streamlines_actor)\nscene.add(cc_ROI_actor)\nscene.add(vol_actor)\nscene.add(vol_actor2)\nwindow.record(scene, n_frames=1,\n              out_path='life_candidates.png',\n              size=(800, 800))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: life_candidates.png\n   :align: center\n\n   **Candidate connectome before life optimization**\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we initialize a LiFE model. We import the ``dipy.tracking.life`` module,\nwhich contains the classes and functions that implement the model:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fiber_model = life.FiberModel(gtab)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since we read the streamlines from a file, already in the voxel space, we do\nnot need to transform them into this space. Otherwise, if the streamline\ncoordinates were in the world space (relative to the scanner iso-center, or\nrelative to the mid-point of the AC-PC-connecting line), we would use this::\n\n   inv_affine = np.linalg.inv(hardi_img.affine)\n\nthe inverse transformation from world space to the voxel space as the affine\nfor the following model fit.\n\nThe next step is to fit the model, producing a ``FiberFit`` class instance,\nthat stores the data, as well as the results of the fitting procedure.\n\nThe LiFE model posits that the signal in the diffusion MRI volume can be\nexplained by the streamlines, by the equation\n\n\\begin{align}y = X\\beta\\end{align}\n\nWhere $y$ is the diffusion MRI signal, $\\beta$ are a set of weights on the\nstreamlines and $X$ is a design matrix. This matrix has the dimensions $m$ by\n$n$, where $m=n_{voxels} \\cdot n_{directions}$, and $n_{voxels}$ is the set of\nvoxels in the ROI that contains the streamlines considered in this model. The\n$i^{th}$ column of the matrix contains the expected contributions of the\n$i^{th}$ streamline (arbitrarily ordered) to each of the voxels. $X$ is a\nsparse matrix, because each streamline traverses only a small percentage of the\nvoxels. The  expected contributions of the streamline are calculated using a\nforward model, where each node of the streamline is modeled as a cylindrical\nfiber compartment with Gaussian diffusion, using the diffusion tensor model.\nSee [Pestilli2014]_ for more detail on the model, and variations of this model.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fiber_fit = fiber_model.fit(data, candidate_sl, affine=np.eye(4))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The ``FiberFit`` class instance holds various properties of the model fit. For\nexample, it has the weights $\\beta$, that are assigned to each streamline. In\nmost cases, a tractography through some region will include redundant\nstreamlines, and these streamlines will have $\\beta_i$ that are 0.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1)\nax.hist(fiber_fit.beta, bins=100, histtype='step')\nax.set_xlabel('Fiber weights')\nax.set_ylabel('# fibers')\nfig.savefig('beta_histogram.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: beta_histogram.png\n   :align: center\n\n   **LiFE streamline weights**\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use $\\beta$ to filter out these redundant streamlines, and generate an\noptimized group of streamlines:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "optimized_sl = [np.row_stack(candidate_sl)[np.where(fiber_fit.beta > 0)[0]]]\nscene = window.Scene()\nscene.add(actor.streamtube(optimized_sl, cmap.line_colors(optimized_sl)))\nscene.add(cc_ROI_actor)\nscene.add(vol_actor)\nwindow.record(scene, n_frames=1, out_path='life_optimized.png',\n              size=(800, 800))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: life_optimized.png\n   :align: center\n\n   **Streamlines selected via LiFE optimization**\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The new set of streamlines should do well in fitting the data, and redundant\nstreamlines have presumably been removed (in this case, about 50% of the\nstreamlines).\n\nBut how well does the model do in explaining the diffusion data? We can\nquantify that: the ``FiberFit`` class instance has a `predict` method, which\ncan be used to invert the model and predict back either the data that was used\nto fit the model, or other unseen data (e.g. in cross-validation, see\n`kfold_xval`).\n\nWithout arguments, the ``.predict()`` method will predict the diffusion signal\nfor the same gradient table that was used in the fit data, but ``gtab`` and\n``S0`` keyword arguments can be used to predict for other acquisition schemes\nand other baseline non-diffusion-weighted signals.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model_predict = fiber_fit.predict()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will focus on the error in prediction of the diffusion-weighted data, and\ncalculate the root of the mean squared error.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model_error = model_predict - fiber_fit.data\nmodel_rmse = np.sqrt(np.mean(model_error[:, 10:] ** 2, -1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As a baseline against which we can compare, we calculate another error term. In\nthis case, we assume that the weight for each streamline is equal\nto zero. This produces the naive prediction of the mean of the signal in each\nvoxel.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "beta_baseline = np.zeros(fiber_fit.beta.shape[0])\npred_weighted = np.reshape(opt.spdot(fiber_fit.life_matrix, beta_baseline),\n                           (fiber_fit.vox_coords.shape[0],\n                            np.sum(~gtab.b0s_mask)))\nmean_pred = np.empty((fiber_fit.vox_coords.shape[0], gtab.bvals.shape[0]))\nS0 = fiber_fit.b0_signal"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since the fitting is done in the demeaned S/S0 domain, we need\nto add back the mean and then multiply by S0 in every voxel:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mean_pred[..., gtab.b0s_mask] = S0[:, None]\nmean_pred[..., ~gtab.b0s_mask] =\\\n    (pred_weighted + fiber_fit.mean_signal[:, None]) * S0[:, None]\nmean_error = mean_pred - fiber_fit.data\nmean_rmse = np.sqrt(np.mean(mean_error ** 2, -1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we can compare the overall distribution of errors between these two\nalternative models of the ROI. We show the distribution of differences in error\n(improvement through model fitting, relative to the baseline model). Here,\npositive values denote an improvement in error with model fit, relative to\nwithout the model fit.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1)\nax.hist(mean_rmse - model_rmse, bins=100, histtype='step')\nax.text(0.2, 0.9, 'Median RMSE, mean model: %.2f' % np.median(mean_rmse),\n        horizontalalignment='left',\n        verticalalignment='center',\n        transform=ax.transAxes)\nax.text(0.2, 0.8, 'Median RMSE, LiFE: %.2f' % np.median(model_rmse),\n        horizontalalignment='left',\n        verticalalignment='center',\n        transform=ax.transAxes)\nax.set_xlabel('RMS Error')\nax.set_ylabel('# voxels')\nfig.savefig('error_histograms.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: error_histograms.png\n   :align: center\n\n   Improvement in error with fitting of the LiFE model.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Second, we can show the spatial distribution of the two error terms,\nand of the improvement with the model fit:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vol_model = np.ones(data.shape[:3]) * np.nan\nvol_model[fiber_fit.vox_coords[:, 0],\n          fiber_fit.vox_coords[:, 1],\n          fiber_fit.vox_coords[:, 2]] = model_rmse\nvol_mean = np.ones(data.shape[:3]) * np.nan\nvol_mean[fiber_fit.vox_coords[:, 0],\n         fiber_fit.vox_coords[:, 1],\n         fiber_fit.vox_coords[:, 2]] = mean_rmse\nvol_improve = np.ones(data.shape[:3]) * np.nan\nvol_improve[fiber_fit.vox_coords[:, 0],\n            fiber_fit.vox_coords[:, 1],\n            fiber_fit.vox_coords[:, 2]] = mean_rmse - model_rmse\nsl_idx = 49\nfig = plt.figure()\nfig.subplots_adjust(left=0.05, right=0.95)\nax = AxesGrid(fig, 111,\n              nrows_ncols=(1, 3),\n              label_mode=\"1\",\n              share_all=True,\n              cbar_location=\"top\",\n              cbar_mode=\"each\",\n              cbar_size=\"10%\",\n              cbar_pad=\"5%\")\nax[0].matshow(np.rot90(t1_data[sl_idx, :, :]), cmap=matplotlib.cm.bone)\nim = ax[0].matshow(np.rot90(vol_model[sl_idx, :, :]), cmap=matplotlib.cm.hot)\nax.cbar_axes[0].colorbar(im)\nax[1].matshow(np.rot90(t1_data[sl_idx, :, :]), cmap=matplotlib.cm.bone)\nim = ax[1].matshow(np.rot90(vol_mean[sl_idx, :, :]), cmap=matplotlib.cm.hot)\nax.cbar_axes[1].colorbar(im)\nax[2].matshow(np.rot90(t1_data[sl_idx, :, :]), cmap=matplotlib.cm.bone)\nim = ax[2].matshow(np.rot90(vol_improve[sl_idx, :, :]),\n                   cmap=matplotlib.cm.RdBu)\nax.cbar_axes[2].colorbar(im)\nfor lax in ax:\n    lax.set_xticks([])\n    lax.set_yticks([])\nfig.savefig(\"spatial_errors.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: spatial_errors.png\n   :align: center\n\n   Spatial distribution of error and improvement.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This image demonstrates that in many places, fitting the LiFE model results in\nsubstantial reduction of the error.\n\nNote that for full-brain tractographies *LiFE* can require large amounts of\nmemory. For detailed memory profiling of the algorithm, based on the\nstreamlines generated in `example_probabilistic_fiber_tracking`, see [this\nIPython notebook](http://nbviewer.ipython.org/gist/arokem/bc29f34ebc97510d9def).\n\nFor the Matlab implementation of LiFE, head over to [Franco Pestilli's github\nwebpage](http://francopestilli.github.io/life/).\n\n## References\n\n.. [Pestilli2014] Pestilli, F., Yeatman, J, Rokem, A. Kay, K. and Wandell B.A.\n   (2014). Validation and statistical inference in living connectomes. Nature\n   Methods 11: 1058-1063. doi:10.1038/nmeth.3098\n\n.. include:: ../links_names.inc\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
}