{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Extracting AFQ tract profiles from segmented bundles\n\nIn this example, we will extract the values of a statistic from a\nvolume, using the coordinates along the length of a bundle. These are called\n`tract profiles`\n\nOne of the challenges of extracting tract profiles is that some of the\nstreamlines in a bundle may diverge significantly from the bundle in some\nlocations. To overcome this challenge, we will use a strategy similar to that\ndescribed in [Yeatman2012]_: We will weight the contribution of each streamline\nto the bundle profile based on how far this streamline is from the mean\ntrajectory of the bundle at that location.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import dipy.stats.analysis as dsa\nimport dipy.tracking.streamline as dts\nfrom dipy.segment.clustering import QuickBundles\nfrom dipy.segment.metricspeed import AveragePointwiseEuclideanMetric\nfrom dipy.segment.featurespeed import ResampleFeature\nfrom dipy.data.fetcher import get_two_hcp842_bundles\nimport dipy.data as dpd\nfrom dipy.io.streamline import load_trk\nfrom dipy.io.image import load_nifti\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport os.path as op" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get started, we will grab the bundles.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bundles_path = dpd.fetch_bundles_2_subjects()\nbundles_folder = bundles_path[1]\n\ncst_l_file = op.join(bundles_folder, \"bundles_2_subjects\", \"subj_2\", \"bundles\",\n \"bundles_cst.left.trk\")\naf_l_file = op.join(bundles_folder, \"bundles_2_subjects\", \"subj_2\", \"bundles\",\n \"bundles_af.left.trk\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Either way, we can use the `dipy.io` API to read in the bundles from file.\n`load_trk` returns both the streamlines, as well as header information.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "cst_l = load_trk(cst_l_file, \"same\", bbox_valid_check=False).streamlines\naf_l = load_trk(af_l_file, \"same\", bbox_valid_check=False).streamlines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next step, we need to make sure that all the streamlines in each bundle\nare oriented the same way. For example, for the CST, we want to make sure that\nall the bundles have their cortical termination at one end of the streamline.\nThis is that when we later extract values from a volume, we won't have\ndifferent streamlines going in opposite directions.\n\nTo orient all the streamlines in each bundles, we will create standard\nstreamlines, by finding the centroids of the left AF and CST bundle models.\n\nThe advantage of using the model bundles is that we can use the same standard\nfor different subjects, which means that we'll get roughly the same orientation\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model_af_l_file, model_cst_l_file = get_two_hcp842_bundles()\n\nmodel_af_l = load_trk(model_af_l_file, \"same\",\n bbox_valid_check=False).streamlines\nmodel_cst_l = load_trk(model_cst_l_file, \"same\",\n bbox_valid_check=False).streamlines\n\n\nfeature = ResampleFeature(nb_points=100)\nmetric = AveragePointwiseEuclideanMetric(feature)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we are going to include all of the streamlines in the single cluster\nfrom the streamlines, we set the threshold to `np.inf`. We pull out the\ncentroid as the standard.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "qb = QuickBundles(np.inf, metric=metric)\n\ncluster_cst_l = qb.cluster(model_cst_l)\nstandard_cst_l = cluster_cst_l.centroids[0]\n\ncluster_af_l = qb.cluster(model_af_l)\nstandard_af_l = cluster_af_l.centroids[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use the centroid streamline for each atlas bundle as the standard to orient\nall of the streamlines in each bundle from the individual subject. Here, the\naffine used is the one from the transform between the atlas and individual\ntractogram. This is so that the orienting is done relative to the space of the\nindividual, and not relative to the atlas space.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "oriented_cst_l = dts.orient_by_streamline(cst_l, standard_cst_l)\noriented_af_l = dts.orient_by_streamline(af_l, standard_af_l)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read volumetric data from an image corresponding to this subject.\n\nFor the purpose of this, we've extracted only the FA within the bundles in\nquestion, but in real use, this is where you would add the FA map of your\nsubject.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "files, folder = dpd.fetch_bundle_fa_hcp()\n\nfa, fa_affine = load_nifti(op.join(folder, \"hcp_bundle_fa.nii.gz\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate weights for each bundle:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "w_cst_l = dsa.gaussian_weights(oriented_cst_l)\nw_af_l = dsa.gaussian_weights(oriented_af_l)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then use the weights to calculate the tract profiles for each bundle\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "profile_cst_l = dsa.afq_profile(fa, oriented_cst_l, fa_affine,\n weights=w_cst_l)\n\nprofile_af_l = dsa.afq_profile(fa, oriented_af_l, fa_affine,\n weights=w_af_l)\n\nfig, (ax1, ax2) = plt.subplots(1, 2)\n\nax1.plot(profile_cst_l)\nax1.set_ylabel(\"Fractional anisotropy\")\nax1.set_xlabel(\"Node along CST\")\nax2.plot(profile_af_l)\nax2.set_xlabel(\"Node along AF\")\nfig.savefig(\"tract_profiles\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: tract_profiles.png\n :align: center\n\n Bundle profiles for the fractional anisotropy in the left CST (left) and\n left AF (right).\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n\n.. [Yeatman2012] Yeatman, Jason D., Robert F. Dougherty, Nathaniel J. Myall,\n Brian A. Wandell, and Heidi M. Feldman. 2012. \"Tract Profiles of White\n Matter Properties: Automating Fiber-Tract Quantification\" PloS One 7 (11):\n e49790.\n\n.. [Garyfallidis17] Garyfallidis et al. Recognition of white matter bundles\n using local and global streamline-based registration and clustering,\n Neuroimage, 2017.\n\n.. [Garyfallidis12] Garyfallidis E. et al., QuickBundles a method for\n tractography simplification, Frontiers in Neuroscience, vol 6, no 175, 2012.\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 }