{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Tracking with the Sparse Fascicle Model\n\nTracking requires a per-voxel model. Here, the model is the Sparse Fascicle\nModel (SFM), described in [Rokem2015]_. This model reconstructs the diffusion\nsignal as a combination of the signals from different fascicles (see also\n`sfm-reconst`).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.core.gradients import gradient_table\nfrom dipy.data import get_sphere, get_fnames\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti, load_nifti_data\nfrom dipy.direction.peaks import peaks_from_model\nfrom dipy.io.streamline import save_trk\nfrom dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.reconst.csdeconv import auto_response_ssst\nfrom dipy.reconst import sfm\nfrom dipy.tracking import utils\nfrom dipy.tracking.local_tracking import LocalTracking\nfrom dipy.tracking.streamline import (select_random_set_of_streamlines,\n transform_streamlines,\n Streamlines)\nfrom dipy.tracking.stopping_criterion import ThresholdStoppingCriterion\nfrom dipy.viz import window, actor, colormap, has_fury\nfrom numpy.linalg import inv\n\n# Enables/disables interactive visualization\ninteractive = False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To begin, we read the Stanford HARDI data set into memory:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hardi_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 data set provides a label map (generated using [FreeSurfer](https://surfer.nmr.mgh.harvard.edu/)), in which the white matter voxels are\nlabeled as either 1 or 2:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "white_matter = (labels == 1) | (labels == 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step in tracking is generating a model from which tracking directions\ncan be extracted in every voxel.\n\nFor the SFM, this requires first that we define a canonical response function\nthat will be used to deconvolve the signal in every voxel\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "response, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We initialize an SFM model object, using this response function and using the\ndefault sphere (362 vertices, symmetrically distributed on the surface of the\nsphere):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sphere = 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": [ "We fit this model to the data in each voxel in the white-matter mask, so that\nwe can use these directions in tracking:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pnm = peaks_from_model(sf_model, data, sphere,\n relative_peak_threshold=.5,\n min_separation_angle=25,\n mask=white_matter,\n parallel=True,\n num_processes=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A ThresholdStoppingCriterion object is used to segment the data to track only\nthrough areas in which the Generalized Fractional Anisotropy (GFA) is\nsufficiently high.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "stopping_criterion = ThresholdStoppingCriterion(pnm.gfa, .25)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tracking will be started from a set of seeds evenly distributed in the white\nmatter:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "seeds = utils.seeds_from_mask(white_matter, affine, density=[2, 2, 2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the sake of brevity, we will take only the first 1000 seeds, generating\nonly 1000 streamlines. Remove this line to track from many more points in all\nof the white matter\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "seeds = seeds[:1000]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have the necessary components to construct a tracking pipeline and\nexecute the tracking\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "streamline_generator = LocalTracking(pnm, stopping_criterion, seeds, affine,\n step_size=.5)\nstreamlines = Streamlines(streamline_generator)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we will create a visualization of these streamlines, relative to this\nsubject's T1-weighted anatomy:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "t1_fname = get_fnames('stanford_t1')\nt1_data, t1_aff = load_nifti(t1_fname)\ncolor = colormap.line_colors(streamlines)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To speed up visualization, we will select a random sub-set of streamlines to\ndisplay. This is particularly important, if you track from seeds throughout the\nentire white matter, generating many streamlines. In this case, for\ndemonstration purposes, we select a subset of 900 streamlines.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plot_streamlines = select_random_set_of_streamlines(streamlines, 900)\n\nif has_fury:\n streamlines_actor = actor.streamtube(\n list(transform_streamlines(plot_streamlines, inv(t1_aff))),\n colormap.line_colors(streamlines), linewidth=0.1)\n\n vol_actor = actor.slicer(t1_data)\n\n vol_actor.display(40, None, None)\n vol_actor2 = vol_actor.copy()\n vol_actor2.display(None, None, 35)\n\n scene = window.Scene()\n scene.add(streamlines_actor)\n scene.add(vol_actor)\n scene.add(vol_actor2)\n\n window.record(scene, out_path='tractogram_sfm.png', size=(800, 800))\n if interactive:\n window.show(scene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: tractogram_sfm.png\n :align: center\n\n **Sparse Fascicle Model tracks**\n\nFinally, we can save these streamlines to a 'trk' file, for use in other\nsoftware, or for further analysis.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_sfm_detr.trk\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n\n.. [Rokem2015] Ariel Rokem, Jason D. Yeatman, Franco Pestilli, Kendrick\n N. Kay, Aviv Mezer, Stefan van der Walt, Brian A. Wandell (2015). Evaluating\n the accuracy of diffusion MRI models in white matter. PLoS ONE 10(4):\n e0123272. doi:10.1371/journal.pone.0123272\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 }