{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Using the free water elimination model to remove DTI free water contamination\n\nAs shown previously (see `example_reconst_dti`), the diffusion tensor model\nis a simple way to characterize diffusion anisotropy. However, in regions near\nthe ventricles and parenchyma, anisotropy can be underestimated by partial\nvolume effects of the cerebral spinal fluid (CSF). This free water contamination\ncan particularly corrupt Diffusion Tensor Imaging analysis of microstructural\nchanges when different groups of subjects show different brain morphology (e.g.\nbrain ventricle enlargement associated with brain tissue atrophy that occurs in\nseveral brain pathologies and aging).\n\nA way to remove this free water influences is to expand the DTI model to take\ninto account an extra compartment representing the contributions of free water\ndiffusion [Pasternak2009]_. The expression of the expanded DTI model is shown\nbelow:\n\n\\begin{align}S(\\mathbf{g}, b) = S_0(1-f)e^{-b\\mathbf{g}^T \\mathbf{D}\n    \\mathbf{g}}+S_0fe^{-b D_{iso}}\\end{align}\n\nwhere $\\mathbf{g}$ and $b$ are diffusion gradient direction and weighted (more\ninformation see `example_reconst_dti`), $S(\\mathbf{g}, b)$ is the\ndiffusion-weighted signal measured, $S_0$ is the signal in a measurement with no\ndiffusion weighting, $\\mathbf{D}$ is the diffusion tensor, $f$ the volume\nfraction of the free water component, and $D_{iso}$ is the isotropic value of\nthe free water diffusion (normally set to $3.0  imes 10^{-3} mm^{2}s^{-1}$).\n\nIn this example, we show how to process a diffusion weighting dataset using an\nadapted version of the free water elimination proposed by [Hoy2014]_.\n\nThe full details of Dipy's free water DTI implementation was published in\n[Henriques2017]_. Please cite this work if you use this algorithm.\n\nLet's start by importing the relevant modules:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport dipy.reconst.fwdti as fwdti\nimport dipy.reconst.dti as dti\nimport matplotlib.pyplot as plt\nfrom dipy.data import fetch_hbn\nimport os.path as op\nimport nibabel as nib\nfrom dipy.core.gradients import gradient_table"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Without spatial constrains the free water elimination model cannot be solved\nin data acquired from one non-zero b-value [Hoy2014]_. Therefore, here we\ndownload a dataset that was acquired with multiple b-values.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = fetch_hbn([\"NDARAA948VFH\"])[1]\ndwi_path = op.join(\n       data_path, \"derivatives\", \"qsiprep\", \"sub-NDARAA948VFH\",\n       \"ses-HBNsiteRU\", \"dwi\")\n\nimg = nib.load(op.join(\n       dwi_path,\n       \"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.nii.gz\"))\n\ngtab = gradient_table(\n       op.join(dwi_path,\n\"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.bval\"),\n       op.join(dwi_path,\n\"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.bvec\"))\n\ndata = np.asarray(img.dataobj)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The free water DTI model can take some minutes to process the full data set.\nThus, we use a brain mask that was calculated during pre-processing, to remove\nthe background of the image and avoid unnecessary calculations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mask_img = nib.load(\n       op.join(dwi_path,\n\"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-brain_mask.nii.gz\"))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Moreover, for illustration purposes we process only one slice of the data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mask = mask_img.get_fdata()\n\ndata_small = data[:, :, 50:51]\nmask_small = mask[:, :, 50:51]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The free water elimination model fit can then be initialized by instantiating\na FreeWaterTensorModel class object:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fwdtimodel = fwdti.FreeWaterTensorModel(gtab)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The data can then be fitted using the ``fit`` function of the defined model\nobject:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fwdtifit = fwdtimodel.fit(data_small, mask=mask_small)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This 2-steps procedure will create a FreeWaterTensorFit object which contains\nall the diffusion tensor statistics free for free water contamination. Below\nwe extract the fractional anisotropy (FA) and the mean diffusivity (MD) of the\nfree water diffusion tensor.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "FA = fwdtifit.fa\nMD = fwdtifit.md"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For comparison we also compute the same standard measures processed by the\nstandard DTI model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dtimodel = dti.TensorModel(gtab)\n\ndtifit = dtimodel.fit(data_small, mask=mask_small)\n\ndti_FA = dtifit.fa\ndti_MD = dtifit.md"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Below the FA values for both free water elimination DTI model and standard DTI\nmodel are plotted in panels A and B, while the respective MD values are plotted\nin panels D and E. For a better visualization of the effect of the free water\ncorrection, the differences between these two metrics are shown in panels C and\nE. In addition to the standard diffusion statistics, the estimated volume\nfraction of the free water contamination is shown on panel G.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig1, ax = plt.subplots(2, 4, figsize=(12, 6),\n                        subplot_kw={'xticks': [], 'yticks': []})\n\nfig1.subplots_adjust(hspace=0.3, wspace=0.05)\nax.flat[0].imshow(FA[:, :, 0].T, origin='lower',\n                  cmap='gray', vmin=0, vmax=1)\nax.flat[0].set_title('A) fwDTI FA')\nax.flat[1].imshow(dti_FA[:, :, 0].T, origin='lower',\n                  cmap='gray', vmin=0, vmax=1)\nax.flat[1].set_title('B) standard DTI FA')\n\nFAdiff = abs(FA[:, :, 0] - dti_FA[:, :, 0])\nax.flat[2].imshow(FAdiff.T, cmap='gray', origin='lower', vmin=0, vmax=1)\nax.flat[2].set_title('C) FA difference')\n\nax.flat[3].axis('off')\n\nax.flat[4].imshow(MD[:, :, 0].T, origin='lower',\n                  cmap='gray', vmin=0, vmax=2.5e-3)\nax.flat[4].set_title('D) fwDTI MD')\nax.flat[5].imshow(dti_MD[:, :, 0].T, origin='lower',\n                  cmap='gray', vmin=0, vmax=2.5e-3)\nax.flat[5].set_title('E) standard DTI MD')\n\nMDdiff = abs(MD[:, :, 0] - dti_MD[:, :, 0])\nax.flat[6].imshow(MDdiff.T, origin='lower', cmap='gray', vmin=0, vmax=2.5e-3)\nax.flat[6].set_title('F) MD difference')\n\nF = fwdtifit.f\n\nax.flat[7].imshow(F[:, :, 0].T, origin='lower',\n                  cmap='gray', vmin=0, vmax=1)\nax.flat[7].set_title('G) free water volume')\n\nplt.show()\nfig1.savefig('In_vivo_free_water_DTI_and_standard_DTI_measures.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: In_vivo_free_water_DTI_and_standard_DTI_measures.png\n   :align: center\n\n   In vivo diffusion measures obtain from the free water DTI and standard\n   DTI. The values of Fractional Anisotropy for the free water DTI model and\n   standard DTI model and their difference are shown in the upper panels (A-C),\n   while respective MD values are shown in the lower panels (D-F). In addition\n   the free water volume fraction estimated from the fwDTI model is shown in\n   panel G.\n\nFrom the figure, one can observe that the free water elimination model\nproduces in general higher values of FA and lower values of MD than the\nstandard DTI model. These differences in FA and MD estimation are expected\ndue to the suppression of the free water isotropic diffusion components.\nUnexpected high amplitudes of FA are however observed in the periventricular\ngray matter. This is a known artefact of regions associated to voxels with high\nwater volume fraction (i.e. voxels containing basically CSF). We are able to\nremove this problematic voxels by excluding all FA values associated with\nmeasured volume fractions above a reasonable threshold of 0.7:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "FA[F > 0.7] = 0\ndti_FA[F > 0.7] = 0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Above we reproduce the plots of the in vivo FA from the two DTI fits and where\nwe can see that the inflated FA values were practically removed:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig1, ax = plt.subplots(1, 3, figsize=(9, 3),\n                        subplot_kw={'xticks': [], 'yticks': []})\n\nfig1.subplots_adjust(hspace=0.3, wspace=0.05)\nax.flat[0].imshow(FA[:, :, 0].T, origin='lower',\n                  cmap='gray', vmin=0, vmax=1)\nax.flat[0].set_title('A) fwDTI FA')\nax.flat[1].imshow(dti_FA[:, :, 0].T, origin='lower',\n                  cmap='gray', vmin=0, vmax=1)\nax.flat[1].set_title('B) standard DTI FA')\n\nFAdiff = abs(FA[:, :, 0] - dti_FA[:, :, 0])\nax.flat[2].imshow(FAdiff.T, cmap='gray', origin='lower', vmin=0, vmax=1)\nax.flat[2].set_title('C) FA difference')\n\nplt.show()\nfig1.savefig('In_vivo_free_water_DTI_and_standard_DTI_corrected.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: In_vivo_free_water_DTI_and_standard_DTI_corrected.png\n   :align: center\n\n   In vivo FA measures obtain from the free water DTI (A) and standard\n   DTI (B) and their difference (C). Problematic inflated FA values of the\n   images were removed by dismissing voxels above a volume fraction threshold\n   of 0.7.\n\n## References\n.. [Pasternak2009] Pasternak, O., Sochen, N., Gur, Y., Intrator, N., Assaf, Y.,\n   2009. Free water elimination and mapping from diffusion MRI. Magn. Reson.\n   Med. 62(3): 717-30. doi: 10.1002/mrm.22055.\n\n.. [Hoy2014] Hoy, A.R., Koay, C.G., Kecskemeti, S.R., Alexander, A.L., 2014.\n   Optimization of a free water elimination two-compartmental model for\n   diffusion tensor imaging. NeuroImage 103, 323-333. doi:\n   10.1016/j.neuroimage.2014.09.053\n\n.. [Henriques2017] Henriques, R.N., Rokem, A., Garyfallidis, E., St-Jean, S.,\n   Peterson E.T., Correia, M.M., 2017. [Re] Optimization of a free water\n   elimination two-compartment model for diffusion tensor imaging.\n   ReScience volume 3, issue 1, article number 2\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
}