{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Calculate DSI-based scalar maps\n\nWe show how to calculate two DSI-based scalar maps: return to origin\nprobability (RTOP) [Descoteaux2011]_ and mean square displacement (MSD)\n[Wu2007]_, [Wu2008]_ on your dataset.\n\nFirst import the necessary modules:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\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\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 = load_nifti(fraw)\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)\n\nprint('data.shape (%d, %d, %d, %d)' % data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instantiate the Model and apply it to the data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dsmodel = DiffusionSpectrumModel(gtab, qgrid_size=35, filter_width=18.5)" ] }, { "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[30:70, 20:80, data.shape[2] // 2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Normalize the signal by the b0\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dataslice = dataslice / (dataslice[..., 0, None]).astype(float)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the return to origin probability on the signal\nthat corresponds to the integral of the signal.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Calculating... rtop_signal')\nrtop_signal = dsmodel.fit(dataslice).rtop_signal()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we calculate the return to origin probability on the propagator, that\ncorresponds to its central value. By default the propagator is divided by its\nsum in order to obtain a properly normalized pdf, however this normalization\nchanges the values of RTOP, therefore in order to compare it with the RTOP\npreviously calculated on the signal we turn the normalized parameter to false.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Calculating... rtop_pdf')\nrtop_pdf = dsmodel.fit(dataslice).rtop_pdf(normalized=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In theory, these two measures must be equal,\nto show that we calculate the mean square error on this two measures.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mse = np.sum((rtop_signal - rtop_pdf) ** 2) / rtop_signal.size\nprint(\"mse = %f\" % mse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "mse = 0.000000\n\nLeaving the normalized parameter to the default changes the values of the\nRTOP but not the contrast between the voxels.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Calculating... rtop_pdf_norm')\nrtop_pdf_norm = dsmodel.fit(dataslice).rtop_pdf()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's calculate the mean square displacement on the normalized propagator.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Calculating... msd_norm')\nmsd_norm = dsmodel.fit(dataslice).msd_discrete()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Turning the normalized parameter to false makes it possible to calculate\nthe mean square displacement on the propagator without normalization.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Calculating... msd')\nmsd = dsmodel.fit(dataslice).msd_discrete(normalized=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show the RTOP images and save them in rtop.png.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = plt.figure(figsize=(6, 6))\nax1 = fig.add_subplot(2, 2, 1, title='rtop_signal')\nax1.set_axis_off()\nind = ax1.imshow(rtop_signal.T, interpolation='nearest', origin='lower')\nplt.colorbar(ind)\nax2 = fig.add_subplot(2, 2, 2, title='rtop_pdf_norm')\nax2.set_axis_off()\nind = ax2.imshow(rtop_pdf_norm.T, interpolation='nearest', origin='lower')\nplt.colorbar(ind)\nax3 = fig.add_subplot(2, 2, 3, title='rtop_pdf')\nax3.set_axis_off()\nind = ax3.imshow(rtop_pdf.T, interpolation='nearest', origin='lower')\nplt.colorbar(ind)\nplt.savefig('rtop.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: rtop.png\n :align: center\n\n Return to origin probability.\n\nShow the MSD images and save them in msd.png.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = plt.figure(figsize=(7, 3))\nax1 = fig.add_subplot(1, 2, 1, title='msd_norm')\nax1.set_axis_off()\nind = ax1.imshow(msd_norm.T, interpolation='nearest', origin='lower')\nplt.colorbar(ind)\nax2 = fig.add_subplot(1, 2, 2, title='msd')\nax2.set_axis_off()\nind = ax2.imshow(msd.T, interpolation='nearest', origin='lower')\nplt.colorbar(ind)\nplt.savefig('msd.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: msd.png\n :align: center\n\n Mean square displacement.\n\n.. [Descoteaux2011] Descoteaux M. et al., \"Multiple q-shell diffusion\n propagator imaging\", Medical Image Analysis, vol 15, no 4, p. 603-621,\n 2011.\n\n.. [Wu2007] Wu Y. et al., \"Hybrid diffusion imaging\", NeuroImage, vol 36,\n p. 617-629, 2007.\n\n.. [Wu2008] Wu Y. et al., \"Computation of Diffusion Function Measures in\n q-Space Using Magnetic Resonance Hybrid Diffusion Imaging\", IEEE\n Transactions on Medical Imaging, vol 27, no 6, p. 858-865, 2008.\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 }