{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Reconstruction with Multi-Shell Multi-Tissue CSD\n\nThis example shows how to use Multi-Shell Multi-Tissue Constrained Spherical\nDeconvolution (MSMT-CSD) introduced by Tournier et al. [Jeurissen2014]_. This\ntutorial goes through the steps involved in implementing the method.\n\nThis method provides improved White Matter(WM), Grey Matter (GM), and\nCerebrospinal fluid (CSF) volume fraction maps, which is otherwise\noverestimated in the standard CSD (SSST-CSD). This is done by using b-value\ndependencies of the different tissue types to estimate ODFs. This method thus\nextends the SSST-CSD introduced in [Tournier2007]_.\n\nThe reconstruction of the fiber orientation distribution function\n(fODF) in MSMT-CSD involves the following steps:\n 1. Generate a mask using Median Otsu (optional step)\n 2. Denoise the data using MP-PCA (optional step)\n 3. Generate Anisotropic Powermap (if T1 unavailable)\n 4. Fit DTI model to the data\n 5. Tissue Classification (needs to be at least two classes of tissues)\n 6. Estimation of the fiber response function\n 7. Use the response function to reconstruct the fODF\n\nFirst, we import all the modules we need from dipy as follows:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport dipy.reconst.shm as shm\nimport dipy.direction.peaks as dp\nimport matplotlib.pyplot as plt\n\nfrom dipy.denoise.localpca import mppca\nfrom dipy.core.gradients import gradient_table, unique_bvals_tolerance\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti\nfrom dipy.segment.mask import median_otsu\nfrom dipy.reconst.mcsd import (auto_response_msmt,\n mask_for_response_msmt,\n response_from_mask_msmt)\nfrom dipy.segment.tissue import TissueClassifierHMRF\nfrom dipy.reconst.mcsd import MultiShellDeconvModel, multi_shell_fiber_response\nfrom dipy.viz import window, actor\n\nfrom dipy.data import get_sphere, get_fnames\nsphere = get_sphere('symmetric724')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this example, we use fetch to download a multi-shell dataset which was\nkindly provided by Hansen and Jespersen (more details about the data are\nprovided in their paper [Hansen2016]_). The total size of the downloaded data\nis 192 MBytes, 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": [ "For the sake of simplicity, we only select two non-zero b-values for this\nexample.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bvals = gtab.bvals\nbvecs = gtab.bvecs\n\nsel_b = np.logical_or(np.logical_or(bvals == 0, bvals == 1000), bvals == 2000)\ndata = data[..., sel_b]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The gradient table is also selected to have the selected b-values (0, 1000 and\n2000)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "gtab = gradient_table(bvals[sel_b], bvecs[sel_b])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We make use of the ``median_otsu`` method to generate the mask for the data as\nfollows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "b0_mask, mask = median_otsu(data, median_radius=2, numpass=1, vol_idx=[0, 1])\n\n\nprint(data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As one can see from its shape, the selected data contains a total of 67\nvolumes of images corresponding to all the diffusion gradient directions\nof the selected b-values and call the ``mppca`` as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "denoised_arr = mppca(data, mask=mask, patch_radius=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will use the denoised array (``denoised_arr``) obtained from ``mppca``\nin the rest of the steps in the tutorial.\n\nAs for the next step, we generate the anisotropic powermap introduced by\n[DellAcqua2014]_. To do so, we make use of the Q-ball Model as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "qball_model = shm.QballModel(gtab, 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We generate the peaks from the ``qball_model`` as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "peaks = dp.peaks_from_model(model=qball_model, data=denoised_arr,\n relative_peak_threshold=.5,\n min_separation_angle=25,\n sphere=sphere, mask=mask)\n\nap = shm.anisotropic_power(peaks.shm_coeff)\n\nplt.matshow(np.rot90(ap[:, :, 10]), cmap=plt.cm.bone)\nplt.savefig(\"anisotropic_power_map.png\")\nplt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: anisotropic_power_map.png\n :align: center\n\n Anisotropic Power Map (Axial Slice)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(ap.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above figure is a visualization of the axial slice of the Anisotropic\nPower Map. It can be treated as a pseudo-T1 for classification purposes\nusing the Hidden Markov Random Fields (HMRF) classifier, if the T1 image\nis not available.\n\nAs we can see from the shape of the Anisotropic Power Map, it is 3D and can be\nused for tissue classification using HMRF. The\nHMRF needs the specification of the number of classes. For the case of MSMT-CSD\nthe ``nclass`` parameter needs to be ``>=2``. In our case, we set it to 3:\nnamely corticospinal fluid (csf), white matter (wm) and gray matter (gm).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nclass = 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, the smoothness factor of the segmentation. Good performance is achieved\nwith values between 0 and 0.5.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "beta = 0.1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then call the ``TissueClassifierHMRF`` with the parameters specified as\nabove:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "hmrf = TissueClassifierHMRF()\ninitial_segmentation, final_segmentation, PVE = hmrf.classify(ap, nclass, beta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we get the tissues segmentation from the final_segmentation.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "csf = np.where(final_segmentation == 1, 1, 0)\ngm = np.where(final_segmentation == 2, 1, 0)\nwm = np.where(final_segmentation == 3, 1, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we want the response function for each of the three tissues and for each\nbvalues. This can be achieved in two different ways. If the case that tissue\nsegmentation is available or that one wants to see the tissue masks used to\ncompute the response functions, a combination of the functions\n``mask_for_response_msmt`` and ``response_from_mask`` is needed.\n\nThe ``mask_for_response_msmt`` function will return a mask of voxels within a\ncuboid ROI and that meet some threshold constraints, for each tissue and bvalue.\nMore precisely, the WM mask must have a FA value above a given threshold. The GM\nmask and CSF mask must have a FA below given thresholds and a MD below other\nthresholds.\n\nNote that for ``mask_for_response_msmt``, the gtab and data should be for\nbvalues under 1200, for optimal tensor fit.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mask_wm, mask_gm, mask_csf = mask_for_response_msmt(gtab, data, roi_radii=10,\n wm_fa_thr=0.7,\n gm_fa_thr=0.3,\n csf_fa_thr=0.15,\n gm_md_thr=0.001,\n csf_md_thr=0.0032)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If one wants to use the previously computed tissue segmentation in addition to\nthe threshold method, it is possible by simply multiplying both masks together.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mask_wm *= wm\nmask_gm *= gm\nmask_csf *= csf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The masks can also be used to calculate the number of voxels for each tissue.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nvoxels_wm = np.sum(mask_wm)\nnvoxels_gm = np.sum(mask_gm)\nnvoxels_csf = np.sum(mask_csf)\n\nprint(nvoxels_wm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, the ``response_from_mask`` function will return the msmt response\nfunctions using precalculated tissue masks.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "response_wm, response_gm, response_csf = response_from_mask_msmt(gtab, data,\n mask_wm,\n mask_gm,\n mask_csf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we can also get directly the response functions by calling the\n``auto_response_msmt`` function, which internally calls\n``mask_for_response_msmt`` followed by ``response_from_mask``. By doing so, we\ndon't have access to the masks and we might have problems with high bvalues\ntensor fit.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "auto_response_wm, auto_response_gm, auto_response_csf = \\\n auto_response_msmt(gtab, data, roi_radii=10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see below, adding the tissue segmentation can change the results\nof the response functions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(\"Responses\")\nprint(response_wm)\nprint(response_gm)\nprint(response_csf)\nprint(\"Auto responses\")\nprint(auto_response_wm)\nprint(auto_response_gm)\nprint(auto_response_csf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point, there are two options on how to use those response functions. We\nwant to create a MultiShellDeconvModel, which takes a response function as\ninput. This response function can either be directly in the current format, or\nit can be a MultiShellResponse format, as produced by the\n``multi_shell_fiber_response`` method. This function assumes a 3 compartments\nmodel (wm, gm, csf) and takes one response function per tissue per bvalue. It is\nimportant to note that the bvalues must be unique for this function.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ubvals = unique_bvals_tolerance(gtab.bvals)\nresponse_mcsd = multi_shell_fiber_response(sh_order=8,\n bvals=ubvals,\n wm_rf=response_wm,\n gm_rf=response_gm,\n csf_rf=response_csf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As mentioned, we can also build the model directly and it will call\n``multi_shell_fiber_response`` internally. Important note here, the function\n``unique_bvals_tolerance`` is used to keep only unique bvalues from the gtab\ngiven to the model, as input for ``multi_shell_fiber_response``. This may\nintroduce differences between the calculated response of each method, depending\non the bvalues given to ``multi_shell_fiber_response`` externally.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "response = np.array([response_wm, response_gm, response_csf])\nmcsd_model_simple_response = MultiShellDeconvModel(gtab, response, sh_order=8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this technique only works for a 3 compartments model (wm, gm, csf). If\none wants more compartments, a custom MultiShellResponse object must be used. It\ncan be inspired by the ``multi_shell_fiber_response`` method.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we build the MSMT-CSD model with the ``response_mcsd`` as input. We then\ncall the ``fit`` function to fit one slice of the 3D data and visualize it.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mcsd_model = MultiShellDeconvModel(gtab, response_mcsd)\nmcsd_fit = mcsd_model.fit(denoised_arr[:, :, 10:11])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The volume fractions of tissues for each voxel are also accessible, as well as\nthe sh coefficients for all tissues. One can also get each sh tissue separately\nusing ``all_shm_coeff`` for each compartment (isotropic) and\n``shm_coeff`` for white matter.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "vf = mcsd_fit.volume_fractions\nsh_coeff = mcsd_fit.all_shm_coeff\ncsf_sh_coeff = sh_coeff[..., 0]\ngm_sh_coeff = sh_coeff[..., 1]\nwm_sh_coeff = mcsd_fit.shm_coeff" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model allows to predict a signal from sh coefficients. There are two ways of\ndoing this.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mcsd_pred = mcsd_fit.predict()\nmcsd_pred = mcsd_model.predict(mcsd_fit.all_shm_coeff)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the fit obtained in the previous step, we generate the ODFs which can be\nvisualized as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mcsd_odf = mcsd_fit.odf(sphere)\n\nprint(\"ODF\")\nprint(mcsd_odf.shape)\nprint(mcsd_odf[40, 40, 0])\n\nfodf_spheres = actor.odf_slicer(mcsd_odf, sphere=sphere, scale=1,\n norm=False, colormap='plasma')\n\ninteractive = False\nscene = window.Scene()\nscene.add(fodf_spheres)\nscene.reset_camera_tight()\n\nprint('Saving illustration as msdodf.png')\nwindow.record(scene, out_path='msdodf.png', size=(600, 600))\n\nif interactive:\n window.show(scene)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: msdodf.png\n :align: center\n\n MSMT-CSD Peaks and ODFs.\n\n## References\n\n.. [Jeurissen2014] B. Jeurissen, et al., \"Multi-tissue constrained spherical\n deconvolution for improved analysis of multi-shell\n diffusion MRI data.\" NeuroImage 103 (2014): 411-426.\n\n.. [Tournier2007] J-D. Tournier, F. Calamante and A. Connelly, \"Robust\n determination of the fibre orientation distribution in\n diffusion MRI: Non-negativity constrained super-resolved\n spherical deconvolution\", Neuroimage, vol. 35, no. 4,\n pp. 1459-1472, (2007).\n\n.. [Hansen2016] B. Hansen and SN. Jespersen, \" Data for evaluation of fast\n kurtosis strategies, b-value optimization and exploration\n of diffusion MRI contrast\", Scientific Data 3: 160072\n doi:10.1038/sdata.2016.72, (2016)\n\n.. [DellAcqua2014] F. Dell'Acqua, et. al., \"Anisotropic Power Maps: A\n diffusion contrast to reveal low anisotropy tissues from\n HARDI data\", Proceedings of International Society for\n Magnetic Resonance in Medicine. Milan, Italy, (2014).\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 }