{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Introduction to Basic Tracking\n\nLocal fiber tracking is an approach used to model white matter fibers by\ncreating streamlines from local directional information. The idea is as\nfollows: if the local directionality of a tract/pathway segment is known, one\ncan integrate along those directions to build a complete representation of that\nstructure. Local fiber tracking is widely used in the field of diffusion MRI\nbecause it is simple and robust.\n\nIn order to perform local fiber tracking, three things are needed: 1) A method\nfor getting directions from a diffusion data set. 2) A method for identifying\nwhen the tracking must stop. 3) A set of seeds from which to\nbegin tracking. This example shows how to combine the 3 parts described above\nto create a tractography reconstruction from a diffusion data set.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To begin, let's load an example HARDI data set from Stanford. If you have\nnot already downloaded this data set, the first time you run this example you\nwill need to be connected to the internet and this dataset will be downloaded\nto your computer.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Enables/disables interactive visualization\ninteractive = False\n\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, load_nifti_data\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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This dataset provides a label map in which all white matter tissues are\nlabeled either 1 or 2. Let's create a white matter mask to restrict tracking to\nthe white matter.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "white_matter = (labels == 1) | (labels == 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. The first thing we need to begin fiber tracking is a way of getting\ndirections from this diffusion data set. In order to do that, we can fit the\ndata to a Constant Solid Angle ODF Model. This model will estimate the\nOrientation Distribution Function (ODF) at each voxel. The ODF is the\ndistribution of water diffusion as a function of direction. The peaks of an ODF\nare good estimates for the orientation of tract segments at a point in the\nimage. Here, we use ``peaks_from_model`` to fit the data and calculate the\nfiber directions in all voxels of the white matter.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.reconst.csdeconv import auto_response_ssst\nfrom dipy.reconst.shm import CsaOdfModel\nfrom dipy.data import default_sphere\nfrom dipy.direction import peaks_from_model\n\nresponse, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)\ncsa_model = CsaOdfModel(gtab, sh_order=6)\ncsa_peaks = peaks_from_model(csa_model, data, default_sphere,\n relative_peak_threshold=.8,\n min_separation_angle=45,\n mask=white_matter)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For quality assurance we can also visualize a slice from the direction field\nwhich we will use as the basis to perform the tracking. The visualization will\nbe done using the ``fury`` python package\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.viz import window, actor, has_fury\n\nif has_fury:\n scene = window.Scene()\n scene.add(actor.peak_slicer(csa_peaks.peak_dirs,\n csa_peaks.peak_values,\n colors=None))\n\n window.record(scene, out_path='csa_direction_field.png', size=(900, 900))\n\n if interactive:\n window.show(scene, size=(800, 800))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: csa_direction_field.png\n :align: center\n\n **Direction Field (peaks)**\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Next we need some way of restricting the fiber tracking to areas with good\ndirectionality information. We've already created the white matter mask,\nbut we can go a step further and restrict fiber tracking to those areas where\nthe ODF shows significant restricted diffusion by thresholding on\nthe generalized fractional anisotropy (GFA).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.tracking.stopping_criterion import ThresholdStoppingCriterion\n\nstopping_criterion = ThresholdStoppingCriterion(csa_peaks.gfa, .25)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, for quality assurance, we can also visualize a slice of the GFA and the\nresulting tracking mask.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n\nsli = csa_peaks.gfa.shape[2] // 2\nplt.figure('GFA')\nplt.subplot(1, 2, 1).set_axis_off()\nplt.imshow(csa_peaks.gfa[:, :, sli].T, cmap='gray', origin='lower')\n\nplt.subplot(1, 2, 2).set_axis_off()\nplt.imshow((csa_peaks.gfa[:, :, sli] > 0.25).T, cmap='gray', origin='lower')\n\nplt.savefig('gfa_tracking_mask.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: gfa_tracking_mask.png\n :align: center\n\n An example of a tracking mask derived from the generalized fractional\n anisotropy (GFA).\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Before we can begin tracking, we need to specify where to \"seed\" (begin)\nthe fiber tracking. Generally, the seeds chosen will depend on the pathways\none is interested in modeling. In this example, we'll use a $2 \\times 2 \\times\n2$ grid of seeds per voxel, in a sagittal slice of the corpus callosum.\nTracking from this region will give us a model of the corpus callosum tract.\nThis slice has label value ``2`` in the label's image.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.tracking import utils\n\nseed_mask = (labels == 2)\nseeds = utils.seeds_from_mask(seed_mask, affine, density=[2, 2, 2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can bring it all together using ``LocalTracking``, using\nthe EuDX algorithm [Garyfallidis12]_. ``EuDX`` [Garyfallidis12]_ is a fast\nalgorithm that we use here to generate streamlines. This algorithm is what is\nused here and the default option when providing the output of peaks directly\nin LocalTracking.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.tracking.local_tracking import LocalTracking\nfrom dipy.tracking.streamline import Streamlines\n\n# Initialization of LocalTracking. The computation happens in the next step.\nstreamlines_generator = LocalTracking(csa_peaks, stopping_criterion, seeds,\n affine=affine, step_size=.5)\n# Generate streamlines object\nstreamlines = Streamlines(streamlines_generator)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will then display the resulting streamlines using the ``fury``\npython package.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.viz import colormap\n\nif has_fury:\n # Prepare the display objects.\n color = colormap.line_colors(streamlines)\n\n streamlines_actor = actor.line(streamlines,\n colormap.line_colors(streamlines))\n\n # Create the 3D display.\n scene = window.Scene()\n scene.add(streamlines_actor)\n\n # Save still images for this static example. Or for interactivity use\n window.record(scene, out_path='tractogram_EuDX.png', size=(800, 800))\n if interactive:\n window.show(scene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: tractogram_EuDX.png\n :align: center\n\n **Corpus Callosum using EuDx**\n\nWe've created a deterministic set of streamlines using the EuDX algorithm. This\nis so called deterministic because if you repeat the fiber tracking (keeping\nall the inputs the same) you will get exactly the same set of streamlines.\nWe can save the streamlines as a Trackvis file so it can be loaded into other\nsoftware for visualization or further analysis.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.io.streamline import save_trk\n\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_EuDX.trk\", streamlines)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n.. [Garyfallidis12] Garyfallidis E., \"Towards an accurate brain tractography\",\nPhD thesis, University of Cambridge, 2012.\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 }