Single- and Multiple-Shell Uniform Sampling Schemes for Diffusion MRI Using Spherical Codes

TL;DR

This paper introduces a novel spherical code approach for designing uniform sampling schemes in diffusion MRI, improving angular separation and rotational invariance over existing methods like EEM and GEEM.

Contribution

It presents the first application of spherical code formulation for single and multi-shell sampling in diffusion MRI, along with five algorithms to optimize sampling schemes.

Findings

SC methods achieve larger angular separation

SC methods show better rotational invariance

Codes and tutorials released in DMRITool

Abstract

In diffusion MRI (dMRI), a good sampling scheme is important for efficient acquisition and robust reconstruction. Diffusion weighted signal is normally acquired on single or multiple shells in q-space. Signal samples are typically distributed uniformly on different shells to make them invariant to the orientation of structures within tissue, or the laboratory coordinate frame. The Electrostatic Energy Minimization (EEM) method, originally proposed for single shell sampling scheme in dMRI, was recently generalized to multi-shell schemes, called Generalized EEM (GEEM). GEEM has been successfully used in the Human Connectome Project (HCP). However, EEM does not directly address the goal of optimal sampling, i.e., achieving large angular separation between sampling points. In this paper, we propose a more natural formulation, called Spherical Code (SC), to directly maximize the minimal angle between different samples in single or multiple shells. We consider not only continuous problems to design single or multiple shell sampling schemes, but also discrete problems to uniformly extract sub-sampled schemes from an existing single or multiple shell scheme, and to order samples in an existing scheme. We propose five algorithms to solve the above problems, including an incremental SC (ISC), a sophisticated greedy algorithm called Iterative Maximum Overlap Construction (IMOC), an 1-Opt greedy method, a Mixed Integer Linear Programming (MILP) method, and a Constrained Non-Linear Optimization (CNLO) method. To our knowledge, this is the first work to use the SC formulation for single or multiple shell sampling schemes in dMRI. Experimental results indicate that SC methods obtain larger angular separation and better rotational invariance than the state-of-the-art EEM and GEEM. The related codes and a tutorial have been released in DMRITool.

Tables

Table 1. TABLE I: Proposed algorithms to solve the problems listed in Section  II-B .

Algorithms	Problems	Initialization
ISC	P-C, P-D, P-O	No
IMOC	P-C, P-D	No
1-Opt	P-C, P-D	Yes (IMOC)
CNLO	P-C	Yes (IMOC+1-Opt)
MILP	P-C, P-D	Not required (IMOC + 1-Opt)

Table 2. TABLE II: Problems listed in Section  II-B and suggested algorithms.

Problems	Suggested algorithms	Initialization
P-C	IMOC + 1-Opt + CNLO	No
P-D	MILP	Not required (IMOC+1-Opt)
P-O	ISC	No

Table 3. TABLE III: Covering radii of multi-shell sampling schemes with 28 × 3 28 3 28\times 3 and 90 × 3 90 3 90\times 3 samples generated by various methods with the corresponding running time.

	Shell 1 (28)	Shell 2 (28)	Shell 3 (28)	Combined ($28\times 3$)	Running Time
SOA [20]	${15.9}^{\circ }$	${15.8}^{\circ }$	${16.7}^{\circ }$	${15.1}^{\circ }$	36 min
Solid Rotation [24]	${25.7}^{\circ }$	${25.7}^{\circ }$	${25.7}^{\circ }$	${7.4}^{\circ }$	not available
IGEEM [18]	${19.2}^{\circ }$	${19.7}^{\circ }$	${19.3}^{\circ }$	${4.7}^{\circ }$	0.3 s
GEEM [18]	${23.6}^{\circ }$	${24.3}^{\circ }$	${23.4}^{\circ }$	${12.3}^{\circ }$	3.7 s
ISC ($N=20481$)	${21.3}^{\circ }$	${19.3}^{\circ }$	${21.1}^{\circ }$	${10.5}^{\circ }$	0.3 s
MILP ($N=321$)	${23.8}^{\circ }$	${23.8}^{\circ }$	${24.3}^{\circ }$	${13.3}^{\circ }$	10 min
IMOC ($N=20481$)	${24.3}^{\circ }$	${24.3}^{\circ }$	${24.3}^{\circ }$	${14.0}^{\circ }$	10.5 s
IMOC + 1-Opt ($N=20481$)	${24.3}^{\circ }$	${24.4}^{\circ }$	${24.3}^{\circ }$	${14.0}^{\circ }$	12.8 s
IMOC + 1-Opt + CNLO	${\text{26.1}}^{\circ }$	${\text{26.3}}^{\circ }$	${\text{26.9}}^{\circ }$	${\text{14.4}}^{\circ }$	39.7 s
EEM (CAMINO) [11, 25]	${25.7}^{\circ }$	${25.7}^{\circ }$	${25.7}^{\circ }$	${15.6}^{\circ }$	not available
GRASSC 5 [34]	${27.8}^{\circ }$	${27.8}^{\circ }$	${27.8}^{\circ }$	${16.2}^{\circ }$	not available
	Shell 1 (90)	Shell 2 (90)	Shell 3 (90)	Combined ($90\times 3$)	Running Time
IGEEM [18]	${10.8}^{\circ }$	${10.3}^{\circ }$	${10.5}^{\circ }$	${2.4}^{\circ }$	0.6 s
GEEM [18]	${12.9}^{\circ }$	${13.3}^{\circ }$	${12.6}^{\circ }$	${7.2}^{\circ }$	11 min
ISC ($N=20481$)	${10.4}^{\circ }$	${9.7}^{\circ }$	${10.4}^{\circ }$	${4.6}^{\circ }$	0.6 s
MILP ($N=321$)	${13.3}^{\circ }$	${13.5}^{\circ }$	${13.3}^{\circ }$	${7.9}^{\circ }$	10 min
IMOC ($N=20481$)	${13.5}^{\circ }$	${13.5}^{\circ }$	${13.5}^{\circ }$	${7.8}^{\circ }$	7.8 s
IMOC + 1-Opt ($N=20481$)	${13.6}^{\circ }$	${13.5}^{\circ }$	${13.5}^{\circ }$	${7.8}^{\circ }$	29.6 s
IMOC + 1-Opt + CNLO	${\text{14.6}}^{\circ }$	${\text{14.6}}^{\circ }$	${\text{14.7}}^{\circ }$	${\text{8.4}}^{\circ }$	31 min
EEM (CAMINO) [11, 25]	${15.1}^{\circ }$	${15.1}^{\circ }$	${15.1}^{\circ }$	not available	not available
GRASSC 5 [34]	${15.7}^{\circ }$	${15.7}^{\circ }$	${15.7}^{\circ }$	not available	not available

Table 4. TABLE IV: Angular differences between estimated directions and ground-truth fiber directions using different schemes by various methods.

	SOA [20]	IGEEM [17, 18]	GEEM [18]	IMOC	IMOC + 1-Opt	IMOC + 1-Opt + CNLO
Angular Difference	${1.41}^{\circ }\pm {0.75}^{\circ }$	${1.67}^{\circ }\pm {0.81}^{\circ }$	${1.23}^{\circ }\pm {0.73}^{\circ }$	${1.43}^{\circ }\pm {0.77}^{\circ }$	${1.31}^{\circ }\pm {0.71}^{\circ }$	${\text{1.16}}^{\circ }\pm {\text{0.64}}^{\circ }$

Table 5. TABLE V: Normalized RMSE of DWI signal reconstruction using tensorial SPFI.

	SOA [20]	IGEEM [17, 18]	GEEM [18]	IMOC	IMOC + 1-Opt	IMOC + 1-Opt + CNLO
$28\times 3$, no noise	$0.46\%\pm 0.048\%$	$0.42\%\pm 0.046\%$	$0.38\%\pm 0.062\%$	$0.38\%\pm 0.0058\%$	$0.36\%\pm \text{0.054}\%$	$\text{0.35\%}\pm 0.055\%$
$90\times 3$, no noise	not available	$0.176\%\pm 0.0044\%$	$0.174\%\pm 0.0022\%$	$0.179\%\pm 0.016\%$	$0.176\%\pm 0.0035$	$\text{0.172\%}\pm \text{0.0021}\%$
$28\times 3$, SNR=$30$	$4.39\%\pm 0.58\%$	$4.33\%\pm 0.60\%$	$\text{4.31}\%\pm 0.58\%$	$4.35\%\pm 0.61\%$	$4.33\%\pm 0.59\%$	$\text{4.31}\%\pm \text{0.56}\%$
$90\times 3$, SNR=$30$	not available	$2.57\%\pm 0.34\%$	$2.55\%\pm \text{0.33}\%$	$2.58\%\pm 0.35\%$	$2.56\%\pm 0.34\%$	$\text{2.54}\%\pm \text{0.33}\%$

Equations

$$\underset{\{{\mathbf{u}}_i{\}}_{i=1}^K}{\min }\sum _{i<j}{f}_2({\mathbf{u}}_i,{\mathbf{u}}_j),f_{\alpha }({\mathbf{u}}_i,{\mathbf{u}}_j)=\frac{1}{\Vert {\mathbf{u}}_i-{\mathbf{u}}_j{\Vert }_2^{\alpha }}+\frac{1}{\Vert {\mathbf{u}}_i+{\mathbf{u}}_j{\Vert }_2^{\alpha }},$$

$$d_i(\{{\mathbf{u}}_l{\}}_{l=1}^K)=\underset{j\ne i,\forall j}{\min }\arccos |{\mathbf{u}}_i^T{\mathbf{u}}_j|.$$

$$d(\{{\mathbf{u}}_l{\}}_{l=1}^K)=\underset{i}{\min }{d}_i(\{{\mathbf{u}}_l{\}}_{l=1}^K)=\underset{j\ne i,\forall i,\forall j}{\min }\arccos |{\mathbf{u}}_i^T{\mathbf{u}}_j|,$$

$$\underset{\{{\mathbf{u}}_i\in \mathbb{D}{\}}_{i=1}^K}{\max }d(\{{\mathbf{u}}_i{\}}_{i=1}^K)$$

$$\underset{\{{\mathbf{u}}_{s,i}\in \mathbb{D}\}}{\max }w{S}^{-1}\sum _{s=1}^{S}d(\{{\mathbf{u}}_{s,i}{\}}_{i=1}^{K_s})+(1-w)d(\{{\mathbf{u}}_{s,i}{\}}_{i=1,\dots ,K_s;s=1,\dots ,S}),$$

$${\theta }^{\text{ub}}(2K)=\arccos \sqrt{4-{\csc }^2\left(\frac{\pi K}{6(K-1)}\right)}$$

$$\underset{\theta ,\{{\mathbf{u}}_i{\}}_{i=1}^K}{\max }\theta$$

$$\arccos (|{\mathbf{u}}_i^T{\mathbf{u}}_j|)\ge \theta ,\forall i<j\le K;$$

$${\mathbf{u}}_i^T{\mathbf{u}}_i=1,\forall i;$$

$$\underset{\theta ,\{{\mathbf{u}}_i{\}}_{i=1}^K}{\max }\theta$$

$$\arccos (|{\mathbf{u}}_i^T{\mathbf{u}}_j|)\ge \theta ,\text{if}|{\text{p}}_i^T{\text{p}}_j|\ge \cos (2{\delta }_0+{\theta }^{\text{ub}}),\forall i<j\le K;$$

$$\arccos (|{\mathbf{u}}_i^T{\text{p}}_i|)\le {\delta }_0,\forall i;$$

$${\mathbf{u}}_i^T{\mathbf{u}}_i=1,\forall i;$$

$$\underset{\{{\theta }_s\},{\theta }_0,\{{\mathbf{u}}_{s,i}\}}{\max }w\frac{1}{S}\sum _{i=1}^S{\theta }_s+(1-w){\theta }_0$$

$$\arccos (|{\mathbf{u}}_{s,i}^T{\mathbf{u}}_{s,j}|)\ge {\theta }_s,\text{if}|{\text{p}}_{s,i}^T{\text{p}}_{s,j}|\ge \cos (2{\delta }_0+{\theta }_s^{\text{ub}}),$$

$$\forall s,\forall i<j\le K_s;$$

$$\arccos (|{\mathbf{u}}_{s,i}^T{\mathbf{u}}_{s^{\prime },j}|)\ge {\theta }_0,\text{if}|{\text{p}}_{s,i}^T{\text{p}}_{s^{\prime },j}|\ge \cos (2{\delta }_0+{\theta }_0^{\text{ub}}),$$

$$\forall s<s^{\prime },\forall i\le K_s,\forall j\le K_{s^{\prime }};$$

$$\arccos (|{\mathbf{u}}_{s,i}^T{\text{p}}_{s,i}|)\le {\delta }_0,\forall s,\forall i;$$

$${\theta }_s\ge {\theta }_0,\forall s;$$

$${\mathbf{u}}_{s,i}^T{\mathbf{u}}_{s,i}=1,\forall s,i;$$

$$\underset{\theta ,\{h_i{\}}_{i=1}^N}{\max }\theta$$

$$\arccos (|{\mathbf{u}}_i^T{\mathbf{u}}_j|)\ge \theta -(2-h_i-h_j)M,\forall i>j$$

$${\theta }^{\text{lb}}\le \theta \le {\theta }^{\text{ub}}(2K)$$

$$\sum _{i=1}^N{h}_i=K;h_i=0,1,\forall i$$

$$\underset{\{{\theta }_s\},\{h_{s,i}\}}{\max }w{S}^{-1}\sum _{s=1}^S{\theta }_s+(1-w){\theta }_0$$

$$\arccos (|{\mathbf{u}}_i^T{\mathbf{u}}_j|)\ge {\theta }_s-(2-h_{s,i}-h_{s,j})M,\forall s,i>j$$

$$\arccos (|{\mathbf{u}}_i^T{\mathbf{u}}_j|)\ge {\theta }_0-(2-h_{s,i}-h_{s^{\prime },j})M,\forall s,s^{\prime },i>j$$

$${\theta }_s^{\text{lb}}\le {\theta }_s\le {\theta }^{\text{ub}}(2K_s),\forall s;{\theta }_0^{\text{lb}}\le {\theta }_0\le {\theta }^{\text{ub}}(2\sum _{s=1}^S{K}_s)$$

$$\sum _{i=1}^N{h}_{s,i}=K_s,\forall s;\sum _{s=1}^S{h}_{s,i}\le 1,\forall i;h_{s,i}=0,1,\forall i,s$$

$$\arccos (|{\mathbf{u}}_{s,i}^T{\mathbf{u}}_{s,j}|)\ge {\theta }_s-(2-h_{s,i}-h_{s,j})M,\forall s,i>j$$

$$\arccos (|{\mathbf{u}}_{s,i}^T{\mathbf{u}}_{s^{\prime },j}|)\ge {\theta }_0-(2-h_{s,i}-h_{s^{\prime },j})M,\forall s,s^{\prime },i>j$$

$$f_2({\mathbf{u}}_i,{\mathbf{u}}_j)$$

$$\in \left[\frac{1}{1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|-2\delta -{\delta }^2{)}^2},\frac{1}{1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|+2\delta +{\delta }^2{)}^2}\right]$$

$$f_2({\text{x}}_i,{\text{x}}_j)=\frac{1}{1-({\text{x}}_i^T{\text{x}}_j{)}^2}\in \left[\frac{1}{1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|-2\delta -{\delta }^2{)}^2},\frac{1}{1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|+2\delta +{\delta }^2{)}^2}\right]$$

$$\frac{1}{1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|+2\delta +{\delta }^2{)}^2}-\frac{1}{1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|-2\delta -{\delta }^2{)}^2}$$

$$=\frac{4|{\mathbf{u}}_i^T{\mathbf{u}}_j|(2\delta +{\delta }^2)}{(1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|+2\delta +{\delta }^2{)}^2)(1-(|{\mathbf{u}}_i^T{\mathbf{u}}_j|-2\delta -{\delta }^2{)}^2)}$$

$$\le \frac{4c_0(2\delta +{\delta }^2)}{(1-(c_0+2\delta +{\delta }^2{)}^2)(1-(c_0-2\delta -{\delta }^2{)}^2)}=L(c_0,\delta )$$

$$\left|\sum _{i<j}{f}_2({\mathbf{u}}_i,{\mathbf{u}}_j)-\sum _{i<j}{f}_2({\text{x}}_i,{\text{x}}_j)\right|\le K(K-1)L(c_0,\delta )$$

$${\text{x}}_1^T{\text{x}}_2$$

$$={\mathbf{u}}_1^T{\mathbf{u}}_2+{\mathbf{u}}_1^T({\mathbf{C}}_1^T+{\mathbf{C}}_2){\mathbf{u}}_2+{\mathbf{u}}_1^T{\mathbf{C}}_1^T{\mathbf{C}}_2{\mathbf{u}}_2$$

$$|{\text{x}}_1^T{\text{x}}_2-{\mathbf{u}}_1^T{\mathbf{u}}_2|\le c$$

Full text

I Introduction

Diffusion MRI (dMRI) is a unique imaging methodology to explore white matter structure and architecture in the human brain by measuring diffusion of water molecules. In dMRI, signal measurements are a limited number of samples of the diffusion signal $E(\mathbf{q})$ in 3D $\mathbf{q}$-space. Reconstruction in dMRI entails fitting the signal samples by using a well-designed model, including Diffusion Tensor Imaging (DTI) [1], and various models in High Angular Resolution Diffusion Imaging (HARDI) [2, 3, 4, 5, 6, 7, 8, 9, 10]. After model fitting, meaningful quantities can be estimated, including diffusion tensor and fractional anisotropy in DTI [1], the ensemble average propagator, diffusion Orientation Distribution Function (dODF), fiber ODF (fODF), and the Generalized Fractional Anisotropy (GFA) in HARDI [2, 3, 4, 5, 6, 7, 8, 9], etc. An appropriate data sampling scheme in $\mathbf{q}$-space is important for all such dMRI acquisition and reconstruction applications in order to recover as much information as possible by using a minimal number of measurements.

The optimal sampling scheme that yields the best reconstruction quality is dependent on the diffusion model, reconstruction algorithm, and the ground truth diffusion signal.

1. A typical DTI sampling scheme is a single shell scheme with about 30 samples and a low b value ($\leq$1200\text{,}\mathrm{s}\mathrm{/}\mathrm{m}\mathrm{m}^{2}) [[11](#bib.bib11)], while multi-shell schemes are preferred (or necessarily needed in some cases) with a higher maximal b value ($\geq$2000\text{\,}\mathrm{s}\mathrm{/}\mathrm{m}\mathrm{m}^{2}) in HARDI methods [5, 4, 10, 9].
2. Even for the same diffusion model, different reconstruction algorithms, including least squares, robust non-linear estimation with various regularizations, etc., may have different preferred schemes [12]. For a given diffusion model, there may be several model fitting (e.g., reconstruction) algorithms. Robust estimation may prefer larger b values due to its robustness to noise [13, 14], and compressed sensing based reconstruction could reduce the number of data samples, compared with traditional least square based reconstruction [15, 16].
3. Different variabilities of DWI signals may prefer different schemes. An efficient scheme should allocate more samples in the domains with larger signal changes than in domains with smaller signal changes. Error propagation analysis showed that the optimal b value for DTI is inversely proportional to the mean diffusivity of the latent signal [11]. Therefore, it is infeasible to develop a general optimal sampling scheme that works best for all signal variabilities and all reconstruction methods. However, a necessary property for an optimal sampling scheme is that the samples should be spherically uniformly distributed with no directional preference, because 1) an optimal sampling scheme should be invariant to the orientation of tissue structures and the laboratory coordinate frame in which the data is scanned;
4. we have to design one single scheme for all different variabilities of signals in human brain. Thus, existing works [11, 17, 18, 19, 20, 21, 22, 23, 24] focus on designing single or multiple shell spherically uniform sampling schemes with input parameters $S,\{K_{s}\}_{s=1}^{S}$, where $S$ is the number of shells and $K_{s}$ is the number of samples in the $s$-th shell. These parameters are user-determined based on the consideration of diffusion signals, models, reconstruction algorithms, noise level, and the scan time limit.

Uniform single shell sampling schemes are widely used in dMRI, where samples in $\mathbf{q}$-space are uniformly distributed on a sphere (i.e., a fixed b-value). Spherical tessellation using spherical polyhedrons can generate uniform single shell schemes with some specific fixed numbers of samples determined by polyhedron types [2]. On the other hand, the Electrostatic Energy Minimization (EEM) method proposed in dMRI by Jones et al., [11] is the most popular way to generate a uniform single shell sampling scheme with an arbitrary number of samples. EEM considers the samples as electrons on sphere, and estimates the sample configuration by minimizing the electrostatic energy based on Coulomb’s law:

[TABLE]

where the second term in $f_{\alpha}(\mathbf{u}_{i},\mathbf{u}_{j})$ is used because antipodal symmetric samples provide redundant information in dMRI data reconstruction. The problem in Eq. (1) is the so-called Thomson problem 222http://en.wikipedia.org/wiki/Thomson_problem. Some best known solutions of EEM have been collected in CAMINO [25]. For the multi-shell case, recently it was reported that staggered samples in different shells are better for reconstruction, compared with repeated spherical samples in different shells [17, 18, 19, 20, 24]. [19] and [20] tried to design a multi-shell scheme with $S$ shells by separating a single shell scheme into $S$ subsets respectively in different ways to minimize the electrostatic energy. [24] proposed to a rotate single shell scheme and use the rotated scheme in different shells. [17] and [18] generalized EEM from the single shell case to the multiple shell case, called Generalized EEM (GEEM), by considering the electrostatic energies both in each individual shell and in the combined shell with all samples. The obtained multi-shell schemes by GEEM showed large angular separations in shells than [19, 20], and have been successfully used in the Human Connectome Project (HCP) [26]. Although EEM and GEEM are widely used, the electrostatic energy formulation does not directly maximize the angular separation between measurements. Samples in $\mathbf{q}$-space are not electrons, and it is still not clear how electrostatic energy relates to dMRI reconstruction.

A good sampling scheme should have large angular separation such that the reconstruction has large angular resolution and good rotational invariance. Thus we propose a mathematically more natural way for sampling scheme design, which directly maximizes the minimal angular difference between sampling points, i.e., covering radius, on the unit sphere. This way to determine sample configuration on sphere is essentially called the Spherical Code (SC) formulation 333http://mathworld.wolfram.com/SphericalCode.html, or Tammes problem 444https://en.wikipedia.org/wiki/Tammes_problem. In many diffusion MRI applications, the information redundancy is considered in local spatial neighborhoods and local angular neighborhoods [27, 28, 29, 30], which is based on a natural assumption that if the samples $\mathbf{u}_{i}$ and $\mathbf{u}_{j}$ are close in the sphere, then the DWI signal samples $E(q\mathbf{u}_{i})$ and $E(q\mathbf{u}_{j})$ are also close. Thus, based on this assumption, maximizing the separation angles between nearest sampling points is related to minimizing the information redundancy in measured DWI samples. The SC formulation directly maximizes the minimal separation angles between the nearest samples. While in the EEM formulation, the electrostatic energy is determined by distances between all sample pairs, not just the nearest samples. Moreover, the covering radius (i.e., minimal separation angle between samples) is a natural and direct measurement of “angular resolution” of spherical sampling schemes, compared with the electrostatic energy. Some pioneer works of HARDI uses finer tessellation gradient table (i.e., larger number of samples with smaller separation angles) to investigate the non-Gaussianity of DWI signals [2, 31]. Wedeen et al., [32] used $\sqrt{\frac{4\pi}{K}}$ to roughly measure the angular resolution of schemes. Caruyer et al., [17, 18] evaluated schemes by different methods using the covering radius. Thus, it is natural to directly maximize the covering radius of sampling schemes.

SC is related with the famous thirteen spheres problem raised by Isaac Newton, and has been well studied in mathematics [33, 34, 35]. Tóth [33] proved an upper bound for the covering radius of $K$ samples on a sphere. There is also a collection of best known solutions for the SC problem with the antipodal symmetry constraint in $\mathbb{S}^{2}$ [34] 555http://neilsloane.com/grass/dim3/. It can be proven that the SC problem is equivalent to EEM problem in Eq. (1) when $\alpha\to\infty$. Thus, Papadakis et al., [36] tried to approximately solve SC problem by iteratively increasing $\alpha$ in Eq. (1). However as $\alpha$ increases, the optimization problem has more local minima and becomes more unstable because the cost function tends to be discontinuous and floating-point overflow occurs [36]. Since the SC formulation directly maximizes the minimal distance between samples, it yields larger angular separation than EEM formulation. However to our knowledge, the SC formulation and the collected best known solutions have not yet been used in the dMRI domain. Although SC is well-studied in the mathematical literature, its current formulation is limited to a continuous single shell, and is not applicable to multiple shells, nor to the discrete problem to uniformly extract subsets from an existing scheme.

In this paper, we propose the SC formulation in diffusion MRI with total five algorithms, not only for continuous problems to design single and multiple shell schemes, but also for discrete problems to uniformly sub-sample an existing scheme, and for ordering an existing scheme. Preliminary parts of this work were reported in two conference papers [37, 38]. In this paper, we categorize sampling scheme design problems in dMRI, propose total 5 algorithms in a unified framework of minimizing the mean covering radius, and provide additional examples, results, derivations, and insights.

The paper is organized as follows. Section II proposes the SC formulation and shows how the existing sampling scheme design problems in dMRI can be formulated as optimization problems using SC. Section III proposes 5 algorithms to solve the problems in Section II. Section IV demonstrates some experimental results of synthetic data and real data.

II Spherical Code (SC) and Sampling Problems in DMRI

II-A Spherical Code (SC) Formulation

For a given set of spherical samples $\{\mathbf{u}_{i}\in\mathbb{S}^{2}\}_{i=1}^{K}$, we define the covering radius 3 for the $i$-th sample as

[TABLE]

where the absolute value operator is used because the antipodal symmetric samples have the same role in dMRI data reconstruction. We also define the covering radius for all samples as the minimal angular distance between all sample pairs, i.e.,

[TABLE]

The SC formulation on a single shell is to find $K$ samples $\{\mathbf{u}_{i}\}_{i=1}^{K}$ such that the covering radius is maximized:

[TABLE]

where $\mathbb{D}\subseteq\mathbb{S}^{2}$ is the solution domain. If $\mathbb{D}=\mathbb{S}^{2}$, Eq. (4) is a continuous SC (CSC) problem for selecting $K$ points from the continuous sphere $\mathbb{S}^{2}$. If $\mathbb{D}=\{\mathbf{u}_{n}\}_{n=1}^{N}$, a set of $N$ predetermined points on $\mathbb{S}^{2}$, then Eq. (4) is a discrete SC (DSC) problem for selecting $K$ from $N$ points. For CSC, [34] proposed to iteratively optimize a continuous cost function that approximates the original cost function in Eq. (4). The authors also released a collection of best known solutions for the SC problem in $\mathbb{S}^{2}$ [34] 5. Note that the original SC problem only considers CSC, i.e., $\mathbb{D}=\mathbb{S}^{2}$, while in this paper we propose both CSC and DSC formulations and generalize them to multi-shell case for designing sampling schemes in dMRI.

For multi-shell sampling, SC is to find a set of points $\{\mathbf{u}_{s,i}\}$ by maximizing the weighted mean of covering radii:

[TABLE]

where $S$ is the number of shells, $K_{s}$ is the number of points on the $s$-th shell, $\mathbf{u}_{s,i}$ is the $i$-th point on the $s$-th shell, and $w$ is a weighting factor between the mean covering radius of the $S$ shells and the covering radius of the combined shell containing all points from the $S$ shells. $w$ is normally set to $0.5$. Because of the second term, the estimated samples in different shell are staggered if $w\neq 1$.

II-B Uniform Sampling Scheme Design Problems in dMRI

Spherical sampling problems in diffusion MRI can be classified into two categories, the continuous category (P-C) which determines samples in the continues sphere, and the discrete category (P-D), i.e., the sub-sampling problem, which determines samples in a discrete set of spherical samples. These two problem categories include the following sampling scheme design problems.

• •

Single-shell continuous problem (P-C-S). Given a number $K$, how does one determine $K$ uniform samples on sphere, such that these samples are separated as far as possible? The problem is for single shell scheme design. P-C-S can be solved by CSC in Eq. (4) with $\mathbb{D}=\mathbb{S}^{2}$.

• •

Multi-shell continuous problem (P-C-M). Given $S$ numbers $\{K_{s}\}_{s=1}^{S}$, how does one determine a scheme with $S$ shells, $K_{s}$ points in the $s$-th shell, such that samples are separated as far as possible not only in each single shell, but also in the combined shell with all samples? The problem is for multi-shell sampling scheme design. P-C-M can be solved by CSC in Eq. (5) with $\mathbb{D}=\mathbb{S}^{2}$.

• •

Single subset from single set problem (P-D-SS). Given $N$ known points $\{\mathbf{u}_{i}\}_{i=1}^{N}$ on sphere, how does one uniformly select $K$ samples from the given $N$ samples, such that these samples are separated as far as possible? The problem is to reduce the number of samples in an existing single shell scheme. It also relates to P-C-S, because after discretizing the continuous sphere using many points, solving P-D-SS is an approximation of solving P-C-S. P-D-SS can be solved by DSC in Eq. (4) with $\mathbb{D}=\{\mathbf{u}_{i}\}_{i=1}^{N}$.

• •

Multiple subsets from multiple sets problem (P-D-MM). Given points $\{\mathbf{u}_{s,i}\}$, $i=1,2,\dots,N_{s},s=1,2,\dots,S$, how does one uniformly select $K_{s}$ from the $N_{s}$ samples for the $s$-th shell, such that samples are separated as far as possible not only in each single shell, but also in the combined shell with all samples? P-D-MM is to reduce the number of samples in an existing multi-shell scheme. It can be solved by DSC in Eq. (5) with $\mathbb{D}=\{\mathbf{u}_{s,i}\}$.

• •

Multiple subsets from single set problem (P-D-MS). Given a single set of points $\{\mathbf{u}_{i}\}_{i=1}^{N}$ on sphere, how does one uniformly select multiple sets of samples from a given single set of samples, such that the samples are separated as far as possible not only in each set, but also in the combined set with all samples? It relates to P-C-M by discretizing the continuous sphere using many points. P-D-MS can be solved by DSC in Eq. (5) with $\mathbb{D}=\{\mathbf{u}_{i}\}_{i=1}^{N}$. [19, 20] actually tried to solve P-D-MS by minimizing electrostatic energy in different heuristic search ways. We propose to solve it by using a better algorithm called Mixed Integer Linear Programming (MILP).

• •

Acquisition ordering problem (P-O). Given a single or multiple shell scheme, how does one optimize the acquisition order of the samples, such that an interruption at any point of the acquisition can still yield a nearly uniform sampling scheme [22, 23]? This problem can be solved by using incremental SC (ISC) as described in III-A.

P-C-S is a special case of P-C-M, and P-D-SS is a special case of P-D-MM (with $S=1$ and $w=1$). In the following context, although P-O belongs to the discrete category, we will use P-C to denote P-C-S and P-C-M, P-D to denote P-D-SS, P-D-MM, and P-D-MS.

III Algorithms for the SC formulation

In this section, we propose five algorithms to solve the above sampling scheme design problems. Table I summarizes the algorithms applicable to the corresponding problems. Note that although one algorithm may work for several problems, only the problems in bold are suggested for good performance. ISC and IMOC do not need an initialization. MILP does not require an initialization, although a good initialization may reduce the computation time. The algorithms in brackets in the third column are suggested for good initializations. The algorithms can also be cascaded, i.e., one algorithm can use the result by another algorithm as the initialization. Table II lists the suggested algorithms for P-O, P-C and P-D, where IMOC + 1-Opt + CNLO means CNLO uses the scheme by IMOC + 1-Opt as its initialization.

III-A Incremental SC (ISC)

The similar incremental strategy used for incremental EEM (IEEM) [39] and incremental GEEM (IGEEM) [17, 18] can be also used to propose a greedy incremental SC estimation, called incremental SC (ISC). ISC works for P-D in both Eq. (4) and Eq. (5) with $\mathbb{D}=\{\mathbf{u}_{i}\}_{i=1}^{N}$. At step $k$, we estimate one point $\mathbf{u}\in\mathbb{D}$ that maximizes the cost function based on the $k-1$ points estimated in previous iterations. ISC can be also applied to P-C, i.e., $\mathbb{D}=\mathbb{S}^{2}$, by approximating $\mathbb{S}^{2}$ using a large number of uniformly distributed points. In practice, we normally use $20481$ points from a 7 order tessellation of the icosahedron. ISC can obtain reasonable schemes in seconds. For P-O to optimize acquisition order of a given scheme $\{\mathbf{u}_{s,i}\}$, we can use ISC with $\mathbb{D}=\{\mathbf{u}_{s,i}\}$. In this way, the first partial set of samples are nearly uniformly separated. It was reported in [40] that optimal ordering reduces the risk of systematic errors due to gradient heating and subject motion.

III-B Iterative Maximum Overlap Construction (IMOC)

Incremental methods (IEEM [39], IGEEM [18] and ISC) are to obtain a scheme with reasonable uniform coverage when the acquisition is terminated and only a subset of samples are used. These incremental methods all have the same limitations.

1. Because they are developed to work reasonably well for any subset with first several samples, the cost function (i.e., electrostatic energy or covering radius) of the full scheme is not optimal.
2. When using a finer uniform sample set, the cost function of the obtained sampling scheme may not be improved, which is beyond our expectation.

We propose a sophisticated greedy algorithm, called Iterative Maximum Overlap Construction (IMOC), to iteratively search the covering radius and sequentially generate samples with maximal spherical overlap on the sphere. IMOC works for both P-C and P-D. See the IMOC in Alg. 1 for designing single shell schemes. IMOC uses MOC in Alg. 2 to verify whether a candidate covering radius $\theta$ is able to construct an acceptable scheme with $K$ samples. We define the coverage set of a point x as $C(\textbf{{x}},\theta)=\{\textbf{{y}}\mid\arccos(|\textbf{{y}}^{T}\textbf{{x}}|)<\theta,\textbf{{y}}\in\mathbb{S}^{2}\}$. With a given candidate $\theta$, MOC generates samples one by one. In the $i$-th step, MOC finds the best point $\mathbf{u}_{i}$ whose coverage $C(\mathbf{u}_{i},\theta)$ has the largest overlap with the existing total coverage set CS, and then adds $C(\mathbf{u}_{i},\theta)$ into CS. See Alg. 2. If all $K$ points can be obtained by MOC using $\theta$, then the best covering radius is no less than $\theta$, which means $\theta$ can be increased. If MOC cannot find $K$ points using $\theta$, then $\theta$ can be decreased. Therefore, IMOC performs a binary search of the covering radius $\theta$ from $(0,\theta^{\text{ub}}(2K))$, where $\theta^{\text{ub}}(2K)$ in Eq. (6) is an upper bound of the covering radius with $K$ points with antipodal symmetry consideration (i.e., total $2K$ points) [33].

[TABLE]

See Fig. 1 (a) for a demonstration of MOC when solving P-C-S in $\mathbb{S}^{1}$ in the 2D case. It is obvious that IMOC yields the globally optimal solution in $\mathbb{S}^{1}$, while ISC and IGEEM cannot.

IMOC and MOC algorithms for multiple shell sampling scheme design are shown in Alg. 3 and 4. For multiple shell case, IMOC needs to search $(\theta_{0},\theta_{1},\dots,\theta_{S})$ in $S+1$ dimensional space, where $\theta_{0}$ is the covering radius of the combined shell with all samples and $\theta_{s}$ ($s>0$) is the covering radius of the $s$-th shell. There is a compromise between $\theta_{0}$ and $\theta_{s}$ ($s>0$). The cost function in Eq. (5) essentially converts a multi-objective optimization problem into a single objective optimization using a weighted summation. There are infinite Pareto optimal covering radii to the multi-objective optimization problem. IMOC only searches the covering radii in one dimensional space determined by the upper bounds $\{\theta_{s}^{\text{ub}}\}_{s=0}^{S}$, which in practice provides good covering radii for both individual shells and the combined shell. See Alg. 3 and Fig. 1 (b) and (c). MOC in Alg. 4 first selects one point $\mathbf{u}_{s,1}$ for each shell. Then it selects the best shell $s$ and the best point $\mathbf{u}_{s,i}$ whose coverage set $C(\mathbf{u}_{s,i},\theta_{s})$ has the largest overlap area with the existing total coverage set $\text{CS}_{s}\cup\text{CS}_{0}$ among all shells and all candidate points outside the existing total coverage set.

Similarly with ISC, for P-C, IMOC can analytically select the best point in each iteration in 2D case in $\mathbb{S}^{1}$. See Fig. 1 (a). However, in $\mathbb{S}^{2}$, IMOC requires a fine uniform discretization of the sphere to find the best point in each step with the largest overlap area by counting the number of points in overlaping sets. For efficient implementation of IMOC, we consider 3 issues.

1. KD-tree is used for efficient nearest neighbor search.
2. The overlap area $A_{s}(\textbf{{x}})$ for a candidate point x in the $s$-th shell can be stored for next step. In the $i+1$-th step, $A_{s}(\textbf{{x}})$ is different from the $i$-th step and needs to be re-calculated, only if the determined point $\mathbf{u}_{s^{\prime},i}$ in the $i$-th step satisfies $\arccos(|\textbf{{x}}^{T}\mathbf{u}_{s^{\prime},i}|)<\theta_{0}$, if $s\neq s^{\prime}$, or $\arccos(|\textbf{{x}}^{T}\mathbf{u}_{s^{\prime},i}|)<\theta_{s}$, if $s=s^{\prime}$.
3. In MOC, for P-C, the candidate points are selected from a small “outside surface” of CS, i.e., $\{\textbf{{x}}\ |\ \textbf{{x}}\notin\text{CS},\exists\textbf{{y}}\in\text{CS},\text{s.t.}\arccos(|\textbf{{x}}^{T}\textbf{{y}}|)\leq\delta\}$, instead of the entire complementary set $(\mathbb{S}^{2}-\text{CS})$. $\delta$ is set as twice of the covering radius of the fine uniform sampling set used in IMOC. It is obvious that in the 2D case, this modification does not change the result because the point in each step of MOC is always in the “outside surface”. See Fig. 1 (a). For P-C-S and P-C-M in $\mathbb{S}^{2}$, the results by IMOC do not change neither in our experiments, although we cannot give a proof. These three issues significantly speed up MOC, because the first one reduces the time for each neighborhood search, and the other two largely reduce the number of neighborhood searches. In our implementation, IMOC only requires a few seconds on an ordinary laptop.

III-C 1-Opt Algorithm

We propose a greedy algorithm called 1-Opt. For both P-C and P-D, with a given initialization scheme, the 1-Opt algorithm at each step finds the best point pair ($\mathbf{u}_{s,i}$ from current result and $\textbf{{x}}\in\mathbb{D}$) such that after setting $\mathbf{u}_{s,i}$ as x the covering radius for $\mathbf{u}_{s,i}$ can be best increased. For single shell case, 1-Opt is to increase $d_{i}$ defined in Eq. (2). For multi-shell case, 1-Opt is to increase both covering radii of $\mathbf{u}_{s,i}$ in the combined shell, i.e., $d_{s,i}=\operatorname*{arg\,max}_{(s,i)\neq(s^{\prime},j),\forall j,\forall s^{\prime}}\arccos|\mathbf{u}_{s,i}^{T}\mathbf{u}_{s^{\prime},j^{\prime}}|$ and in the $s$-th shell, i.e., $d_{i}=\operatorname*{arg\,max}_{j\neq i,\forall j}\arccos|\mathbf{u}_{s,i}^{T}\mathbf{u}_{s,j}|$. See Alg. 5. 1-Opt is to increase the covering radius for every individual sample. 1-Opt requires an initialization, and a good initialization normally results in a good result. In practice, we use the result by IMOC as the initialization of 1-Opt.

Note that if IMOC is used to design a scheme with many samples for P-C, then there sometimes is a hole area in the result of IMOC. This is because of the discretization used in IMOC for P-C. By applying 1-Opt to the IMOC result, even though 1-Opt sometimes cannot increase the overall covering radius for the whole scheme, it can better arrange the points to fix the hole. See Fig. 3.

III-D Constrained Non-Linear Optimization (CNLO)

For P-C, SC problem in Eq. (4) can be converted into a continuous optimization, called Constrained Non-Linear Optimization (CNLO):

[TABLE]

where Eq. (7b) means all separation angles are larger than $\theta$, and Eq. (7c) means the estimated samples are on the unit sphere. After the problem is solved, the solution $\theta^{*}$ is the estimated covering radius and there exists $i$ and $j$ such that $\arccos|(\mathbf{u}_{i}^{*})^{T}\mathbf{u}_{j}^{*}|=\theta^{*}$. Note that both cost function and constraints are continuous.

CNLO in Eq. (7) has total $\frac{K(K-1)}{2}$ inequality constraints. Thus for P-C-S or P-C-M with a large number of samples, CNLO may be very slow due to large number of constraints. CNLO is a local optimization method which requires an initialization. Thus for a given initialization $\{\textbf{{p}}_{i}\}$ we propose a sped up version of CNLO to reduce the number of inequality constraints via

[TABLE]

where $\delta_{0}$ controls the difference between the final results $\{\mathbf{u}_{i}\}$ and the initialization $\{\textbf{{p}}_{i}\}$, and the inequality constraints which have two initialization points far from $2\delta_{0}+\theta^{\text{ub}}$ do not need to be considered in Eq. (8b), because $\theta^{\text{ub}}$ is an upper bound. If $\delta_{0}=0$, then the result is just the initialization. If $\delta_{0}=\pi/2$, then the result of Eq. (8) is the same as the result of Eq. (7). When setting $\delta_{0}$ as a small float number, the number of constraints can be largely reduced and CNLO can be sped up.

Similarly, with a given initialization $\{\textbf{{p}}_{s,i}\}$, P-C-M in Eq. (5) can be solved using CNLO in Eq. (9):

[TABLE]

$\delta_{0}$ controls the difference between the result and the initialization. Eq. (9b) means that the separation angles of samples in the $s$-th shell are larger than its corresponding covering radius $\theta_{s}$. Eq. (9c) means that the separation angles of samples in two different shells are larger than the covering radius $\theta_{0}$ for the combined shell with all samples. Eq. (9e) means that the covering radii for all single shells are always larger than the covering radius for the combined shell, in other words, Eq. (9e) means that the separation angles of samples in the same shell are larger than $\theta_{0}$.

CNLO in Eq. (8) and Eq. (9) are continuous non-convex constrained optimization problems. We use sequential quadratic programming (SQP) to solve CNLO. In each step, SQP solves a quadratic programming problem which locally approximates the original non-linear problem. In practice, we use the SQP solver implemented in NLOPT library [41] (NLOPT_LD_SLSQP). CNLO yields a locally optimal solution with a given initialization. We suggest to use the result by IMOC + 1-Opt as the initialization. In practice, we set $\delta_{0}=0.1$, and keep performing CNLO by updating the initialization $\{\textbf{{p}}_{i}\}$ (or $\{\textbf{{p}}_{s,i}\}$) using the previously obtained result.

III-E Mixed Integer Linear Programming (MILP)

We propose a Mixed Integer Linear Programming (MILP) algorithm to solve P-D. For $\mathbb{D}=\{\mathbf{u}_{n}\}_{n=1}^{N}$,  Eq. (4) can be solved using MILP:

[TABLE]

where $h_{i}=1$ means that $\mathbf{u}_{i}$ is selected as one of the $K$ points, $\theta^{\text{lb}}$ is a lower bound of the covering radius, which can be set to [math] or the covering radius from an existing sampling scheme, $\theta^{\text{ub}}(2K)$ is an upper bound of the covering radius in Eq. (6), and $M$ is the difference between the maximal distance (denoted as $d_{\text{max}}$) and the minimal distance (denoted as $d_{\text{min}}$) of any two points $\mathbf{u}_{i},\mathbf{u}_{j}\in\mathbb{D}$, $i\neq j$. After solving MILP, the solution of $\theta$, denoted as $\theta^{*}$, is the covering radius of the selected $K$ samples. The constraint in Eq. (10b) only takes effect when two points are both selected, i.e., $h_{i}=h_{j}=1$, and is automatically satisfied when $h_{i}=0$ or $h_{j}=0$, because we set $M$ large enough such that $\arccos(|\mathbf{u}_{i}^{T}\mathbf{u}_{j}|)\geq d_{\text{min}}\geq\theta^{*}-(2-h_{i}-h_{j})M$, when $h_{i}+h_{j}\leq 1$.

Similarly,  Eq. (5) with $\mathbb{D}=\{\mathbf{u}_{s,i}\}$ can be solved using MILP:

[TABLE]

where $h_{s,i}=1$ means that $\mathbf{u}_{i}$ is selected as one of the $K_{s}$ points on the $s$-th shell, $\theta^{\text{lb}}_{s}$ and $\theta^{\text{ub}}(2K_{s})$ are the lower and upper bounds of the covering radius on the $s$-th shell, and $\theta^{\text{lb}}_{0}$ and $\theta^{\text{ub}}(2\sum_{i=1}^{S}K_{s})$ are the lower and upper bounds for the combined shell with all points. Similarly with Eq. (10b), constraints in Eq. (11b) and Eq. (11c) are respectively for the first and second terms in Eq. (11a). $\sum_{s=1}^{S}h_{s,i}\leq 1$ makes $\mathbf{u}_{i}$ to be selected at most in one shell such that the estimated samples are staggered across shells. Note that the above MILP in Eq. (11) is for P-D-MS. For P-D-MM, Eq. (11b) and Eq. (11c) need to be replaced by Eq. (12a) and Eq. (12b)

[TABLE]

where $h_{s,i}=1$ means $\mathbf{u}_{s,i}$ is selected as one of the $K_{s}$ points on the $s$-th shell, and $\sum_{s=1}^{S}h_{s,i}\leq 1$ can be removed from Eq. (11e).

Note that although in general the MILP is NP hard, it can be solved with a globally optimal solution using branch and bound method which iteratively solves the relaxed LP program. In our implementation, we solve Eq. (10) and Eq. (11) using GUROBI [42], which can yield the globally optimal solution or at least a reasonable solution within 10 minutes for P-D. In practice, we progressively increase the lower bound $\theta^{\text{lb}}$ based on the solutions estimated in previous iterations to find a better feasible solution within 10 minutes which is good enough in experiments.

IV Experiments

IV-A Effect of Discretization in Incremental Greedy Methods

Incremental greedy methods (i.e., IEEM [39], IGEEM [18], ISC, and IMOC), all require a uniform sample set as a discretization of $\mathbb{S}^{2}$ to determine the best sample in each step. It is obvious that a finer uniform set may not yield better results for ISC in 2D case. We test the discretization effect in $\mathbb{S}^{2}$ for both single shell and multi-shell scheme design. We consider two uniform sample sets respectively with $81$ samples and $20482$ samples from spherical tessellation. With these two uniform sets, IEEM, ISC, IMOC, and IMOC + 1-Opt are performed to generate single shell scheme with $K\in[5,80]$, where IMOC + 1-Opt means 1-Opt with the initialization by IMOC. The left subfigure of Fig. 2 depicts the covering radii of schemes obtained by different methods with different uniform sets. It also shows the best known single shell schemes using the SC concept collected by Dr. Sloane 5, denoted as GRASSC, and the best known single shell schemes using EEM collected in CAMINO [25]. When $K\in[50,80]$ is close to $81$ (the number of uniform samples in the coarse uniform set), covering radii by IEEM and ISC are not improved using the finer set, and when $K<50$ is far from $81$, IEEM and ISC obtain larger covering radii using finer discretization. When using the finer uniform set, covering radii by IMOC and IMOC + 1-Opt are improved, close to the best known solutions in GRASSC, and are much larger than covering radii of the schemes in CAMINO. ISC, IMOC, IMOC + 1-Opt are performed to generate schemes with 3 shells, $K\in[5,25]$ per shell by using these two uniform sets. The right subfigure of Fig. 2 shows mean of covering radii of three shells and the covering radii of the combined shell. It also demonstrates that when $K\times 3$ is close to $81$, covering radii of the combined shell by ISC are not improved using the finer set, while IMOC has no such issue. The experiment shows that IMOC and IMOC + 1-Opt generally obtain better results than ISC and IEEM for P-C. CNLO can further improve the results by IMOC + 1-Opt, and the covering radii by IMOC + 1-Opt + CNLO obtains the best results, very closed to the solutions in GRASSC in single shell case. Thus, as shown in Table II, ISC is suggested for P-O, and IMOC + 1-Opt + CNLO is suggested for P-C.

IV-B Covering Radii in Sampling Schemes

We evaluate the proposed method by generating multi-shell schemes with $K\times 3$ samples. We test two cases with $K=28$ and $K=90$, considering these two cases were also used for evaluation in [18], and the multi-shell scheme with $90\times 3$ samples by IGEEM was also used in HCP [18, 26].

Table III shows the covering radii of the schemes by different methods, where IMOC + 1-Opt + CNLO means CNLO using the result by IMOC + 1-Opt as the initialization. Sparse and Optimal Acquisition (SOA) in [20], which separates a uniform scheme into several sets by using a greedy algorithm, is performed using the software released by the author 666https://sites.google.com/site/hispeedpackets/Home/soa-design. The result by SOA for the scheme with $90\times 3$ samples is not shown, because it did not terminate after running for 2 days. The solid rotation method in [24] starts from the single shell scheme with $28$ samples by EEM, and finds the optimal three rotations for the 3 shells, such that the electrostatic energy in the combined shell is minimized. The scheme by the solid rotation method is the 3 rotated EEM schemes, and the covering radii of the scheme were extracted from [18]. Note that we run GEEM by using the latest software released by the authors 777https://github.com/ecaruyer/qspace, which obtains better covering radii for individual shells but a slightly worse covering radius for the combined shell than the original results shown in [18]. As references, the table also lists the best known results for single shell with $K$ and $K\times 3$ samples using EEM from CAMINO [25] and the SC formulation from 5 [34]. Note that the precomputed schemes by EEM and GRASSC are only for the single shell case, not for the multi-shell case. If this single shell scheme (with $28$ or $90$ samples) is repeated for 3 shells (i.e., with 3 different b values), then the covering radius for the combined shell (with $84$ samples or $270$ samples) is actually $0^{\circ}$. The scheme with $84$ samples by EEM is not related with the scheme with 28 samples by EEM. IMOC obtains better covering radii than incremental methods, i.e., IGEEM, ISC, and it is even better than MILP, which is because the globally optimal schemes by MILP are determined only from a coarse uniform set with $321$ samples, not from a fine uniform set. It is impractical for MILP to work with $20481$ samples, because MILP is NP hard. For both $28\times 3$ and $90\times 3$ cases, IMOC + 1-Opt + CNLO obtains best schemes with largest covering radii in both 3 individual shells and the combined shell, and the covering radii by IMOC + 1-Opt + CNLO in these 3 shells are even better than best known scheme collected in CAMINO using EEM [11, 25], similarly with the left subfigure in Fig. 2.

Table III also shows the running time for different methods on an ordinary laptop. Incremental methods (IGEEM, ISC, IMOC) are all very fast. IGEEM and ISC have similar running time, because they can be implemented in the same efficient way to maintain the table of increased cost functions between the selected sample set and each unselected sample, and update the table when each new sample is selected. A similar acceleration strategy can be used for 1-Opt. MILP is terminated after 10 minutes, because we found that the cost function of MILP changes very slow after 10 minutes. Compared with the $28\times 3$ scheme, the running time of estimating the $90\times 3$ scheme by all methods except IMOC increases largely. The running time of GEEM increases from $3.7$ seconds to $11$ minites, and the running time of CNLO increases from $39.7$ seconds to $31$ minutes, because of the larger number of constraints. The running time of IMOC decreases for the $90\times 3$ scheme, probably because the iteration number of MOC decreases for the $90\times 3$ scheme, compared with the $28\times 3$ scheme. It should be noted that the running time of scheme design methods is actually not very important for practical usage, because the schemes are optimized only once for scanners.

Fig. 3 shows the schemes by IMOC, IMOC + 1-Opt, and IMOC + 1-Opt + CNLO. It shows that

1. IMOC + 1-Opt can fix the hole in the scheme by IMOC;
2. IMOC + 1-Opt increases the covering radii of individual shells of the scheme by IMOC, although the covering radius of the combined shell keeps the same ($7.78^{\circ}$);
3. CNLO can refine the scheme by IMOC + 1-Opt. Thus for P-C, we do not use IMOC as a stand-alone method. Instead, IMOC + 1-Opt + CNLO is suggested for P-C, where two greedy methods (i.e., IMOC + 1-Opt) are used to set initialization of gradient based optimization method CNLO.

IV-C Rotational Invariance in Reconstruction

We test rotational invariance of the schemes with $28\times 3$ samples by different methods in Table III. We generate synthetic diffusion signal from a mixture of cylinder model [43], where two cylinders cross with $60^{\circ}$ and share the cylinder parameters of length $L=$5\text{,}\mathrm{m}\mathrm{m}, radius $\rho=$5\text{\,}\mathrm{\SIUnitSymbolMicro m}, diffusivity $D_{0}=$2.02\text{\times}{10}^{-3}\text{,}\mathrm{m}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$. With each tested scheme, we rotate the model and generate signals $20481$ times with the rotation angles determined by the uniform sample set with $20481$ samples. Then we perform Non-Negative Spherical Deconvolution (NNSD) [9] with the spherical harmonic (SH) order $8$ to estimate fODFs, detect the peaks, and calculate the mean angular differences by comparing the detected peaks with the ground-truth fiber directions in these $20481$ tests. A good sampling scheme should give low mean angular differences with low standard deviations. Table IV lists the mean and standard deviation of angular differences obtained by NNSD using different schemes. IMOC + 1-Opt +CNLO yields the lowest angular differences with the lowest deviation, and IMOC has the lower mean and deviation than IGEEM.

We compare the $28\times 3$ schemes and $90\times 3$ schemes in DWI signal reconstruction by using tensorial Spherical Polar Fourier Imaging (SPFI) [16]. The DWI signals are synthesized using a mixture of tensor model with two symmetric tensors with eigenvalues $[1.7,0.2,0.2]\times 10^{-3}\,\text{mm}^{2}/\text{s}$ with a crossing angle of $60^{\circ}$. We rotate the model and generate signals $20481$ times with the rotation angles determined by the uniform sample set with $20481$ samples. Tensorial SPFI with spherical order of $8$ and radial order of $4$ is used to fit the data and re-sample the DWI signals at continues 3 shells with a fine tessellation of $20481$ samples. Then normalized mean Root-Mean-Squared-Error (RMSE) of DWI signals in 3 shells is calculated. We also add Rician noise with SNR=$30$ (i.e., standard deviation of $1/30$) for 1000 trails, and re-calculate mean normalized RMSE in reconstruction using noisy samples. Table V shows the mean and standard deviation of normalized RMSE in reconstruction using noise-free samples and noisy samples. It shows that IMOC + 1-Opt + CNLO obtains lower mean and std of normalized RMSE values in reconstruction, especially when the number of samples is small and SNR is high. However, when using a larger number of noisy samples in SPFI (e.g., $90\times 3$ scheme with noise, SNR=30), the normalized RMSE results by IMOC + 1-Opt + CNLO and GEEM are similar. It also shows that the $90\times 3$ scheme by IMOC has much larger standard deviation (i.e., rotational invariance) of normalized RMSE than schemes by other methods, which is because of the hole area effect in IMOC. IMOC + 1-Opt and IMOC + 1-Opt + CNLO reduce the rotational invariance of IMOC.

IV-D MILP for Separation of Sampling Schemes

We evaluate the effectiveness of the proposed MILP in Eq. (11) with $w=1$ to solve P-D-MS, by verifying whether it can separate a set of samples into several subsets, keeping the points in each subset as uniform as possible. We compared MILP with the subsetpoints program in CAMINO [11, 25] and dirsplit in MRtrix 888http://mrtrix.readthedocs.io/en/latest/reference/commands/dirsplit.html for the same task. For this evaluation, we randomly mix two sets of uniform points, where one set has $81$ points generated by spherical tessellation (with covering radius $15.9^{\circ}$) and the other set has $60$ points (with covering radius $18.3^{\circ}$) from CAMINO generated by EEM. Separation of these $141$ points into subsets respectively with $81$ and $60$ samples should ideally give results that match the original uniform point sets. MILP obtained results that exactly match the original point sets within 1 second. subsetpoints in CAMINO, which minimizes electrostatic energy by using simulated annealing, saves the result every hour when it runs. It obtains two incorrect points in each subset after running for 2 hours, 7 incorrect points after 8 hours, and the correct result after 9 hours. Although the correct result has been obtained, the program continues to run for hours until the simulated annealing temperature is finally small enough. dirsplit in MRtrix initializes two sets and then randomly swaps 2 directions in each step to reduce the electrostatic energy. dirsplit can only separate the scheme into 2 schemes with nearly equal numbers 71 and 70. With default parameters and multi-threading, it finished running after 16 minutes. However, the covering radii in these two sets are $6.4^{\circ}$ and $5.1^{\circ}$ respectively, different from the original two uniform sets.

IV-E MILP for sub-sampling of multi-shell HCP scheme

In order to reduce the scanning time and reuse the existing data, we may want to uniformly select a sub-sampled scheme from an existing scheme. IGEEM has been successfully used in multi-shell diffusion data acquisition in HCP [17, 18, 26]. The scheme has 3 shells with b values $1000$, $2000$, and $3000$ $\mathrm{s}\mathrm{/}\mathrm{m}\mathrm{m}^{2}$, 90 samples per shell. See top left subfigure in Fig. 4. MILP and a random selection are performed to obtain two sub-sampled schemes with $45\times 3$ samples. It is a problem of P-D-MM, which to our knowledge can not be solved by any existing released routine (e.g., in CAMINO or MRtrix). MILP yields the covering radius $7.0^{\circ}$ in the combined shell, while the random selection yields $2.4^{\circ}$, which is the same as the original HCP scheme. See first row in Fig. 4. We also perform NNSD with SH order $8$ on a subject with the full scheme of $90\times 3$ samples and two sub-sampled schemes with $45\times 3$ samples from MILP and random selection. The corresponding fODF fields with detected peaks are displayed in Fig. 4 with GFA as backgrounds. Peaks are detected from the fODFs with GFA higher than $0.5$. Fig. 4 shows that the estimated fODF and peaks from sub-sampled scheme by MILP are much closer to the results from the full scheme, compared with the random selection. Note that fODFs near grey matter areas by the scheme from random selection are very noisy, compared with the results using the full scheme. That is because that random selection cannnot keep the uniformity in the sub-sampled scheme. In the region visualized in Fig. 4, the mean angular difference of the detected peaks is $4.2^{\circ}$ between the results from the full scheme and the sub-sampled scheme by MILP, and $5.9^{\circ}$ between the results from the full scheme and the sub-sampled scheme by random selection. The normalized RMSE of the estimated GFA of fODFs is $7.5\%$ for the sub-sampled scheme by MILP and is $18.9\%$ for the scheme by random selection.

V Discussion

V-A Sampling scheme design without antipodal symmetry

The SC formulation and EEM formulation [11, 17, 18, 19, 20, 21, 22, 23] both are to design sampling schemes with the antipodal symmetry constraint, based on prior knowledge that the diffusion signal is antipodal symmetric, i.e., $E(\mathbf{q})=E(-\mathbf{q})$. However, some works proposed to model asymmetry of diffusion signals [44] or fiber ODFs [45], which may prefer sampling schemes without the antipodal symmetry constraint. Moreover general uniform single or multiple shell sampling schemes (i.e., without the antipodal symmetry constraint) are also be useful in some applications other than diffusion MRI.

Although the proposed SC formulation and 5 algorithms focus on sampling scheme design with antipodal symmetry, they can be easily performed to design general sampling schemes without the antipodal symmetry constraint. Actually the original SC formulation in mathematics [33, 35] does not consider the antipodal symmetry constraint. If the absolute value operator is removed from the definition of covering radius in Eq. (3), then Eq. (4) and Eq. (5) are SC formulations respectively for general single and multiple shell scheme design, where Eq. (4) is the original SC formulation in [33, 35]. ISC and 1-Opt can be used for general uniform sampling scheme design straightforwardly, based on the revised definition of covering radius. For IMOC, if we revise the definition of the coverage set of a point x as $C(\textbf{{x}},\theta)=\{\textbf{{y}}\mid\arccos(\textbf{{y}}^{T}\textbf{{x}})<\theta,\textbf{{y}}\in\mathbb{S}^{2}\}$, without the absolute value operator, then IMOC can be used for general uniform scheme design. In Fig. 1 (a), with the revised definition, the covering set does not include the yellow arc generated by the dashed samples in the antipodal positions. For CNLO and MILP, if we remove the absolute operator in Eq. (7), (8), (9), (10) and (11), then CNLO and MILP can be used for general uniform sampling scheme design without the antipodal symmetry constraint. For each proposed algorithm, the strategy and procedure remain the same if the antipodal symmetry constraint is removed.

V-B Sampling scheme design with other considerations

V-B1 Diffusion models and model fitting

As discussed in the introduction part, both the SC formulation and EEM formulation are to design uniform sampling schemes in single or multiple shells with parameters $S$ and $\{K_{s}\}_{s=1}^{S}$. This is based on the fact that diffusion signals in voxels have various preferred orientations, while we need to design one scheme for all signals in human brain. The number of b values $S$, and the number of samples in each shell $\{K_{s}\}_{s=1}^{S}$, are both user determined. Different diffusion models, e.g., tensor model and various HARDI models, may prefer different distributions of b values in the radial part of the $\mathbf{q}$-space, and different number of samples in low or b values. Moreover, even for the same model, different model fitting algorithms may prefer different schemes [12]. The samples with high b values capture more orientational information, while have lower single-to-noise ratio (SNR), than the samples with low b values. Thus robust estimation algorithms [14, 9, 16], which are less sensitive to noise, may prefer larger number of samples with high b values.

V-B2 Eddy current correction

Recent works on eddy current correction [46, 29, 47] showed that three kinds of acquisitions are useful for eddy current correction, i.e., 1) the blip-up blip-down acquisition, which scans the the data twice with flipped phase-encoding directions; 2) the acquisition with gradients on the whole sphere that includes both $\{\mathbf{u}_{i}\}$ and $\{-\mathbf{u}_{i}\}$; 3) the acquisition with gradients on the whole sphere without the antipodal symmetry constraint. Schemes by SC can be directly used in the blip-up blip-down acquisition with flipped phase-encoding directions. and can be used in the second kind of acquisition by incorporating the flipped gradients $\{-\mathbf{u}_{i}\}$. For the third acquisition on the whole sphere without the antipodal symmetry constraint, we can perform the proposed optimization methods without the constraint, as described in Section V-A.

V-B3 Gradient reorientation

The uniform schemes generated by SC or EEM are used to scan DWI data. Some preprocessing steps, including motion correction, eddy current correction, registration of DWIs, etc., may require gradient reorientation [48, 49]. Gradient reorientation changes each individual gradient, and also the separation angles between gradients. Thus even though the designed schemes by SC for scanning data are optimized with large angular separation, the reorientated scheme for model fitting may be sub-optimal. This issue also exists for schemes by EEM, where the electrostatic energy is sub-optimal after reorientation. The gradient reorientation issue has not been considered in existing works on uniform scheme design [11, 17, 18, 19, 20, 21, 22, 23], because the transform matrix for each gradient vector is dependent on DWI images and also the image registration algorithm used in the preprocessing steps. In this paper, we show that in both SC and EEM formulations, uniform sampling scheme design before gradient reorientation is an approximation of uniform sampling scheme design after gradient reorientation, in an approximation algorithm point of view 999https://en.wikipedia.org/wiki/Approximation_algorithm.

Intuitively, if the transform matrix is close to the identity matrix, then the reoriented scheme are closed to the original scheme which has been well optimized. This intuition can be quantitatively described in Proposition 1, which shows that if the transform matrix is inside a ball of the identity matrix with radius $\delta$, then the inner product of two reoriented gradient vectors $\textbf{{x}}_{1}^{T}\textbf{{x}}_{2}$ is inside a range centered at the inner product of the original two gradient vectors $\mathbf{u}_{1}^{T}\mathbf{u}_{2}$ with a radius of $2\delta+\delta^{2}$. The proof of this proposition is shown in Appendix A.

Proposition 1

Angular separation after gradient reorientation. Let $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ be two gradient vectors in $\mathbb{S}^{2}$. Denote the two gradient vectors after reorientation as $\textbf{{x}}_{1}$ and $\textbf{{x}}_{2}$, and the corresponding transform matrices as $\mathbf{T}_{1}$ and $\mathbf{T_{2}}$ such that $\textbf{{x}}_{i}=\mathbf{T}_{i}\mathbf{u}_{i}$, $i=1,2$. If $\|\mathbf{T}_{i}-\mathbf{I}\|_{2}\leq\delta$, where $\mathbf{I}$ is the identity matrix, then the reoriented gradient vectors satisfies $|\textbf{{x}}_{1}^{T}\textbf{{x}}_{2}-\mathbf{u}_{1}^{T}\mathbf{u}_{2}|\leq 2\delta+\delta^{2}$.

Note that Proposition 1 works for gradient reorientation using rotation matrices and also general affine transforms [48], where if $\mathbf{A}_{i}$ is the affine matrix used to reorient $\mathbf{u}_{i}$ such that $\textbf{{x}}_{i}=\frac{\mathbf{A}_{i}\mathbf{u}_{i}}{\|\mathbf{A}_{i}\mathbf{u}_{i}\|_{2}}$, then $\mathbf{T}_{i}=\frac{\mathbf{A}_{i}}{\|\mathbf{A}_{i}\mathbf{u}_{i}\|_{2}}$. Proposition 1 means that the separation angle of two reoriented gradients $\arccos(\textbf{{x}}_{1}^{T}\textbf{{x}}_{2})$ is in $[\arccos(\mathbf{u}_{1}^{T}\mathbf{u}_{2}+2\delta+\delta^{2}),\arccos(\mathbf{u}_{1}^{T}\mathbf{u}_{2}-2\delta-\delta^{2})]$. And the minimal separation angle of any two reoriented gradients is also in a range defined by the separation angle of two original gradient vectors. Thus, in an approximation algorithm point of view 9, maximizing the angular separation of the original scheme for scanner is to approximately maximize the angular separation of the reorientated scheme for model fitting. In practice, the transform matrix between a reoriented gradient vector and an original gradient vector in motion correction or eddy current correction is close to the identity matrix, i.e., $\|\mathbf{T}-\mathbf{I}\|_{2}<\delta$, and $\delta$ is a small value. If the transform matrix $\mathbf{T}$ is the rotation about $z$-axis with the degree of $5^{\circ}$, then $\epsilon=0.123$ that is about $12\%$ of the norm of $\mathbf{I}$. In practice, the rotation angle of the transform matrix is smaller than $5^{\circ}$. Thus the change of the separation angle after gradient reorientation is small.

For the EEM formulation, as shown in Eq. (1), the electrostatic energy between $\mathbf{u}_{i}$ and $\mathbf{u}_{j}$ is

[TABLE]

Based on Proposition 1, after gradient reorientation, the electrostatic energy between $\textbf{{x}}_{i}$ and $\textbf{{x}}_{j}$ is

[TABLE]

Note that here we assume $|\mathbf{u}_{i}^{T}\mathbf{u}_{j}|+2\delta+\delta^{2}<1$. The length of the range is

[TABLE]

where $c_{0}=\max_{i<j}|\mathbf{u}_{i}^{T}\mathbf{u}_{j}|<1-2\delta-\delta^{2}$. Thus, the electrostatic energy of the reoriented sampling scheme $\sum_{i<j}f_{2}(\mathbf{u}_{i},\mathbf{u}_{j})$ , and the electrostatic energy of the original scheme $\sum_{i<j}f_{2}(\textbf{{x}}_{i},\textbf{{x}}_{j})$, are both inside a range, and we have

[TABLE]

where $K$ is the number of samples in the single shell scheme. Thus, minimizing the electrostatic energy of the original scheme is to approximately minimize the energy of the reoriented scheme. When $\delta$ is small, i.e., the transform matrix is close to the identity matrix, the range length in Eq. (16) is small and the electrostatic energies of the original scheme and the reoriented scheme are close.

We have shown that uniform sampling scheme design before gradient reorientation is an approximation of uniform sampling scheme design after gradient reorientation. However, it should be noted that if the transform matrix due to gradient reorientation is far from the identity matrix, the approximation error can be large. As shown in experiments, the proposed SC method yields larger angular separation (i.e., covering radius) between samples than the state-of-the-art EEM and GEEM method. If we assume that gradient reorientation reduces the covering radius of all schemes by the same amount ($<5^{\circ}$), then, after gradient reorientation, the scheme by SC still has larger covering radius than the scheme by EEM.

V-B4 Minimizing gradient heating

The SC formulation in this paper is to design uniform sampling scheme in shells, similarly with the works based on EEM [11, 17, 18, 19, 20, 21, 22, 23], which have not yet considered gradient heating specifically. We have proposed to use ISC for P-O, where the first partial set of samples are nearly uniformly separated. It was reported in [40] that optimal ordering not only ensures that usable data are available if the scan needs to be cut short, but also reduces the risk of systematic errors due to gradient heating and subject motion. Minimizing gradient heating is beyond of the content of this paper, and gradient heating may be reduced in different ways [50, 51]. It can be a future work to investigate which way of reducing gradient heating is more appropriate to be incorporated into current SC framework in this paper.

VI Conclusion

In this paper, we proposed to use the Spherical Code concept, which directly maximizes angular separations between samples, to design single- and multi-shell uniform sampling schemes. We showed that various continuous and discrete sampling scheme design problems can be formulated using SC. We proposed five novel algorithms (i.e., ISC, IMOC, 1-Opt, MILP, CNLO) based on SC for uniform sampling scheme design. See Table I. ISC is mainly to order samples in an existing scheme. IMOC and 1-Opt are to quickly obtain schemes in deterministic greedy ways. Although IMOC is a greedy incremental method, it obtains globally optimal solution in $\mathbb{S}^{1}$. For the sphere $\mathbb{S}^{2}$ in dMRI, covering radii of schemes by IMOC are very close to the best known solutions in single shell case 5, and are much larger than EEM in CAMINO. CNLO is a local optimization method which refines a given initialization (normally from IMOC+1-Opt). MILP is used to perform sub-sampling in an existing scheme. It can obtain globally optimal solution or a good result within 10 minutes. The experiments showed that compared with the state-of-the-art methods [18, 20] based on electrostatic energy which is used in HCP, the multi-shell scheme by IMOC + 1-Opt + CNLO yields larger covering radii in general, and better rotational invariance in dMRI reconstruction when the SNR is high and the number of samples is small. However, when a large number of noisy samples are used in DWI reconstruction, the sampling schemes by different methods may have similar reconstruction results. The proposed SC methods are to maximize the angular separation between samples, independent of diffusion models and diffusion reconstruction. The optimal sampling scheme that yields the best reconstruction quality is dependent on the diffusion model, reconstruction algorithm, and the ground truth diffusion signal. It is a future work to optimize the sampling schemes for a given specific diffusion model and a given reconstruction algorithm, so that the final result of the diffusion reconstruction is optimized. We have released some demos and codes in DMRITool 1.

Appendix A Proof of Proposition 1

Denote $\mathbf{T}_{i}=\mathbf{I}+\mathbf{C}_{i}$, then we have $\|\mathbf{C}_{i}\|_{2}\leq\delta$, and

[TABLE]

Let $c=|\mathbf{u}_{1}^{T}(\mathbf{C}_{1}^{T}+\mathbf{C}_{2})\mathbf{u}_{2}|+|\mathbf{u}_{1}^{T}\mathbf{C}_{1}^{T}\mathbf{C}_{2}\mathbf{u}_{2}|$, then

[TABLE]

Because $|\mathbf{u}_{1}^{T}(\mathbf{C}_{1}^{T}+\mathbf{C}_{2})\mathbf{u}_{2}|\leq\|\mathbf{C}_{1}^{T}+\mathbf{C}_{2}\|_{2}\leq 2\delta$, and $|\mathbf{u}_{1}^{T}\mathbf{C}_{1}^{T}\mathbf{C}_{2}\mathbf{u}_{2}|\leq\|\mathbf{C}_{1}^{T}\|_{2}\|\mathbf{C}_{2}\|_{2}=\delta^{2}$, we have $c\leq 2\delta+\delta^{2}$. Then

[TABLE]

Acknowledgments

The authors would like to thank Dr. Ruiliang Bai and Dr. Alexandru Avram for useful discussions. This work was supported in part by funds provided by the Intramural Research Program of the Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD) (ZIA-HD000266), and NIH grants (NS093842 and EB022880). The data were provided in part by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.