{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Intravoxel incoherent motion\nThe intravoxel incoherent motion (IVIM) model describes diffusion\nand perfusion in the signal acquired with a diffusion MRI sequence\nthat contains multiple low b-values. The IVIM model can be understood\nas an adaptation of the work of Stejskal and Tanner [Stejskal65]_\nin biological tissue, and was proposed by Le Bihan [LeBihan84]_.\nThe model assumes two compartments: a slow moving compartment,\nwhere particles diffuse in a Brownian fashion as a consequence of thermal\nenergy, and a fast moving compartment (the vascular compartment), where\nblood moves as a consequence of a pressure gradient. In the first compartment,\nthe diffusion coefficient is $\\mathbf{D}$ while in the second compartment, a\npseudo diffusion term $\\mathbf{D^*}$ is introduced that describes the\ndisplacement of the blood elements in an assumed randomly laid out vascular\nnetwork, at the macroscopic level. According to [LeBihan84]_,\n$\\mathbf{D^*}$ is greater than $\\mathbf{D}$.\n\nThe IVIM model expresses the MRI signal as follows:\n\n .. math::\n S(b)=S_0(fe^{-bD^*}+(1-f)e^{-bD})\n\nwhere $\\mathbf{b}$ is the diffusion gradient weighing value (which is dependent\non the measurement parameters), $\\mathbf{S_{0}}$ is the signal in the absence\nof diffusion gradient sensitization, $\\mathbf{f}$ is the perfusion\nfraction, $\\mathbf{D}$ is the diffusion coefficient and $\\mathbf{D^*}$ is\nthe pseudo-diffusion constant, due to vascular contributions.\n\nIn the following example we show how to fit the IVIM model on a\ndiffusion-weighted dataset and visualize the diffusion and pseudo-diffusion\ncoefficients. First, we import all relevant modules:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n\nfrom dipy.reconst.ivim import IvimModel\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_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We get an IVIM dataset using DIPY_'s data fetcher ``read_ivim``.\nThis dataset was acquired with 21 b-values in 3 different directions.\nVolumes corresponding to different directions were registered to each\nother, and averaged across directions. Thus, this dataset has 4 dimensions,\nwith the length of the last dimension corresponding to the number\nof b-values. In order to use this model the data should contain signals\nmeasured at 0 bvalue.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fraw, fbval, fbvec = get_fnames('ivim')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The gtab contains a GradientTable object (information about the gradients e.g.\nb-values and b-vectors). We get the data from the file using\n``load_nifti_data``.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data = load_nifti_data(fraw)\nbvals, bvecs = read_bvals_bvecs(fbval, fbvec)\ngtab = gradient_table(bvals, bvecs, b0_threshold=0)\nprint('data.shape (%d, %d, %d, %d)' % data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data has 54 slices, with 256-by-256 voxels in each slice. The fourth\ndimension corresponds to the b-values in the gtab. Let us visualize the data\nby taking a slice midway(z=33) at $\\mathbf{b} = 0$.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "z = 33\nb = 0\n\nplt.imshow(data[:, :, z, b].T, origin='lower', cmap='gray',\n interpolation='nearest')\nplt.axhline(y=100)\nplt.axvline(x=170)\nplt.savefig(\"ivim_data_slice.png\")\nplt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: ivim_data_slice.png\n :align: center\n\n Heat map of a slice of data\n\nThe region around the intersection of the cross-hairs in the figure\ncontains cerebral spinal fluid (CSF), so it should have a very high\n$\\mathbf{f}$ and $\\mathbf{D^*}$, the area just medial to that is white matter\nso that should be lower, and the region more laterally contains a mixture of\ngray matter and CSF. That should give us some contrast to see the\nvalues varying across the regions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x1, x2 = 90, 155\ny1, y2 = 90, 170\ndata_slice = data[x1:x2, y1:y2, z, :]\n\nplt.imshow(data[x1:x2, y1:y2, z, b].T, origin='lower',\n cmap=\"gray\", interpolation='nearest')\nplt.savefig(\"CSF_slice.png\")\nplt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: CSF_slice.png\n :align: center\n\n Heat map of the CSF slice selected.\n\nNow that we have prepared the datasets we can go forward with\nthe ivim fit. We provide two methods of fitting the parameters of the IVIM\nmulti-exponential model explained above. We first fit the model with a simple\nfitting approach by passing the option `fit_method='trr'`. This method uses\na two-stage approach: first, a linear fit used to get quick initial guesses\nfor the parameters $\\mathbf{S_{0}}$ and $\\mathbf{D}$ by considering b-values\ngreater than ``split_b_D`` (default: 400))and assuming a mono-exponential\nsignal. This is based on the assumption that at high b-values the signal can be\napproximated as a mono exponential decay and by taking the logarithm of the\nsignal values a linear fit can be obtained. Another linear fit for ``S0``\n(bvals < ``split_b_S0`` (default: 200)) follows and ``f`` is estimated using\n$1 - S0_{prime}/S0$. Then a non-linear least-squares fitting is performed to\nfit ``D_star`` and ``f``. If the ``two_stage`` flag is set to ``True`` while\ninitializing the model, a final non-linear least squares fitting is performed\nfor all the parameters. All initializations for the model such as ``split_b_D``\nare passed while creating the ``IvimModel``. If you are using Scipy 0.17, you\ncan also set bounds by setting ``bounds=([0., 0., 0.,0.], [np.inf, 1., 1., 1.]))``\nwhile initializing the ``IvimModel``.\n\nFor brevity, we focus on a small section of the slice as selected above,\nto fit the IVIM model. First, we instantiate the IvimModel object.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ivimmodel = IvimModel(gtab, fit_method='trr')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To fit the model, call the `fit` method and pass the data for fitting.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ivimfit = ivimmodel.fit(data_slice)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fit method creates a IvimFit object which contains the\nparameters of the model obtained after fitting. These are accessible\nthrough the `model_params` attribute of the IvimFit object.\nThe parameters are arranged as a 4D array, corresponding to the spatial\ndimensions of the data, and the last dimension (of length 4)\ncorresponding to the model parameters according to the following\norder : $\\mathbf{S_{0}, f, D^*, D}$.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ivimparams = ivimfit.model_params\nprint(\"ivimparams.shape : {}\".format(ivimparams.shape))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we see, we have a 20x20 slice at the height z = 33. Thus we\nhave 400 voxels. We will now plot the parameters obtained from the\nfit for a voxel and also various maps for the entire slice.\nThis will give us an idea about the diffusion and perfusion in\nthat section. Let(i, j) denote the coordinate of the voxel. We have\nalready fixed the z component as 33 and hence we will get a slice\nwhich is 33 units above.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "i, j = 10, 10\nestimated_params = ivimfit.model_params[i, j, :]\nprint(estimated_params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can map the perfusion and diffusion maps for the slice. We\nwill plot a heatmap showing the values using a colormap. It will be\nuseful to define a plotting function for the heatmap here since we\nwill use it to plot for all the IVIM parameters. We will need to specify\nthe lower and upper limits for our data. For example, the perfusion\nfractions should be in the range (0,1). Similarly, the diffusion and\npseudo-diffusion constants are much smaller than 1. We pass an argument\ncalled ``variable`` to out plotting function which gives the label for\nthe plot.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plot_map(raw_data, variable, limits, filename):\n fig, ax = plt.subplots(1)\n lower, upper = limits\n ax.set_title('Map for {}'.format(variable))\n im = ax.imshow(raw_data.T, origin='lower', clim=(lower, upper),\n cmap=\"gray\", interpolation='nearest')\n fig.colorbar(im)\n fig.savefig(filename)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us get the various plots with `fit_method = 'trr'` so that we can visualize\nthem in one page\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plot_map(ivimfit.S0_predicted, \"Predicted S0\", (0, 10000), \"predicted_S0.png\")\nplot_map(data_slice[:, :, 0], \"Measured S0\", (0, 10000), \"measured_S0.png\")\nplot_map(ivimfit.perfusion_fraction, \"f\", (0, 1), \"perfusion_fraction.png\")\nplot_map(ivimfit.D_star, \"D*\", (0, 0.01), \"perfusion_coeff.png\")\nplot_map(ivimfit.D, \"D\", (0, 0.001), \"diffusion_coeff.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we will fit the same model with a more refined optimization process with\n`fit_method='VarPro'` (for \"Variable Projection\"). The VarPro computes the IVIM\nparameters using the MIX approach [Farooq16]_. This algorithm uses three\ndifferent optimizers. It starts with a differential evolution algorithm and\nfits the parameters in the power of exponentials. Then the fitted parameters in\nthe first step are utilized to make a linear convex problem. Using a convex\noptimization, the volume fractions are determined. The last step is non-linear\nleast-squares fitting on all the parameters. The results of the first and\nsecond optimizers are utilized as the initial values for the last step of the\nalgorithm.\n\nAs opposed to the `'trr'` fitting method, this approach does not need to set\nany thresholds on the bvals to differentiate between the perfusion\n(pseudo-diffusion) and diffusion portions and fits the parameters\nsimultaneously. Making use of the three step optimization mentioned above\nincreases the convergence basin for fitting the multi-exponential functions of\nmicrostructure models. This method has been described in further detail in\n[Fadnavis19]_ and [Farooq16]_.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ivimmodel_vp = IvimModel(gtab, fit_method='VarPro')\nivimfit_vp = ivimmodel_vp.fit(data_slice)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just like the `'trr'` fit method, `'VarPro'` creates a IvimFit object which\ncontains the parameters of the model obtained after fitting. These are\naccessible through the `model_params` attribute of the IvimFit object.\nThe parameters are arranged as a 4D array, corresponding to the spatial\ndimensions of the data, and the last dimension (of length 4)\ncorresponding to the model parameters according to the following\norder : $\\mathbf{S_{0}, f, D^*, D}$.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "i, j = 10, 10\nestimated_params = ivimfit_vp.model_params[i, j, :]\nprint(estimated_params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compare the fit using `fit_method='VarPro'` and `fit_method='trr'`, we can\nplot one voxel's signal and the model fit using both methods.\n\nWe will use the `predict` method of the IvimFit objects, to get the predicted\nsignal, based on each one of the model fit methods.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1)\n\nax.scatter(gtab.bvals, data_slice[i, j, :],\n color=\"green\", label=\"Measured signal\")\n\nivim_trr_predict = ivimfit.predict(gtab)[i, j, :]\n\nax.plot(gtab.bvals, ivim_trr_predict, label=\"trr prediction\")\n\nS0_est, f_est, D_star_est, D_est = ivimfit.model_params[i, j, :]\n\ntext_fit = \"\"\"trr param estimates: \\n S0={:06.3f} f={:06.4f}\\n\n D*={:06.5f} D={:06.5f}\"\"\".format(S0_est, f_est, D_star_est, D_est)\n\nax.text(0.65, 0.80, text_fit, horizontalalignment='center',\n verticalalignment='center', transform=plt.gca().transAxes)\n\nivim_predict_vp = ivimfit_vp.predict(gtab)[i, j, :]\nax.plot(gtab.bvals, ivim_predict_vp, label=\"VarPro prediction\")\n\nax.set_xlabel(\"bvalues\")\nax.set_ylabel(\"Signals\")\n\nS0_est, f_est, D_star_est, D_est = ivimfit_vp.model_params[i, j, :]\n\ntext_fit = \"\"\"VarPro param estimates: \\n S0={:06.3f} f={:06.4f}\\n\n D*={:06.5f} D={:06.5f}\"\"\".format(S0_est, f_est, D_star_est, D_est)\n\nax.text(0.65, 0.50, text_fit, horizontalalignment='center',\n verticalalignment='center', transform=plt.gca().transAxes)\n\nfig.legend(loc='upper right')\nfig.savefig(\"ivim_voxel_plot.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: ivim_voxel_plot.png\n :align: center\n\n Plot of the signal from one voxel.\n\nLet us get the various plots with `fit_method = 'VarPro'` so that we can\nvisualize them in one page\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nplot_map(ivimfit_vp.S0_predicted, \"Predicted S0\", (0, 10000),\n \"predicted_S0.png\")\nplot_map(data_slice[..., 0], \"Measured S0\", (0, 10000), \"measured_S0.png\")\nplot_map(ivimfit_vp.perfusion_fraction, \"f\", (0, 1), \"perfusion_fraction.png\")\nplot_map(ivimfit_vp.D_star, \"D*\", (0, 0.01), \"perfusion_coeff.png\")\nplot_map(ivimfit_vp.D, \"D\", (0, 0.001), \"diffusion_coeff.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: predicted_S0.png\n :align: center\n\n Heatmap of S0 predicted from the fit\n\n.. figure:: measured_S0.png\n :align: center\n\n Heatmap of measured signal at bvalue = 0.\n\n.. figure:: perfusion_fraction.png\n :align: center\n\n Heatmap of perfusion fraction values predicted from the fit\n\n.. figure:: perfusion_coeff.png\n :align: center\n\n Heatmap of perfusion coefficients predicted from the fit.\n\n.. figure:: diffusion_coeff.png\n :align: center\n\n Heatmap of diffusion coefficients predicted from the fit\n\nReferences:\n\n.. [Stejskal65] Stejskal, E. O.; Tanner, J. E. (1 January 1965).\n \"Spin Diffusion Measurements: Spin Echoes in the Presence\n of a Time-Dependent Field Gradient\". The Journal of Chemical\n Physics 42 (1): 288. Bibcode: 1965JChPh..42..288S.\n doi:10.1063/1.1695690.\n\n.. [LeBihan84] Le Bihan, Denis, et al. \"Separation of diffusion\n and perfusion in intravoxel incoherent motion MR\n imaging.\" Radiology 168.2 (1988): 497-505.\n\n.. [Fadnavis19] Fadnavis, Shreyas et.al. \"MicroLearn: Framework for machine\n learning, reconstruction, optimization and microstructure\n modeling, Proceedings of: International Society of Magnetic\n Resonance in Medicine (ISMRM), Montreal, Canada, 2019.\n\n.. [Farooq16] Farooq, Hamza, et al. \"Microstructure Imaging of Crossing (MIX)\n White Matter Fibers from diffusion MRI.\" Scientific reports 6\n (2016).\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 }