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