{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Reconstruction with the Sparse Fascicle Model\n\nIn this example, we will use the Sparse Fascicle Model (SFM) [Rokem2015]_, to\nreconstruct the fiber Orientation Distribution Function (fODF) in every voxel.\n\nFirst, we import the modules we will use in this example:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import dipy.reconst.sfm as sfm\nimport dipy.data as dpd\nimport dipy.direction.peaks as dpp\nfrom dipy.io.image import load_nifti\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.core.gradients import gradient_table\nfrom dipy.viz import window, actor"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the purpose of this example, we will use the Stanford HARDI dataset (150\ndirections, single b-value of 2000 $s/mm^2$) that can be automatically\ndownloaded. If you have not yet downloaded this data-set in one of the other\nexamples, you will need to be connected to the internet the first time you run\nthis example. The data will be stored for subsequent runs, and for use with\nother examples.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "hardi_fname, hardi_bval_fname, hardi_bvec_fname = dpd.get_fnames('stanford_hardi')\ndata, affine = load_nifti(hardi_fname)\n\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\n# Enables/disables interactive visualization\ninteractive = False"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reconstruction of the fiber ODF in each voxel guides subsequent tracking\nsteps. Here, the model is the Sparse Fascicle Model, described in\n[Rokem2014]_. This model reconstructs the diffusion signal as a combination of\nthe signals from different fascicles. This model can be written as:\n\n\\begin{align}y = X\\beta\\end{align}\n\nWhere $y$ is the signal and $\\beta$ are weights on different points in the\nsphere. The columns of the design matrix, $X$ are the signals in each point in\nthe measurement that would be predicted if there was a fascicle oriented in the\ndirection represented by that column. Typically, the signal used for this\nkernel will be a prolate tensor with axial diffusivity 3-5 times higher than\nits radial diffusivity. The exact numbers can also be estimated from examining\nparts of the brain in which there is known to be only one fascicle (e.g. in\ncorpus callosum).\n\nSparsity constraints on the fiber ODF ($\\beta$) are set through the Elastic Net\nalgorithm [Zou2005]_.\n\nElastic Net optimizes the following cost function:\n\n\\begin{align}\\sum_{i=1}^{n}{(y_i - \\hat{y}_i)^2} + \\alpha (\\lambda \\sum_{j=1}^{m}{w_j}+(1-\\lambda) \\sum_{j=1}^{m}{w^2_j}\\end{align}\n\nwhere $\\hat{y}$ is the signal predicted for a particular setting of $\\beta$,\nsuch that the left part of this expression is the squared loss function;\n$\\alpha$ is a parameter that sets the balance between the squared loss on\nthe data, and the regularization constraints. The regularization parameter\n$\\lambda$ sets the `l1_ratio`, which controls the balance between L1-sparsity\n(low sum of weights), and low L2-sparsity (low sum-of-squares of the weights).\n\nJust like Constrained Spherical Deconvolution (see `reconst-csd`), the SFM\nrequires the definition of a response function. We'll take advantage of the\nautomated algorithm in the :mod:`csdeconv` module to find this response\nfunction:\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.reconst.csdeconv import auto_response_ssst\nresponse, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The ``response`` return value contains two entries. The first is an array with\nthe eigenvalues of the response function and the second is the average S0 for\nthis response.\n\nIt is a very good practice to always validate the result of\n``auto_response_ssst``. For, this purpose we can print it and have a look\nat its values.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(response)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "(array([ 0.0014,  0.00029,  0.00029]), 416.206)\n\nWe initialize an SFM model object, using these values. We will use the default\nsphere (362  vertices, symmetrically distributed on the surface of the sphere),\nas a set of putative fascicle directions that are considered in the model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sphere = dpd.get_sphere()\nsf_model = sfm.SparseFascicleModel(gtab, sphere=sphere,\n                                   l1_ratio=0.5, alpha=0.001,\n                                   response=response[0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the purpose of the example, we will consider a small volume of data\ncontaining parts of the corpus callosum and of the centrum semiovale\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_small = data[20:50, 55:85, 38:39]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Fitting the model to this small volume of data, we calculate the ODF of this\nmodel on the sphere, and plot it.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sf_fit = sf_model.fit(data_small)\nsf_odf = sf_fit.odf(sphere)\n\nfodf_spheres = actor.odf_slicer(sf_odf, sphere=sphere, scale=0.8,\n                                colormap='plasma')\n\nscene = window.Scene()\nscene.add(fodf_spheres)\n\nprint('Saving illustration as sf_odfs.png')\nwindow.record(scene, out_path='sf_odfs.png', size=(1000, 1000))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can extract the peaks from the ODF, and plot these as well\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sf_peaks = dpp.peaks_from_model(sf_model,\n                                data_small,\n                                sphere,\n                                relative_peak_threshold=.5,\n                                min_separation_angle=25,\n                                return_sh=False)\n\n\nscene.clear()\nfodf_peaks = actor.peak_slicer(sf_peaks.peak_dirs, sf_peaks.peak_values)\nscene.add(fodf_peaks)\n\nprint('Saving illustration as sf_peaks.png')\nwindow.record(scene, out_path='sf_peaks.png', size=(1000, 1000))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we plot both the peaks and the ODFs, overlaid:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fodf_spheres.GetProperty().SetOpacity(0.4)\nscene.add(fodf_spheres)\n\nprint('Saving illustration as sf_both.png')\nwindow.record(scene, out_path='sf_both.png', size=(1000, 1000))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: sf_both.png\n   :align: center\n\n   SFM Peaks and ODFs.\n\nTo see how to use this information in tracking, proceed to `sfm-track`.\n\n## References\n\n.. [Rokem2015] Ariel Rokem, Jason D. Yeatman, Franco Pestilli, Kendrick\n   N. Kay, Aviv Mezer, Stefan van der Walt, Brian A. Wandell\n   (2015). Evaluating the accuracy of diffusion MRI models in white\n   matter. PLoS ONE 10(4): e0123272. doi:10.1371/journal.pone.0123272\n\n.. [Zou2005] Zou H, Hastie T (2005). Regularization and variable\n   selection via the elastic net. J R Stat Soc B:301-320\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
}