{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# DKI MultiTensor Simulation\n\nIn this example we show how to simulate the Diffusion Kurtosis Imaging (DKI)\ndata of a single voxel. DKI captures information about the non-Gaussian\nproperties of water diffusion which is a consequence of the existence of tissue\nbarriers and compartments. In these simulations compartmental heterogeneity is\ntaken into account by modeling different compartments for the intra- and\nextra-cellular media of two populations of fibers. These simulations are\nperformed according to [RNH2015]_.\n\nWe first import all relevant modules.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom dipy.sims.voxel import (multi_tensor_dki, single_tensor)\nfrom dipy.data import get_fnames\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.core.gradients import gradient_table\nfrom dipy.reconst.dti import (decompose_tensor, from_lower_triangular)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the simulation we will need a GradientTable with the b-values and\nb-vectors. Here we use the GradientTable of the sample DIPY_ dataset\n``small_64D``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fimg, fbvals, fbvecs = get_fnames('small_64D')\nbvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "DKI requires data from more than one non-zero b-value. Since the dataset\n``small_64D`` was acquired with one non-zero b-value we artificially produce a\nsecond non-zero b-value.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bvals = np.concatenate((bvals, bvals * 2), axis=0)\nbvecs = np.concatenate((bvecs, bvecs), axis=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The b-values and gradient directions are then converted to DIPY's\n``GradientTable`` format.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gtab = gradient_table(bvals, bvecs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In ``mevals`` we save the eigenvalues of each tensor. To simulate crossing\nfibers with two different media (representing intra and extra-cellular media),\na total of four components have to be taken in to account (i.e. the first two\ncompartments correspond to the intra and extra cellular media for the first\nfiber population while the others correspond to the media of the second fiber\npopulation)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mevals = np.array([[0.00099, 0, 0],\n                   [0.00226, 0.00087, 0.00087],\n                   [0.00099, 0, 0],\n                   [0.00226, 0.00087, 0.00087]])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In ``angles`` we save in polar coordinates ($\\theta, \\phi$) the principal\naxis of each compartment tensor. To simulate crossing fibers at 70$^{\\circ}$\nthe compartments of the first fiber are aligned to the X-axis while the\ncompartments of the second fiber are aligned to the X-Z plane with an angular\ndeviation of 70$^{\\circ}$ from the first one.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "angles = [(90, 0), (90, 0), (20, 0), (20, 0)]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In ``fractions`` we save the percentage of the contribution of each\ncompartment, which is computed by multiplying the percentage of contribution\nof each fiber population and the water fraction of each different medium\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fie = 0.49  # intra-axonal water fraction\nfractions = [fie*50, (1 - fie)*50, fie*50, (1 - fie)*50]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Having defined the parameters for all tissue compartments, the elements of the\ndiffusion tensor (DT), the elements of the kurtosis tensor (KT) and the DW\nsignals simulated from the DKI model can be obtain using the function\n``multi_tensor_dki``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "signal_dki, dt, kt = multi_tensor_dki(gtab, mevals, S0=200, angles=angles,\n                                      fractions=fractions, snr=None)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also add Rician noise with a specific SNR.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "signal_noisy, dt, kt = multi_tensor_dki(gtab, mevals, S0=200,\n                                        angles=angles, fractions=fractions,\n                                        snr=10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For comparison purposes, we also compute the DW signal if only the diffusion\ntensor components are taken into account. For this we use DIPY's function\n``single_tensor`` which requires that dt is decomposed into its eigenvalues and\neigenvectors.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt_evals, dt_evecs = decompose_tensor(from_lower_triangular(dt))\nsignal_dti = single_tensor(gtab, S0=200, evals=dt_evals, evecs=dt_evecs,\n                           snr=None)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we can visualize the values of the different version of simulated\nsignals for all assumed gradient directions and bvalues.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(signal_dti, label='noiseless dti')\nplt.plot(signal_dki, label='noiseless dki')\nplt.plot(signal_noisy, label='with noise')\nplt.legend()\nplt.show()\nplt.savefig('simulated_dki_signal.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: simulated_dki_signal.png\n   :align: center\n\n   Simulated signals obtain from the DTI and DKI models.\n\nNon-Gaussian diffusion properties in tissues are responsible to smaller signal\nattenuation for larger bvalues when compared to signal attenuation from free\nGaussian water diffusion. This feature can be shown from the figure above,\nsince signals simulated from the DKI models reveals larger DW signal\nintensities than the signals obtained only from the diffusion tensor\ncomponents.\n\n## References\n\n.. [RNH2015] R. Neto Henriques et al., \"Exploring the 3D geometry of the\n   diffusion kurtosis tensor - Impact on the development of robust tractography\n   procedures and novel biomarkers\", NeuroImage (2015) 111, 85-99.\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
}