Super-Resolution Asymmetry Spectrum Imaging (SR-ASI)

Abstract

We introduce a technique for super-resolution reconstruction of diffusion MRI, harnessing fiber-continuity (FC) as a constraint in a global whole-brain optimization framework. FC is a biologically-motivated constraint that relates orientation information between neighboring voxels. We show that it can be used to effectively constrain the inverse problem of recovering high-resolution data from low-resolution data. Since voxels are inter-related by FC, we devise a global optimization framework that allows solutions pertaining to all voxels to be solved simultaneously.

Methods

The low resolution (LR) DW volumes are considered as down-sampled versions of high resolution (HR) DW volumes that need to be estimated. For this purpose, we let \(\mathbf{Y} \in \mathbb{R}^{N_{v} \times M_g}\) denote a set of LR DW volumes with \(N_{v}\) voxels per volume and \(M_{g}\) gradient directions across different shells. The LR volumes \(\mathbf{Y}\) are related to the HR volumes \(\mathbf{S} \in \mathbb{R}^{M_{v} \times M_{g}}\) by \(\mathbf{Y} = \mathbf{DS} + \boldsymbol{\mu}\), where \(\mathbf{D}\) is a down-sampling matrix that averages neighboring voxels and \(\boldsymbol{\mu}\) denotes noise. The sampling of \(\mathbf{Y}\) can be encoded by a binary matrix \(\mathbf{B}\) such that \(\mathbf{B} \odot \mathbf{Y} = \mathbf{B} \odot (\mathbf{DS} + \boldsymbol{\mu})\), where \(\odot\) is the element-wise multiplication operator. The model is general and can account for different spatio-angular subsampling schemes determined by \(\mathbf{B}\) and \(\mathbf{D}\).

Our method uses fiber continuity across voxels to regularize the inverse problem of recovering the HR volumes. Fiber continuity is a natural constraint supporting the fact that orientations should be consistent along fiber tracts. A model based on asymmetric fiber orientation distribution functions (AFODFs) is used to represent the sub-voxel configurations. In DMRI, the signal profile is typically assumed to be antipodal symmetric, implying that orientation in a positive hemisphere is always identical to its counterpart in the negative hemisphere. Symmetric FODFs are not sufficient in regions with fanning and bending fiber trajectories. To improve sub-voxel fiber continuity, we estimate AFODFs by incorporating information from neighboring voxels, allowing FODF asymmetry to better represent complex configurations {Cetin_Karayumak2018-om,Bastiani2017-hm,Reisert2012-xo}. The discontinuity of AFODFs, \(\boldsymbol{\mathcal{F}}_{\text{aniso}}\), over all voxels and directions (\(\mathbf{u}\in\mathbb{S}^2\)) is measured as \(\Phi(\boldsymbol{\mathcal{F}}_{\text{aniso}})= \int_{\mathbf{u}\in\mathbb{S}^2}\left\| \boldsymbol{\mathcal{F}}_{\text{aniso}}(\mathbf{u})\otimes \mathbf{W}(\mathbf{u})- \boldsymbol{\mathcal{F}}_{\text{aniso}}(\mathbf{u})\right\|_{\text{F}}^2 d\mathbf{u}\), where \(\mathbf{W}\) is a normalized directional probability distribution function as defined in {Wu2020-ox} and \(\otimes\) is a column-wise multiplication operator.

We use a multi-tissue spherical convolution model to represent the diffusion signal at each voxel of \(\mathbf{S}\). Specifically, \(\mathbf{S}\) can be decomposed as \(\mathbf{S}=\boldsymbol{\mathcal{RX}}\), with \(\boldsymbol{\mathcal{R}}=\left[\mathbf{R}_{\text{aniso}}, \mathbf{R}_{\text{aniso}}, \mathbf{R}_{\text{iso},1},\cdots,{\mathbf{R}_{\text{iso},n}} \right]\) being the multi-tissue axially symmetric response function (RF) consisting of an anisotropic component \(\mathbf{R}_{\text{aniso}}\) and \(n\) isotropic compartments \(\{\mathbf{R}_{\text{iso},i}\}_{i=1}^n\) as defined in {Wu2019-pe}. Each column of \(\boldsymbol{\mathcal{X}} = \left[ \mathbf{X}_{\text{aniso}}^{+}, \mathbf{X}_{\text{aniso}}^{-}, \mathbf{X}_{\text{iso},1},\cdots,\mathbf{X}_{\text{iso},n} \right]^{\top}\) represents the spherical harmonics (SH) coefficients of AFODF at each voxel location. The AFODF is represented by two sets of even-order spherical harmonics {Wu2020-ox}, capturing potentially asymmetric sub-voxel fiber configurations. That is, \(\boldsymbol{\mathcal{F}}_{\text{aniso}}(\mathbf{u})= \begin{cases} \mathbf{A}(\mathbf{u})\cdot\boldsymbol{\mathbf{X}}_{\text{aniso}}^{+}, & \mathbf{u}\in\mathbb{S}_{+}^2,\\ \mathbf{A}(\mathbf{u})\cdot\boldsymbol{\mathbf{X}}_{\text{aniso}}^{-}, & \mathbf{u}\in\mathbb{S}_{-}^2,\\ \end{cases}\) where \(\mathbf{A}(\mathbf{u})\) is the real even-order SH basis function sampled in direction \(\mathbf{u}\), and the column vectors of \(\boldsymbol{\mathbf{X}}_{\text{aniso}}^{+}\) and \(\boldsymbol{\mathbf{X}}_{\text{aniso}}^{-}\) are the corresponding SH coefficients for the positive and negative hemispheres (\(\mathbb{S}_{+}^2\) and \(\mathbb{S}_{-}^2\)).

Based on Eq.eqref{eq:sr}, we estimate \(\mathbf{S}\) and \(\boldsymbol{\mathcal{X}}\) by solving \(\min_{\mathbf{S},\boldsymbol{\mathcal{X}}} \frac{1}{2}\left\| \mathbf{B} \odot ( \mathbf{Y} - \mathbf{DS} )\right \|^2 + \lambda \left\| \boldsymbol{\mathcal{L}}\boldsymbol{\mathcal{X}} \right\|^2, \\ ~\text{s.t.}~ \left\{ \begin{array}{lr} \mathbf{S} = \boldsymbol{\mathcal{R}}\boldsymbol{\mathcal{X}},\\ \boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{X}} \succeq 0,\\ \Phi(\boldsymbol{\mathcal{F}_{\text{aniso}}})\le \epsilon,\\ \mathbf{S} \succeq 0, \end{array} \right.\) where the Laplace-Beltrami regularization {Tournier2007-kg} is realized via \(\boldsymbol{\mathcal{L}}=\text{diag}\left[\mathbf{L},\mathbf{L},0,\cdots,0\right]\) where \(\mathbf{L}_{jj}= \ell_j(\ell_j+1)\) with \(\ell_j\) being the order associated with the \(j\)-th coefficient. \(\lambda\) is the regularization parameter and the matrix \(\boldsymbol{\mathcal{A}}= \text{diag}\left[\mathbf{A},\mathbf{A},a_0,\cdots,a_0\right]\) maps the SH coefficients \(\boldsymbol{\mathcal{X}}\) to the AFODF amplitudes for imposing AFODF non-negativity. \(a_0=\sqrt{\frac{1}{4\pi}}\) is the \(0\)-th order SH basis function. Problem eqref{eq:model2} is solved using alternating direction method of multipliers (ADMM). Specifically, the problem can be cast as a strictly convex quadratic programming (QP) problem, which can be solved by an operator splitting solver for quadratic programs (OSQP){Stellato2020-fp}. Following common super-resolution frameworks {Shi2015-dv,Ning2016-kk,Lum2015-wk}, we used linearly interpolated data to initialize ADMM.

The WM RF is estimated for each b-shell from the average of the reoriented diffusion signal profiles of voxels with high anisotropy in WM. Also, isotropic RFs are set with diffusivity values ranging from \(0.1 \times 10^{-3}\,\text{mm}^2/\text{s}\) to \(3 \times 10^{-3}\,\text{mm}^2/\text{s}\) and step size \(0.2 \times 10^{-3}\,\text{mm}^2/\text{s}\). This will allow better characterization of tissue microstructure in infant brain.

We used k-fold cross-validation to determine the optimal regularization parameter \(\lambda\) in~eqref{eq:model2}. Two thousand random WM voxels from an BCP DMRI dataset and 61 candidate values for \(\lambda\), logarithmically spanning the interval between 0.01 and 10, were considered {Wu2020-vl,Wu2020-jq,Wu2020-vl,Bastiani2017-hm}. We found that the optimal value was 0.01. This value was applied to all experiments. The method is publicly available as part of the Computational Medical Imaging Toolkit (url{https://www.comedi.io}).

References

1. Wu, Y., Lin, W., Shen, D., Yap, P.-T., Consortium, U.B.C.P., others, 2019d. Asymmetry Spectrum Imaging for Baby Diffusion Tractography, in: International Conference on Information Processing in Medical Imaging. Springer, pp. 319–331.

2. Wu, Y., Lin, W., Shen, D., Yap, P.-T., 2019b. Improving Tractography in Baby Diffusion MRI via Asymmetric Spectrum Imaging, in: Proceedings of the International Society of Magnetic Resonance in Medicine (ISMRM).

3. Wu, Y., Hong, Y., Ahmad, S., Chang, W., Lin, W., Shen, D., Yap, P,-T. 2020. Globally Optimized Super-Resolution of Diffusion MRI Data via Fiber Continuity, in: International Conference on Information Processing in Medical Imaging.