Motion correction, eddy current correction and susceptibility distortion correction were applied to each volume of the DW images before averaging over the four repetitions using the toolbox in FSL (FMRIB Software Library, University of Oxford, http://fsl.fmrib.ox.ac.uk/fsl/). Motion and eddy current distortions were corrected using a linear registration to the b0 image (eddy_correct, FSL) for each volume within each repetition. Susceptibility distortion was corrected by calculating

Fig. 1 The top panel shows simulated ODFs of a single fiber orientation (0°, 30°, 60°, 90°). The ODFs were simulated with a single tensor with axial diffusivity of 1200 mm^{2}/s and radial diffusivity of 250 mm^{2}/s at b = 2500 s/mm^{2}. The middle panel shows magnitude components of SH coefficients. The bottom panel shows phase components of SH coefficients. l and m are orders of Legendre function in the spherical harmonic bases. For illustrating purpose, SH coefficients of l = 0 to 4 and m = — l to l are shown here

the geometric distortion and signal loss from the field map and was compensated in the DW images (fudge, FSL). A final motion correction was applied to all four repetitions by rigidly registering the b0 images from each repetition before averaging.

In the diffusion space, the averaged DW images were then used to calculate the structural ODF profiles of the amygdala using in-house MATLAB programs [22] with a Q-ball Imaging (QBI) algorithm [14]. Spherical Harmonics (SH) coefficients of each ODF profile were extracted up to an order of 6, i.e., l_{max} = 6. As a symmetric ODF was assumed and odd orders contain only noise information, only coefficients of even orders were kept [19,20]. A total of 28 SH coefficient pairs (i.e., magnitude and phase) that represent the shape and orientation of ODF in each voxel were entered into the subsequent spectral clustering. Figure 1 shows simulated ODFs of a single fiber orientation (0°, 30°, 60°, 90°) and their SH coefficients. Figure 2 shows simulated ODFs for crossing fibers with a rotating 2nd fiber (0°, 30°, 60°, 90°) and their SH coefficients. Consistent with observations described in [19], the shape (e.g., number of crossing fibers) and the orientation (e.g., rotation angle) of the ODFs are described by the combination of magnitude and phase components of the SH coefficients.

Fig. 2 The top panel shows simulated ODFs for crossing fibers with a rotating 2nd fiber (0°, 30°, 60°, 90°). The ODFs were simulated with two tensors at b = 2500 s/mm^{2}. Each tensor has axial diffusivity of 1200 mm^{2}/s and radial diffusivity of 250 mm^{2}/s. The middle panel shows the magnitude components of the SH coefficients. The bottom panel shows the phase components of the SH coefficients. l and m are orders of Legendre function in the spherical harmonic bases; and the SH coefficients of l = 0 to 4 and m = — l to l are shown here