Multi-shot diffusion imaging holds great potential for enabling high spatial resolution diffusion imaging. The technique can also achieve lower TE to examine structures that are proximal to inhomogeneous fields/metal implants and for studies at higher field strengths. However, reconstruction of MS-DW data is a non-trivial process. The diffusion gradients designed to encode microscopic motion will also encode macroscopic motion between shots rendering different phases to different shot images. Conventional MS-DW reconstruction methods estimate the motion-induced phase associated with the individual shots as the first step which is then either removed using a conjugate-phase approach or modeled while combining the images from the shots. Phase estimation typically involves reconstructing a low-resolution version of the individual shot images either from navigator data or using highly regularized reconstructions of the individual shot data. The quality of MS-DW reconstruction greatly depends on the estimated phase. If the phase estimates are corrupted either due to high noise or under-sampling, the reconstruction will suffer. Here we develop a method to reconstruct MS-DW data that do not depend on motion-induced phase estimates.

This work exploits the recent theoretical developments introduced in a class of methods that uses annihilating filters in k-space to constrain k-space data estimation^1-5.The recovery of missing k-space data is posed as a structured matrix completion problem when the matrix qualifies for low-rankness; the low-rankedness arise when high linear dependencies exist between the k-space data points in a neighborhood. We derive here an annihilating filter formulation for the multi-shot multi-channel diffusion data. Denoting the DWI from $$$i^{th}$$$ coil and $$$l^{th}$$$ shot as $$$m_{il}(x)$$$, we can write the following two relations:

(i)For a given coil $$$i$$$, the shots $$$l,j$$$ differ by a phase $$$\Phi(x)$$$ in image domain: $$$m_{il}(x)=m_{ij}(x)\Phi_l(x),$$$ which can be equivalently expressed as a convolution in k-space. Since $$$\Phi(x)$$$ is a smooth function, it can be approximated using a filter ($$$\hat{\Phi}(k)$$$) of finite support in k-space. The convolution with this filter can be implemented as multiplication using Toeplitz matrices: $$${\cal{T}}(\hat{m}_{ij}(k))\cdot\hat{\Phi}_l(k)={\cal{T}}(\hat{m}_{il}(k)), \forall j,l$$$. Here, $$$\hat{(.)}$$$ denotes variables in k-space after 2-D Fourier transform (FT). Since this relation holds true for all the coils, we get a set of conditions which can be written in matrix form as: \begin{equation}\underbrace{\begin{bmatrix} {\cal{T}}(\hat{m}_{11}) & {\cal{T}}(\hat{m}_{12}) ... & {\cal{T}}(\hat{m}_{1N_s}) \\ \\{\cal{T}}(\hat{m}_{21}) & {\cal{T}}(\hat{m}_{22}) ... & {\cal{T}}(\hat{m}_{2N_s}) \\ \vdots \\{\cal{T}}(\hat{m}_{N_c1}) & {\cal{T}}(\hat{m}_{N_c2}) ... & {\cal{T}}(\hat{m}_{N_cN_s}) \end{bmatrix}}_{\bf{{\cal{T}}}_1(\hat{m})}\underbrace{\left[ \begin{array}{c} \hat \Phi_2 \\ -I\\ \vdots \\ 0 \end{array} \begin{array}{c} 0\\\hat \Phi_3 \\ -I\\ \vdots \end{array} \begin{array}{c} \cdots \end{array}\right]}_{\bf{\hat{\Phi}}}=0 \Rightarrow \bf{{\cal{T}}}_1(\hat{m})\bf{\hat{\Phi}} = 0. \end{equation}The above annihilating filter bank relation holds true for all the coils which implies that the matrix $$$\bf{{\cal{T}}}_1(\hat{m})$$$ is highly low-rank. Thus, even if we don't have an estimate of $$$\bf{\hat{\Phi}}$$$, we can impose a low-rank penalty on the matrix $$$\bf{{\cal{T}}}_1(\hat{m})$$$.

(ii)for a given shot $$$l, m_{il}=p(x)S_i(x); i=1:N_c$$$, where $$$S_i(x)$$$ is the coil sensitivity of the $$$i^{th}$$$ coil and $$$N_c$$$ is the number of coils. The coil sensitivities can be computed from the non-diffusion-weighted image for all slices. Combining the two relations, the reconstruction of the multi-shot k-space data can be written as the following reconstruction problem: $$\tilde {\bf m}= \text{argmin}_{{\bf \hat{m}}} {||\bf{{\cal{A}}\hat{m}-y}}||^2_{\ell_2} + \lambda||{\bf{\cal{T}}_2(\hat{m})}||_* ,$$ where $$$\bf{y}$$$ is the matrix of measured multi-channel multi-shot k-space data of dimension $$$N_1xN_2xN_cxN_{s}$$$. The operator $$$\bf{{\cal{A}}}$$$ is given by $$$\bf{{\cal{F}}}\circ\bf{{\cal{S}}}\circ\bf{{\cal{F}}^{-1}}$$$ where $$$\bf{{\cal{F}}}$$$ and $$$\bf{{\cal{F}}^{-1}}$$$ represent the FT and the inverse FT operation and $$$\bf{{\cal{S}}}$$$ represents multiplication by coil sensitivities. $$$\lambda$$$ is the regularization parameter and $$${\bf{\cal{T}}_2}$$$ is a block Toeplitz matrix consisting of the Toeplitz matrices of the multi-shot k-space data from all the channels. The above nuclear norm minimization problem can be solved using singular value thresholding.

Diffusion data was collected on GEMR9507T scanner at the University of Iowa equipped with a 32-channel Tx/Rx coil using a 4-shot Stajeskal-Tanner sequence (b-value=1000s/mm2; FOV=220mmx220mm; matrix size=128x128, slice thickness=1.7mm, TE=69.6ms). We compare the reconstruction of the MS-DW data to the standard MUSE method⁶ that relies on phase estimates for the recovery of the multi-shot data.