1""" 2============================================================================= 3Using the free water elimination model to remove DTI free water contamination 4============================================================================= 5 6As shown previously (see :ref:`example_reconst_dti`), the diffusion tensor 7model is a simple way to characterize diffusion anisotropy. However, in regions 8near the cerebral ventricle and parenchyma can be underestimated by partial 9volume effects of the cerebral spinal fluid (CSF). This free water 10contamination can particularly corrupt Diffusion Tensor Imaging analysis of 11microstructural changes when different groups of subjects show different brain 12morphology (e.g. brain ventricle enlargement associated with brain tissue 13atrophy that occurs in several brain pathologies and ageing). 14 15A way to remove this free water influences is to expand the DTI model to take 16into account an extra compartment representing the contributions of free water 17diffusion [Pasternak2009]_. The expression of the expanded DTI model is shown 18below: 19 20.. math:: 21 22 S(\mathbf{g}, b) = S_0(1-f)e^{-b\mathbf{g}^T \mathbf{D} 23 \mathbf{g}}+S_0fe^{-b D_{iso}} 24 25where $\mathbf{g}$ and $b$ are diffusion gradient direction and weighted 26(more information see :ref:`example_reconst_dti`), $S(\mathbf{g}, b)$ is the 27diffusion-weighted signal measured, $S_0$ is the signal in a measurement with 28no diffusion weighting, $\mathbf{D}$ is the diffusion tensor, $f$ the volume 29fraction of the free water component, and $D_{iso}$ is the isotropic value of 30the free water diffusion (normally set to $3.0 \times 10^{-3} mm^{2}s^{-1}$). 31 32In this example, we show how to process a diffusion weighting dataset using 33an adapted version of the free water elimination proposed by [Hoy2014]_. 34 35The full details of Dipy's free water DTI implementation was published in 36[Henriques2017]_. Please cite this work if you use this algorithm. 37 38Let's start by importing the relevant modules: 39""" 40 41import numpy as np 42import dipy.reconst.fwdti as fwdti 43import dipy.reconst.dti as dti 44import matplotlib.pyplot as plt 45from dipy.data import fetch_cenir_multib 46from dipy.data import read_cenir_multib 47from dipy.segment.mask import median_otsu 48 49""" 50Without spatial constrains the free water elimination model cannot be solved 51in data acquired from one non-zero b-value [Hoy2014]_. Therefore, here we 52download a dataset that was required from multiple b-values. 53""" 54 55fetch_cenir_multib(with_raw=False) 56 57""" 58From the downloaded data, we read only the data acquired with b-values up to 592000 $s/mm^2$ to decrease the influence of non-Gaussian diffusion 60effects of the tissue which are not taken into account by the free water 61elimination model [Hoy2014]_. 62""" 63 64bvals = [200, 400, 1000, 2000] 65 66img, gtab = read_cenir_multib(bvals) 67 68data = np.asarray(img.dataobj) 69 70affine = img.affine 71 72""" 73The free water DTI model can take some minutes to process the full data set. 74Thus, we remove the background of the image to avoid unnecessary calculations. 75""" 76 77maskdata, mask = median_otsu(data, vol_idx=[0, 1], median_radius=4, numpass=2, 78 autocrop=False, dilate=1) 79 80""" 81Moreover, for illustration purposes we process only an axial slice of the 82data. 83""" 84 85axial_slice = 40 86 87mask_roi = np.zeros(data.shape[:-1], dtype=bool) 88mask_roi[:, :, axial_slice] = mask[:, :, axial_slice] 89 90""" 91The free water elimination model fit can then be initialized by instantiating 92a FreeWaterTensorModel class object: 93""" 94 95fwdtimodel = fwdti.FreeWaterTensorModel(gtab) 96 97""" 98The data can then be fitted using the ``fit`` function of the defined model 99object: 100""" 101 102fwdtifit = fwdtimodel.fit(data, mask=mask_roi) 103 104 105""" 106This 2-steps procedure will create a FreeWaterTensorFit object which contains 107all the diffusion tensor statistics free for free water contaminations. Below 108we extract the fractional anisotropy (FA) and the mean diffusivity (MD) of the 109free water diffusion tensor.""" 110 111FA = fwdtifit.fa 112MD = fwdtifit.md 113 114""" 115For comparison we also compute the same standard measures processed by the 116standard DTI model 117""" 118 119dtimodel = dti.TensorModel(gtab) 120 121dtifit = dtimodel.fit(data, mask=mask_roi) 122 123dti_FA = dtifit.fa 124dti_MD = dtifit.md 125 126""" 127Below the FA values for both free water elimination DTI model and standard DTI 128model are plotted in panels A and B, while the repective MD values are ploted 129in panels D and E. For a better visualization of the effect of the free water 130correction, the differences between these two metrics are shown in panels C and 131E. In addition to the standard diffusion statistics, the estimated volume 132fraction of the free water contamination is shown on panel G. 133""" 134 135fig1, ax = plt.subplots(2, 4, figsize=(12, 6), 136 subplot_kw={'xticks': [], 'yticks': []}) 137 138fig1.subplots_adjust(hspace=0.3, wspace=0.05) 139ax.flat[0].imshow(FA[:, :, axial_slice].T, origin='lower', 140 cmap='gray', vmin=0, vmax=1) 141ax.flat[0].set_title('A) fwDTI FA') 142ax.flat[1].imshow(dti_FA[:, :, axial_slice].T, origin='lower', 143 cmap='gray', vmin=0, vmax=1) 144ax.flat[1].set_title('B) standard DTI FA') 145 146FAdiff = abs(FA[:, :, axial_slice] - dti_FA[:, :, axial_slice]) 147ax.flat[2].imshow(FAdiff.T, cmap='gray', origin='lower', vmin=0, vmax=1) 148ax.flat[2].set_title('C) FA difference') 149 150ax.flat[3].axis('off') 151 152ax.flat[4].imshow(MD[:, :, axial_slice].T, origin='lower', 153 cmap='gray', vmin=0, vmax=2.5e-3) 154ax.flat[4].set_title('D) fwDTI MD') 155ax.flat[5].imshow(dti_MD[:, :, axial_slice].T, origin='lower', 156 cmap='gray', vmin=0, vmax=2.5e-3) 157ax.flat[5].set_title('E) standard DTI MD') 158 159MDdiff = abs(MD[:, :, axial_slice] - dti_MD[:, :, axial_slice]) 160ax.flat[6].imshow(MDdiff.T, origin='lower', cmap='gray', vmin=0, vmax=2.5e-3) 161ax.flat[6].set_title('F) MD difference') 162 163F = fwdtifit.f 164 165ax.flat[7].imshow(F[:, :, axial_slice].T, origin='lower', 166 cmap='gray', vmin=0, vmax=1) 167ax.flat[7].set_title('G) free water volume') 168 169plt.show() 170fig1.savefig('In_vivo_free_water_DTI_and_standard_DTI_measures.png') 171 172""" 173 174.. figure:: In_vivo_free_water_DTI_and_standard_DTI_measures.png 175 :align: center 176 177 In vivo diffusion measures obtain from the free water DTI and standard 178 DTI. The values of Fractional Anisotropy for the free water DTI model and 179 standard DTI model and their difference are shown in the upper panels (A-C), 180 while respective MD values are shown in the lower panels (D-F). In addition 181 the free water volume fraction estimated from the fwDTI model is shown in 182 panel G. 183 184From the figure, one can observe that the free water elimination model 185produces in general higher values of FA and lower values of MD than the 186standard DTI model. These differences in FA and MD estimation are expected 187due to the suppression of the free water isotropic diffusion components. 188Unexpected high amplitudes of FA are however observed in the periventricular 189gray matter. This is a known artefact of regions associated to voxels with high 190water volume fraction (i.e. voxels containing basically CSF). We are able to 191remove this problematic voxels by excluding all FA values associated with 192measured volume fractions above a reasonable threshold of 0.7: 193""" 194 195FA[F > 0.7] = 0 196dti_FA[F > 0.7] = 0 197 198""" 199Above we reproduce the plots of the in vivo FA from the two DTI fits and where 200we can see that the inflated FA values were practically removed: 201""" 202 203fig1, ax = plt.subplots(1, 3, figsize=(9, 3), 204 subplot_kw={'xticks': [], 'yticks': []}) 205 206fig1.subplots_adjust(hspace=0.3, wspace=0.05) 207ax.flat[0].imshow(FA[:, :, axial_slice].T, origin='lower', 208 cmap='gray', vmin=0, vmax=1) 209ax.flat[0].set_title('A) fwDTI FA') 210ax.flat[1].imshow(dti_FA[:, :, axial_slice].T, origin='lower', 211 cmap='gray', vmin=0, vmax=1) 212ax.flat[1].set_title('B) standard DTI FA') 213 214FAdiff = abs(FA[:, :, axial_slice] - dti_FA[:, :, axial_slice]) 215ax.flat[2].imshow(FAdiff.T, cmap='gray', origin='lower', vmin=0, vmax=1) 216ax.flat[2].set_title('C) FA difference') 217 218plt.show() 219fig1.savefig('In_vivo_free_water_DTI_and_standard_DTI_corrected.png') 220 221""" 222 223.. figure:: In_vivo_free_water_DTI_and_standard_DTI_corrected.png 224 :align: center 225 226 In vivo FA measures obtain from the free water DTI (A) and standard 227 DTI (B) and their difference (C). Problematic inflated FA values of the 228 images were removed by dismissing voxels above a volume fraction threshold 229 of 0.7. 230 231References 232---------- 233.. [Pasternak2009] Pasternak, O., Sochen, N., Gur, Y., Intrator, N., Assaf, Y., 234 2009. Free water elimination and mapping from diffusion MRI. Magn. Reson. 235 Med. 62(3): 717-30. doi: 10.1002/mrm.22055. 236.. [Hoy2014] Hoy, A.R., Koay, C.G., Kecskemeti, S.R., Alexander, A.L., 2014. 237 Optimization of a free water elimination two-compartmental model for 238 diffusion tensor imaging. NeuroImage 103, 323-333. doi: 239 10.1016/j.neuroimage.2014.09.053 240.. [Henriques2017] Henriques, R.N., Rokem, A., Garyfallidis, E., St-Jean, S., 241 Peterson E.T., Correia, M.M., 2017. [Re] Optimization of a free water 242 elimination two-compartment model for diffusion tensor imaging. 243 ReScience volume 3, issue 1, article number 2 244 245""" 246