Overview

Echo planar imaging (EPI)¹ is commonly used clinically for its speed, but is sensitive to system field inhomogeneity² and subject-specific magnetic susceptibility.³ Off-resonance is typically neglected in image reconstruction, resulting in geometric distortion along the phase-encoding dimension and signal dropout in regions of susceptibility variation. Multi-shot techniques encoding an auxiliary point-spread-function (PSF) dimension overcome these limitations,^4–6 showing potential in diffusion and anatomic applications.^7,8 However, multi-shot encoding increases scan time, motivating acquisition and reconstruction techniques to push acceleration. PSF-EPI anatomic scans may also exhibit residual flow-induced ghosting artifacts.⁸ This work proposes a model-based iterative reconstruction (MBIR) framework with low-rank constraints for anatomic PSF-encoded EPI scans. Circumventing modular unfolding operations on intermediary datasets, we reconstruct undistorted images directly from undersampled data. Images generated with MBIR have reduced residual ghosting and improved signal recovery in regions of susceptibility, demonstrated with brain MRI data acquired on a compact 3T system.⁹

Methods

Signal Model & Image Reconstruction
The PSF-encoded EPI signal measured during readout $$$u\in [0,U)$$$ of phase-encoding line $$$v\in [0,V)$$$ within an auxiliary pass $$$p \in [0,P)$$$ has triple Fourier encoding respectively in $$$k_x$$$, $$$k_y$$$, and $$$k_a$$$. The discretized¹⁰ signal can be modeled as
$$
\mathbf{G}[u,v,p,c]
=
\sum_{i=0}^{M-1}
\sum_{l=0}^{N-1}
s_c[i,l]
f[i,l]
e^{-\jmath \Delta \omega [i,l] \tau[v]}
e^{-\jmath (k_x[u] x[i] + k_y[v] y[l] + k_a[p] y[l])}
+
\mathbf{\epsilon}[u,v,p,c]
\tag{1}
$$
where we neglect T2* decay, assume row-flipping and pre-correction of eddy current and gradient delays,^11,12 presume linear gradients,¹³ and model off-resonance effects manifesting only along the phase-encoding dimension during high-bandwidth EPI readouts.¹⁴ Here, $$$f$$$ is the target signal in pixel indices $$$i\in[0,M)$$$ and $$$l\in[0,N)$$$, $$$\Delta \omega$$$ is the B₀ field map, $$$s_c$$$ is the sensitivity profile for coil $$$c\in [0,C)$$$, $$$\tau[v]=vT$$$ for echo spacing $$$T$$$, and $$$\mathbf{\epsilon}$$$ is Gaussian noise.

Without access to a field map, we target a composite unknown $$$\mathbf{\Theta}[i,l,q]\overset{\Delta}{=} \rho[i,l] e^{-\jmath \Delta \omega[i,l] \gamma_q T}$$$, i.e., the product of magnetization and the off-resonance field, where $$$\rho(x,y) \equiv f(x,y) e^{-\jmath \Delta \omega(x,y)t_0}$$$ for time $$$t_0$$$ between RF isodelay and initial readout. We partition the measurement window into $$$Q$$$ time segments,¹⁵ each with $$$q$$$^th segment center $$$\gamma_q$$$, assuming constant temporal field inhomogeneity evolution in each segment,¹⁶ enabling FFTs and implicitly presuming low-rank properties. This yields a 3-way tensor with third unfolding $$$\mathbf{\Theta}_{[3]} = \sum_{i,l,q} \delta_{i+lM} \mathbf{\Theta}[i,l,q] \delta_q^T$$$, where $$$\delta$$$ is Kronecker's delta; scaled by a time-bandpass filter $$$\mathbf{\Phi}_q$$$ segmenting the acquisition window into $$$Q$$$ time segments and $$$Q+1$$$ breakpoints:
$$
\mathbf{G}[u,v,p,c]
=
\sum_{q=0}^{Q}
\mathbf{\Phi}_q[v]
\sum_{i=0}^{M-1}
\sum_{l=0}^{N-1}
s_c[i,l]
\rho[i,l]
e^{-\jmath \Delta \omega [i,l]\gamma_q T}
e^{-\jmath (k_x[u] x[i] + k_y[v] y[l] + k_a[p] y[l])}
+
\mathbf{\epsilon}[u,v,p,c]
\tag{2}
$$
This abstracts to
$$
\mathbf{G}_{[3]}
=
\sum_{q=0}^{Q}
\sum_{c=0}^{C-1}
\big(
\delta_c \otimes \mathbf{\Lambda}_q\text{diag}(\mathbf{S}_{[3]} \delta_c )
\big)
\mathbf{\Theta}
\mathbf{\Phi}
+
\mathbf{\epsilon},
\tag{3}
$$
where $$$\mathbf{S}_{[3]}$$$ is the $$$MN\times C$$$ unfolding of sensitivity maps, with Kronecker's product $$$\otimes$$$ and triple Fourier encoding operator
$$
\mathbf{\Lambda}_q
\overset{\Delta}{=}
(\mathbf{\Xi}_\text{a} \mathbf{F}_\text{a} \otimes I)
\big( \text{diag} (\mathbf{\Xi}_\text{p} \mathbf{F}_\text{p})^T \otimes I \big)
(I \otimes \mathbf{F}_{\text{ro}})
\tag{4}
$$
$$$\mathbf{F}_\text{a}$$$, $$$\mathbf{F}_\text{p}$$$, and $$$\mathbf{F}_{\text{ro}}$$$ respectively denote auxiliary, phase-encoding, and readout Fourier encoding; and $$$\mathbf{\Xi}$$$ are undersampling operators. The target $$$\mathbf{\Theta}$$$ can be estimated with least-squares optimization:
$$
\hat{\mathbf{\Theta}}
=
\underset{\mathbf{\Theta}}{\text{arg min}}
\bigg\{
\Omega(\mathbf{\Theta})\overset{\Delta}{=}
\bigg\lVert
\sum_{q=0}^{Q}
\sum_{c=0}^{C-1}
\big( \delta_c \otimes \mathbf{\Lambda}_q \text{diag} (\mathbf{S}_{[3]} \delta_c) \big)
\mathbf{\Theta}
\mathbf{\Phi}
-
\mathbf{G}_{[3]}
\bigg\rVert_F^2
\bigg\}
\tag{5}
$$
where $$$\lVert \cdot \rVert_F$$$ is the Frobenius norm. We solve (5) with conjugate gradient iteration.

Data Acquisition & Processing
Under an IRB-approved protocol, a healthy volunteer was scanned on a compact 3T MRI system.⁹ T2-weighted, PSF-encoded EPI (T2-DIADEM)^7,8 was completed according to Table 1, with in-plane acceleration factor 3 and reduced PSF encodings for a reduced field of view (rFOV) factor of 8. Preceding MBIR or standard reconstruction, row flipping, ramp sampling, and phase correction in k-space were completed, and sensitivity profiles estimated using vendor-provided tools. In a standard reconstruction, PSF-encoded images were unfolded in the phase-encoding dimension using SENSE;^17,18 then resulting data were unfolded in the PSF-encoding dimension before calculating distortion-free and distorted images.⁶ Instead, MBIR was completed with $$$Q=20$$$ time segments and 10 conjugate gradient iterations. Reconstruction and processing were completed in Matlab (The MathWorks, Inc., Natick, MA, USA). Post-processing included zero-padding, Fermi window apodization, intensity correction, and image-domain gradient nonlinearity correction.

Results

Images reconstructed with MBIR showed reduced ghosting artifact and modest advantages in recovering signal in regions of susceptibility. Figures 1 and 2 show comparisons between standard images and images generated with MBIR. Ghosting artifact is seen anterior to the basilar artery and near the optic nerve in standard images, but is reduced for MBIR. Figure 3 illustrates signal recovery from a region near the sinuses with line profiles, with increased recovery for MBIR compared to standard and sharper transition between air and anatomical signal.

Discussion & Conclusion

Images reconstructed with the proposed framework have reduced ghosting and improved signal recovery in regions adjacent to susceptibility. These preliminary findings suggest MBIR may enhance the capabilities of T2-weighted anatomic imaging with PSF-encoded EPI. With improved image quality at fixed acquisition time and reconstruction directly from undersampled data, MBIR could enable further acceleration, a key limitation to broad adoption of PSF-EPI techniques. Future work will encompass regularization to explore further potential for artifact suppression and signal recovery, and adaptation for diffusion applications.

Acknowledgements

This work is supported by National Institutes of Health (NIH) U01 EB024450, NIH U01 EB026979, Mayo Clinic Graduate School of Biomedical Sciences, and the National Science Foundation (NSF) Graduate Research Fellowship Program (GRFP).