Merry Mani^{1}, Hemant Kumar Aggarwal^{1}, Vincent Magnotta^{1}, and Mathews Jacob^{1}

Multi-shot diffusion-weighted imaging reconstructions are challenged by the inter-shot phase inconsistency that exists between the data from different shots. The MUSSELS algorithm enabled the direct reconstruction of the multi-shot k-space data by posing it as a low-rank based matrix recovery problem. The iterative algorithm has been shown to successfully recover the missing k-space samples in accelerated and non-accelerated acquisitions. However, the reconstruction time increases as the number of shots\acceleration increases. We propose a new formulation based on iterative re-weighted least squares that increase the computational efficiency of the matrix completion by several folds to speed up the recovery of multi-shot data.

Multi-shot echo-planar imaging (EPI) provides a straightforward means to improve the spatial resolution of diffusion-weighted images (DWIs). The multiple short EPI readouts limit blurring and distortion in the DWIs and enables short echo-times suitable for high-resolution imaging when combined with partial Fourier (pF) acceleration. Since the k-space of a given DWI is sampled over multiple TRs using the above scheme, a phase inconsistency exists between the DWIs corresponding to the different shots. A typical IFFT-based reconstruction will show significant ghosting and aliasing in the DWIs if the k-space data are combined without phase compensation. The recently proposed MUSSELS method^{1} enables a novel approach to combine the k-space data from the multiple shots without needing phase compensation. The method aims to recover the missing data in each k-space shot using matrix completion algorithms. The DWI is then recovered as the sum-of-squares of the different shot images; the images will be free of artifacts since the phase may be discarded.

The MUSSELS method typically involve several iterations to achieve reasonably stable recovery of the missing k-space samples especially for accelerated acquisitions such as pF, where only 55-65% of the samples per shot are known. In such cases, additional properties such as virtual shots^{1,2} are utilized to improve the reconstruction. However, the subsequent doubling of the shots and the matrix completion further increase the reconstruction time. We propose to speed up the k-space data recovery using an iterative re-weighted least squares (IRLS) strategy^{3,4}, which provides an estimate of the FIR filters at each step to make the reconstruction faster.

The key idea in the MUSSELS method is the low-rank property of the block-Hankel matrices which are derived from the annihilation relations that exist between the k-space shots^{1}. These relations are summarized as \begin{equation}~\underbrace{\begin{bmatrix}{\bf{\cal{H}}}(\widehat{m_1})~~{\bf{\cal{H}}}(\widehat{m_2})~~...~~{\bf{\cal{H}}}(\widehat{m_{N_s}})\end{bmatrix}}_{{{\bf{H}}}_1(\bf\widehat{m})}\underbrace{\begin{bmatrix} \hspace{1mm}\widehat{\boldsymbol\phi_2}\hspace{5mm}0\hspace{20mm}0\hspace{6mm}\widehat{\boldsymbol\phi_3}\hspace{6mm}\\~\hspace{-5mm}-\widehat{\boldsymbol\phi_1}\hspace{4mm}\widehat{\boldsymbol\phi_3}~\hspace{14mm}0\hspace{6mm}0\hspace{9mm}\\~\hspace{1mm}0\hspace{2mm}~-\widehat{\boldsymbol\phi_2}~\hspace{3mm}...~\hspace{3mm}\vdots~\hspace{2mm}-\widehat{\boldsymbol\phi_1}~\hspace{2mm}...\\~\hspace{1mm}0\hspace{8mm}~0~\hspace{20mm}0~\hspace{6mm}0~\hspace{8mm}\\~\hspace{1mm}\vdots\hspace{8mm}~\vdots~\hspace{18mm}~\widehat{\boldsymbol\phi_{Ns}}\hspace{3mm}\vdots~\hspace{8mm}\\~ \hspace{1mm}0\hspace{7mm}~0~\hspace{14mm}~-\widehat{\boldsymbol\phi_{Ns}}\hspace{2mm}0~\hspace{8mm}\\\end{bmatrix}}_{\bf{\hat{\Phi}}}=\begin{bmatrix}0~~0~~...~~0~~0~~...~\end{bmatrix}.\hspace{20mm}[1]\end{equation}where $$${\bf{m}}_i$$$ and $$${\bf{\phi}}_i$$$ are the multi-shot image data and the phase of the $$$i^{th}$$$ shot, $$$\widehat{{\bf{m}}_i}$$$ and $$$\widehat{{\bf{\phi}}_i}$$$ are their Fourier samples respectively. The above relation establishes the low-rank property of $$${{{\bf{H}}}_1(\bf\widehat{m})}$$$ which is utilized in the MUSSELS formulation for msDW k-space data recovery as:\begin{equation}\label{use_Hankel}\widehat{\tilde{\bf{m}}}=argmin_{{\bf{\widehat{m}}}}~\underbrace{||{\mathcal{A}}\left({\bf{\widehat{m}}}\right)-{\bf{\widehat{y}}}||^2_{\ell_2}}_{\mbox{data~consistency}}~+~\lambda~\underbrace{||{\bf{\bf{H_1}}}({\bf{\widehat{{m}}}})||_*,}_{\mbox{low-rank~penalty}}\hspace{20mm}[2]\end{equation}without the knowledge of the annihilation filters $$$\widehat{\bf{\Phi}}$$$. To improve the reconstruction of pF data, $$${{{\bf{H}}}_1(\bf\widehat{m})}$$$ can be modified to account for virtual shots using the conjugate symmetry property of the k-space data as\begin{equation}{{{\bf{H}}}_1(\bf\widehat{m})}=\begin{bmatrix}{\bf{\cal{H}}}(\widehat{m_1})~~~{\bf{\cal{H}}}(\widehat{m_2})~~...~~{\bf{\cal{H}}}(\widehat{m_{N_s}})~~~{\bf{\cal{H}}}(\widehat{m_1}^\dagger)~~{\bf{\cal{H}}}(\widehat{m_2}^\dagger)~~...~~{\bf{\cal{H}}}(\widehat{m_{N_s}}^\dagger)\end{bmatrix}\hspace{20mm}[3]\end{equation}where $$$\dagger$$$ represents the flipped conjugate of the k-space data. The resulting cost function is minimized using an alternating minimization scheme that alternates between a data consistency enforcement and matrix completion using singular value shrinkage.

The IRLS formulation allows one to re-write the nuclear norm term as \begin{equation}||{{{\bf{H}}}_1(\bf\widehat{m})}||_* \le ||{{{\bf{H}}}_1(\bf\widehat{m})}W^{1/2}||_F^2,\hspace{20mm}[4]\end{equation}where $$$W=[{{{\bf{H}}}_1^*(\bf\widehat{m})}{{{\bf{H}}}_1(\bf\widehat{m})}+~\epsilon~I]^{-1/2}.$$$ Using this property, we re-write Eq. 2 as\begin{equation}\widehat{\tilde{\bf{m}}}=~argmin_{{\bf{\widehat{m}}}}~||{\mathcal{A}}\left({\bf{\widehat{m}}}\right)-{\bf{\widehat{y}}}||^2_{\ell_2}~+~\lambda~||{\bf{H_1}}({\bf{\widehat{{m}}}})W^{1/2}||_F^2,\hspace{20mm}[5]\end{equation}where $$$_F$$$ represents the Frobenius norm. The above IRLS formulation alternates between a data consistency enforcement and a denoiser, which projects the k-space data to a constraint set. During each iteration, $$$W$$$ provides an estimate of annihilation filters $$$\widehat{\bf{\Phi}}$$$ which leads to faster recovery of the k-space samples. Moreover, the computational complexity of the singular value decomposition (SVD) is also reduced. If $$$U\Sigma~V^*$$$ represents the SVD of $$${{{\bf{H}}}_1(\bf\widehat{m})}$$$, the SVD of the gram matrix $$${{{\bf{H}}}_1^*(\bf\widehat{m})}{{{\bf{H}}}_1(\bf\widehat{m})}$$$ is given by $$$U\Sigma^2U^*$$$. Hence $$$W$$$ can be computed from the SVD of the gram matrix as $$$U\Sigma^{-1}$$$. Note that the SVD of $$${{{\bf{H}}}_1^*(\bf\widehat{m})}{{{\bf{H}}}_1(\bf\widehat{m})}$$$ involves less computational complexity than that of $$${{{\bf{H}}}_1(\bf\widehat{m})}$$$ even as the number of shots increase for the virtual shot strategy.

We tested the new IRLS-MUSSELS on the recovery of high-resolution DWI acquired using pF schemes. A 4-shot dual-spin-echo diffusion sequence was used for data acquisition with b-value=700s/mm2, 60 diffusion directions and one non-DWI using a 32-channel head coil. The imaging parameters were: FOV=210mmx210mm, matrix size=256x152, slice thickness=4mm, TE=84ms, NEX=1.

Figure 1 shows the reconstruction of a DWI using various implementations of MUSSELS. Here, the reconstruction accommodating virtual shots show recovery of fine structural details. The proposed IRLS implementation accommodating virtual shots achieve the same image quality, but reduce the reconstruction time by a factor of 7. Figure 2 shows the result of a diffusion tensor fit from the 60 DWIs reconstructed using IRLS MUSSELS without and with virtual shots. The improved recovery of fine structural details using virtual shots is clearly evident. Similar improvements are observed in other slices also (see fig 3).

1. Mani, M., Jacob, M., Kelley, D. and Magnotta, V. Multi-shot sensitivity-encoded diffusion data recovery using structured low-rank matrix completion (MUSSELS). Magn Reson Med. doi: 10.1002/mrm.26382, 2016

2. Blaimer, M., Gutberlet, M., Kellman, P., Breuer, F., Kostler, H., and Griswold, M.Virtual Coil Concept for Improved Parallel MRI Employing Conjugate Symmetric Signals. Magnetic Resonance in Medicine 61:93-102 (2009)

3. Chartrand, R.; Yin, W. "Iteratively reweighted algorithms for compressive sensing". IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2008. pp. 3869-3872.

4. Gregory Ongie, Sampurna Biswas, and Mathews Jacob, "Convex recovery of continuous domain piecewise constant images from non-uniform Fourier samples", IEEE Trans. Signal Process., vol. 66, no. 1, pp. 236-250, 2017.

Figure 1: Comparison of MUSSELS reconstruction without virtual shots (a) and with virtual shots (b) for partial Fourier data recovery using cost function [2]. The reconstruction employing virtual shots (b) show fine details in the images as highlighted by the arrows. (c) shows the reconstruction with virtual shots implemented using the proposed IRLS scheme given in [5]. The IRLS scheme achieves the same reconstruction quality and reduces the reconstruction time by a factor of 7. The reconstruction is performed in Matlab on Intel Xeon 3.5GHz CPU with 32GB RAM.

Figure 2: Comparison of IRLS based MUSSELS reconstruction without (top row) and with virtual shots (bottom row). (a),
(b) and (c) shows the FA map, the color-coded FA maps and the primary eigen vector maps. The zoomed region from the
color-coded FA maps highlight the improved reconstruction of the fine structural details.

Figure 3: Comparison of IRLS based MUSSELS reconstruction without (top row) and with virtual shots (bottom row). (a) shows a diffusion-weighted image and (b) shows the color-coded the FA maps computed from all 60 directions. The improved recovery of the fine structural detail is observed in regions highlighted by the arrows.