1""" 2========================================================== 3Signal Reconstruction Using Spherical Harmonics 4========================================================== 5 6This example shows how you can use a spherical harmonics (SH) function to 7reconstruct any spherical function using DIPY_. In order to generate a 8signal, we will use the sphere created in :ref:`example_gradients_spheres`. 9""" 10 11import numpy as np 12from gradients_spheres import sph 13 14""" 15We now need to create our initial signal. To do so, we will use our sphere's 16vertices as the sampled points of our spherical function (SF). We will 17use ``multi_tensor_odf`` to simulate an ODF. For more information on how to use 18DIPY_ to simulate a signal and ODF, see :ref:`example_simulate_multi_tensor`. 19""" 20 21from dipy.sims.voxel import multi_tensor_odf 22 23mevals = np.array([[0.0015, 0.00015, 0.00015], 24 [0.0015, 0.00015, 0.00015]]) 25angles = [(0, 0), (60, 0)] 26odf = multi_tensor_odf(sph.vertices, mevals, angles, [50, 50]) 27 28 29from dipy.viz import window, actor 30 31# Enables/disables interactive visualization 32interactive = False 33 34scene = window.Scene() 35scene.SetBackground(1, 1, 1) 36 37odf_actor = actor.odf_slicer(odf[None, None, None, :], sphere=sph) 38odf_actor.RotateX(90) 39scene.add(odf_actor) 40 41print('Saving illustration as symm_signal.png') 42window.record(scene, out_path='symm_signal.png', size=(300, 300)) 43if interactive: 44 window.show(scene) 45 46""" 47.. figure:: symm_signal.png 48 :align: center 49 50 Illustration of the simulated signal sampled on a sphere of 64 points 51 per hemisphere 52 53We can now express this signal as a series of SH coefficients using 54``sf_to_sh``. This function converts a series of SF coefficients in a series of 55SH coefficients. For more information on SH basis, see :ref:`sh-basis`. For 56this example, we will use the ``descoteaux07`` basis up to a maximum SH order 57of 8. 58""" 59 60from dipy.reconst.shm import sf_to_sh 61 62# Change this value to try out other bases 63sh_basis = 'descoteaux07' 64# Change this value to try other maximum orders 65sh_order = 8 66 67sh_coeffs = sf_to_sh(odf, sph, sh_order, sh_basis) 68 69""" 70``sh_coeffs`` is an array containing the SH coefficients multiplying the SH 71functions of the chosen basis. We can use it as input of ``sh_to_sf`` to 72reconstruct our original signal. We will now reproject our signal on a high 73resolution sphere using ``sh_to_sf``. 74""" 75 76from dipy.data import get_sphere 77from dipy.reconst.shm import sh_to_sf 78 79high_res_sph = get_sphere('symmetric724').subdivide(2) 80reconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order, sh_basis) 81 82scene.clear() 83odf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph) 84odf_actor.RotateX(90) 85scene.add(odf_actor) 86 87print('Saving output as symm_reconst.png') 88window.record(scene, out_path='symm_reconst.png', size=(300, 300)) 89if interactive: 90 window.show(scene) 91 92""" 93.. figure:: symm_reconst.png 94 :align: center 95 96 Reconstruction of a symmetric signal on a high resolution sphere using a 97 symmetric basis 98 99While a symmetric SH basis works well for reconstructing symmetric SF, it fails 100to do so on asymmetric signals. We will now create such a signal by using a 101different ODF for each hemisphere of our sphere. 102""" 103 104mevals = np.array([[0.0015, 0.0003, 0.0003]]) 105angles = [(0, 0)] 106odf2 = multi_tensor_odf(sph.vertices, mevals, angles, [100]) 107 108n_pts_hemisphere = int(sph.vertices.shape[0] / 2) 109asym_odf = np.append(odf[:n_pts_hemisphere], odf2[n_pts_hemisphere:]) 110 111scene.clear() 112odf_actor = actor.odf_slicer(asym_odf[None, None, None, :], sphere=sph) 113odf_actor.RotateX(90) 114scene.add(odf_actor) 115 116print('Saving output as asym_signal.png') 117window.record(scene, out_path='asym_signal.png', size=(300, 300)) 118if interactive: 119 window.show(scene) 120 121""" 122.. figure:: asym_signal.png 123 :align: center 124 125 Illustration of an asymmetric signal sampled on a sphere of 64 126 points per hemisphere 127 128Let's try to reconstruct this SF using a symmetric SH basis. 129""" 130 131sh_coeffs = sf_to_sh(asym_odf, sph, sh_order, sh_basis) 132reconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order, sh_basis) 133 134scene.clear() 135odf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph) 136odf_actor.RotateX(90) 137scene.add(odf_actor) 138 139print('Saving output as asym_reconst.png') 140window.record(scene, out_path='asym_reconst.png', size=(300, 300)) 141if interactive: 142 window.show(scene) 143 144""" 145.. figure:: asym_reconst.png 146 :align: center 147 148 Reconstruction of an asymmetric signal using a symmetric SH basis 149 150As we can see, a symmetric basis fails to properly represent asymmetric SF. 151Fortunately, DIPY_ also implements full SH bases, which can deal with symmetric 152as well as asymmetric signals. For this tutorial, we will demonstrate it using 153the ``descoteaux07`` full SH basis by setting ``full_basis=true``. 154""" 155 156sh_coeffs = sf_to_sh(asym_odf, sph, sh_order, 157 sh_basis, full_basis=True) 158reconst = sh_to_sf(sh_coeffs, high_res_sph, sh_order, 159 sh_basis, full_basis=True) 160 161scene.clear() 162odf_actor = actor.odf_slicer(reconst[None, None, None, :], sphere=high_res_sph) 163odf_actor.RotateX(90) 164scene.add(odf_actor) 165 166print('Saving output as asym_reconst_full.png') 167window.record(scene, out_path='asym_reconst_full.png', size=(300, 300)) 168if interactive: 169 window.show(scene) 170 171""" 172.. figure:: asym_reconst_full.png 173 :align: center 174 175 Reconstruction of an asymmetric signal using a full SH basis 176 177As we can see, a full SH basis properly reconstruct asymmetric signal. 178 179.. include:: ../links_names.inc 180""" 181