{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Using the RESTORE algorithm for robust tensor fitting\n\nThe diffusion tensor model takes into account certain kinds of noise (thermal),\nbut not other kinds, such as \"physiological\" noise. For example, if a subject\nmoves during acquisition of one of the diffusion-weighted samples, this\nmight have a substantial effect on the parameters of the tensor fit calculated\nin all voxels in the brain for that subject. One of the pernicious consequences\nof this is that it can lead to wrong interpretation of group differences. For\nexample, some groups of participants (e.g. young children, patient groups,\netc.) are particularly prone to motion and differences in tensor parameters and\nderived statistics (such as FA) due to motion would be confounded with actual\ndifferences in the physical properties of the white matter. An example of this\nwas shown in a paper by Yendiki et al. [Yendiki2013]_.\n\nOne of the strategies to deal with this problem is to apply an automatic method\nfor detecting outliers in the data, excluding these outliers and refitting the\nmodel without the presence of these outliers. This is often referred to as\n\"robust model fitting\". One of the common algorithms for robust tensor fitting\nis called RESTORE, and was first proposed by Chang et al. [Chang2005]_.\n\nIn the following example, we will demonstrate how to use RESTORE on a simulated\ndataset, which we will corrupt by adding intermittent noise.\n\nWe start by importing a few of the libraries we will use.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The module ``dipy.reconst.dti`` contains the implementation of tensor fitting,\nincluding an implementation of the RESTORE algorithm.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import dipy.reconst.dti as dti" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``dipy.data`` is used for small datasets that we use in tests and examples.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import dipy.data as dpd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``dipy.io.image`` is for loading / saving imaging datasets\n``dipy.io.gradients`` is for loading / saving our bvals and bvecs\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.io.image import load_nifti\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.core.gradients import gradient_table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``dipy.viz`` package is used for 3D visualization and matplotlib for 2D\nvisualizations:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.viz import window, actor\nimport matplotlib.pyplot as plt\n\nimport dipy.denoise.noise_estimate as ne\n\n# Enables/disables interactive visualization\ninteractive = False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If needed, the ``fetch_stanford_hardi`` function will download the raw dMRI\ndataset of a single subject. The size of this dataset is 87 MBytes. You only\nneed to fetch once.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hardi_fname, hardi_bval_fname, hardi_bvec_fname = dpd.get_fnames('stanford_hardi')\ndata, affine = load_nifti(hardi_fname)\n\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We initialize a DTI model class instance using the gradient table used in the\nmeasurement. By default, ``dti.TensorModel`` will use a weighted least-squares\nalgorithm (described in [Chang2005]_) to fit the parameters of the model. We\ninitialize this model as a baseline for comparison of noise-corrupted models:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dti_wls = dti.TensorModel(gtab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the purpose of this example, we will focus on the data from a region of\ninterest (ROI) surrounding the Corpus Callosum. We define that ROI as the\nfollowing indices:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "roi_idx = (slice(20, 50), slice(55, 85), slice(38, 39))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And use them to index into the data:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data = data[roi_idx]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This dataset is not very noisy, so we will artificially corrupt it to simulate\nthe effects of \"physiological\" noise, such as subject motion. But first, let's\nestablish a baseline, using the data as it is:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fit_wls = dti_wls.fit(data)\n\nfa1 = fit_wls.fa\nevals1 = fit_wls.evals\nevecs1 = fit_wls.evecs\ncfa1 = dti.color_fa(fa1, evecs1)\nsphere = dpd.default_sphere" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We visualize the ODFs in the ROI using ``dipy.viz`` module:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "scene = window.Scene()\nscene.add(actor.tensor_slicer(evals1, evecs1, scalar_colors=cfa1,\n sphere=sphere, scale=0.3))\nprint('Saving illustration as tensor_ellipsoids_wls.png')\nwindow.record(scene, out_path='tensor_ellipsoids_wls.png', size=(600, 600))\nif interactive:\n window.show(scene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: tensor_ellipsoids_wls.png\n :align: center\n\n Tensor Ellipsoids.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "scene.clear()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we corrupt the data with some noise. To simulate a subject that moves\nintermittently, we will replace a few of the images with a very low signal\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "noisy_data = np.copy(data)\nnoisy_idx = slice(-10, None) # The last 10 volumes are corrupted\nnoisy_data[..., noisy_idx] = 1.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use the same model to fit this noisy data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fit_wls_noisy = dti_wls.fit(noisy_data)\nfa2 = fit_wls_noisy.fa\nevals2 = fit_wls_noisy.evals\nevecs2 = fit_wls_noisy.evecs\ncfa2 = dti.color_fa(fa2, evecs2)\n\nscene = window.Scene()\nscene.add(actor.tensor_slicer(evals2, evecs2, scalar_colors=cfa2,\n sphere=sphere, scale=0.3))\nprint('Saving illustration as tensor_ellipsoids_wls_noisy.png')\nwindow.record(scene, out_path='tensor_ellipsoids_wls_noisy.png',\n size=(600, 600))\nif interactive:\n window.show(scene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In places where the tensor model is particularly sensitive to noise, the\nresulting tensor field will be distorted\n\n.. figure:: tensor_ellipsoids_wls_noisy.png\n :align: center\n\n Tensor Ellipsoids from noisy data.\n\nTo estimate the parameters from the noisy data using RESTORE, we need to\nestimate what would be a reasonable amount of noise to expect in the\nmeasurement. To do that, we use the ``dipy.denoise.noise_estimate`` module:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sigma = ne.estimate_sigma(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This estimate of the standard deviation will be used by the RESTORE algorithm\nto identify the outliers in each voxel and is given as an input when\ninitializing the TensorModel object:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dti_restore = dti.TensorModel(gtab, fit_method='RESTORE', sigma=sigma)\nfit_restore_noisy = dti_restore.fit(noisy_data)\nfa3 = fit_restore_noisy.fa\nevals3 = fit_restore_noisy.evals\nevecs3 = fit_restore_noisy.evecs\ncfa3 = dti.color_fa(fa3, evecs3)\n\nscene = window.Scene()\nscene.add(actor.tensor_slicer(evals3, evecs3, scalar_colors=cfa3,\n sphere=sphere, scale=0.3))\nprint('Saving illustration as tensor_ellipsoids_restore_noisy.png')\nwindow.record(scene, out_path='tensor_ellipsoids_restore_noisy.png',\n size=(600, 600))\nif interactive:\n window.show(scene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: tensor_ellipsoids_restore_noisy.png\n :align: center\n\n Tensor Ellipsoids from noisy data recovered with RESTORE.\n\nThe tensor field looks rather restored to its noiseless state in this\nimage, but to convince ourselves further that this did the right thing, we will\ncompare the distribution of FA in this region relative to the baseline, using\nthe RESTORE estimate and the WLS estimate [Chung2006]_.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig_hist, ax = plt.subplots(1)\nax.hist(np.ravel(fa2), color='b', histtype='step', label='WLS')\nax.hist(np.ravel(fa3), color='r', histtype='step', label='RESTORE')\nax.hist(np.ravel(fa1), color='g', histtype='step', label='Original')\nax.set_xlabel('Fractional Anisotropy')\nax.set_ylabel('Count')\nplt.legend()\nfig_hist.savefig('dti_fa_distributions.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: dti_fa_distributions.png\n :align: center\n\n\nThis demonstrates that RESTORE can recover a distribution of FA that more\nclosely resembles the baseline distribution of the noiseless signal, and\ndemonstrates the utility of the method to data with intermittent\nnoise. Importantly, this method assumes that the tensor is a good\nrepresentation of the diffusion signal in the data. If you have reason to\nbelieve this is not the case (for example, you have data with very high b\nvalues and you are particularly interested in locations in the brain in which\nfibers cross), you might want to use a different method to fit your data.\n\n\n## References\n\n.. [Yendiki2013] Yendiki, A, Koldewynb, K, Kakunooria, S, Kanwisher, N, and\n Fischl, B. (2013). Spurious group differences due to head motion in a\n diffusion MRI study. Neuroimage.\n\n.. [Chang2005] Chang, L-C, Jones, DK and Pierpaoli, C (2005). RESTORE: robust\n estimation of tensors by outlier rejection. MRM, 53: 1088-95.\n\n.. [Chung2006] Chung, SW, Lu, Y, Henry, R-G, (2006). Comparison of bootstrap\n approaches for estimation of uncertainties of DTI parameters. NeuroImage 33,\n 531-541.\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 }