Gradients and Spheres

This example shows how you can create gradient tables and sphere objects using DIPY.

Usually, as we saw in Getting started with DIPY, you load your b-values and b-vectors from disk and then you can create your own gradient table. But this time let’s say that you are an MR physicist and you want to design a new gradient scheme or you are a scientist who wants to simulate many different gradient schemes.

Now let’s assume that you are interested in creating a multi-shell acquisition with 2-shells, one at b=1000 \(s/mm^2\) and one at b=2500 \(s/mm^2\). For both shells let’s say that we want a specific number of gradients (64) and we want to have the points on the sphere evenly distributed.

This is possible using the disperse_charges which is an implementation of electrostatic repulsion Jones et al.[1] .

Let’s start by importing the necessary modules.

import numpy as np

from dipy.core.gradients import gradient_table
from dipy.core.sphere import HemiSphere, Sphere, disperse_charges
from dipy.viz import actor, window

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

rng = np.random.default_rng()
n_pts = 64
theta = np.pi * rng.random(n_pts)
phi = 2 * np.pi * rng.random(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.

hsph_updated, potential = disperse_charges(hsph_initial, 5000)

In hsph_updated we have the updated HemiSphere with the points nicely distributed on the hemisphere. Let’s visualize them.

# Enables/disables interactive visualization
interactive = False

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

scene.add(actor.point(hsph_initial.vertices, window.colors.red, point_radius=0.05))
scene.add(actor.point(hsph_updated.vertices, window.colors.green, point_radius=0.05))

window.record(scene=scene, out_path="initial_vs_updated.png", size=(300, 300))
if interactive:
    window.show(scene)

Illustration of electrostatic repulsion of red points which become green points.

We can also create a sphere from the hemisphere and show it in the following way.

sph = Sphere(xyz=np.vstack((hsph_updated.vertices, -hsph_updated.vertices)))

scene.clear()
scene.add(actor.point(sph.vertices, window.colors.green, point_radius=0.05))

window.record(scene=scene, out_path="full_sphere.png", size=(300, 300))
if interactive:
    window.show(scene)

Full sphere.

It is time to create the Gradients. For this purpose we will use the function gradient_table and fill it with the hsph_updated vectors that we created above.

vertices = hsph_updated.vertices
values = np.ones(vertices.shape[0])

We need two stacks of vertices, one for every shell, and we need two sets of b-values, one at 1000 \(s/mm^2\), and one at 2500 \(s/mm^2\), as we discussed previously.

bvecs = np.vstack((vertices, vertices))
bvals = np.hstack((1000 * values, 2500 * values))

We can also add some b0s. Let’s add one at the beginning and one at the end.

bvecs = np.insert(bvecs, (0, bvecs.shape[0]), np.array([0, 0, 0]), axis=0)
bvals = np.insert(bvals, (0, bvals.shape[0]), 0)

print(bvals)
print(bvecs)

[   0. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000.
 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000.
 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000.
 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000.
 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000.
 1000. 1000. 1000. 1000. 1000. 2500. 2500. 2500. 2500. 2500. 2500. 2500.
 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500.
 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500.
 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500.
 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500.
 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500. 2500.    0.]
[[ 0.          0.          0.        ]
 [ 0.37618308  0.9219296   0.09236939]
 [ 0.84851297 -0.03638374  0.52792231]
 [ 0.74825058  0.22093446  0.62554699]
 [-0.53927213  0.61071147  0.57984228]
 [ 0.98013998 -0.04149057  0.19391791]
 [ 0.02461784  0.17667659  0.98396105]
 [ 0.38988413  0.64317566  0.65902612]
 [ 0.034099   -0.56509772  0.82431901]
 [-0.53657134 -0.42997535  0.72609393]
 [ 0.64074659  0.75712022  0.12732942]
 [ 0.33356492 -0.59797103  0.72881074]
 [-0.87196661 -0.438395    0.21790836]
 [-0.57971596 -0.66583658  0.46967121]
 [ 0.05332005 -0.92972909  0.36436354]
 [ 0.9362261  -0.34243955  0.07884059]
 [-0.06798288  0.88064643  0.46887119]
 [-0.22158404  0.96259869  0.15589831]
 [-0.11817505 -0.98566508  0.12041184]
 [ 0.07775453  0.97909906  0.18793423]
 [ 0.59215688 -0.58883211  0.55011543]
 [-0.76124197 -0.39515796  0.51416034]
 [ 0.20528509  0.44720638  0.87055412]
 [-0.26686263 -0.86280497  0.42936223]
 [ 0.91424223  0.24593995  0.32198554]
 [-0.39083134  0.83127594  0.3952609 ]
 [-0.98460361  0.09323201  0.14786318]
 [ 0.59754846 -0.02731328  0.80136747]
 [ 0.26676707 -0.95252422  0.14674105]
 [-0.33213127  0.1657259   0.92856004]
 [ 0.54151196  0.39620269  0.74147706]
 [-0.23737174  0.69135686  0.68240776]
 [-0.9754604  -0.21994822  0.00998977]
 [-0.75217015  0.39591338  0.52677573]
 [-0.4738356  -0.13250632  0.8705871 ]
 [ 0.3828242   0.18446696  0.90521687]
 [ 0.06287377 -0.30268959  0.95101309]
 [ 0.40736892 -0.32483307  0.85354206]
 [ 0.00835214 -0.77899518  0.62697429]
 [-0.63389085  0.17834519  0.75257916]
 [ 0.33187429 -0.80924439  0.48475042]
 [ 0.26555688  0.86491999  0.42590252]
 [-0.11559352  0.45365613  0.88364826]
 [-0.46552182 -0.86076035  0.20586657]
 [ 0.54765401 -0.79829491  0.25059992]
 [ 0.67233139 -0.30977608  0.67231636]
 [-0.30804697 -0.66739083  0.67800926]
 [-0.92255043 -0.15295486  0.35426757]
 [-0.68129686  0.66739787  0.30069032]
 [-0.72819261 -0.11716411  0.67528371]
 [ 0.76793885 -0.56892483  0.29426937]
 [-0.71596847 -0.67856042  0.16414904]
 [-0.86827268  0.11905662  0.48158912]
 [ 0.86433006 -0.30318534  0.40126324]
 [ 0.85016541  0.50980957  0.13157878]
 [ 0.73969083  0.49829401  0.45228371]
 [-0.77208605  0.63549788  0.00505749]
 [ 0.08328025  0.70678329  0.70251105]
 [ 0.24893417 -0.07181093  0.96585453]
 [-0.25658977 -0.40553316  0.87732807]
 [-0.16116844 -0.08904087  0.98290206]
 [ 0.54870029  0.72975071  0.40791163]
 [-0.52079595  0.85085764  0.06937471]
 [-0.88742445  0.39675529  0.23465525]
 [-0.43039872  0.42534922  0.79613754]
 [ 0.37618308  0.9219296   0.09236939]
 [ 0.84851297 -0.03638374  0.52792231]
 [ 0.74825058  0.22093446  0.62554699]
 [-0.53927213  0.61071147  0.57984228]
 [ 0.98013998 -0.04149057  0.19391791]
 [ 0.02461784  0.17667659  0.98396105]
 [ 0.38988413  0.64317566  0.65902612]
 [ 0.034099   -0.56509772  0.82431901]
 [-0.53657134 -0.42997535  0.72609393]
 [ 0.64074659  0.75712022  0.12732942]
 [ 0.33356492 -0.59797103  0.72881074]
 [-0.87196661 -0.438395    0.21790836]
 [-0.57971596 -0.66583658  0.46967121]
 [ 0.05332005 -0.92972909  0.36436354]
 [ 0.9362261  -0.34243955  0.07884059]
 [-0.06798288  0.88064643  0.46887119]
 [-0.22158404  0.96259869  0.15589831]
 [-0.11817505 -0.98566508  0.12041184]
 [ 0.07775453  0.97909906  0.18793423]
 [ 0.59215688 -0.58883211  0.55011543]
 [-0.76124197 -0.39515796  0.51416034]
 [ 0.20528509  0.44720638  0.87055412]
 [-0.26686263 -0.86280497  0.42936223]
 [ 0.91424223  0.24593995  0.32198554]
 [-0.39083134  0.83127594  0.3952609 ]
 [-0.98460361  0.09323201  0.14786318]
 [ 0.59754846 -0.02731328  0.80136747]
 [ 0.26676707 -0.95252422  0.14674105]
 [-0.33213127  0.1657259   0.92856004]
 [ 0.54151196  0.39620269  0.74147706]
 [-0.23737174  0.69135686  0.68240776]
 [-0.9754604  -0.21994822  0.00998977]
 [-0.75217015  0.39591338  0.52677573]
 [-0.4738356  -0.13250632  0.8705871 ]
 [ 0.3828242   0.18446696  0.90521687]
 [ 0.06287377 -0.30268959  0.95101309]
 [ 0.40736892 -0.32483307  0.85354206]
 [ 0.00835214 -0.77899518  0.62697429]
 [-0.63389085  0.17834519  0.75257916]
 [ 0.33187429 -0.80924439  0.48475042]
 [ 0.26555688  0.86491999  0.42590252]
 [-0.11559352  0.45365613  0.88364826]
 [-0.46552182 -0.86076035  0.20586657]
 [ 0.54765401 -0.79829491  0.25059992]
 [ 0.67233139 -0.30977608  0.67231636]
 [-0.30804697 -0.66739083  0.67800926]
 [-0.92255043 -0.15295486  0.35426757]
 [-0.68129686  0.66739787  0.30069032]
 [-0.72819261 -0.11716411  0.67528371]
 [ 0.76793885 -0.56892483  0.29426937]
 [-0.71596847 -0.67856042  0.16414904]
 [-0.86827268  0.11905662  0.48158912]
 [ 0.86433006 -0.30318534  0.40126324]
 [ 0.85016541  0.50980957  0.13157878]
 [ 0.73969083  0.49829401  0.45228371]
 [-0.77208605  0.63549788  0.00505749]
 [ 0.08328025  0.70678329  0.70251105]
 [ 0.24893417 -0.07181093  0.96585453]
 [-0.25658977 -0.40553316  0.87732807]
 [-0.16116844 -0.08904087  0.98290206]
 [ 0.54870029  0.72975071  0.40791163]
 [-0.52079595  0.85085764  0.06937471]
 [-0.88742445  0.39675529  0.23465525]
 [-0.43039872  0.42534922  0.79613754]
 [ 0.          0.          0.        ]]

Both b-values and b-vectors look correct. Let’s now create the GradientTable.

gtab = gradient_table(bvals, bvecs=bvecs)

scene.clear()

We can also visualize the gradients. Let’s color the first shell blue and the second shell cyan.

colors_b1000 = window.colors.blue * np.ones(vertices.shape)
colors_b2500 = window.colors.cyan * np.ones(vertices.shape)
colors = np.vstack((colors_b1000, colors_b2500))
colors = np.insert(colors, (0, colors.shape[0]), np.array([0, 0, 0]), axis=0)
colors = np.ascontiguousarray(colors)

scene.add(actor.point(gtab.gradients, colors, point_radius=100))

window.record(scene=scene, out_path="gradients.png", size=(300, 300))
if interactive:
    window.show(scene)

Diffusion gradients.

References

[1]

Derek K. Jones, Mark A. Horsfield, and Andrew Simmons. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magnetic Resonance in Medicine, 42(3):515–25, 1999. URL: http://www.ncbi.nlm.nih.gov/pubmed/10467296, doi:10.1002/(SICI)1522-2594(199909)42:3<515::AID-MRM14>3.0.CO;2-Q.