This study aims to parcellate brain white matter (WM) into supervoxels with homogeneous diffusion property measured by Bregman divergence between orientation distribution functions (ODF) derived from q-ball imaging (QBI). Current WM parcellation method on large affinity matrix needs enormous computation resources¹, or voxels in result supervoxel are disconnected². In this study, we propose a fast and robust flow-based WM supervoxel parcellation method. Parcellating the white matter into supervoxels with homogeneous diffusion properties can benefit functional magnetic resonance imaging (fMRI) analysis³, WM connectivity network construction⁴, tractography, and region-of-interest based analysis⁵.

The overview of the proposed method can be seen in Fig. 2. Our method was performed on WM graph, based on a flow-based supervoxel parcellation algorithm as described below.

(1) WM graph construction: We used FreeSurfer⁶ to segment T1 image, and selected WM tissues to construct a WM binary image. Voxels in WM binary image were then transferred to a node in WM graph. An edge was generated to connect each pair of adjacent voxels.

(2) Spherical harmonics (SH) expansion: We reconstructed ODF from diffusion weighted imaging (DWI) according to Descoteaux’s method⁷ using DSI-studio⁸. ODF is a spherical function and can be represented as a linear combination of SH. Because ODF is real and centro-symmetric, we performed an up-to 4th even order real SH expansion:

$${\bf D} = \sum_{{\it l}=0,2,4} \sum_{m=-{\it l}}^{{\it l}}f_{lm}Y_{lm}$$

where $$$Y_{lm}$$$ is m-order l-degree real Laplace spherical harmonic, and $$$f_{lm}$$$ is corresponding expansion coefficient. We also performed SH expansion to logarithm of ODF values. Then, two $$$15\times1$$$ vectors of SH coefficients were attached to each node of WM graph.

(3) Lattice seeding: We performed lattice seeding on whole WM image space to uniformly split image into grids. Since given lattice resolution is not always integral multiple of voxel size, we used an accumulator to sum residual grid resolution. As compensation, grid was enlarged when residual grid resolution was not less than a voxel size. For expected supervoxel number, a gradient descent search was used to find best grid resolution. We selected WM voxels nearest to grid centers as seeds to initialize supervoxel segmentation.

(4) Flow assignment: Starting from seeds, a flow assigned values to all nodes along edges. If flow from seed $$$s$$$ meet node $$$n$$$, the value to be assigned was defined as integral of ODF differences along trajectory of the flow:

$$v_n=\int_{\Gamma} d{\bf D}$$

If the to-be-assigned value was less than predefined value, the flow assigned new value to $$$n$$$, and labeled it as $$$s$$$.

We used functional Bregman divergence⁹ to define difference between ODFs. For ODF $$$g$$$, let

$$\phi\left[g\right]=\int_{S}g\ln(g)d\Omega$$

According to definition of functional Bregman divergence $$$d_{\phi}\left[f,g\right]=\phi\left[f\right]-\phi\left[g\right]-\delta\phi\left[g;f-g\right]$$$, relative entropy (also Kullback-Leibler divergence) between ODF $$$f$$$ and $$$g$$$ is then defined as:

$$d_{KL}\left[f,g\right]=\int_{S}{f\ln\frac{f}{g}-\left(f-g\right)d\Omega}$$

To obtain a symmetric non-negative measure, we used total average of $$$d_{KL}$$$ (also Jensen-Shannon divergence):

$$d_{JS}\left(f,g\right)=\frac{1}{2}\left(d_{KL}\left(f,g\right)+d_{KL}\left(g,f\right)\right)=\frac{1}{2}\int_{S}f\ln{f}-f\ln{g}+g\ln{g}-g\ln{f}{d\Omega}$$

Let $$$C^f$$$, $$$C^g$$$ denote SH coefficients of ODF $$$f$$$ and $$$g$$$, and $$$C^{\ln{f}}$$$, $$$C^{\ln{g}}$$$ denote SH coefficients of logarithms $$$ln{f}$$$ and $$$ln{g}$$$. Due to orthogonality of SH, $$$d_{JS}$$$ can be represented as:

$$d_{JS}\left(f,g\right)=\frac{1}{2}\sum_{{\it l}=0,2,4}\sum_{m=-{\it l}}^{\it l} C_{l,m}^{f}C_{l,m}^{\ln{f}}-C_{l,m}^{f}C_{l,m}^{\ln{g}}+C_{l,m}^{g}C_{l,m}^{\ln{g}}-C_{l,m}^{g}C_{l,m}^{\ln{f}}$$

(5) Reflow: Assignment flows could intersect each other, and hence some nodes were departed from corresponding seeds. We called a reflow process to fix this issue. Firstly, a flow searched connected nodes from same seed, and labeled them as “connected”. For “unconnected” nodes, we erased values and labels, and pushed their “connected” neighbors to call another flow assignment. The reflow process was called several times until all nodes were labeled as “connected”, to guarantee integrity of supervoxels.

(6) Update seeds: A gradient descent method was used to found new seeds (supervoxel centers) from original seeds. A new seed was a node within supervoxel that minimizes summed values. Then, a new flow assignment was called from new seeds. We terminated the assignment-update loop when seeds remained unchanged, or oscillated.

We applied the proposed method to high-quality data in Human Connetome Project (HCP) (http://www.humanconnectomeproject.org). We parcellated WM into 5000 parcels. Fig. 2 illustrates the connectivity network with superovoxels as nodes. It can be seen that local fiber orientations within a supervoxel are homogeneous, as shown in Fig. 3. Fig. 4 illustrates the usage of WM supervoxels to extract corticospinal tracts, cingulum of cingulate gyrus, forceps major, and forceps minor.In this work, we propose a WM supervoxel parcellation method in QBI. This method is a feasible and efficient approach to generate WM parcels with homogenous diffusion properties. It could benefit future relevant neuroscience researches, like WM connectivity network construction, ROI-based analysis, tractography, etc.