{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# MultiTensor Simulation\n\nIn this example we show how someone can simulate the signal and the ODF of a\nsingle voxel using a MultiTensor.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nfrom dipy.core.sphere import disperse_charges, HemiSphere\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_sphere\nfrom dipy.sims.voxel import multi_tensor, multi_tensor_odf\nfrom dipy.viz import window, actor"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the simulation we will need a GradientTable with the b-values and\nb-vectors. To create one, we can first create some random points on a\n``HemiSphere`` using spherical polar coordinates.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_pts = 64\ntheta = np.pi * np.random.rand(n_pts)\nphi = 2 * np.pi * np.random.rand(n_pts)\nhsph_initial = HemiSphere(theta=theta, phi=phi)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we call ``disperse_charges`` which will iteratively move the points so\nthat the electrostatic potential energy is minimized.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "hsph_updated, potential = disperse_charges(hsph_initial, 5000)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We need two stacks of ``vertices``, one for every shell, and we need two sets\nof b-values, one at 1000 $s/mm^2$, and one at 2500 $s/mm^2$, as we discussed\npreviously.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vertices = hsph_updated.vertices\nvalues = np.ones(vertices.shape[0])\n\nbvecs = np.vstack((vertices, vertices))\nbvals = np.hstack((1000 * values, 2500 * values))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also add some b0s. Let's add one at the beginning and one at the end.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bvecs = np.insert(bvecs, (0, bvecs.shape[0]), np.array([0, 0, 0]), axis=0)\nbvals = np.insert(bvals, (0, bvals.shape[0]), 0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's now create the ``GradientTable``.\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.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mevals = np.array([[0.0015, 0.0003, 0.0003],\n                   [0.0015, 0.0003, 0.0003]])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In ``angles`` we save in polar coordinates ($\\theta, \\phi$) the principal\naxis of each tensor.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "angles = [(0, 0), (60, 0)]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In ``fractions`` we save the percentage of the contribution of each tensor.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fractions = [50, 50]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The function ``multi_tensor`` will return the simulated signal and an array\nwith the principal axes of the tensors in cartesian coordinates.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "signal, sticks = multi_tensor(gtab, mevals, S0=100, 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, sticks = multi_tensor(gtab, mevals, S0=100, angles=angles,\n                                    fractions=fractions, snr=20)\n\nplt.plot(signal, label='noiseless')\n\nplt.plot(signal_noisy, label='with noise')\nplt.legend()\n# plt.show()\nplt.savefig('simulated_signal.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: simulated_signal.png\n   :align: center\n\n   **Simulated MultiTensor signal**\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the ODF simulation we will need a sphere. Because we are interested in a\nsimulation of only a single voxel, we can use a sphere with very high\nresolution. We generate that by subdividing the triangles of one of DIPY_'s\ncached spheres, which we can read in the following way.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sphere = get_sphere('repulsion724')\nsphere = sphere.subdivide(2)\n\nodf = multi_tensor_odf(sphere.vertices, mevals, angles, fractions)\n\n# Enables/disables interactive visualization\ninteractive = False\n\nscene = window.Scene()\n\nodf_actor = actor.odf_slicer(odf[None, None, None, :], sphere=sphere,\n                             colormap='plasma')\nodf_actor.RotateX(90)\n\nscene.add(odf_actor)\n\nprint('Saving illustration as multi_tensor_simulation')\nwindow.record(scene, out_path='multi_tensor_simulation.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: multi_tensor_simulation.png\n   :align: center\n\n   Simulating a MultiTensor ODF.\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
}