Angular Upsampling in Infant Diffusion MRI Using Neighborhood Matching in x-q Space

Diffusion MRI requires sufficient coverage of the diffusion wavevector space, also known as the q-space, to adequately capture the pattern of water diffusion in various directions and scales. As a result, the acquisition time can be prohibitive for individuals who are unable to stay still in the scanner for an extensive period of time, such as infants. To address this problem, in this paper we harness non-local self-similar information in the x-q space of diffusion MRI data for q-space upsampling. Specifically, we first perform neighborhood matching to establish the relationships of signals in x-q space. The signal relationships are then used to regularize an ill-posed inverse problem related to the estimation of high angular resolution diffusion MRI data from its low-resolution counterpart. Our framework allows information from curved white matter structures to be used for effective regularization of the otherwise ill-posed problem. Extensive evaluations using synthetic and infant diffusion MRI data demonstrate the effectiveness of our method. Compared with the widely adopted interpolation methods using spherical radial basis functions and spherical harmonics, our method is able to produce high angular resolution diffusion MRI data with greater quality, both qualitatively and quantitatively.

Diffusion MRI requires sufficient coverage of the diffusion wavevector space, also known as the q-space, to adequately capture the pattern of water diffusion in various directions and scales. As a result, the acquisition time can be prohibitive for individuals who are unable to stay still in the scanner for an extensive period of time, such as infants.
To address this problem, in this paper we harness non-local self-similar information in the x-q space of diffusion MRI data for q-space upsampling. Specifically, we first perform neighborhood matching to establish the relationships of signals in x-q space. The signal relationships are then used to regularize an ill-posed inverse problem related to the estimation of high angular resolution diffusion MRI data from its low-resolution counterpart. Our framework allows information from curved white matter structures to be used for effective regularization of the otherwise ill-posed problem. Extensive evaluations using synthetic and infant diffusion MRI data demonstrate the effectiveness of our method. Compared with the widely adopted interpolation methods using spherical radial basis functions and spherical harmonics, our method is able to produce high angular resolution diffusion MRI data with greater quality, both qualitatively and quantitatively.

INTRODUCTION
Infant brain development involves complex cerebral growth and maturation with the white matter (WM) undergoing rapid myelination and synaptogenesis (Qiu et al., 2015). Diffusion MRI (DMRI) has been widely employed to study this developmental process in vivo (Yap et al., 2011;Huang et al., 2013;Dubois et al., 2014;Qiu et al., 2015). For instance, using diffusion tensor imaging (DTI), researchers have observed an increase in the fractional anisotropy (FA) during the first few years of life (Dubois et al., 2014;Qiu et al., 2015), implying more restriction on water movement owing to the ensheathment of oligodendrocytes around the axons. Mean diffusivity (Dubois et al., 2014;Qiu et al., 2015) and structural connectivity (Yap et al., 2011;Huang et al., 2013) have also been used to study early brain development.
Existing brain development studies mainly rely on DTI. However, DTI utilizes a simple but insufficient tensor model to characterize local fiber configurations. More sophisticated diffusion models are needed to better characterize realistic configurations, such as crossing, fanning, and bending, and to tease out inter-and extra-cellular compartments (Ning et al., 2015;Yap et al., 2016b;Ye et al., 2016). However, unlike the tensor model that only requires 6 diffusion-weighted (DW) images and one non-DW image, these advanced models requires high angular resolution (HAR) DMRI data. The angular resolution of DMRI data is determined by the number of gradients used in data acquisition. Each gradient corresponds to a point in q-space. A larger number of DW images allows the use of advanced diffusion models but prolong the acquisition time, which is prohibitive in clinical settings.
In practice, the window of opportunity for imaging infants is short. To put this in perspective, in the Human Connectome Project (HCP) (Van Essen et al., 2012) each individual was allotted a DMRI scan time of about an hour. However, in the Baby Connectome Project (BCP) (Fallik, 2016;Cao et al., 2017;Howell et al., 2018), the tolerable scan time is well below 15 min. Infants are typically scanned without sedation while they are asleep. The scanning may also need to be terminated prematurely if the infant is awakened by the loud acoustic noise and sudden vibrations caused by the rapid switching of gradient amplitude and polarity (McJury and Shellock, 2000;Hutter et al., 2017). The short acquisition time precludes a denser coverage of q-space, limiting studies to simpler models such as the diffusion tensor.
To increase the angular resolution without inducing additional acquisition time, post-acquisition angular upsampling methods have been proposed. Tuch (2004) proposed to interpolate q-space using spherical radial basis functions (SRBFs). This method interpolates one point in q-space by weighted averaging of angularly neighboring q-space measurements. An alternative method is interpolation using spherical harmonics (SHs) (Descoteaux et al., 2007). For each voxel, this method first decomposes the DMRI measurements into an SH coefficient vectors. DMRI measurements are then reconstructed with the coefficients and the corresponding SH basis. Despite the promising performance of these methods, one major limitation is that only q-space information is considered during interpolation, and valuable x-space information is overlooked.
To overcome this limitation, in this paper, we propose to harness joint x-q space information for angular upsampling of DMRI data. For this purpose, we first establish signal correspondences in x-q space using a robust neighborhood matching technique described in Chen et al. (2016aChen et al. ( , 2017a. We then recover the HAR DMRI data in a regularization framework based on the signal correspondences. Angular upsampling is achieved by using the signal relationships in the joint x-q space. Extensive experiments on synthetic and infant DMRI data demonstrate that our method is able to recover HAR DMRI data with a remarkably improved quality. Part of this work has been reported in our recent conference paper (Chen et al., 2017b). Herein, we present (1) New insights into the application of our method to infant DMRI data; (2) More detailed descriptions of the proposed method; (3) More complete mathematical derivation details; (4) Comparison with the stateof-the-art method, SH interpolation; (5) New experimental results on synthetic data and infant data acquired at different time points. None of these is part of the conference publication.
The rest of the paper is organized as follows. In Section 2, we give a detailed description of the proposed method. In Section 3, we demonstrate the effectiveness of the method with both synthetic and real infant DMRI data. In Section 4, we provide further discussion on this work. Finally, in Section 5, we conclude this work.

METHODS
Each signal in the x-q space is associated with a voxel location in the x-space and a set of diffusion gradient direction and strength in the q-space. The signals are typically collected via a set of DW images, each corresponding to a point in q-space. In this section, we will first detail in how to establish the signal relationships in x-q space and then clarify how to utilize these relationships to regularize the ill-posed inverse problem associated with angular upsampling of DMRI data.

Signal Correspondences in x-q Space
We utilize a neighborhood matching technique to determine signal correspondences in x-q space. Neighborhood matching techniques have a wide range of applications in medical image analysis, including statistical group comparison (Chen et al., 2015), atlas building (Yang et al., 2017;Saghafi et al., 2017;Kim et al., 2017), fiber orientation estimation (Chen et al., 2016b), and denoising (Chen et al., 2016a,c). Instead of using the conventional methods designed mainly for x-space (Buades et al., 2005;Chen et al., 2016c), we employ the method described in Chen et al. (2016aChen et al. ( , 2017a for accurately establishing signal correspondences in x-q space. The signal relationships in x-q space is quantified using similarity weights (Chen et al., 2016a(Chen et al., , 2017a. A large weight is assigned to two matching signals, while the weight of two mismatching signals is low. To compute the similarity weight, robust rotation-invariant features are first computed for each node in x-q space. As in Chen et al. (2016aChen et al. ( , 2017a, we represent the q-space using a graph and employ a technique, named graph framelet transforms (GFTs), for the feature computation. The resulting features are then employed to compute the similarity weights, essentially establishing the signal correspondence in x-q space. Figure 1 illustrates an overview of this method. We flesh out these steps in the following.

Graph Representation for q-space
As illustrated in Figure 1, we use a graph G to represent the q-space. For this purpose, we compute the adjacency weight between each two nodes in q-space using two kernels for diffusion gradient directions and strengths. In this way, the geometric relationships in q-space are encoded in the graph G where a large edge weight indicates two nodes sharing similar gradient directions and diffusion weightings. Since q-space is now represented as a graph, we can view the signals at each voxel as a function f defined on the graph.

Feature Computation
We then compute the multi-scale features using GFTs. The key idea of GFT is to slice the frequency spectrum of f in a multiscale fashion by using a set of masks. This involves how to perform convolutions on a graph. As clarified in Dong (2017), the eigenvectors and eigenvalues of graph Laplacian are the Fourier spectrum and basis of one graph. Therefore, the convolution can be defined in the transform domain by slicing the Fourier spectrum using a mask. In this way, the graph framelet analysis transform W up to level L and mask R can be defined as where B L,R := {(1, 1), (1, 2), . . . , (1, R), (2, 1), . . . , (L, R)} ∪ {(L, 0)} and W l,r is a matrix representing the GFT operator associated with l and r. We then have the feature vector for graph node k as φ[k] : = {α l,r [k] : (l, r) ∈ B L,R }. GFT has three major advantages: (1) Due to the multi-scale nature, it provides rich characterizations of the signal information; (2) As proven in Yap et al. (2016a); Dong (2017), the GFT feature is rotation-invariant, allowing the neighborhood matching in the curved domain. (3) As shown in Equation (1), GFT is linear and thus computationally efficient.

Neighborhood Matching in x-q Space
The resulting GFT features are then employed for neighborhood matching. We define the similarity weight as a Gaussian function of GFT feature distance, i.e., (2) For the x-q space neighborhood matching in DMRI data, we extend Equation (2) with the consideration of spatial location and diffusion-weighting variations. The new form of Equation (2) is where i and j are the indices of two spatial locations; b k is the b-value for the gradient with index k. Through neighborhood matching, similarity weights are assign to paired nodes in x-q space, ultimately establishing dense node-to-node correspondences.

Angular Upsampling
We use the x-q space data relationships determined in the previous section to guide data upsampling in the q-space. In general, recovering the high angular resolution (HAR) data from the low angular resolution (LAR) data is an ill-posed inverse problem. Solution can however be feasible by imposing structure via harnessing prior information to reduce the dimensionality of the problem. Our approach uses the signal correlation in the x-q space to help reduce the complexity of the problem by imposing that the reconstructed signal should be smooth in a non-local sense. That is, neighboring points in the product space of the signal space and the x-q space should be reconstructed using similar values. Unlike the commonly used x-space regularization (Coupé et al., 2013), which takes into account spatial correlation, x-q space regularization allows correlation across DW images collected using different gradient directions and strengths to be considered. This is fitting when considering the fact that WM structures might be curved and hence causing rapid changes within a DW image. We formulate the angular upsampling problem in a regularization framework and define our objective function as where x and y are two vectors representing the HAR and LAR data, respectively; A is a q-space downsampling operator; R i,k is an operator that extracts the diffusion signal associated with index (i, k). The penalty function consists of a data fidelity term and a regularization term based on x-q space neighborhood matching. The data fidelity term well preserves the information of original LAR data, while the regularization term encourages that each signal to be represented by its matching signals, essentially establishing the non-local smoothness (Protter et al., 2009).

Optimization
To minimize (4), we compute the derivative and equate it to zero: Based on the facts that the neighborhood is symmetric (i.e., if (j, l) ∈ V(i, k), then (i, k) ∈ V(j, l)) and the weights are symmetric (i.e., w[i, k; j, l] = w[j, l; i, k]) (Protter et al., 2009), we have Equation (5) can be simplified using (6), giving Equation (7) can be solved directly but involves the inversion of a very large matrix, therefore we choose instead to use fixedpoint iteration to solve the problem, as suggested in Protter et al. (2009). If we let x n be the solution at iteration n, the following can be proven to be convergent (Protter et al., 2009): The solution x can be obtained iteratively using is a diagonal matrix, therefore the matrix inversion in (9) can be done effectively.

Initialization
The data are transformed so that the noise is Gaussian distributed as described in Koay et al. (2009). The algorithm is then initialized using an upsampled version of y, which is obtained via interpolation using SHs.

Neighborhood Matching
Neighborhood matching is performed based on the upsampled version of y. The resulting weights remain unchanged until a solution x is obtained. In principle, we can use x to re-estimate the weights and rerun the algorithm to obtain a refined solution. However, our experimental results indicate that the benefit of doing so is minimal. Therefore, we will only show results without weight re-estimation.

Stopping Criterion
We stop the algorithm when the mean absolute difference (MAD) between the outcomes of two iterations, i.e., x n−1 and x n , is less than a constant tol. We define tol = βσ G , where σ G is the standard derivation of the Gaussian noise and β is a constant.

EXPERIMENTS
The proposed angular upsampling method was evaluated using both synthetic and real data. Through grid search, we found that λ = 100 and β = 10 −3 give the best results. We compared the proposed method with two baselines, SRBF interpolation and SH interpolation.

Synthetic Data
For quantitative evaluation, we generated a set of synthetic data using phantomαs (Caruyer et al., 2014) and the fiber geometric model of ISBI 2013 HARDI challenge 1 . We utilized 321 gradient directions, uniformly distributed on the surface of a sphere, to simulate the HAR DMRI data. The number of gradient directions was reduced to 81 for the simulation of corresponding LAR counterpart. In the data simulation, we used three b-values, including 1, 000, 2, 000, 3, 000 s/mm 2 . After obtaining the noisefree LAR data, we added four levels (SNR = 15,20,25,30) of 32-channel noncentral chi (nc-χ) noise to the data to simulate noise disturbances. The noise-free HAR DMRI data was used as the ground truth for quantitative evaluations.

Real Data
DMRI data were acquired for three infants at three different time points: 0 month, 6 months, and 12 months. All enrolled subjects had written informed consent provided by parent/guardian. The experimental protocols were approved by the Institutional Review Board of the University of North Carolina (UNC) School of Medicine. The study was carried out in accordance with the approved guidelines. All the data are acquired using a Siemens 3T Magnetom Prisma MR scanner and a standard imaging protocol: 140 × 140 imaging matrix, 1.5 × 1.5 × 1.5 mm 3 resolution, TE = 88 ms, TR = 2,365 ms, 32-channel receiver coil, b = 700, 1500, 3000 s/mm 2 , and 144 non-collinear gradient directions. We uniformly selected 72 gradient directions to generate the LAR data for evaluation.

Evaluation Methods
Quantitative and qualitative evaluations were performed as described in the following: 1. RMSE maps: We computed voxel-wise RMSE value between two sets of DMRI datasets to measure their similarity locally. 2. FA images: We computed the FA images using the iterative weighted tensor fitting method presented in Salvador et al. (2005). 3. Absolute difference (AD) maps: We compute the AD map between one FA image and the ground truth FA image to evaluate the performance of the algorithm locally. 4. Mean normalized absolute difference (MNAD): We computed MNAD value by 1) computing AD map, 2) normalizing the AD map using the ground truth FA image voxel-wisely, and 3) computing the mean value of the normalized map within the brain region. 5. Peak signal-to-noise ratio (PSNR): PSNR is used for quantitative evaluation of FA images, and is defined as where MAX is the maximum FA value, which is 1 in our case. 6. Fiber ODFs: We compute fiber orientation distribution functions (ODFs) using the method presented in Yap et al. (2016b), which caters to multiple tissue types using multi-shell data.

Quantitative Comparison -Synthetic Data
Using the FA image of noise-free HAR data as ground truth, we evaluated the quality of the upsampled data using MNAD and PSNR. The results, shown in Figures 2 and 3, indicate that the proposed method outperforms SRBF interpolation and SH interpolation for all noise levels. The largest improvement over the second best method, SH interpolation, is 0.017 in term of MNAD when SNR = 25. The corresponding PSNR improvement is 5.74 dB.

DW Images -Synthetic Data
The full views and close-up views of DW images, shown in the top two rows of Figure 4, indicate that the proposed method results in better structural contrast. The close-up views of RMSE maps, shown in the bottom row of Figure 4, indicate that our method gives lower RMSE than SRBF interpolation and SH interpolation, which demonstrates that the upsampled data given by our method is closer to ground truth.

FA Images -Synthetic Data
The top row of Figure 5 shows the FA images given by ground truth data and upsampled data. We use warm colors to represent large FA values for a better visualization. Compared with SRBF interpolation and SH interpolation, our method produces an FA image closer to the ground truth. This  observation is further confirmed by the AD maps of the FA images, shown in the bottom row of Figure 5. Our method yields the lowest MAD value, indicating the best performance.

Fiber ODFs -Synthetic Data
Accurate ODF estimation relies on sufficient angular samples. The ODFs, shown in Figure 6, indicate that our method gives clean and coherent ODFs that are close to the ground truth. In contrast, spurious peaks are introduced by SRBF interpolation and SH interpolation.

Diffusion Signal Profiles -Synthetic Data
For a more direct visualization of the upsampled data, we rendered the signal values on a sphere. The results, shown in Figure 7, indicate that our method gives values that are close to the ground truth.

MNAD Comparison -Real Data
We also computed MNAD values for the quantitative evaluation of real data experimental results. Figure 8 shows the MNAD between the FA images given by the upsampled data and the original HAR data. For all time points, our method outperforms SRBF interpolation and SH interpolation, with a largest MNAD     reduction of 0.026 over the second best method at 12-months time point.

DW Images -Real Data
The observations from Figures 9 and 10 for real data are consistent with that in Figure 4. The DW image given by our method shows more subtle structural details and is closer to the one in the original HAR data.

FA Images -Real Data
Figures 11 and 12 further confirm our observation in Figure 5. Our method produces low AD values and gives an FA FIGURE 10 | Close-Up Views of DW Images -Real Data. Regional close-up views of the DW images of a 12-month infant subject.
FIGURE 11 | FA Images -Real Data. Comparison of FA images using infant data. The color of FA images represents the FA value, e.g., warmer color means a larger FA value.
Frontiers in Neuroinformatics | www.frontiersin.org 9 September 2018 | Volume 12 | Article 57 image that is close to that given by the original HAR data. Figure 13 indicates that our method gives clean and coherent ODFs that are very similar to those given by the original HAR data. In contrast, SRBF interpolation and SH interpolation produce ODFs with a large number of spurious peaks.

DISCUSSIONS
Our method is effective because it preserves the sharpness of signal profiles in q-space during upsampling. Utilizing non-local smoothness as prior, it takes into account signal similarity in xq space and avoids the pitfall of averaging over disparate signals. For the curved white matter structures, considering only signal correlation in x-space (i.e., a fixed point in q-space) is problematic because the signal changes rapidly across space. On the other hand, considering only signal correlation in q-space (i.e., a fixed point in x-space) causes smoothing of anisotropic signal profiles due to sharp changes across q-space measurements. Our method harnesses the fact that the signal is smooth in the joint x-q space, even for highly curved structures.
There are some recent works on using compressive sensing (CS) (Baraniuk, 2007) to recover high-resolution (HR) DMRI data from under-sampled k-q space data. For instance, Mani et al. (2015) performed HR DMRI data reconstruction by imposing sparsity on the coefficients of ODFs and by reducing the total variation (TV) of the HR DMRI data. Cheng et al. (2015) proposed a method, called 6D-CS-DMRI, to recover the ensemble average propagator (EAP) and HR DMRI data simultaneously in a CS framework. The associated ill-posed inverse problem was regularized by the sparsity of the coefficients of EAP, the DW image smoothness enforced by TV, and the sparsity of the wavelet coefficients of DW images. Despite achieving promising performance in recovering HR DMRI data, these methods need dedicated imaging protocols, restricting their widespread application. In contrast, our method enhances the angular resolution post-acquisition and thus avoids special imaging protocols.
In practice, our method is limited by its large memory requirement. The memory issue is mainly caused by the need to store the matching weights used in the fixed-point iteration algorithm. To put this in perspective, if we have n x-q space points, the number of matching weights then becomes n × m, where m is the size of x-q space search volume. This indicates that the size of matching weights is about m times larger than the DMRI data, causing memory issues. A straightforward solution to this problem is to perform angular upsampling in overlapped blocks and then combine the results to form the HAR DMRI data. Another solution is to reduce the number of matching weights to save memory cost. For instance, we can select a certain number of top weights, instead of using all the weights given by neighborhood matching.

CONCLUSION
We have presented a regularization framework for q-space upsampling. The relationships of signals in x-q space are used to regularize the inverse problem associated with recovering the HAR DMRI data. Extensive experiments on synthetic and infant DMRI data indicate that our method is able to produce HAR DMRI data with significantly improved quality. Future research effort will be directed to extending the current framework for resolution enhancement in joint x-q space.

AUTHOR CONTRIBUTIONS
GC and P-TY implemented the code and designed the experiments. GC drafted the manuscript. P-TY revised the manuscript. WL provided the clinical data. BD, YZ, and DS participated in idea discussion and reviewed the manuscript.