{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Signal Reconstruction Using Spherical Harmonics\n\nThis example shows how you can use a spherical harmonics (SH) function to\nreconstruct any spherical function using DIPY_. In order to generate a\nsignal, we will need to generate an evenly distributed sphere.\nLet's import some standard packages.\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\npolar 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. In ``hsph_updated`` we\nhave the updated ``HemiSphere`` with the points nicely distributed on the\nhemisphere.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "hsph_updated, potential = disperse_charges(hsph_initial, 5000)\nsphere = Sphere(xyz=np.vstack((hsph_updated.vertices, -hsph_updated.vertices)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now need to create our initial signal. To do so, we will use our sphere's\nvertices as the sampled points of our spherical function (SF). We will\nuse ``multi_tensor_odf`` to simulate an ODF. For more information on how to use\nDIPY_ to simulate a signal and ODF, see `example_simulate_multi_tensor`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.sims.voxel import multi_tensor_odf\n\nmevals = np.array([[0.0015, 0.00015, 0.00015],\n                   [0.0015, 0.00015, 0.00015]])\nangles = [(0, 0), (60, 0)]\nodf = multi_tensor_odf(sphere.vertices, mevals, angles, [50, 50])\n\n\nfrom dipy.viz import window, actor\n\n# Enables/disables interactive visualization\ninteractive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\n\nodf_actor = actor.odf_slicer(odf[None, None, None, :], sphere=sphere)\nodf_actor.RotateX(90)\nscene.add(odf_actor)\n\nprint('Saving illustration as symm_signal.png')\nwindow.record(scene, out_path='symm_signal.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: symm_signal.png\n   :align: center\n\n   Illustration of the simulated signal sampled on a sphere of 64 points\n   per hemisphere\n\nWe can now express this signal as a series of SH coefficients using\n``sf_to_sh``. This function converts a series of SF coefficients in a series of\nSH coefficients. For more information on SH basis, see `sh-basis`. For\nthis example, we will use the ``descoteaux07`` basis up to a maximum SH order\nof 8.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.reconst.shm import sf_to_sh\n\n# Change this value to try out other bases\nsh_basis = 'descoteaux07'\n# Change this value to try other maximum orders\nsh_order = 8\n\nsh_coeffs = sf_to_sh(odf, sphere, sh_order, sh_basis)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "``sh_coeffs`` is an array containing the SH coefficients multiplying the SH\nfunctions of the chosen basis. We can use it as input of ``sh_to_sf`` to\nreconstruct our original signal. We will now reproject our signal on a high\nresolution sphere using ``sh_to_sf``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.data import get_sphere\nfrom dipy.reconst.shm import sh_to_sf\n\nhigh_res_sph = get_sphere('symmetric724').subdivide(2)\nreconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order, sh_basis)\n\nscene.clear()\nodf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph)\nodf_actor.RotateX(90)\nscene.add(odf_actor)\n\nprint('Saving output as symm_reconst.png')\nwindow.record(scene, out_path='symm_reconst.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: symm_reconst.png\n   :align: center\n\n   Reconstruction of a symmetric signal on a high resolution sphere using a\n   symmetric basis\n\nWhile a symmetric SH basis works well for reconstructing symmetric SF, it fails\nto do so on asymmetric signals. We will now create such a signal by using a\ndifferent ODF for each hemisphere of our sphere.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mevals = np.array([[0.0015, 0.0003, 0.0003]])\nangles = [(0, 0)]\nodf2 = multi_tensor_odf(sphere.vertices, mevals, angles, [100])\n\nn_pts_hemisphere = int(sphere.vertices.shape[0] / 2)\nasym_odf = np.append(odf[:n_pts_hemisphere], odf2[n_pts_hemisphere:])\n\nscene.clear()\nodf_actor = actor.odf_slicer(asym_odf[None, None, None, :], sphere=sphere)\nodf_actor.RotateX(90)\nscene.add(odf_actor)\n\nprint('Saving output as asym_signal.png')\nwindow.record(scene, out_path='asym_signal.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: asym_signal.png\n   :align: center\n\n   Illustration of an asymmetric signal sampled on a sphere of 64\n   points per hemisphere\n\nLet's try to reconstruct this SF using a symmetric SH basis.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sh_coeffs = sf_to_sh(asym_odf, sphere, sh_order, sh_basis)\nreconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order, sh_basis)\n\nscene.clear()\nodf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph)\nodf_actor.RotateX(90)\nscene.add(odf_actor)\n\nprint('Saving output as asym_reconst.png')\nwindow.record(scene, out_path='asym_reconst.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: asym_reconst.png\n   :align: center\n\n   Reconstruction of an asymmetric signal using a symmetric SH basis\n\nAs we can see, a symmetric basis fails to properly represent asymmetric SF.\nFortunately, DIPY_ also implements full SH bases, which can deal with symmetric\nas well as asymmetric signals. For this tutorial, we will demonstrate it using\nthe ``descoteaux07`` full SH basis by setting ``full_basis=true``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sh_coeffs = sf_to_sh(asym_odf, sphere, sh_order,\n                     sh_basis, full_basis=True)\nreconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order,\n                   sh_basis, full_basis=True)\n\nscene.clear()\nodf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph)\nodf_actor.RotateX(90)\nscene.add(odf_actor)\n\nprint('Saving output as asym_reconst_full.png')\nwindow.record(scene, out_path='asym_reconst_full.png', size=(300, 300))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: asym_reconst_full.png\n    :align: center\n\n    Reconstruction of an asymmetric signal using a full SH basis\n\nAs we can see, a full SH basis properly reconstruct asymmetric signal.\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
}