{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# An introduction to the Probabilistic Direction Getter\n\nProbabilistic fiber tracking is a way of reconstructing white matter\nconnections using diffusion MR imaging. Like deterministic fiber tracking, the\nprobabilistic approach follows the trajectory of a possible pathway step by\nstep starting at a seed, however, unlike deterministic tracking, the tracking\ndirection at each point along the path is chosen at random from a distribution.\nThe distribution at each point is different and depends on the observed\ndiffusion data at that point. The distribution of tracking directions at each\npoint can be represented as a probability mass function (PMF) if the possible\ntracking directions are restricted to discrete numbers of well distributed\npoints on a sphere.\n\nThis example is an extension of the `example_tracking_introduction_eudx`\nexample. We'll begin by repeating a few steps from that example, loading the\ndata and fitting a Constrained Spherical Deconvolution (CSD) model.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from 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, load_nifti_data\nfrom dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,\n                                   auto_response_ssst)\nfrom dipy.tracking import utils\nfrom dipy.tracking.local_tracking import LocalTracking\nfrom dipy.tracking.streamline import Streamlines\nfrom dipy.tracking.stopping_criterion import ThresholdStoppingCriterion\nfrom dipy.viz import window, actor, colormap, has_fury\n\n\n# Enables/disables interactive visualization\ninteractive = False\n\nhardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\nlabel_fname = get_fnames('stanford_labels')\n\ndata, affine, hardi_img = load_nifti(hardi_fname, return_img=True)\nlabels = load_nifti_data(label_fname)\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\nseed_mask = (labels == 2)\nwhite_matter = (labels == 1) | (labels == 2)\nseeds = utils.seeds_from_mask(seed_mask, affine, density=1)\n\nresponse, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)\ncsd_model = ConstrainedSphericalDeconvModel(gtab, response, sh_order=6)\ncsd_fit = csd_model.fit(data, mask=white_matter)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use the GFA of the CSA model to build a stopping criterion.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.reconst.shm import CsaOdfModel\n\ncsa_model = CsaOdfModel(gtab, sh_order=6)\ngfa = csa_model.fit(data, mask=white_matter).gfa\nstopping_criterion = ThresholdStoppingCriterion(gfa, .25)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The Fiber Orientation Distribution (FOD) of the CSD model estimates the\ndistribution of small fiber bundles within each voxel. We can use this\ndistribution for probabilistic fiber tracking. One way to do this is to\nrepresent the FOD using a discrete sphere. This discrete FOD can be used by the\n``ProbabilisticDirectionGetter`` as a PMF for sampling tracking directions. We\nneed to clip the FOD to use it as a PMF because the latter cannot have negative\nvalues. Ideally, the FOD should be strictly positive, but because of noise\nand/or model failures sometimes it can have negative values.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.direction import ProbabilisticDirectionGetter\nfrom dipy.data import small_sphere\nfrom dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.io.streamline import save_trk\n\nfod = csd_fit.odf(small_sphere)\npmf = fod.clip(min=0)\nprob_dg = ProbabilisticDirectionGetter.from_pmf(pmf, max_angle=30.,\n                                                sphere=small_sphere)\nstreamline_generator = LocalTracking(prob_dg, stopping_criterion, seeds,\n                                     affine, step_size=.5)\nstreamlines = Streamlines(streamline_generator)\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_probabilistic_dg_pmf.trk\")\n\nif has_fury:\n    scene = window.Scene()\n    scene.add(actor.line(streamlines, colormap.line_colors(streamlines)))\n    window.record(scene, out_path='tractogram_probabilistic_dg_pmf.png',\n                  size=(800, 800))\n    if interactive:\n        window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_probabilistic_dg_pmf.png\n   :align: center\n\n   **Corpus Callosum using probabilistic direction getter from PMF**\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "One disadvantage of using a discrete PMF to represent possible tracking\ndirections is that it tends to take up a lot of memory (RAM). The size of the\nPMF, the FOD in this case, must be equal to the number of possible tracking\ndirections on the hemisphere, and every voxel has a unique PMF. In this case\nthe data is ``(81, 106, 76)`` and ``small_sphere`` has 181 directions so the\nFOD is ``(81, 106, 76, 181)``. One way to avoid sampling the PMF and holding it\nin memory is to build the direction getter directly from the spherical harmonic\n(SH) representation of the FOD. By using this approach, we can also use a\nlarger sphere, like ``default_sphere`` which has 362 directions on the\nhemisphere, without having to worry about memory limitations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.data import default_sphere\n\nprob_dg = ProbabilisticDirectionGetter.from_shcoeff(csd_fit.shm_coeff,\n                                                    max_angle=30.,\n                                                    sphere=default_sphere)\nstreamline_generator = LocalTracking(prob_dg, stopping_criterion, seeds,\n                                     affine, step_size=.5)\nstreamlines = Streamlines(streamline_generator)\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_probabilistic_dg_sh.trk\")\n\nif has_fury:\n    scene = window.Scene()\n    scene.add(actor.line(streamlines, colormap.line_colors(streamlines)))\n    window.record(scene, out_path='tractogram_probabilistic_dg_sh.png',\n                  size=(800, 800))\n    if interactive:\n        window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_probabilistic_dg_sh.png\n   :align: center\n\n   **Corpus Callosum using probabilistic direction getter from SH**\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Not all model fits have the ``shm_coeff`` attribute because not all models use\nthis basis to represent the data internally. However we can fit the ODF of any\nmodel to the spherical harmonic basis using the ``peaks_from_model`` function.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.direction import peaks_from_model\n\npeaks = peaks_from_model(csd_model, data, default_sphere, .5, 25,\n                         mask=white_matter, return_sh=True, parallel=True,\n                         num_processes=2)\nfod_coeff = peaks.shm_coeff\n\nprob_dg = ProbabilisticDirectionGetter.from_shcoeff(fod_coeff, max_angle=30.,\n                                                    sphere=default_sphere)\nstreamline_generator = LocalTracking(prob_dg, stopping_criterion, seeds,\n                                     affine, step_size=.5)\nstreamlines = Streamlines(streamline_generator)\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_probabilistic_dg_sh_pfm.trk\")\n\nif has_fury:\n    scene = window.Scene()\n    scene.add(actor.line(streamlines, colormap.line_colors(streamlines)))\n    window.record(scene, out_path='tractogram_probabilistic_dg_sh_pfm.png',\n                  size=(800, 800))\n    if interactive:\n        window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_probabilistic_dg_sh_pfm.png\n   :align: center\n\n   **Corpus Callosum using probabilistic direction getter from SH (\n   peaks_from_model)**\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. 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
}