{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Gradients and Spheres\n\nThis example shows how you can create gradient tables and sphere objects using\nDIPY_.\n\nUsually, as we saw in `example_quick_start`, you load your b-values and\nb-vectors from disk and then you can create your own gradient table. But\nthis time let's say that you are an MR physicist and you want to design a new\ngradient scheme or you are a scientist who wants to simulate many different\ngradient schemes.\n\nNow let's assume that you are interested in creating a multi-shell\nacquisition with 2-shells, one at b=1000 $s/mm^2$ and one at b=2500 $s/mm^2$.\nFor both shells let's say that we want a specific number of gradients (64) and\nwe want to have the points on the sphere evenly distributed.\n\nThis is possible using the ``disperse_charges`` which is an implementation of\nelectrostatic repulsion [Jones1999]_.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom dipy.core.sphere import disperse_charges, Sphere, HemiSphere"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can first create some random points on a ``HemiSphere`` using spherical polar\ncoordinates.\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 that\nthe 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": [
        "In ``hsph_updated`` we have the updated ``HemiSphere`` with the points nicely\ndistributed on the hemisphere. Let's visualize them.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.viz import window, actor\n\n# Enables/disables interactive visualization\ninteractive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\n\nscene.add(actor.point(hsph_initial.vertices, window.colors.red,\n                      point_radius=0.05))\nscene.add(actor.point(hsph_updated.vertices, window.colors.green,\n                      point_radius=0.05))\n\nprint('Saving illustration as initial_vs_updated.png')\nwindow.record(scene, out_path='initial_vs_updated.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: initial_vs_updated.png\n   :align: center\n\n   Illustration of electrostatic repulsion of red points which become\n   green points.\n\nWe can also create a sphere from the hemisphere and show it in the\nfollowing way.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sph = Sphere(xyz=np.vstack((hsph_updated.vertices, -hsph_updated.vertices)))\n\nscene.clear()\nscene.add(actor.point(sph.vertices, window.colors.green, point_radius=0.05))\n\nprint('Saving illustration as full_sphere.png')\nwindow.record(scene, out_path='full_sphere.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: full_sphere.png\n   :align: center\n\n   Full sphere.\n\nIt is time to create the Gradients. For this purpose we will use the\nfunction ``gradient_table`` and fill it with the ``hsph_updated`` vectors that\nwe created above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.core.gradients import gradient_table\n\nvertices = hsph_updated.vertices\nvalues = np.ones(vertices.shape[0])"
      ]
    },
    {
      "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": [
        "bvecs = 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)\n\nprint(bvals)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "::\n\n    [    0.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.\n      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.\n      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.\n      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.\n      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.\n      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.\n      1000.  1000.  1000.  1000.  1000.  2500.  2500.  2500.  2500.  2500.\n      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.\n      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.\n      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.\n      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.\n      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.\n      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.     0.]\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(bvecs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "::\n\n    [[ 0.          0.          0.        ]\n     [-0.80451777 -0.16877559  0.56944355]\n     [ 0.32822557 -0.94355999  0.04430036]\n     [-0.23584135 -0.96241331  0.13468285]\n     [-0.39207424 -0.73505312  0.55314981]\n     [-0.32539386 -0.16751384  0.93062235]\n     [-0.82043195 -0.39411534  0.41420347]\n     [ 0.65741493  0.74947875  0.07802061]\n     [ 0.88853765  0.45303621  0.07251925]\n     [ 0.39638642 -0.15185138  0.90543855]\n                     ...\n     [ 0.10175269  0.08197111  0.99142681]\n     [ 0.50577702 -0.37862345  0.77513476]\n     [ 0.42845026  0.40155296  0.80943535]\n     [ 0.26939707  0.81103868  0.51927014]\n     [-0.48938584 -0.43780086  0.75420946]\n     [ 0.          0.          0.        ]]\n\nBoth b-values and b-vectors look correct. Let's now create the\n``GradientTable``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gtab = gradient_table(bvals, bvecs)\n\nscene.clear()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also visualize the gradients. Let's color the first shell blue and\nthe second shell cyan.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "colors_b1000 = window.colors.blue * np.ones(vertices.shape)\ncolors_b2500 = window.colors.cyan * np.ones(vertices.shape)\ncolors = np.vstack((colors_b1000, colors_b2500))\ncolors = np.insert(colors, (0, colors.shape[0]), np.array([0, 0, 0]), axis=0)\ncolors = np.ascontiguousarray(colors)\n\nscene.add(actor.point(gtab.gradients, colors, point_radius=100))\n\nprint('Saving illustration as gradients.png')\nwindow.record(scene, out_path='gradients.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: gradients.png\n   :align: center\n\n   Diffusion gradients.\n\n## References\n\n.. [Jones1999] Jones, DK. et al. Optimal strategies for measuring diffusion in\n   anisotropic systems by magnetic resonance imaging, Magnetic Resonance in\n   Medicine, vol 42, no 3, 515-525, 1999.\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
}