{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Reconstruct with Diffusion Spectrum Imaging\n\nWe show how to apply Diffusion Spectrum Imaging [Wedeen08]_ to\ndiffusion MRI datasets of Cartesian keyhole diffusion gradients.\n\nFirst import the necessary modules:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames, get_sphere\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti\nfrom dipy.reconst.dsi import DiffusionSpectrumModel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Download and get the data filenames for this tutorial.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fraw, fbval, fbvec = get_fnames('taiwan_ntu_dsi')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "img contains a nibabel Nifti1Image object (data) and gtab contains a\nGradientTable object (gradient information e.g. b-values). For example to read\nthe b-values it is possible to write print(gtab.bvals).\n\nLoad the raw diffusion data and the affine.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data, affine, voxel_size = load_nifti(fraw, return_voxsize=True)\nbvals, bvecs = read_bvals_bvecs(fbval, fbvec)\nbvecs[1:] = (bvecs[1:] /\n np.sqrt(np.sum(bvecs[1:] * bvecs[1:], axis=1))[:, None])\ngtab = gradient_table(bvals, bvecs)\nprint('data.shape (%d, %d, %d, %d)' % data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "data.shape ``(96, 96, 60, 203)``\n\nThis dataset has anisotropic voxel sizes, therefore reslicing is necessary.\n\nInstantiate the Model and apply it to the data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dsmodel = DiffusionSpectrumModel(gtab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's just use one slice only from the data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dataslice = data[:, :, data.shape[2] // 2]\n\ndsfit = dsmodel.fit(dataslice)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load an odf reconstruction sphere\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sphere = get_sphere('repulsion724')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the ODFs with this specific sphere\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ODF = dsfit.odf(sphere)\n\nprint('ODF.shape (%d, %d, %d)' % ODF.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ODF.shape ``(96, 96, 724)``\n\nIn a similar fashion it is possible to calculate the PDFs of all voxels\nin one call with the following way\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "PDF = dsfit.pdf()\n\nprint('PDF.shape (%d, %d, %d, %d, %d)' % PDF.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "PDF.shape ``(96, 96, 17, 17, 17)``\n\nWe see that even for a single slice this PDF array is close to 345 MBytes so we\nreally have to be careful with memory usage when use this function with a full\ndataset.\n\nThe simple solution is to generate/analyze the ODFs/PDFs by iterating through\neach voxel and not store them in memory if that is not necessary.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.core.ndindex import ndindex\n\nfor index in ndindex(dataslice.shape[:2]):\n pdf = dsmodel.fit(dataslice[index]).pdf()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you really want to save the PDFs of a full dataset on the disc we recommend\nusing memory maps (``numpy.memmap``) but still have in mind that even if you do\nthat for example for a dataset of volume size ``(96, 96, 60)`` you will need\nabout 2.5 GBytes which can take less space when reasonable spheres\n(with < 1000 vertices) are used.\n\nLet's now calculate a map of Generalized Fractional Anisotropy (GFA) [Tuch04]_\nusing the DSI ODFs.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.reconst.odf import gfa\n\nGFA = gfa(ODF)\n\nimport matplotlib.pyplot as plt\n\nfig_hist, ax = plt.subplots(1)\nax.set_axis_off()\nplt.imshow(GFA.T)\nplt.savefig('dsi_gfa.png', bbox_inches='tight')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: dsi_gfa.png\n :align: center\n\nSee also `example_reconst_dsi_metrics` for calculating different types\nof DSI maps.\n\n\n.. [Wedeen08] Wedeen et al., Diffusion spectrum magnetic resonance imaging\n (DSI) tractography of crossing fibers, Neuroimage, vol 41, no 4,\n 1267-1277, 2008.\n\n.. [Tuch04] Tuch, D.S, Q-ball imaging, MRM, vol 52, no 6, 1358-1372, 2004.\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 }