1"""
2=====================
3Gradients and Spheres
4=====================
5
6This example shows how you can create gradient tables and sphere objects using
7DIPY_.
8
9Usually, as we saw in :ref:`example_quick_start`, you load your b-values and
10b-vectors from disk and then you can create your own gradient table. But
11this time let's say that you are an MR physicist and you want to design a new
12gradient scheme or you are a scientist who wants to simulate many different
13gradient schemes.
14
15Now let's assume that you are interested in creating a multi-shell
16acquisition with 2-shells, one at b=1000 $s/mm^2$ and one at b=2500 $s/mm^2$.
17For both shells let's say that we want a specific number of gradients (64) and
18we want to have the points on the sphere evenly distributed.
19
20This is possible using the ``disperse_charges`` which is an implementation of
21electrostatic repulsion [Jones1999]_.
22"""
23
24import numpy as np
25from dipy.core.sphere import disperse_charges, Sphere, HemiSphere
26
27"""
28We can first create some random points on a ``HemiSphere`` using spherical polar
29coordinates.
30"""
31
32n_pts = 64
33theta = np.pi * np.random.rand(n_pts)
34phi = 2 * np.pi * np.random.rand(n_pts)
35hsph_initial = HemiSphere(theta=theta, phi=phi)
36
37"""
38Next, we call ``disperse_charges`` which will iteratively move the points so that
39the electrostatic potential energy is minimized.
40"""
41
42hsph_updated, potential = disperse_charges(hsph_initial, 5000)
43
44"""
45In ``hsph_updated`` we have the updated ``HemiSphere`` with the points nicely
46distributed on the hemisphere. Let's visualize them.
47"""
48
49from dipy.viz import window, actor
50
51# Enables/disables interactive visualization
52interactive = False
53
54scene = window.Scene()
55scene.SetBackground(1, 1, 1)
56
57scene.add(actor.point(hsph_initial.vertices, window.colors.red,
58                      point_radius=0.05))
59scene.add(actor.point(hsph_updated.vertices, window.colors.green,
60                      point_radius=0.05))
61
62print('Saving illustration as initial_vs_updated.png')
63window.record(scene, out_path='initial_vs_updated.png', size=(300, 300))
64if interactive:
65    window.show(scene)
66
67"""
68.. figure:: initial_vs_updated.png
69   :align: center
70
71   Illustration of electrostatic repulsion of red points which become
72   green points.
73
74We can also create a sphere from the hemisphere and show it in the
75following way.
76
77"""
78
79sph = Sphere(xyz=np.vstack((hsph_updated.vertices, -hsph_updated.vertices)))
80
81scene.clear()
82scene.add(actor.point(sph.vertices, window.colors.green, point_radius=0.05))
83
84print('Saving illustration as full_sphere.png')
85window.record(scene, out_path='full_sphere.png', size=(300, 300))
86if interactive:
87    window.show(scene)
88
89"""
90.. figure:: full_sphere.png
91   :align: center
92
93   Full sphere.
94
95It is time to create the Gradients. For this purpose we will use the
96function ``gradient_table`` and fill it with the ``hsph_updated`` vectors that
97we created above.
98"""
99
100from dipy.core.gradients import gradient_table
101
102vertices = hsph_updated.vertices
103values = np.ones(vertices.shape[0])
104
105"""
106We need two stacks of ``vertices``, one for every shell, and we need two sets
107of b-values, one at 1000 $s/mm^2$, and one at 2500 $s/mm^2$, as we discussed
108previously.
109"""
110
111bvecs = np.vstack((vertices, vertices))
112bvals = np.hstack((1000 * values, 2500 * values))
113
114"""
115We can also add some b0s. Let's add one at the beginning and one at the end.
116"""
117
118bvecs = np.insert(bvecs, (0, bvecs.shape[0]), np.array([0, 0, 0]), axis=0)
119bvals = np.insert(bvals, (0, bvals.shape[0]), 0)
120
121print(bvals)
122
123"""
124
125::
126
127    [    0.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
128      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
129      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
130      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
131      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
132      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
133      1000.  1000.  1000.  1000.  1000.  2500.  2500.  2500.  2500.  2500.
134      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.
135      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.
136      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.
137      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.
138      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.
139      2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.  2500.     0.]
140
141"""
142
143print(bvecs)
144
145"""
146
147::
148
149    [[ 0.          0.          0.        ]
150     [-0.80451777 -0.16877559  0.56944355]
151     [ 0.32822557 -0.94355999  0.04430036]
152     [-0.23584135 -0.96241331  0.13468285]
153     [-0.39207424 -0.73505312  0.55314981]
154     [-0.32539386 -0.16751384  0.93062235]
155     [-0.82043195 -0.39411534  0.41420347]
156     [ 0.65741493  0.74947875  0.07802061]
157     [ 0.88853765  0.45303621  0.07251925]
158     [ 0.39638642 -0.15185138  0.90543855]
159                     ...
160     [ 0.10175269  0.08197111  0.99142681]
161     [ 0.50577702 -0.37862345  0.77513476]
162     [ 0.42845026  0.40155296  0.80943535]
163     [ 0.26939707  0.81103868  0.51927014]
164     [-0.48938584 -0.43780086  0.75420946]
165     [ 0.          0.          0.        ]]
166
167Both b-values and b-vectors look correct. Let's now create the
168``GradientTable``.
169"""
170
171gtab = gradient_table(bvals, bvecs)
172
173scene.clear()
174
175"""
176We can also visualize the gradients. Let's color the first shell blue and
177the second shell cyan.
178"""
179
180colors_b1000 = window.colors.blue * np.ones(vertices.shape)
181colors_b2500 = window.colors.cyan * np.ones(vertices.shape)
182colors = np.vstack((colors_b1000, colors_b2500))
183colors = np.insert(colors, (0, colors.shape[0]), np.array([0, 0, 0]), axis=0)
184colors = np.ascontiguousarray(colors)
185
186scene.add(actor.point(gtab.gradients, colors, point_radius=100))
187
188print('Saving illustration as gradients.png')
189window.record(scene, out_path='gradients.png', size=(300, 300))
190if interactive:
191    window.show(scene)
192
193"""
194.. figure:: gradients.png
195   :align: center
196
197   Diffusion gradients.
198
199References
200----------
201
202.. [Jones1999] Jones, DK. et al. Optimal strategies for measuring diffusion in
203   anisotropic systems by magnetic resonance imaging, Magnetic Resonance in
204   Medicine, vol 42, no 3, 515-525, 1999.
205
206.. include:: ../links_names.inc
207
208"""
209