"""
==========================================================
Signal Reconstruction Using Spherical Harmonics
==========================================================

This example shows how you can use a spherical harmonics (SH) function to
reconstruct any spherical function using DIPY_. In order to generate a
signal, we will need to generate an evenly distributed sphere.
Let's import some standard packages.
"""

import numpy as np
from dipy.core.sphere import disperse_charges, Sphere, HemiSphere

"""
We can first create some random points on a ``HemiSphere`` using spherical
polar coordinates.
"""

n_pts = 64
theta = np.pi * np.random.rand(n_pts)
phi = 2 * np.pi * np.random.rand(n_pts)
hsph_initial = HemiSphere(theta=theta, phi=phi)

"""
Next, we call ``disperse_charges`` which will iteratively move the points so
that the electrostatic potential energy is minimized. In ``hsph_updated`` we
have the updated ``HemiSphere`` with the points nicely distributed on the
hemisphere.
"""

hsph_updated, potential = disperse_charges(hsph_initial, 5000)
sphere = Sphere(xyz=np.vstack((hsph_updated.vertices, -hsph_updated.vertices)))

"""
We now need to create our initial signal. To do so, we will use our sphere's
vertices as the sampled points of our spherical function (SF). We will
use ``multi_tensor_odf`` to simulate an ODF. For more information on how to use
DIPY_ to simulate a signal and ODF, see :ref:`example_simulate_multi_tensor`.
"""

from dipy.sims.voxel import multi_tensor_odf

mevals = np.array([[0.0015, 0.00015, 0.00015],
                   [0.0015, 0.00015, 0.00015]])
angles = [(0, 0), (60, 0)]
odf = multi_tensor_odf(sphere.vertices, mevals, angles, [50, 50])


from dipy.viz import window, actor

# Enables/disables interactive visualization
interactive = False

scene = window.Scene()
scene.SetBackground(1, 1, 1)

odf_actor = actor.odf_slicer(odf[None, None, None, :], sphere=sphere)
odf_actor.RotateX(90)
scene.add(odf_actor)

print('Saving illustration as symm_signal.png')
window.record(scene, out_path='symm_signal.png', size=(300, 300))
if interactive:
    window.show(scene)

"""
.. figure:: symm_signal.png
   :align: center

   Illustration of the simulated signal sampled on a sphere of 64 points
   per hemisphere

We can now express this signal as a series of SH coefficients using
``sf_to_sh``. This function converts a series of SF coefficients in a series of
SH coefficients. For more information on SH basis, see :ref:`sh-basis`. For
this example, we will use the ``descoteaux07`` basis up to a maximum SH order
of 8.
"""

from dipy.reconst.shm import sf_to_sh

# Change this value to try out other bases
sh_basis = 'descoteaux07'
# Change this value to try other maximum orders
sh_order = 8

sh_coeffs = sf_to_sh(odf, sphere, sh_order, sh_basis)

"""
``sh_coeffs`` is an array containing the SH coefficients multiplying the SH
functions of the chosen basis. We can use it as input of ``sh_to_sf`` to
reconstruct our original signal. We will now reproject our signal on a high
resolution sphere using ``sh_to_sf``.
"""

from dipy.data import get_sphere
from dipy.reconst.shm import sh_to_sf

high_res_sph = get_sphere('symmetric724').subdivide(2)
reconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order, sh_basis)

scene.clear()
odf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph)
odf_actor.RotateX(90)
scene.add(odf_actor)

print('Saving output as symm_reconst.png')
window.record(scene, out_path='symm_reconst.png', size=(300, 300))
if interactive:
    window.show(scene)

"""
.. figure:: symm_reconst.png
   :align: center

   Reconstruction of a symmetric signal on a high resolution sphere using a
   symmetric basis

While a symmetric SH basis works well for reconstructing symmetric SF, it fails
to do so on asymmetric signals. We will now create such a signal by using a
different ODF for each hemisphere of our sphere.
"""

mevals = np.array([[0.0015, 0.0003, 0.0003]])
angles = [(0, 0)]
odf2 = multi_tensor_odf(sphere.vertices, mevals, angles, [100])

n_pts_hemisphere = int(sphere.vertices.shape[0] / 2)
asym_odf = np.append(odf[:n_pts_hemisphere], odf2[n_pts_hemisphere:])

scene.clear()
odf_actor = actor.odf_slicer(asym_odf[None, None, None, :], sphere=sphere)
odf_actor.RotateX(90)
scene.add(odf_actor)

print('Saving output as asym_signal.png')
window.record(scene, out_path='asym_signal.png', size=(300, 300))
if interactive:
    window.show(scene)

"""
.. figure:: asym_signal.png
   :align: center

   Illustration of an asymmetric signal sampled on a sphere of 64
   points per hemisphere

Let's try to reconstruct this SF using a symmetric SH basis.
"""

sh_coeffs = sf_to_sh(asym_odf, sphere, sh_order, sh_basis)
reconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order, sh_basis)

scene.clear()
odf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph)
odf_actor.RotateX(90)
scene.add(odf_actor)

print('Saving output as asym_reconst.png')
window.record(scene, out_path='asym_reconst.png', size=(300, 300))
if interactive:
    window.show(scene)

"""
.. figure:: asym_reconst.png
   :align: center

   Reconstruction of an asymmetric signal using a symmetric SH basis

As we can see, a symmetric basis fails to properly represent asymmetric SF.
Fortunately, DIPY_ also implements full SH bases, which can deal with symmetric
as well as asymmetric signals. For this tutorial, we will demonstrate it using
the ``descoteaux07`` full SH basis by setting ``full_basis=true``.
"""

sh_coeffs = sf_to_sh(asym_odf, sphere, sh_order,
                     sh_basis, full_basis=True)
reconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order,
                   sh_basis, full_basis=True)

scene.clear()
odf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph)
odf_actor.RotateX(90)
scene.add(odf_actor)

print('Saving output as asym_reconst_full.png')
window.record(scene, out_path='asym_reconst_full.png', size=(300, 300))
if interactive:
    window.show(scene)

"""
.. figure:: asym_reconst_full.png
    :align: center

    Reconstruction of an asymmetric signal using a full SH basis

As we can see, a full SH basis properly reconstruct asymmetric signal.

.. include:: ../links_names.inc

"""