{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Using Various Stopping Criterion for Tractography\nThe stopping criterion determines if the tracking stops or continues at each\ntracking position. The tracking stops when it reaches an ending region\n(e.g. low FA, gray matter or corticospinal fluid regions) or exits the image\nboundaries. The tracking also stops if the direction getter has no direction\nto follow.\n\nEach stopping criterion determines if the stopping is 'valid' or\n'invalid'. A streamline is 'valid' when the stopping criterion determines if\nthe streamline stops in a position classified as 'ENDPOINT' or 'OUTSIDEIMAGE'.\nA streamline is 'invalid' when it stops in a position classified as\n'TRACKPOINT' or 'INVALIDPOINT'. These conditions are described below. The\n'LocalTracking' generator can be set to output all generated streamlines\nor only the 'valid' ones. See Girard et al. (2004) [Girard2014]_ and Smith et\nal.(2012) [Smith2012]_ for more details on these methods.\n\nThis example is an extension of the\n`example_tracking_deterministic` example. We begin by loading the\ndata, creating a seeding mask from white matter voxels of the corpus callosum,\nfitting a Constrained Spherical Deconvolution (CSD) reconstruction\nmodel and creating the maximum deterministic direction getter.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames, default_sphere\nfrom dipy.direction import DeterministicMaximumDirectionGetter\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti, load_nifti_data\nfrom dipy.io.streamline import save_trk\nfrom dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,\n                                   auto_response_ssst)\nfrom dipy.reconst.dti import fractional_anisotropy, TensorModel\nfrom dipy.tracking import utils\nfrom dipy.tracking.local_tracking import LocalTracking\nfrom dipy.tracking.streamline import Streamlines\nfrom dipy.tracking.stopping_criterion import (ActStoppingCriterion,\n                                              BinaryStoppingCriterion,\n                                              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_, _, f_pve_wm = get_fnames('stanford_pve_maps')\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\nwhite_matter = load_nifti_data(f_pve_wm)\n\nseed_mask = (labels == 2)\nseed_mask[white_matter < 0.5] = 0\nseeds = utils.seeds_from_mask(seed_mask, affine, density=2)\n\nresponse, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)\ncsd_model = ConstrainedSphericalDeconvModel(gtab, response)\ncsd_fit = csd_model.fit(data, mask=white_matter)\n\ndg = DeterministicMaximumDirectionGetter.from_shcoeff(csd_fit.shm_coeff,\n                                                      max_angle=30.,\n                                                      sphere=default_sphere)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Threshold Stopping Criterion\nA scalar map can be used to define where the tracking stops. The threshold\nstopping criterion uses a scalar map to stop the tracking whenever the\ninterpolated scalar value is lower than a fixed threshold. Here, we show\nan example using the fractional anisotropy (FA) map of the DTI model.\nThe threshold stopping criterion uses a trilinear interpolation at the\ntracking position.\n\n**Parameters**\n\n- metric_map: numpy array [:, :, :]\n- threshold: float\n\n**Stopping States**\n\n- 'ENDPOINT': stops at a position where metric_map < threshold; the streamline\nreached the target stopping area.\n- 'OUTSIDEIMAGE': stops at a position outside of metric_map; the streamline\nreached an area outside the image where no direction data is available.\n- 'TRACKPOINT': stops at a position because no direction is available; the\nstreamline is stopping where metric_map >= threshold, but there is no valid\ndirection to follow.\n- 'INVALIDPOINT': N/A.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tensor_model = TensorModel(gtab)\ntenfit = tensor_model.fit(data, mask=labels > 0)\nFA = fractional_anisotropy(tenfit.evals)\n\nthreshold_criterion = ThresholdStoppingCriterion(FA, .2)\n\nfig = plt.figure()\nmask_fa = FA.copy()\nmask_fa[mask_fa < 0.2] = 0\nplt.xticks([])\nplt.yticks([])\nplt.imshow(mask_fa[:, :, data.shape[2] // 2].T, cmap='gray', origin='lower',\n           interpolation='nearest')\nfig.tight_layout()\nfig.savefig('threshold_fa.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: threshold_fa.png\n :align: center\n\n **Thresholded fractional anisotropy map.**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "streamline_generator = LocalTracking(dg,\n                                     threshold_criterion,\n                                     seeds,\n                                     affine,\n                                     step_size=.5,\n                                     return_all=True)\nstreamlines = Streamlines(streamline_generator)\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_probabilistic_thresh_all.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_deterministic_thresh_all.png',\n                  size=(800, 800))\n    if interactive:\n        window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_deterministic_thresh_all.png\n :align: center\n\n **Corpus Callosum using deterministic tractography with a thresholded\n fractional anisotropy mask.**\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Binary Stopping Criterion\nA binary mask can be used to define where the tracking stops. The binary\nstopping criterion stops the tracking whenever the tracking position is outside\nthe mask. Here, we show how to obtain the binary stopping criterion from\nthe white matter mask defined above. The binary stopping criterion uses a\nnearest-neighborhood interpolation at the tracking position.\n\n**Parameters**\n\n- mask: numpy array [:, :, :]\n\n**Stopping States**\n\n- 'ENDPOINT': stops at a position where mask = 0; the streamline\nreached the target stopping area.\n- 'OUTSIDEIMAGE': stops at a position outside of metric_map; the streamline\nreached an area outside the image where no direction data is available.\n- 'TRACKPOINT': stops at a position because no direction is available; the\nstreamline is stopping where mask > 0, but there is no valid direction to\nfollow.\n- 'INVALIDPOINT': N/A.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "binary_criterion = BinaryStoppingCriterion(white_matter == 1)\n\nfig = plt.figure()\nplt.xticks([])\nplt.yticks([])\nfig.tight_layout()\nplt.imshow(white_matter[:, :, data.shape[2] // 2].T, cmap='gray',\n           origin='lower', interpolation='nearest')\n\nfig.savefig('white_matter_mask.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: white_matter_mask.png\n :align: center\n\n **White matter binary mask.**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "streamline_generator = LocalTracking(dg,\n                                     binary_criterion,\n                                     seeds,\n                                     affine,\n                                     step_size=.5,\n                                     return_all=True)\nstreamlines = Streamlines(streamline_generator)\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_deterministic_binary_all.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_deterministic_binary_all.png',\n                  size=(800, 800))\n    if interactive:\n        window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_deterministic_binary_all.png\n :align: center\n\n **Corpus Callosum using deterministic tractography with a binary white\n matter mask.**\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## ACT Stopping Criterion\nAnatomically-constrained tractography (ACT) [Smith2012]_ uses information from\nanatomical images to determine when the tractography stops. The ``include_map``\ndefines when the streamline reached a 'valid' stopping region (e.g. gray\nmatter partial volume estimation (PVE) map) and the ``exclude_map`` defines\nwhen the streamline reached an 'invalid' stopping region (e.g. corticospinal\nfluid PVE map). The background of the anatomical image should be added to the\n``include_map`` to keep streamlines exiting the brain (e.g. through the\nbrain stem). The ACT stopping criterion uses a trilinear interpolation\nat the tracking position.\n\n**Parameters**\n\n- ``include_map``: numpy array ``[:, :, :]``,\n- ``exclude_map``: numpy array ``[:, :, :]``,\n\n**Stopping States**\n\n- 'ENDPOINT': stops at a position where ``include_map`` > 0.5; the streamline\nreached the target stopping area.\n- 'OUTSIDEIMAGE': stops at a position outside of ``include_map`` or\n``exclude_map``; the streamline reached an area outside the image where no\ndirection data is available.\n- 'TRACKPOINT': stops at a position because no direction is available; the\nstreamline is stopping where ``include_map`` < 0.5 and ``exclude_map`` < 0.5,\nbut there is no valid direction to follow.\n- 'INVALIDPOINT': ``exclude_map`` > 0.5; the streamline reach a position which\nis anatomically not plausible.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "f_pve_csf, f_pve_gm, f_pve_wm = get_fnames('stanford_pve_maps')\npve_csf_data = load_nifti_data(f_pve_csf)\npve_gm_data = load_nifti_data(f_pve_gm)\npve_wm_data = load_nifti_data(f_pve_wm)\n\nbackground = np.ones(pve_gm_data.shape)\nbackground[(pve_gm_data + pve_wm_data + pve_csf_data) > 0] = 0\n\ninclude_map = pve_gm_data\ninclude_map[background > 0] = 1\nexclude_map = pve_csf_data\n\nact_criterion = ActStoppingCriterion(include_map, exclude_map)\n\nfig = plt.figure()\nplt.subplot(121)\nplt.xticks([])\nplt.yticks([])\nplt.imshow(include_map[:, :, data.shape[2] // 2].T, cmap='gray',\n           origin='lower', interpolation='nearest')\n\nplt.subplot(122)\nplt.xticks([])\nplt.yticks([])\nplt.imshow(exclude_map[:, :, data.shape[2] // 2].T, cmap='gray',\n           origin='lower', interpolation='nearest')\n\nfig.tight_layout()\nfig.savefig('act_maps.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: act_maps.png\n :align: center\n\n **Include (left) and exclude (right) maps for ACT.**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "streamline_generator = LocalTracking(dg,\n                                     act_criterion,\n                                     seeds,\n                                     affine,\n                                     step_size=.5,\n                                     return_all=True)\nstreamlines = Streamlines(streamline_generator)\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_deterministic_act_all.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_deterministic_act_all.png',\n                  size=(800, 800))\n    if interactive:\n        window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_deterministic_act_all.png\n :align: center\n\n **Corpus Callosum using deterministic tractography with ACT stopping\n criterion.**\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "streamline_generator = LocalTracking(dg,\n                                     act_criterion,\n                                     seeds,\n                                     affine,\n                                     step_size=.5,\n                                     return_all=False)\nstreamlines = Streamlines(streamline_generator)\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_deterministic_act_valid.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_deterministic_act_valid.png',\n                  size=(800, 800))\n    if interactive:\n        window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_deterministic_act_valid.png\n :align: center\n\n **Corpus Callosum using deterministic tractography with ACT stopping\n criterion. Streamlines ending in gray matter region only.**\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The threshold and binary stopping criterion use respectively a scalar map and a\nbinary mask to stop the tracking. The ACT stopping criterion use partial volume\nfraction (PVE) maps from an anatomical image to stop the tracking.\nAdditionally, the ACT stopping criterion determines if the tracking stopped in\nexpected regions (e.g. gray matter) and allows the user to get only\nstreamlines stopping in those regions.\n\n### Notes\nCurrently,the proposed method that cuts streamlines going through\nsubcortical gray matter regions is not implemented. The\nbacktracking technique for streamlines reaching INVALIDPOINT is not\nimplemented either [Smith2012]_.\n\n\n### References\n\n.. [Smith2012] Smith, R. E., Tournier, J.-D., Calamante, F., & Connelly, A.\n    Anatomically-constrained tractography: Improved diffusion MRI\n    streamlines tractography through effective use of anatomical\n    information. NeuroImage, 63(3), 1924-1938, 2012.\n\n.. [Girard2014] Girard, G., Whittingstall, K., Deriche, R., & Descoteaux, M.\n    Towards quantitative connectivity analysis: reducing tractography biases.\n    NeuroImage, 98, 266-278, 2014.\n\n.. include:: ../links_names.inc\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
}