This post describes a reimplementation of the method described by Caruyer et al. (2013). If you already know what this is about, you can skip the information below and simply download the code from GitHub: multishell.py. If you don’t know, a good place to start is Tuch et al. (2004).

Background

When we plan a magnetic resonance imaging protocol, we often have to balance what we would like to have in terms of quantity and quality of images, versus what can realistically be obtained within limited time slots (e.g., 30 minutes or 1 hour) while minimizing discomfort to the research participants, most of whom would rather not spend more than an hour in a somewhat uncomfortable MRI exam table. For diffusion weighted imaging, we want as many directions in q-space as possible, and that the directions are uniformly sampled across all possible directions, which help subsequent modelling of various diffusion parameters, fiber orientation and fiber crossings.

In practice, we need to establish a set of directions that uniformly covers the surface of a sphere. Presets available within the scanner workstation, provided by the manufacturer, are sometimes helpful but these are not entirely customizable. For example, for Siemens scanners, the number of directions can be selected as 6, 10, 12, 20, 30, 64, 128, or 256. But what if the number of directions for your experiment that makes best use of the time available is none of these? The problem becomes more complicated if there are multiple b-values (e.g., in a multishell acquisition or, e.g., with diffusion spectrum imaging). How to find the directions that uniformly sample all possible directions, without repetitions across shells, and while taking into account the multiple b-values?

Note that this is not simply a matter of finding random points uniformly distributed on the surface of a sphere. That could be done by simply generating independent Gaussian noise in the 3D space, then scaling the distance of each point to the center of the coordinate system to a certain radius, e.g., 1 (that is, treating each point as a vector and scaling it to unit norm). Instead, what we want are directions that are as far apart from each other as possible, so that the resulting image for each provides maximally unique information.

Moreover, having an arbitrary number of points equally separated on the surface of a sphere is impossible. Only five configurations exist, with at most 20 points; these configurations correspond to vertex locations of the five Platonic polyhedra. We need, therefore, to be content with values that are as spaced as possible on the surface, already knowing they can’t be all equally spaced.

Cost function

For just one b-value (one shell), a possible solution is to treat each direction as a charged particle in the surface of a sphere centered at the center of the coordinate system. The positions of the particles can then be established such that the electrostatic repulsion forces are minimized (Conturo et al., 1996; Jones et al., 1999). Caruyer et al. (2013) extended the idea to multiple b-values (multiple shells). Using the same notation as in the original paper, consider $v$ the repulsion force between particles represented by their positions (vectors) $\mathbf{u}$ and $\mathbf{v}$ (Equation 3 in the paper):

$v(\mathbf{u}, \mathbf{v}) = \dfrac{1}{|\mathbf{u} - \mathbf{v}|^2} + \dfrac{1}{|\mathbf{u} + \mathbf{v}|^2}$

The above takes into account the direct repulsion between the two particles (first term, with a difference in the denominator) as well as repulsion between their antipodal representations (second term, with a sum in the denominator). Caruyer et al. (2013) then go on to propose the following cost functions for intra-shell and inter-shell electrostatic forces:

Intra-shell uniformity: Minimizing intra-shell electrostatic repulsion forces ensures that the points (directions) within each shell are evenly spread (Equation 2 if the paper):

$V_1 = \dfrac{1}{S} \sum_{s=1}^S \dfrac{1}{K_s^2} \sum_{i<j} v(\mathbf{u}_{s,i}, \mathbf{u}_{s,j})$

where $S$ is the number of shells, and $K_s$ is the number of points (directions) on shell $s$ .

Inter-shell uniformity: Minimizing inter-shell electrostatic repulsion forces ensures that points across all shells, when projected onto a unit sphere, are evenly spread (Equation 4 in the paper):

$V_2 = \dfrac{1}{K^2} \sum_{s,i} \sum_{t,j} v(\mathbf{u}_{s,i}, \mathbf{u}_{t,j})$

where $K = \sum_s K_s$ is the total number of points across all shells.

We have therefore two cost-functions that can be combined with a suitable factor that gives more or less weight for repulsions within- or between-shells:

$V = \alpha V_1 + (1-\alpha) V_2$

under the constraint that all points need to lie on the surface of a sphere of radius 1.

Gradient

The gradient of the cost function $V$ in the Cartesian space is provided in Equation 7 of the paper:

$\dfrac{\partial V}{\partial x_\mathbf{u}} = -2 \left(\dfrac{x_\mathbf{u}-x_\mathbf{v}}{|\mathbf{u}-\mathbf{v}|^4} + \dfrac{x_\mathbf{u}+x_\mathbf{v}}{|\mathbf{u}+\mathbf{v}|^4}\right)$

where $\partial x_\mathbf{u}$ are the axes directions in the Cartesian coordinate system, i.e., $x_\mathbf{u} \in \{x, y, z\}$ , and where again the first term within brackets considers the repulsion forces of the particles themselves, whereas the second considers their antipodal representations.

Antipodal points

The inclusion of the antipodal representations of the points in the cost function is necessary otherwise whenever an even number of directions is considered, they would always lineup perfectly, pairwise, in antipodal sides of the sphere. Such configuration would not be optimal for diffusion imaging as the scanner would sample twice along the same direction and, despite opposite signs, result in the same diffusion weighted images, and thus provide repeated information (Jones et al., 1999). Consider, for example, a minimal case of only two directions: the best configuration (lowest repulsion force, lowest energy) is the one in which the two points are as far away from each other as possible, thus being on opposite sides of the sphere, and both ultimately leading to the same image. By including the antipodal points in the cost function (the second term in the equation for $v(\mathbf{u}, \mathbf{v})$ , as well as the corresponding term in the gradient function) as new charged particles, we avoid this problem: for the minimal case of two particles, the two directions will be orthogonal one to another, thus each one providing maximally unique information.

The alpha parameter

The weight parameter $\alpha$ ranges between 0 and 1 and controls how much we want to emphasize intra-shell repulsion forces (greater $\alpha$ , towards 1) or inter-shell repulsion forces (smaller $\alpha$ , towards 0). While it may seem natural to choose $\alpha$ = 0.50, that not necessarily yields optimal results. In the implementation provided here, the default is set as $\alpha$ = 0.75, and the user can change its value until results are satisfactory.

Effect of varying the parameter alpha for S = 2 shells and K = 100 directions (50 per shell). Low values of alpha cause the vectors not to be well spaced within shell, but well separated between shells; large values of alpha tend to have the opposite effect. Although alpha = 0.50 may seem a reasonable tradeoff, the best result seem to depend on the number of directions, shells, and the number of directions per shell. In the above example, among the 5 configurations, alpha = 0.75 appears to provide the most uniform coverage within and across shells. In the left column, the directions are represented as arrows, in the middle column, the directions are presented as blobs on the surface of each shell, and in the third column, the directions are reprojected to be shown in a single sphere (as if all shells had the same radius). The figure shows only one hemisphere for clarity, but approximately half of the directions shown are in fact in the opposite hemisphere. Click **here** for a larger figure.

Simultaneous optimization

For a number $K$ of desired directions, the best solution can be found by starting with random vectors in the unit sphere, spread them across $S$ shells, and allow them to move freely on each iteration on the surface of the sphere until the energy (repulsion force) is minimized, under the constraint that all directions lie on the unit sphere. Optimization can use the Nelder-Mead method (Jones et al., 1999), albeit the constraint would need to be imposed manually at each iteration, or SLS-QP (Caruyer et al., 2013). The implementation provided here uses SLS-QP from SciPy.

Incremental optimization

The original paper also discusses an alternative optimization scheme whereby only one new direction is added to an existing set, while retaining all previous ones fixed (page 1537 of the paper, 1st column, just before the Results section, and references therein). By iteratively estimating each new direction, one can obtain a quasi-optimal solution that has as main advantage the fact that, if scanning is interrupted before all images have been collected, the collected images are still nearly optimal for later processing. In the implementation provided here, the option to optimize incrementally is provided. A new direction is added to the shell that currently has the highest relative deficit of directions in relation to the desired final number of directions per shell. Although SLS-QP could be used, it is simpler to bypass the constraint by operating with spherical coordinates and using L-BFGS-B instead; this also allows faster reordering, as discussed next.

Optimal reordering of directions

While simultaneous optimization guarantees directions as spread as possible on the surface of each shell, and also when all shells are considered, incremental optimization ensures that informative data can be obtained if the imaging acquisition needs to be interrupted. In principle, you would have to choose between the two. Instead, in the implementation provided here, an option that combines of both schemes is provided: first, all directions are found via simultaneous optimization, then, while keeping the first direction fixed the second direction is the one which, among those already identified, that has the lowest repulsion force in relation to the first; the process is repeated for the third direction, finding the one, among those previously identified and not yet already selected, that has the lowest repulsion force in relation to the previous two; the process is iterated until all directions have been reordered. This can be used also with existing directions (i.e., without any new optimization). To use this feature, use the option –reorder.

In practice

The reimplementation described in this article is called multishell.py. It’s one single file that can be made executable and run as a command directly from the shell (e.g., from bash or zsh). Alternatively, it can be imported into your Python scripts, such that you can invoke the functions within it directly (e.g., the optimizer functions).

The following Python modules are necessary:

NumPy
SciPy
matplotlib if you want to use the –plot option
mpl_toolkits, specifically, the mplot3d toolkit, if you want to plot the directions as blobs (default) or rings

Assuming the file has been downloaded, made executable, and is in the current directory, then usage/syntax information is available with:

./multishell.py --help

File formats

In this reimplementation, the directions can be saved in files that can be imported directly by Siemens (*.dvs), General Electric (*.dat), or Philips scanners (*.txt). The directions can also be saved in FSL format (*.bvec/bval, also used by BIDS), in MRtrix3 format (*.b), or in the tabular format used by two earlier implementations provided by Emmanuel Caruyer (this and this). Conversion between any of these formats is also possible.

Some file formats require the specification of the b-values for each shell before they can be saved (e.g., Siemens, GE, FSL, MRtrix3). Some file formats require the specification of the maximum b-value so they can be correctly be read (e.g., Siemens, GE). The tabular formats used by Caruyer’s implementations do not specify b-values, only shell indices (they can start at 0 or 1 depending on the implementation). Hence, conversion from/to certain formats may require additional options, namely –bvalues and –bmax (more on this below).

Creating and saving a new scheme

The options –new and –save create and save a set of directions, respectively. Assuming that the multishell.py file has been downloaded, made executable, and is in the current directory, the options to create a new sampling scheme are:

./multishell.py --new <method> <K> <S> <distr> ...

./multishell.py --new <method> <Ks> ...

where method can be simultaneous or incremental, K is the total number of directions, S is the number of shells, distr is how the directions should be distributed across shells (evenly, linearly, quadratically), and Ks is a list (provided between single quotes) of directions per shell. To save, use:

./multishell.py --save <prefix> <format> ...

where prefix is the prefix of the file(s) that will be created, and format is the format (siemens, ge, philips, fsl, mrtrix3, caruyer).

For example, say you want to specify a number of directions $K$ = 100, a number of shells $S$ = 3, that the number of directions should be spread across the shells linearly, that they should be optimized simultaneously to be optimally spaced, with weighting parameter $\alpha$ = 0.75, and with b-values for each shell 1000, 2000, and 3000 s/mm², then save the new scheme to a file for use in a Siemens scanner, then use:

./multishell.py --new simultaneous 100 3 linearly --alpha 0.50 --bvalues '[1000,2000,3000]' --save myprefix siemens

Or say you want 4 shells, the first with 15 directions, the second with 29, the third with 47, and the fourth with 68 directions, that they should be optimized incrementally to allow optimal data use in case the acquisition is interrupted, that the b-values should be 800, 1600, 2400, and 3200 s/mm², and want to save as a file for use in a GE scanner, use:

./multishell.py --new incremental '[15,29,47,68]' --bvalues '[800,1600,2400,3200]' --save myprefix ge

Or maybe you want as many shells as directions, with diffusion spectrum imaging (DSI) in mind:

./multishell.py --new simultaneous 250 250 uniform --bvalues 'list(range(1000,3500,10))' --save myprefix siemens

Loading an existing scheme

Files with existing directions can be loaded for conversion, plotting, or manipulation with the option –load:

./multishell.py --load <filename> <format> ...

where filename is the name of the file to be loaded (with extension) and format is their format (siemens, ge, philips, fsl, mrtrix3, caruyer). For the FSL format, either the bvec or bval file can be specified, but both must exist.

Format conversions

Say you have a pair of bvec/bval files from data collected elsewhere, and want to use the same directions in a new acquisition in a Siemens scanner, then use:

./multishell.py --load mydirs.bvec fsl --save myprefix siemens

In another example, suppose you have a set of directions from a Siemens scanner (which does not store actual b-values, just their relative magnitudes in relation to the maximum specified in the MRI parameter card), and want to use them in a Philips scanner (a format that stores the bvalues). When the sequence is run in the Siemens scanner, the bvalue is set as 2500 s/mm². Then use:

./multishell.py --load mysiemens.dvs siemens --bmax 2500 --save myprefix philips

Optimally reordering directions

Perhaps in the last conversion, you also wanted to reorder existing directions so that they allow optimal use of the data if the acquisition is interrupted. Then simply add the option –reorder:

./multishell.py --load mysiemens.dvs siemens --bmax 2500 --reorder --save myprefix philips

Newly created directions that were found by simultaneous optimization can also be reordered. Returning to some of the above examples:

./multishell.py --new simultaneous 100 3 linearly --bvalues '[1000,2000,3000]' --reorder --save myprefix siemens

./multishell.py --new simultaneous 250 250 uniform --bvalues 'list(range(1000,3500,10))' --reorder --save myprefix siemens

The option –alpha can be used when reordering directions if you do not want to stay with the default.

Adding or removing b0

If your existing file includes acquisitions with b-values = 0 (i.e., without the application of diffusion-weighted gradients), you can remove them with the option –removeb0.

./multishell.py --removeb0 ...

You can also add some b-values = 0 at the beginning, intersperse them regularly throughout the acquisition, and/or add some to the end with the option –addb0. Say you want to add 5 b0 volumes at the beginning and 2 at the end:

./multishell.py --addb0 '[5,0,2]' ...

Or maybe you want a few, say 8, interspersed throughout the acquisition:

./multishell.py --addb0 '[5,8,2]' ...

For most formats, the directions associated with b-values = 0 are simply (0,0,0), but for Philips scanners, these need to be unique, so random directions are used instead. Also note that for Philips, the first b-value must always be 0, so with this format, regardless whether –addb0 is used or not, if the first direction does not have a b-value = 0 at the time the file is saved, one will be added at the beginning.

Plotting

Basic plotting of the directions is available with the option –plot. When invoking from the shell, you can specify a file format for saving (e.g., pdf, eps, svg, png); four plots will be produced: with and without reprojecting all shells to unit sphere, and with the directions represented as arrows or as blobs.

./multishell.py --plot svg ...

If instead of a file format you provide the keyword interactive, it will open the plots in a window to allow interactive visualization, manipulation of the figure, and saving.

./multishell.py --plot interactive ...

Since antipodal points are included in the cost function, all directions can, for visualization purposes, be represented in only one hemisphere. However, that is only for display; the vectors themselves have equal likelihood of being in one hemisphere or another, thus covering the whole surface of the sphere.

Order of arguments and of operations

The order the arguments are provided is not relevant. If the same argument is repeated, it is the last one that counts. Inside the program, the operations happen in the following order:

The options –new or –load are executed, thus either creating or loading a set of directions.
If a file in a format that does not include b-values was loaded, or if a new set of directions was just created, the option –bvalues is parsed, thus defining the b-values for each shell.
The option –reorder is executed if provided, thus optimally reordering the directions.
The option –removeb0 is executed if provided, thus removing directions that have been assigned a b-value = 0.
The option –addb0 is executed if provided, thus adding b-values = 0 and corresponding directions.
The option –save is executed if provided, thus saving the set of directions to a file in the specified format.
The option –plot is executed if provided, thus showing the directions (interactive mode) or saved to a file in the specified format (png, pdf, svg, eps, and various others supported by matplotlib).

Integration with your Python code

If you prefer not to use the function from the system shell and prefer instead to call it from your Python code, you can do that too. Place the file multishell.py where you can import it, then call the functions as needed. A simple example is:

import multishell as ms
import numpy as np

# Number of directions
K = 120

# Number of shells
S = 3

# b-values of each shell
blist = [1000, 2000, 3000]

# Distribute the number of directions across shells
Ks = ms.directions_per_shell(K, S, 'quadratically')

# Simultaneous optimization of all directions
vectors, shells = ms.simultaneous_optimization(Ks, alpha=0.75, maxiter=100)

# Optimally reorder the directions
vectors, shells, idx = ms.optimal_reordering(vectors, shells)

# b-values for each corresponding direction
bvalues = np.array(blist)[shells[:,0]][:,None]

# Save in FSL format
fileprefix = 'mydirs'
ms.write_fsl(fileprefix, vectors, bvalues)

# Plot
ms.plot_directions(vectors, bvalues=bvalues, colorby='acquisition', style='quiver', reproject=False)

Why a new implementation?

Two helpful implementations provided by Caruyer are available. One is a webtool that allows up to 5 shells and up to 500 directions be downloaded after estimation via incremental optimization, but it doesn’t allow specification of the weighting parameter $\alpha$ , and only saves the directions in a tabular format that requires further manipulation before it can be used in the scanner. The other is a Python version that uses simultaneous optimization and provides finer control, but the shells are indexed differently from the webtool, and doesn’t include incremental optimization; it also doesn’t save directly in a format that can be used in common scanners. It allows plotting, but in my system figures appear incorrectly cropped. Conversion tools for Siemens and GE are available, but they require additional installation of packages for format changes that are relatively simple. For the webtool, there is always a risk that it may no longer be available in the future; as a concrete example, when the directions were downloaded for use in the Human Connectome Project (Sotiropoulos et al., 2013), it was located at a different address, hosted by the INRIA, now no longer accessible (broken link). Additionally, I could not understand Figure 2 in the original paper (compare with the figures above), and I had to implement it to understand what was going on. Finally, an all-in-one tool was much needed by our group as we plan on projects and need to find best use of scanner time. As this new implementation is helpful to us, it is hoped it will be helpful to you as well.

References

The most important reference is the original paper by Caruyer and colleagues, and which should be cited in case you use this reimplementation in your research:

Caruyer E, Lenglet C, Sapiro G, Deriche R. Design of multishell sampling schemes with uniform coverage in diffusion MRI. Magn Reson Med. 2013 Jun;69(6):1534-40. doi: 10.1002/mrm.24736. Epub 2013 Apr 26. PMID: 23625329; PMCID: PMC5381389.

Other relevant work cited above are:

Conturo TE, McKinstry RC, Akbudak E, Robinson BH. Encoding of anisotropic diffusion with tetrahedral gradients: a general mathematical diffusion formalism and experimental results. Magn Reson Med. 1996 Mar;35(3):399-412. doi: 10.1002/mrm.1910350319. PMID: 8699953.
Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn Reson Med. 1999 Sep;42(3):515-25. PMID: 10467296.
Sotiropoulos SN, Jbabdi S, Xu J, Andersson JL, Moeller S, Auerbach EJ, Glasser MF, Hernandez M, Sapiro G, Jenkinson M, Feinberg DA, Yacoub E, Lenglet C, Van Essen DC, Ugurbil K, Behrens TE; WU-Minn HCP Consortium. Advances in diffusion MRI acquisition and processing in the Human Connectome Project. Neuroimage. 2013 Oct 15;80:125-43. doi: 10.1016/j.neuroimage.2013.05.057. Epub 2013 May 20. PMID: 23702418; PMCID: PMC3720790.
Tuch DS. Q-ball imaging. Magn Reson Med. 2004 Dec;52(6):1358-72. doi: 10.1002/mrm.20279. PMID: 15562495.