1""" 2 3.. _reconst-mcsd: 4 5================================================ 6Reconstruction with Multi-Shell Multi-Tissue CSD 7================================================ 8 9This example shows how to use Multi-Shell Multi-Tissue Constrained Spherical 10Deconvolution (MSMT-CSD) introduced by Tournier et al. [Jeurissen2014]_. This 11tutorial goes through the steps involved in implementing the method. 12 13This method provides improved White Matter(WM), Grey Matter (GM), and 14Cerebrospinal fluid (CSF) volume fraction maps, which is otherwise 15overestimated in the standard CSD (SSST-CSD). This is done by using b-value 16dependencies of the different tissue types to estimate ODFs. This method thus 17extends the SSST-CSD introduced in [Tournier2007]_. 18 19The reconstruction of the fiber orientation distribution function 20(fODF) in MSMT-CSD involves the following steps: 21 1. Generate a mask using Median Otsu (optional step) 22 2. Denoise the data using MP-PCA (optional step) 23 3. Generate Anisotropic Powermap (if T1 unavailable) 24 4. Fit DTI model to the data 25 5. Tissue Classification (needs to be at least two classes of tissues) 26 6. Estimation of the fiber response function 27 7. Use the response function to reconstruct the fODF 28 29First, we import all the modules we need from dipy as follows: 30""" 31 32import numpy as np 33import dipy.reconst.shm as shm 34import dipy.direction.peaks as dp 35import dipy.reconst.dti as dti 36import matplotlib.pyplot as plt 37 38from dipy.denoise.localpca import mppca 39from dipy.core.gradients import gradient_table, unique_bvals_tolerance 40from dipy.io.gradients import read_bvals_bvecs 41from dipy.io.image import load_nifti 42from dipy.segment.mask import median_otsu 43from dipy.reconst.csdeconv import auto_response_ssst 44from dipy.reconst.mcsd import (auto_response_msmt, 45 mask_for_response_msmt, 46 response_from_mask_msmt) 47from dipy.segment.tissue import TissueClassifierHMRF 48from dipy.reconst.mcsd import MultiShellDeconvModel, multi_shell_fiber_response 49from dipy.viz import window, actor 50 51from dipy.direction.peaks import peaks_from_model 52 53from dipy.data import get_sphere, get_fnames 54sphere = get_sphere('symmetric724') 55 56""" 57For this example, we use fetch to download a multi-shell dataset which was 58kindly provided by Hansen and Jespersen (more details about the data are 59provided in their paper [Hansen2016]_). The total size of the downloaded data 60is 192 MBytes, however you only need to fetch it once. 61""" 62 63fraw, fbval, fbvec, t1_fname = get_fnames('cfin_multib') 64 65data, affine = load_nifti(fraw) 66bvals, bvecs = read_bvals_bvecs(fbval, fbvec) 67gtab = gradient_table(bvals, bvecs) 68 69""" 70For the sake of simplicity, we only select two non-zero b-values for this 71example. 72""" 73 74bvals = gtab.bvals 75bvecs = gtab.bvecs 76 77sel_b = np.logical_or(np.logical_or(bvals == 0, bvals == 1000), bvals == 2000) 78data = data[..., sel_b] 79 80""" 81The gradient table is also selected to have the selected b-values (0, 1000 and 822000) 83""" 84 85gtab = gradient_table(bvals[sel_b], bvecs[sel_b]) 86 87""" 88We make use of the ``median_otsu`` method to generate the mask for the data as 89follows: 90""" 91 92b0_mask, mask = median_otsu(data, median_radius=2, numpass=1, vol_idx=[0, 1]) 93 94 95print(data.shape) 96 97""" 98As one can see from its shape, the selected data contains a total of 67 99volumes of images corresponding to all the diffusion gradient directions 100of the selected b-values and call the ``mppca`` as follows: 101""" 102 103denoised_arr = mppca(data, mask=mask, patch_radius=2) 104 105""" 106Now we will use the denoised array (``denoised_arr``) obtained from ``mppca`` 107in the rest of the steps in the tutorial. 108 109As for the next step, we generate the anisotropic powermap introduced by 110[DellAcqua2014]_. To do so, we make use of the Q-ball Model as follows: 111""" 112 113qball_model = shm.QballModel(gtab, 8) 114 115""" 116We generate the peaks from the ``qball_model`` as follows: 117""" 118 119peaks = dp.peaks_from_model(model=qball_model, data=denoised_arr, 120 relative_peak_threshold=.5, 121 min_separation_angle=25, 122 sphere=sphere, mask=mask) 123 124ap = shm.anisotropic_power(peaks.shm_coeff) 125 126plt.matshow(np.rot90(ap[:, :, 10]), cmap=plt.cm.bone) 127plt.savefig("anisotropic_power_map.png") 128plt.close() 129 130""" 131.. figure:: anisotropic_power_map.png 132 :align: center 133 134 Anisotropic Power Map (Axial Slice) 135""" 136 137print(ap.shape) 138 139""" 140The above figure is a visualization of the axial slice of the Anisotropic 141Power Map. It can be treated as a pseudo-T1 for classification purposes 142using the Hidden Markov Random Fields (HMRF) classifier, if the T1 image 143is not available. 144 145As we can see from the shape of the Anisotropic Power Map, it is 3D and can be 146used for tissue classification using HMRF. The 147HMRF needs the specification of the number of classes. For the case of MSMT-CSD 148the ``nclass`` parameter needs to be ``>=2``. In our case, we set it to 3: 149namely corticospinal fluid (csf), white matter (wm) and gray matter (gm). 150""" 151 152nclass = 3 153 154""" 155Then, the smoothness factor of the segmentation. Good performance is achieved 156with values between 0 and 0.5. 157""" 158 159 160beta = 0.1 161 162 163""" 164We then call the ``TissueClassifierHMRF`` with the parameters specified as 165above: 166""" 167 168hmrf = TissueClassifierHMRF() 169initial_segmentation, final_segmentation, PVE = hmrf.classify(ap, nclass, beta) 170 171 172""" 173Then, we get the tissues segmentation from the final_segmentation. 174""" 175 176csf = np.where(final_segmentation == 1, 1, 0) 177gm = np.where(final_segmentation == 2, 1, 0) 178wm = np.where(final_segmentation == 3, 1, 0) 179 180 181""" 182Now, we want the response function for each of the three tissues and for each 183bvalues. This can be achieved in two different ways. If the case that tissue 184segmentation is available or that one wants to see the tissue masks used to 185compute the response functions, a combination of the functions 186``mask_for_response_msmt`` and ``response_from_mask`` is needed. 187 188The ``mask_for_response_msmt`` function will return a mask of voxels within a 189cuboid ROI and that meet some threshold constraints, for each tissue and bvalue. 190More precisely, the WM mask must have a FA value above a given threshold. The GM 191mask and CSF mask must have a FA below given thresholds and a MD below other 192thresholds. 193 194Note that for ``mask_for_response_msmt``, the gtab and data should be for 195bvalues under 1200, for optimal tensor fit. 196""" 197 198mask_wm, mask_gm, mask_csf = mask_for_response_msmt(gtab, data, roi_radii=10, 199 wm_fa_thr=0.7, 200 gm_fa_thr=0.3, 201 csf_fa_thr=0.15, 202 gm_md_thr=0.001, 203 csf_md_thr=0.0032) 204 205""" 206If one wants to use the previously computed tissue segmentation in addition to 207the threshold method, it is possible by simply multiplying both masks together. 208""" 209 210mask_wm *= wm 211mask_gm *= gm 212mask_csf *= csf 213 214""" 215The masks can also be used to calculate the number of voxels for each tissue. 216""" 217 218nvoxels_wm = np.sum(mask_wm) 219nvoxels_gm = np.sum(mask_gm) 220nvoxels_csf = np.sum(mask_csf) 221 222print(nvoxels_wm) 223 224""" 225Then, the ``response_from_mask`` function will return the msmt response 226functions using precalculated tissue masks. 227""" 228 229response_wm, response_gm, response_csf = response_from_mask_msmt(gtab, data, 230 mask_wm, 231 mask_gm, 232 mask_csf) 233 234""" 235Note that we can also get directly the response functions by calling the 236``auto_response_msmt`` function, which internally calls 237``mask_for_response_msmt`` followed by ``response_from_mask``. By doing so, we 238don't have access to the masks and we might have problems with high bvalues 239tensor fit. 240""" 241 242auto_response_wm, auto_response_gm, auto_response_csf = \ 243 auto_response_msmt(gtab, data, roi_radii=10) 244 245""" 246As we can see below, adding the tissue segmentation can change the results 247of the response functions. 248""" 249 250print("Responses") 251print(response_wm) 252print(response_gm) 253print(response_csf) 254print("Auto responses") 255print(auto_response_wm) 256print(auto_response_gm) 257print(auto_response_csf) 258 259 260""" 261At this point, there are two options on how to use those response functions. We 262want to create a MultiShellDeconvModel, which takes a response function as 263input. This response function can either be directly in the current format, or 264it can be a MultiShellResponse format, as produced by the 265``multi_shell_fiber_response`` method. This function assumes a 3 compartments 266model (wm, gm, csf) and takes one response function per tissue per bvalue. It is 267important to note that the bvalues must be unique for this function. 268""" 269 270ubvals = unique_bvals_tolerance(gtab.bvals) 271response_mcsd = multi_shell_fiber_response(sh_order=8, 272 bvals=ubvals, 273 wm_rf=response_wm, 274 gm_rf=response_gm, 275 csf_rf=response_csf) 276 277""" 278As mentionned, we can also build the model directly and it will call 279``multi_shell_fiber_response`` internally. Important note here, the function 280``unique_bvals_tolerance`` is used to keep only unique bvalues from the gtab 281given to the model, as input for ``multi_shell_fiber_response``. This may 282introduce differences between the calculted response of each method, depending 283on the bvalues given to ``multi_shell_fiber_response`` externally. 284""" 285 286response = np.array([response_wm, response_gm, response_csf]) 287mcsd_model_simple_response = MultiShellDeconvModel(gtab, response, sh_order=8) 288 289""" 290Note that this technique only works for a 3 compartments model (wm, gm, csf). If 291one wants more compartments, a custom MultiShellResponse object must be used. It 292can be inspired by the ``multi_shell_fiber_response`` method. 293""" 294 295""" 296Now we build the MSMT-CSD model with the ``response_mcsd`` as input. We then 297call the ``fit`` function to fit one slice of the 3D data and visualize it. 298""" 299 300mcsd_model = MultiShellDeconvModel(gtab, response_mcsd) 301mcsd_fit = mcsd_model.fit(denoised_arr[:, :, 10:11]) 302 303""" 304The volume fractions of tissues for each voxel are also accessible, as well as 305the sh coefficients for all tissues. One can also get each sh tissue separately 306using ``all_shm_coeff`` for each compartment (isotropic) and 307``shm_coeff`` for white matter. 308""" 309 310vf = mcsd_fit.volume_fractions 311sh_coeff = mcsd_fit.all_shm_coeff 312csf_sh_coeff = sh_coeff[..., 0] 313gm_sh_coeff = sh_coeff[..., 1] 314wm_sh_coeff = mcsd_fit.shm_coeff 315 316""" 317The model allows to predict a signal from sh coefficients. There are two ways of 318doing this. 319""" 320 321mcsd_pred = mcsd_fit.predict() 322mcsd_pred = mcsd_model.predict(mcsd_fit.all_shm_coeff) 323 324""" 325From the fit obtained in the previous step, we generate the ODFs which can be 326visualized as follows: 327""" 328 329mcsd_odf = mcsd_fit.odf(sphere) 330 331print("ODF") 332print(mcsd_odf.shape) 333print(mcsd_odf[40, 40, 0]) 334 335fodf_spheres = actor.odf_slicer(mcsd_odf, sphere=sphere, scale=1, 336 norm=False, colormap='plasma') 337 338interactive = False 339scene = window.Scene() 340scene.add(fodf_spheres) 341scene.reset_camera_tight() 342 343print('Saving illustration as msdodf.png') 344window.record(scene, out_path='msdodf.png', size=(600, 600)) 345 346if interactive: 347 window.show(scene) 348 349""" 350.. figure:: msdodf.png 351 :align: center 352 353 MSMT-CSD Peaks and ODFs. 354 355References 356---------- 357 358.. [Jeurissen2014] B. Jeurissen, et al., "Multi-tissue constrained spherical 359 deconvolution for improved analysis of multi-shell 360 diffusion MRI data." NeuroImage 103 (2014): 411-426. 361 362.. [Tournier2007] J-D. Tournier, F. Calamante and A. Connelly, "Robust 363 determination of the fibre orientation distribution in 364 diffusion MRI: Non-negativity constrained super-resolved 365 spherical deconvolution", Neuroimage, vol. 35, no. 4, 366 pp. 1459-1472, (2007). 367 368.. [Hansen2016] B. Hansen and SN. Jespersen, " Data for evaluation of fast 369 kurtosis strategies, b-value optimization and exploration 370 of diffusion MRI contrast", Scientific Data 3: 160072 371 doi:10.1038/sdata.2016.72, (2016) 372 373.. [DellAcqua2014] F. Dell'Acqua, et. al., "Anisotropic Power Maps: A 374 diffusion contrast to reveal low anisotropy tissues from 375 HARDI data", Proceedings of International Society for 376 Magnetic Resonance in Medicine. Milan, Italy, (2014). 377 378.. include:: ../links_names.inc 379 380""" 381