{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Continuous and analytical diffusion signal modelling with MAP-MRI\n\nWe show how to model the diffusion signal as a linear combination of continuous\nfunctions from the MAP-MRI basis [Ozarslan2013]_. This continuous representation\nallows for the computation of many properties of both the signal and diffusion\npropagator.\n\nWe show how to estimate the analytical orientation distribution function (ODF)\nand a variety of scalar indices. These include rotationally invariant quantities\nsuch as the mean squared displacement (MSD), Q-space inverse variance (QIV),\nteturn-to-origin probability (RTOP) and non-Gaussianity (NG). Interestingly, the\nMAP-MRI basis also allows for the computation of directional indices, such as\nthe return-to-axis probability (RTAP), the return-to-plane probability\n(RTPP), and the parallel and perpendicular non-Gaussianity.\n\nThe estimation of these properties from noisy and sparse DWIs requires that the\nfitting of the MAP-MRI basis is constrained and/or regularized. This can be done\nby constraining the diffusion propagator to positive values [Ozarslan2013]_\n[DelaHaije2020]_, and through analytic Laplacian regularization (MAPL)\n[Fick2016a]_.\n\nFirst import the necessary modules:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.reconst import mapmri\nfrom dipy.viz import window, actor\nfrom dipy.data import get_fnames, get_sphere\nfrom dipy.core.gradients import gradient_table\nfrom dipy.io.image import load_nifti\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.viz.plotting import compare_maps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Download and read the data for this tutorial.\n\nMAP-MRI requires multi-shell data, to properly fit the radial part of the basis.\nThe total size of the downloaded data is 187.66 MBytes, however you only need\nto 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')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``data`` contains the voxel data and ``gtab`` contains a ``GradientTable``\nobject (gradient information e.g. b-values). For example, to show the b-values\nit is possible to write::\n\n print(gtab.bvals)\n\nFor the values of the q-space indices to make sense it is necessary to\nexplicitly state the ``big_delta`` and ``small_delta`` parameters in the\ngradient table.\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)\ngtab = gradient_table(bvals, bvecs)\n\nbig_delta = 0.0365 # seconds\nsmall_delta = 0.0157 # seconds\ngtab = gradient_table(bvals=gtab.bvals, bvecs=gtab.bvecs,\n big_delta=big_delta,\n small_delta=small_delta)\n\ndata_small = data[40:65, 50:51]\n\nprint('data.shape (%d, %d, %d, %d)' % data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The MAP-MRI Model can now be instantiated. The ``radial_order`` determines the\nexpansion order of the basis, i.e., how many basis functions are used to\napproximate the signal.\n\nFirst, we must decide to use the anisotropic or isotropic MAP-MRI basis. As was\nshown in [Fick2016a]_, the isotropic basis is best used for tractography\npurposes, as the anisotropic basis has a bias towards smaller crossing angles\nin the ODF. For signal fitting and estimation of scalar quantities the\nanisotropic basis is suggested. The choice can be made by setting\n``anisotropic_scaling=True`` or ``anisotropic_scaling=False``.\n\nFirst, we must select the method of regularization and/or constraining the\nbasis fitting.\n\n- ``laplacian_regularization=True`` makes it use Laplacian regularization\n [Fick2016a]_. This method essentially reduces spurious oscillations in the\n reconstruction by minimizing the Laplacian of the fitted signal. Several\n options can be given to select the regularization weight:\n\n - ``regularization_weighting=GCV`` uses generalized cross-validation\n [Craven1978]_ to find an optimal weight.\n - ``regularization_weighting=0.2`` for example would omit the GCV and\n just set it to 0.2 (found to be reasonable in HCP data [Fick2016a]_).\n - ``regularization_weighting=np.array(weights)`` would make the GCV use\n a custom range to find an optimal weight.\n\n- ``positivity_constraint=True`` makes it use a positivity constraint on the\n diffusion propagator. This method constrains the final solution of the\n diffusion propagator to be positive either globally [DelaHaije2020] or at a\n set of discrete points [Ozarslan2013]_, since negative values should not\n exist.\n\n - Setting ``global_constraints=True`` enforces positivity everywhere.\n With the setting ``global_constraints=False`` positivity is enforced on\n a grid determined by ``pos_grid`` and ``pos_radius``.\n\nBoth methods do a good job of producing viable solutions to the signal fitting\nin practice. The difference is that the Laplacian regularization imposes\nsmoothness over the entire signal, including the extrapolation beyond the\nmeasured signal. In practice this may result in, but does not guarantee,\npositive solutions of the diffusion propagator. The positivity constraint\nguarantees a positive solution which in general results in smooth solutions, but\ndoes not guarantee it.\n\nA suggested strategy is to use a low Laplacian weight together with a positivity\nconstraint. In this way both desired properties are guaranteed in final\nsolution. Higher Laplacian weights are recommended when the data is very noisy.\n\nWe use the package CVXPY (https://www.cvxpy.org/) to solve convex optimization\nproblems when ``positivity_constraint=True``, so we need to first install CVXPY.\nWhen using ``global_constraints=True`` to ensure global positivity, it is\nrecommended to use the MOSEK solver (https://www.mosek.com/) together with CVXPY\nby setting ``cvxpy_solver='MOSEK'``. Different solvers can differ greatly in\nterms of runtime and solution accuracy, and in some cases solvers may show\nwarnings about convergence or recommended option settings.\n\nFor now we will generate the anisotropic models for different combinations.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "radial_order = 6\n\nmap_model_laplacian_aniso = mapmri.MapmriModel(gtab, radial_order=radial_order,\n laplacian_regularization=True,\n laplacian_weighting=.2)\n\nmap_model_positivity_aniso = mapmri.MapmriModel(gtab,\n radial_order=radial_order,\n laplacian_regularization=False,\n positivity_constraint=True)\n\nmap_model_both_aniso = mapmri.MapmriModel(gtab, radial_order=radial_order,\n laplacian_regularization=True,\n laplacian_weighting=.05,\n positivity_constraint=True)\n\nmap_model_plus_aniso = mapmri.MapmriModel(gtab,\n radial_order=radial_order,\n laplacian_regularization=False,\n positivity_constraint=True,\n global_constraints=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that when we use only Laplacian regularization, the ``GCV`` option may\nselect very low regularization weights in very anisotropic white matter such\nas the corpus callosum, resulting in corrupted scalar indices. In this example\nwe therefore choose a preset value of 0.2, which was shown to be quite robust\nand also faster in practice [Fick2016a]_.\n\nWe can then fit the MAP-MRI model to the data:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mapfit_laplacian_aniso = map_model_laplacian_aniso.fit(data_small)\nmapfit_positivity_aniso = map_model_positivity_aniso.fit(data_small)\nmapfit_both_aniso = map_model_both_aniso.fit(data_small)\nmapfit_plus_aniso = map_model_plus_aniso.fit(data_small)\n\nfits = [mapfit_laplacian_aniso, mapfit_positivity_aniso, mapfit_both_aniso,\n mapfit_plus_aniso]\nfit_labels = ['MAPL', 'CMAP', 'CMAPL', 'MAP+']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the fitted models we will first illustrate the estimation of q-space\nindices. We will compare the estimation using only Laplacian regularization\n(MAPL), using only the global positivity constraint (MAP+), using only\npositivity in discrete points (CMAP), or using both Laplacian regularization and\npositivity in discrete points (CMAPL). We first show the RTOP [Ozarslan2013]_.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "compare_maps(fits, maps=['rtop'], fit_labels=fit_labels, map_labels=['RTOP'],\n filename='MAPMRI_rtop.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: MAPMRI_rtop.png\n :align: center\n\nIt can be seen that all maps appear quite smooth and similar, although it is\nclear that the global positivity constraints provide smoother maps than the\ndiscrete constraints. Higher Laplacian weights also produce smoother maps, but\ntend to suppress the estimated RTOP values. The similarity and differences in\nreconstruction can be further illustrated by visualizing the analytic norm of\nthe Laplacian of the fitted signal.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "compare_maps(fits, maps=['norm_of_laplacian_signal'], fit_labels=fit_labels,\n map_labels=['Norm of Laplacian'],\n map_kwargs={'vmin': 0, 'vmax': 3},\n filename='MAPMRI_norm_laplacian.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: MAPMRI_norm_laplacian.png\n :align: center\n\nA high Laplacian norm indicates that the gradient in the three-dimensional\nsignal reconstruction changes a lot -- something that may indicate spurious\noscillations. In the Laplacian reconstruction (left) we see that there are some\nisolated voxels that have a higher norm than the rest. In the positivity\nconstrained reconstruction the norm is already smoother. When both methods are\nused together the overall norm gets smoother still, since both smoothness of the\nsignal and positivity of the propagator are imposed.\n\nFrom now on we will just use the combined approach and the globally constrained\napproach, show all maps we can generate, and explain their significance.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fits = fits[2:]\nfit_labels = fit_labels[2:]\n\ncompare_maps(fits, maps=['msd', 'qiv', 'rtop', 'rtap', 'rtpp'],\n fit_labels=fit_labels,\n map_labels=['MSD', 'QIV', 'RTOP', 'RTAP', 'RTPP'],\n filename='MAPMRI_maps.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: MAPMRI_maps.png\n :align: center\n\nFrom left to right:\n\n- Mean squared displacement (MSD) is a measure for how far protons are able to\n diffuse. a decrease in MSD indicates protons are hindered/restricted more,\n as can be seen by the high MSD in the CSF, but low in the white matter.\n- Q-space Inverse Variance (QIV) is a measure of variance in the signal, which\n is said to have higher contrast to white matter than the MSD\n [Hosseinbor2013]_. QIV has also been shown to have high sensitivity to tissue\n composition change in a simulation study [Fick2016b]_.\n- Return-to-origin probability (RTOP) quantifies the probability that a proton\n will be at the same position at the first and second diffusion gradient\n pulse. A higher RTOP indicates that the volume a spin is inside of is\n smaller, meaning more overall restriction. This is illustrated by the low\n values in CSF and high values in white matter.\n- Return-to-axis probability (RTAP) is a directional index that quantifies\n the probability that a proton will be along the axis of the main eigenvector\n of a diffusion tensor during both diffusion gradient pulses. RTAP has been\n related to the apparent axon diameter [Ozarslan2013]_ [Fick2016a]_ under\n several strong assumptions on the tissue composition and acquisition\n protocol.\n- Return-to-plane probability (RTPP) is a directional index that quantifies the\n probability that a proton will be on the plane perpendicular to the main\n eigenvector of a diffusion tensor during both gradient pulses. RTPP is\n related to the length of a pore [Ozarslan2013]_ but in practice should be\n similar to that of Gaussian diffusion.\n\nIt is also possible to estimate the amount of non-Gaussian diffusion in every\nvoxel [Ozarslan2013]_. This quantity is estimated through the ratio between\nGaussian volume (MAP-MRI's first basis function) and the non-Gaussian volume\n(all other basis functions) of the fitted signal. For this value to be\nphysically meaningful we must use a b-value threshold in the MAP-MRI model. This\nthreshold makes the scale estimation in MAP-MRI use only samples that\nrealistically describe Gaussian diffusion, i.e., at low b-values.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "map_model_both_ng = mapmri.MapmriModel(gtab, radial_order=radial_order,\n laplacian_regularization=True,\n laplacian_weighting=.05,\n positivity_constraint=True,\n bval_threshold=2000)\n\nmap_model_plus_ng = mapmri.MapmriModel(gtab, radial_order=radial_order,\n positivity_constraint=True,\n global_constraints=True,\n bval_threshold=2000)\n\nmapfit_both_ng = map_model_both_ng.fit(data_small)\nmapfit_plus_ng = map_model_plus_ng.fit(data_small)\n\nfits = [mapfit_both_ng, mapfit_plus_ng]\nfit_labels = ['CMAPL', 'MAP+']\n\ncompare_maps(fits, maps=['ng', 'ng_perpendicular', 'ng_parallel'],\n fit_labels=fit_labels,\n map_labels=['NG', 'NG perpendicular', 'NG parallel'],\n filename='MAPMRI_ng.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: MAPMRI_ng.png\n :align: center\n\nOn the left we see the overall NG and on the right the directional perpendicular\nNG and parallel NG. The NG ranges from 1 (completely non-Gaussian) to 0\n(completely Gaussian). The overall NG of a voxel is always higher or equal than\neach of its components. It can be seen that NG has low values in the CSF and\nhigher in the white matter.\n\nIncreases or decreases in these values do not point to a specific\nmicrostructural change, but can indicate clues as to what is happening, similar\nto fractional anisotropy. An initial simulation study that quantifies the added\nvalue of q-space indices over DTI indices is given in [Fick2016b]_.\n\nThe MAP-MRI framework also allows for the estimation of orientation distribution\nfunctions (ODFs). We recommend to use the isotropic implementation for this\npurpose, as the anisotropic implementation has a bias towards smaller crossing\nangles.\n\nFor the isotropic basis we recommend to use a ``radial_order`` of 8, as the\nbasis needs more generic and needs more basis functions to approximate the\nsignal.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "radial_order = 8\nmap_model_both_iso = mapmri.MapmriModel(gtab, radial_order=radial_order,\n laplacian_regularization=True,\n laplacian_weighting=.1,\n positivity_constraint=True,\n anisotropic_scaling=False)\n\nmapfit_both_iso = map_model_both_iso.fit(data_small)" ] }, { "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": [ "Compute the ODFs.\n\nThe radial order ``s`` can be increased to sharpen the results, but it might\nalso make the ODFs noisier. Always check the results visually.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "odf = mapfit_both_iso.odf(sphere, s=2)\nprint('odf.shape (%d, %d, %d, %d)' % odf.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Display the ODFs.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Enables/disables interactive visualization\ninteractive = False\n\nscene = window.Scene()\nsfu = actor.odf_slicer(odf, sphere=sphere, colormap='plasma', scale=0.5)\nsfu.display(y=0)\nsfu.RotateX(-90)\nscene.add(sfu)\nwindow.record(scene, out_path='odfs.png', size=(600, 600))\nif interactive:\n window.show(scene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: odfs.png\n :align: center\n\n Orientation distribution functions (ODFs).\n\n## References\n\n.. [Ozarslan2013] Ozarslan E. et al., \"Mean apparent propagator (MAP) MRI: A\n novel diffusion imaging method for mapping tissue microstructure\",\n NeuroImage, 2013.\n\n.. [Fick2016a] Fick, Rutger HJ, et al. \"MAPL: Tissue microstructure estimation\n using Laplacian-regularized MAP-MRI and its application to HCP data.\"\n NeuroImage (2016).\n\n.. [Craven1978] Craven et al. \"Smoothing Noisy Data with Spline Functions.\"\n NUMER MATH 31.4 (1978): 377-403.\n\n.. [Hosseinbor2013] Hosseinbor et al. \"Bessel fourier orientation\n reconstruction (bfor): an analytical diffusion propagator reconstruction\n for hybrid diffusion imaging and computation of q-space indices. NeuroImage\n 64, 650-670.\n\n.. [Fick2016b] Fick et al. \"A sensitivity analysis of Q-space indices with\n respect to changes in axonal diameter, dispersion and tissue composition.\n ISBI 2016.\n\n.. [DelaHaije2020] 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\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 }