{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Calculate SHORE scalar maps\n\nWe show how to calculate two SHORE-based scalar maps: return to origin\nprobability (RTOP) [Descoteaux2011]_ and mean square displacement (MSD)\n[Wu2007]_, [Wu2008]_ on your data. SHORE can be used with any multiple b-value\ndataset like multi-shell or DSI.\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.shore import ShoreModel" ] }, { "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\nread the 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)\nprint('data.shape (%d, %d, %d, %d)' % data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instantiate the Model.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "asm = ShoreModel(gtab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's just use only 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": [ "Fit the signal with the model and calculate the SHORE coefficients.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "asmfit = asm.fit(dataslice)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the analytical RTOP 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 = asmfit.rtop_signal()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we calculate the analytical RTOP on the propagator,\nthat corresponds to its central value.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Calculating... rtop_pdf')\nrtop_pdf = asmfit.rtop_pdf()" ] }, { "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\nLet's calculate the analytical mean square displacement on the propagator.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('Calculating... msd')\nmsd = asmfit.msd()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show the maps and save them to a file.\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')\nax2.set_axis_off()\nind = ax2.imshow(rtop_pdf.T, interpolation='nearest', origin='lower')\nplt.colorbar(ind)\nax3 = fig.add_subplot(2, 2, 3, title='msd')\nax3.set_axis_off()\nind = ax3.imshow(msd.T, interpolation='nearest', origin='lower', vmin=0)\nplt.colorbar(ind)\nplt.savefig('SHORE_maps.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: SHORE_maps.png\n :align: center\n\n RTOP and MSD calculated using the SHORE model.\n\n## References\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, p.\n 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 }