{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Reconstruction of the diffusion signal with the kurtosis tensor model\n\nThe diffusion kurtosis model is an expansion of the diffusion tensor model\n(see `example_reconst_dti`). In addition to the diffusion tensor (DT), the\ndiffusion kurtosis model quantifies the degree to which water diffusion in\nbiological tissues is non-Gaussian using the kurtosis tensor (KT)\n[Jensen2005]_.\n\nMeasurements of non-Gaussian diffusion from the diffusion kurtosis model are of\ninterest because they can be used to characterize tissue microstructural\nheterogeneity [Jensen2010]_. Moreover, DKI can be used to: 1) derive concrete\nbiophysical parameters, such as the density of axonal fibers and diffusion\ntortuosity [Fierem2011]_ (see `example_reconst_dki_micro`); and 2)\nresolve crossing fibers in tractography and to obtain invariant rotational\nmeasures not limited to well-aligned fiber populations [NetoHe2015]_.\n\nThe diffusion kurtosis model expresses the diffusion-weighted signal as:\n\n\\begin{align}S(n,b)=S_{0}e^{-bD(n)+\frac{1}{6}b^{2}D(n)^{2}K(n)}\\end{align}\n\nwhere $\\mathbf{b}$ is the applied diffusion weighting (which is dependent on\nthe measurement parameters), $S_0$ is the signal in the absence of diffusion\ngradient sensitization, $\\mathbf{D(n)}$ is the value of diffusion along\ndirection $\\mathbf{n}$, and $\\mathbf{K(n)}$ is the value of kurtosis along\ndirection $\\mathbf{n}$. The directional diffusion $\\mathbf{D(n)}$ and kurtosis\n$\\mathbf{K(n)}$ can be related to the diffusion tensor (DT) and kurtosis tensor\n(KT) using the following equations:\n\n\\begin{align}D(n)=\\sum_{i=1}^{3}\\sum_{j=1}^{3}n_{i}n_{j}D_{ij}\\end{align}\n\nand\n\n\\begin{align}K(n)=\frac{MD^{2}}{D(n)^{2}}\\sum_{i=1}^{3}\\sum_{j=1}^{3}\\sum_{k=1}^{3}\n \\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl}\\end{align}\n\nwhere $D_{ij}$ are the elements of the second-order DT, and $W_{ijkl}$ the\nelements of the fourth-order KT and $MD$ is the mean diffusivity. As the DT,\nKT has antipodal symmetry and thus only 15 Wijkl elements are needed to fully\ncharacterize the KT:\n\n\\begin{align}\begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz}\n & ... \\\n & W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy}\n & ... \\\n & W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz}\n & & )\\end{matrix}\\end{align}\n\nIn the following example we show how to fit the diffusion kurtosis model on\ndiffusion-weighted multi-shell datasets and how to estimate diffusion kurtosis\nbased statistics.\n\nFirst, we import all relevant modules:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport dipy.reconst.dki as dki\nimport dipy.reconst.dti as dti\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.segment.mask import median_otsu\nfrom dipy.viz.plotting import compare_maps\nfrom scipy.ndimage import gaussian_filter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "DKI requires multi-shell data, i.e. data acquired from more than one non-zero\nb-value. Here, we use fetch to download a multi-shell dataset which was kindly\nprovided by Hansen and Jespersen (more details about the data are provided in\ntheir paper [Hansen2016]_). The total size of the downloaded data is 192\nMBytes, however you only need to fetch it once.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fraw, fbval, fbvec, t1_fname = get_fnames('cfin_multib')\n\ndata, affine = load_nifti(fraw)\nbvals, bvecs = read_bvals_bvecs(fbval, fbvec)\ngtab = gradient_table(bvals, bvecs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Function ``get_fnames`` downloads and outputs the paths of the data,\n``load_nifti`` returns the data as a nibabel Nifti1Image object, and\n``read_bvals_bvecs`` loads the arrays containing the information about the\nb-values and b-vectors. These later arrays are converted to the GradientTable\nobject required for Dipy_'s data reconstruction.\n\nBefore fitting the data, we perform some data pre-processing. We first compute\na brain mask to avoid unnecessary calculations on the background of the image.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "maskdata, mask = median_otsu(data, vol_idx=[0, 1], median_radius=4, numpass=2,\n autocrop=False, dilate=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the diffusion kurtosis models involves the estimation of a large number\nof parameters [TaxCMW2015]_ and since the non-Gaussian components of the\ndiffusion signal are more sensitive to artefacts [NetoHe2012]_, it might be\nfavorable to suppress the effects of noise and artefacts before diffusion\nkurtosis fitting. In this example the effects of noise and artefacts are\nsuppress by using 3D Gaussian smoothing (with a Gaussian kernel with\nfwhm=1.25) as suggested by pioneer DKI studies (e.g. [Jensen2005]_,\n[NetoHe2012]_). Although here the Gaussian smoothing is used so that results\nare comparable to these studies, it is important to note that more advanced\nnoise and artifact suppression algorithms are available in DIPY_, e.g. the\nMarcenko-Pastur PCA denoising algorithm (`example-denoise-mppca`) and\nthe Gibbs artefact suppression algorithm (`example-denoise-gibbs`).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fwhm = 1.25\ngauss_std = fwhm / np.sqrt(8 * np.log(2)) # converting fwhm to Gaussian std\ndata_smooth = np.zeros(data.shape)\nfor v in range(data.shape[-1]):\n data_smooth[..., v] = gaussian_filter(data[..., v], sigma=gauss_std)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have loaded and pre-processed the data we can go forward\nwith DKI fitting. For this, the DKI model is first defined for the data's\nGradientTable object by instantiating the DiffusionKurtosisModel object in the\nfollowing way:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dkimodel = dki.DiffusionKurtosisModel(gtab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To fit the data using the defined model object, we call the ``fit`` function of\nthis object. For the purpose of this example, we will only fit a single slice of\nthe data:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data_smooth = data_smooth[:, :, 9:10]\nmask = mask[:, :, 9:10]\ndkifit = dkimodel.fit(data_smooth, mask=mask)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fit method creates a DiffusionKurtosisFit object, which contains all the\ndiffusion and kurtosis fitting parameters and other DKI attributes. For\ninstance, since the diffusion kurtosis model estimates the diffusion tensor,\nall standard diffusion tensor statistics can be computed from the\nDiffusionKurtosisFit instance. For example, we can extract the fractional\nanisotropy (FA), the mean diffusivity (MD), the axial diffusivity (AD) and the\nradial diffusivity (RD) from the DiffusionKurtosisiFit instance. Of course,\nthese measures can also be computed from DIPY's ``TensorModel`` fit, and should\nbe analogous; however, theoretically, the diffusion statistics from the kurtosis\nmodel are expected to have better accuracy, since DKI's diffusion tensor are\ndecoupled from higher order terms effects [Veraar2011]_, [NetoHe2012]_. Below we\ncompare the FA, MD, AD, and RD, computed from both DTI and DKI.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tenmodel = dti.TensorModel(gtab)\ntenfit = tenmodel.fit(data_smooth, mask=mask)\n\nfits = [tenfit, dkifit]\nmaps = ['fa', 'md', 'ad', 'rd']\nfit_labels = ['DTI', 'DKI']\nmap_kwargs = [{'vmax': 0.7}, {'vmax': 2e-3}, {'vmax': 2e-3}, {'vmax': 2e-3}]\ncompare_maps(fits, maps, fit_labels=fit_labels, map_kwargs=map_kwargs,\n filename='Diffusion_tensor_measures_from_DTI_and_DKI.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: Diffusion_tensor_measures_from_DTI_and_DKI.png\n :align: center\n\n Diffusion tensor measures obtained from the diffusion tensor estimated\n from DKI (upper panels) and DTI (lower panels).\n\nDTI's diffusion estimates present lower values than DKI's estimates,\nshowing that DTI's diffusion measurements are underestimated by higher order\neffects.\n\nIn addition to the standard diffusion statistics, the DiffusionKurtosisFit\ninstance can be used to estimate the non-Gaussian measures of mean kurtosis\n(MK), the axial kurtosis (AK) and the radial kurtosis (RK).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "maps = ['mk', 'ak', 'rk']\ncompare_maps([dkifit], maps, fit_labels=['DKI'],\n map_kwargs={'vmin': 0, 'vmax': 1.5},\n filename='Kurtosis_tensor_standard_measures.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: Kurtosis_tensor_standard_measures.png\n :align: center\n\n DKI standard kurtosis measures.\n\nThe non-Gaussian behaviour of the diffusion signal is expected to be higher\nwhen tissue water is confined by multiple compartments. MK is, therefore,\nhigher in white matter since it is highly compartmentalized by myelin sheaths.\nThese water diffusion compartmentalization is expected to be more pronounced\nperpendicularly to white matter fibers and thus the RK map presents higher\namplitudes than the AK map.\n\nIt is important to note that kurtosis estimates might present negative estimates\nin deep white matter regions (e.g. the band of dark voxels in the RK map above).\nThese negative kurtosis values are artefactual and might be induced by:\n1) low radial diffusivities of aligned white matter - since it is very hard\nto capture non-Gaussian information in radial direction due to it's low\ndiffusion decays, radial kurtosis estimates (and consequently the mean\nkurtosis estimates) might have low robustness and tendency to exhibit negative\nvalues [NetoHe2012]_;\n2) Gibbs artefacts - MRI images might be corrupted by signal oscillation\nartefact between tissue's edges if an inadequate number of high frequencies of\nthe k-space is sampled. These oscillations might have different signs on\nimages acquired with different diffusion-weighted and inducing negative biases\nin kurtosis parametric maps [Perron2015]_, [NetoHe2018]_.\n\nOne can try to suppress this issue by using the more advance noise and artefact\nsuppression algorithms, e.g., as mentioned above, the MP-PCA denoising\n(`example-denoise-mppca`) and Gibbs Unringing\n(`example-denoise-gibbs`) algorithms. Alternatively, one can overcome this\nartefact by computing the kurtosis values from powder-averaged\ndiffusion-weighted signals. The details on how to compute the kurtosis from\npowder-average signals in dipy are described in follow the tutorial\n(`example-reconst-msdki`). Finally, one can use constrained optimization to\nensure that the fitted parameters are physically plausible [DelaHa2020]_, as we\nwill illustrate in the next section. Ideally though, artefacts such as Gibbs\nringing should be corrected for as well as possible before using constrained\noptimization.\n\n## Constrained optimization for DKI\n\nWhen instantiating the DiffusionKurtosisModel, the model can be set up to use\nconstraints with the option `fit_method='CLS'` (for ordinary least squares) or\nwith `fit_method='CWLS'` (for weighted least squares). Constrained fitting takes\nmore time than unconstrained fitting, but is generally recommended to prevent\nphysically unplausible parameter estimates [DelaHa2020]_. For performance\npurposes it is recommended to use the MOSEK solver (https://www.mosek.com/) by\nsetting ``cvxpy_solver='MOSEK'``. Different solvers can differ greatly in terms\nof runtime and solution accuracy, and in some cases solvers may show warnings\nabout convergence or recommended option settings.\n\n

Note

In certain atypical scenarios, the DKI+ constraints could potentially be\n too restrictive. Always check the results of a constrained fit with their\n unconstrained counterpart to verify that there are no unexpected\n qualitative differences.

\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dkimodel_plus = dki.DiffusionKurtosisModel(gtab, fit_method='CLS')\ndkifit_plus = dkimodel_plus.fit(data_smooth, mask=mask)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now compare the kurtosis measures obtained with the constrained fit to\nthe measures obtained before, where we see that many of the artefactual voxels\nhave now been corrected. In particular outliers caused by pure noise -- instead\nof for example acquisition artefacts -- can be corrected with this method.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "compare_maps([dkifit_plus], ['mk', 'ak', 'rk'], fit_labels=['DKI+'],\n filename='Kurtosis_tensor_standard_measures_plus.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: Kurtosis_tensor_standard_measures_plus.png\n :align: center\n\n DKI standard kurtosis measures obtained with constrained optimization.\n\nWhen using constrained optimization, the expected range of the kurtosis measures\nis also naturally constrained, and so does not typically require additional\nclipping.\n\nFinally, constrained optimization obviates the need for smoothing in many cases:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dkifit_noisy = dkimodel.fit(data[:, :, 9:10], mask=mask)\ndkifit_noisy_plus = dkimodel_plus.fit(data[:, :, 9:10], mask=mask)\n\ncompare_maps([dkifit_noisy, dkifit_noisy_plus], ['mk', 'ak', 'rk'],\n fit_labels=['DKI', 'DKI+'], map_kwargs={'vmin': 0, 'vmax': 1.5},\n filename='Kurtosis_tensor_standard_measures_noisy.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: Kurtosis_tensor_standard_measures_noisy.png\n :align: center\n\n DKI standard kurtosis measures obtained on unsmoothed data with constrained\n optimization.\n\n## Mean kurtosis tensor and kurtosis fractional anisotropy\n\nAs pointed by previous studies [Hansen2013]_, axial, radial and mean kurtosis\ndepends on the information of both diffusion and kurtosis tensor. DKI measures\nthat only depend on the kurtosis tensor include the mean of the kurtosis tensor\n[Hansen2013]_, and the kurtosis fractional anisotropy [GlennR2015]_. These\nmeasures are computed and illustrated below:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "compare_maps([dkifit_plus], ['mkt', 'kfa'], fit_labels=['DKI+'],\n map_kwargs=[{'vmin': 0, 'vmax': 1.5}, {'vmin': 0, 'vmax': 1}],\n filename='Measures_from_kurtosis_tensor_only.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: Measures_from_kurtosis_tensor_only.png\n :align: center\n\n DKI measures obtained from the kurtosis tensor only.\n\nAs reported by [Hansen2013]_, the mean of the kurtosis tensor (MKT) produces\nsimilar maps than the standard mean kurtosis (MK). On the other hand,\nthe kurtosis fractional anisotropy (KFA) maps shows that the kurtosis tensor\nhave different degrees of anisotropy than the FA measures from the diffusion\ntensor.\n\n### References\n.. [Jensen2005] Jensen JH, Helpern JA, Ramani A, Lu H, Kaczynski K (2005).\n Diffusional Kurtosis Imaging: The Quantification of\n Non_Gaussian Water Diffusion by Means of Magnetic Resonance\n Imaging. Magnetic Resonance in Medicine 53: 1432-1440\n.. [Jensen2010] Jensen JH, Helpern JA (2010). MRI quantification of\n non-Gaussian water diffusion by kurtosis analysis. NMR in\n Biomedicine 23(7): 698-710\n.. [Fierem2011] Fieremans E, Jensen JH, Helpern JA (2011). White matter\n characterization with diffusion kurtosis imaging. NeuroImage\n 58: 177-188\n.. [Veraar2011] Veraart J, Poot DH, Van Hecke W, Blockx I, Van der Linden A,\n Verhoye M, Sijbers J (2011). More Accurate Estimation of\n Diffusion Tensor Parameters Using Diffusion Kurtosis Imaging.\n Magnetic Resonance in Medicine 65(1): 138-145\n.. [NetoHe2012] Neto Henriques R, Ferreira H, Correia M, (2012). Diffusion\n kurtosis imaging of the healthy human brain. Master\n Dissertation Bachelor and Master Programin Biomedical\n Engineering and Biophysics, Faculty of Sciences.\n http://repositorio.ul.pt/bitstream/10451/8511/1/ulfc104137_tm_Rafael_Henriques.pdf\n.. [Hansen2013] Hansen B, Lund TE, Sangill R, and Jespersen SN (2013).\n Experimentally and computationally393fast method for estimation\n of a mean kurtosis. Magnetic Resonance in Medicine 69,\n 1754\u20131760.394doi:10.1002/mrm.24743\n.. [GlennR2015] Glenn GR, Helpern JA, Tabesh A, Jensen JH (2015).\n Quantitative assessment of diffusional387kurtosis anisotropy.\n NMR in Biomedicine28, 448\u2013459. doi:10.1002/nbm.3271\n.. [NetoHe2015] Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015).\n Exploring the 3D geometry of the diffusion kurtosis tensor -\n Impact on the development of robust tractography procedures and\n novel biomarkers, NeuroImage 111: 85-99\n.. [Perron2015] Perrone D, Aelterman J, Piz\u030curica A, Jeurissen B, Philips W,\n Leemans A, (2015). The effect of Gibbs ringing artifacts on\n measures derived from diffusion MRI. Neuroimage 120, 441-455.\n https://doi.org/10.1016/j.neuroimage.2015.06.068.\n.. [TaxCMW2015] Tax CMW, Otte WM, Viergever MA, Dijkhuizen RM, Leemans A\n (2014). REKINDLE: Robust extraction of kurtosis INDices with\n linear estimation. Magnetic Resonance in Medicine 73(2):\n 794-808.\n.. [Hansen2016] Hansen, B, Jespersen, SN (2016). Data for evaluation of fast\n kurtosis strategies, b-value optimization and exploration of\n diffusion MRI contrast. Scientific Data 3: 160072\n doi:10.1038/sdata.2016.72\n.. [NetoHe2018] Neto Henriques R (2018). Advanced Methods for Diffusion MRI\n Data Analysis and their Application to the Healthy Ageing Brain\n (Doctoral thesis). https://doi.org/10.17863/CAM.29356\n.. [DelaHa2020] Dela Haije et al. \"Enforcing necessary non-negativity\n constraints for common diffusion MRI models using sum of squares\n programming\". NeuroImage 209, 2020, 116405.\n\n.. include:: ../links_names.inc\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 }