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.90284067 -0.00537135  0.42994171]
 [ 0.30650366  0.39270992  0.86708386]
 [ 0.44292628  0.66189678  0.60473875]
 [ 0.25576943 -0.24342225  0.93558944]
 [-0.59102209  0.13376898  0.79548649]
 [-0.39932533  0.73247928  0.5513741 ]
 [-0.7724451   0.28255915  0.56876084]
 [-0.43101605 -0.87701786  0.21228481]
 [ 0.69641646 -0.56889939  0.4374444 ]
 [-0.91482454  0.26691231  0.30307404]
 [ 0.84912684  0.50428771  0.15709077]
 [-0.52636985  0.47472788  0.70538516]
 [-0.39173709 -0.76043681  0.51795551]
 [-0.61509287 -0.22142987  0.75672292]
 [-0.06903793 -0.18374275  0.98054697]
 [ 0.47898624  0.12108571  0.8694311 ]
 [ 0.54834976  0.79085073  0.27178606]
 [ 0.38765779 -0.72460467  0.56979778]
 [-0.02223442  0.82864694  0.55932985]
 [-0.95826803 -0.26841143  0.09837522]
 [-0.19586481  0.62265621  0.75758578]
 [ 0.96245902  0.18457228  0.19901183]
 [-0.650921   -0.66349878  0.36887832]
 [-0.3344662   0.31958473  0.88656526]
 [ 0.5413797  -0.18705062  0.81970732]
 [-0.51740282  0.81282074  0.26761308]
 [ 0.93233099 -0.25637744  0.2550089 ]
 [ 0.19497009  0.07438314  0.97798457]
 [-0.18790634 -0.66962971  0.71853132]
 [ 0.54940601 -0.46031236  0.69732744]
 [-0.6958735   0.56458087  0.44385642]
 [ 0.77928537 -0.60474684  0.16430329]
 [ 0.48961598 -0.81585091  0.3076743 ]
 [ 0.55595414 -0.83077027  0.02712493]
 [-0.66756998 -0.74442803  0.01331285]
 [ 0.1911333  -0.93204753  0.30782374]
 [-0.94039447 -0.07189224  0.33239998]
 [ 0.16442763  0.65431461  0.73813003]
 [ 0.78716367 -0.28371669  0.54761135]
 [-0.99780287  0.06085034  0.02620423]
 [-0.13794433 -0.90579574  0.40063131]
 [ 0.72608121  0.04440718  0.6861735 ]
 [-0.37543831 -0.05939797  0.92494214]
 [-0.77199191  0.61505562  0.16042157]
 [ 0.34695168  0.93468067  0.07743758]
 [-0.00463038  0.41881976  0.90805758]
 [-0.92410305  0.38177296  0.01682171]
 [-0.83446271 -0.52042196  0.18119869]
 [ 0.82143452  0.30483911  0.48199424]
 [ 0.26653566  0.86662547  0.42180449]
 [ 0.25866064 -0.54526366  0.79735953]
 [-0.80347351 -0.03446184  0.59434224]
 [-0.80210165 -0.35913718  0.47713041]
 [ 0.08293625 -0.79149227  0.60552586]
 [-0.5518957  -0.51576826  0.6552818 ]
 [ 0.59282933  0.38881183  0.70525084]
 [-0.10860558  0.12985734  0.98556679]
 [-0.22496065  0.91427787  0.33688082]
 [-0.01798683 -0.47532041  0.87962889]
 [-0.3352992  -0.38717381  0.85887769]
 [ 0.6979726   0.57781673  0.4230391 ]
 [ 0.04732149  0.9770314   0.2077747 ]
 [-0.07472268 -0.99122149  0.10907096]
 [-0.25511264  0.96684941  0.01094316]
 [ 0.90284067 -0.00537135  0.42994171]
 [ 0.30650366  0.39270992  0.86708386]
 [ 0.44292628  0.66189678  0.60473875]
 [ 0.25576943 -0.24342225  0.93558944]
 [-0.59102209  0.13376898  0.79548649]
 [-0.39932533  0.73247928  0.5513741 ]
 [-0.7724451   0.28255915  0.56876084]
 [-0.43101605 -0.87701786  0.21228481]
 [ 0.69641646 -0.56889939  0.4374444 ]
 [-0.91482454  0.26691231  0.30307404]
 [ 0.84912684  0.50428771  0.15709077]
 [-0.52636985  0.47472788  0.70538516]
 [-0.39173709 -0.76043681  0.51795551]
 [-0.61509287 -0.22142987  0.75672292]
 [-0.06903793 -0.18374275  0.98054697]
 [ 0.47898624  0.12108571  0.8694311 ]
 [ 0.54834976  0.79085073  0.27178606]
 [ 0.38765779 -0.72460467  0.56979778]
 [-0.02223442  0.82864694  0.55932985]
 [-0.95826803 -0.26841143  0.09837522]
 [-0.19586481  0.62265621  0.75758578]
 [ 0.96245902  0.18457228  0.19901183]
 [-0.650921   -0.66349878  0.36887832]
 [-0.3344662   0.31958473  0.88656526]
 [ 0.5413797  -0.18705062  0.81970732]
 [-0.51740282  0.81282074  0.26761308]
 [ 0.93233099 -0.25637744  0.2550089 ]
 [ 0.19497009  0.07438314  0.97798457]
 [-0.18790634 -0.66962971  0.71853132]
 [ 0.54940601 -0.46031236  0.69732744]
 [-0.6958735   0.56458087  0.44385642]
 [ 0.77928537 -0.60474684  0.16430329]
 [ 0.48961598 -0.81585091  0.3076743 ]
 [ 0.55595414 -0.83077027  0.02712493]
 [-0.66756998 -0.74442803  0.01331285]
 [ 0.1911333  -0.93204753  0.30782374]
 [-0.94039447 -0.07189224  0.33239998]
 [ 0.16442763  0.65431461  0.73813003]
 [ 0.78716367 -0.28371669  0.54761135]
 [-0.99780287  0.06085034  0.02620423]
 [-0.13794433 -0.90579574  0.40063131]
 [ 0.72608121  0.04440718  0.6861735 ]
 [-0.37543831 -0.05939797  0.92494214]
 [-0.77199191  0.61505562  0.16042157]
 [ 0.34695168  0.93468067  0.07743758]
 [-0.00463038  0.41881976  0.90805758]
 [-0.92410305  0.38177296  0.01682171]
 [-0.83446271 -0.52042196  0.18119869]
 [ 0.82143452  0.30483911  0.48199424]
 [ 0.26653566  0.86662547  0.42180449]
 [ 0.25866064 -0.54526366  0.79735953]
 [-0.80347351 -0.03446184  0.59434224]
 [-0.80210165 -0.35913718  0.47713041]
 [ 0.08293625 -0.79149227  0.60552586]
 [-0.5518957  -0.51576826  0.6552818 ]
 [ 0.59282933  0.38881183  0.70525084]
 [-0.10860558  0.12985734  0.98556679]
 [-0.22496065  0.91427787  0.33688082]
 [-0.01798683 -0.47532041  0.87962889]
 [-0.3352992  -0.38717381  0.85887769]
 [ 0.6979726   0.57781673  0.4230391 ]
 [ 0.04732149  0.9770314   0.2077747 ]
 [-0.07472268 -0.99122149  0.10907096]
 [-0.25511264  0.96684941  0.01094316]
 [ 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.