Deconvolution‐based distortion correction of EPI using analytic single‐voxel point‐spread functions

To develop a postprocessing algorithm that corrects geometric distortions due to spatial variations of the static magnetic field amplitude, B0, and effects from relaxation during signal acquisition in EPI.


| INTRODUCTION
Echo-planar imaging 1 is the most widely used fast MRI technique, especially in functional MRI (fMRI), diffusion-weighted (dw) MRI or studies of tissue perfusion. Such experiments typically use single-shot acquisitions with a relatively long echo train and a concomitantly low effective bandwidth along the phase-encoding (PE) direction. Consequently, the in-plane resolution may be degraded by blurring related to transverse relaxation with an effective rate R * 2 = 1∕T * 2 during image acquisition. 2 Even more prominent are geometric distortions due to spatially varying offsets of the Larmor frequency, Δ 0 = 2 Δ 0 . Such offsets may result from the presence of compounds of different magnetic susceptibility in the imaged object or from eddy currents induced by switching magnetic field gradients. 2,3 As R * 2 and Δ 0 increase with the amplitude of the static magnetic field, B 0 , the related image artifacts become increasingly prominent at higher field. They can be mitigated with parallel imaging techniques. 4,5 However, applications may be limited by the performance of the receive coil, such as in fMRI at laminar resolution, where the FOV may be of the order of the size of the coil elements. [6][7][8] A panoply of approaches has been developed to correct distortions in EPI. One family of algorithms involves nonrigid registration of the EPI data to geometrically correct structural scans, which are typically available in parallel to fMRI or dw-MRI acquisitions. 9 A physics-based procedure includes the calculation of a B 0 "field map" 10 from the signal phase in a separate multiple gradient-echo (GE) experiment and permits a quantitative a priori characterization of magnetic field inhomogeneities. 11 With the "conjugate phase method," the EPI k-space data can then be corrected by demodulation of the offset-induced phase error. 12 An accelerated version of the conjugate phase method is the so-called "frequency-segmented approach" (or "multifrequency reconstruction" [MFR]). 13,14 These techniques are capable of correcting geometrical shifts for any k-space trajectory. However, they may lack an adequate intensity correction in areas of large B 0 variations. Alternatively, "k-space energy spectrum analysis" has been proposed for dynamic field mapping and distortion correction directly from the EPI data. 15 "Reversed gradient approaches" use information from two EPI acquisitions with opposite polarity of the phase blips. [16][17][18] It is then assumed that both acquisitions yield identical signal intensity but geometric distortions in opposite directions, and that they can be corrected to the same undistorted template. A popular algorithm of this class is TOPUP, 19,20 implemented in FSL. 21 Recent simulations suggest that multiple-PE methods, such as TOPUP, are more accurate than multiple-echo (ME) B 0 mapping, while registration-based approaches yielded the poorest corrections. 22 Point-spread function (PSF) mapping techniques include a preparation scan to obtain a measured PSF n for each image voxel n. [23][24][25][26] They were proposed to outperform field-map techniques if equal time is allowed for PSF n mapping and field mapping. 24 An extended version includes the combination of two distortion-corrected EPI data sets of opposite blip polarity for additional image correction. 27 Most of these techniques inherently assume that the subject (and hence, the distortions) remains static during the experiment. This does not capture local B 0 fluctuations due to head motion, respiration, or scanner drifts. Dynamic distortion correction methods may address this limitation in fMRI time series. 15,[28][29][30][31] They often rely on ME acquisitions or multichannel phase information.
In the current work, we introduce a PSF method for deconvolution-based distortion correction (referred to as "DecoDisCo"). It uses information from a field map obtained from either ME or multiple-PE acquisitions. Unlike previous PSF-based corrections, it relies on prior knowledge of an analytic PSF n that considers the particular k-space trajectory. 32 From this input, a deconvolution operator is devised for a regularized computation of distortion-corrected and intensity-corrected data. Apart from frequency offsets, effects from T * 2 relaxation are also considered in the mathematical description, for potential deblurring. Furthermore, an algorithm is presented for the combination of images recorded with opposite blip polarity. Results from applications to in vivo imaging of human brain are compared to those obtained with existing methods including MFR and TOPUP.

| THEORY
We restrict our consideration to 2D-EPI variants and assume that frequency encoding is performed along the x-axis in N x steps and PE along the y-axis in N y steps. Next, we define matrices A and Ã of dimension N x × N y containing complex-valued signal intensities in an idealized undistorted object and in the measured image, respectively. The image is a convolution of a spatially varying complex PSF and the object. 33 For EPI, only distortions in PE direction are relevant, whereas the frequency-encoding bandwidth is typically large enough to keep distortions below the voxel size. 2 Instead of the entire object A (and image Ã ), it is therefore sufficient to consider an arbitrary object column a in PE dimension. Each voxel n of this column (0 ≤ n ≤ N y − 1) is characterized by its specific object intensity a n and an individual PSF n , depending on the distribution of Larmor frequencies. The PSF n is a function of the k-space trajectory, the TE, the local frequency offset Δ 0,n , and the local relaxation time T * 2,n , whereby the latter two reflect the mean and width of the intravoxel frequency distribution, respectively. 34 Due to the convolution, | 2447 PATZIG eT Al. PSF n contributes intensity to any voxel m of the corresponding column ã of image intensities (0 ≤ m ≤ N y − 1).

| Analytic point-spread functions
The PSF n is obtained by an inverse Fourier transform (ℱ −1 ) of the modulation transfer function of the pulse sequence in k-space. Consideration of Δ 0,n and T * 2,n and a truncated acquisition window leads to the following k-space weighting function for GE-EPI at the location of voxel n 35 : where t is the sampling time of k-space coordinate k y (with t = 0 defined by the onset of the acquisition window); and is the imaginary unit. The k-space lines are separated by increments of Δk y = 1∕FOV y and cover a range from k + y,min = − 2f p − 1 N y Δk y ∕2 to k + y,max = +N y Δk y ∕2 and from k − y,min = −N y Δk y ∕2 to k − y,max = 2f p − 1 N y Δk y ∕2 for bottom-up (indicated by "+") and top-down (indicated by "-") trajectories, respectively. The fraction f p ∈ 1∕2, 1 considers a potential partial Fourier GE (pGE) acquisition in PE direction. 36 Truncation of the acquisition window is accounted for by the indicator function 1 [a,b]  In partial-Fourier acquisitions, missing samples in one tile of k-space are computed as complex conjugates of the (measured) corresponding lines in the opposite tile ("conjugate filling"), which has not been considered above. The PSF n is then obtained as the sum of the result from Equation 3 for the acquired tile plus an equivalent expression (with time reversal) for the computed tile, as follows: An EPI variant referred to as double-shot EPI with centerout trajectories and intrinsic navigation (DEPICTING) 37 acquires two k-space segments in one TR. Both trajectories start at the central line and progresses outward in opposite directions by using opposite phase blip polarities. The sampling time of these trajectories can be written as This leads to a PSF n that is also obtainable from Equation 6 by assuming pGE-EPI with conjugate filling of one half of k-space ( f p = 1∕2): These derivations of PSF n can be adapted for other k-space trajectories. Parallel imaging 4,5 with an acceleration factor R is straightforwardly accounted for by rewriting Δk y → Δk y ⋅ R, FOV y → FOV y ∕R and N y → N y ∕R, whereas Δy = FOV y ∕N y remains unchanged. The corresponding PSF n is therefore obtained by replacing N y by N y ∕R in Equations 2, 4, and 7.

| Distortion correction by deconvolution
The image intensity ã m of a voxel m is the sum of all object intensities a n weighted by the value of PSF n at a position y mn = (m − n) ⋅ Δy according to the distance between voxels n and m: A compact expression for an entire image column (along the PE direction) is where PSF is a matrix of size N y × N y . Figure 1A shows an experimentally obtained Δ 0,n -column. The magnitude of the corresponding PSF + fGE matrix is given in Figure 1B. Nonzero PSF n y mn ⋅ a n .
(10) a = PSF ⋅ a, Δ 0,n values lead to apparent off-diagonal shifts of the center peak in the PSF + fGE matrix. An experimental T * 2,n -column may be constructed in a similar fashion (not shown), with shorter T * 2,n values leading to apparent broadening of the center peak. The least-squares solution of Equation 10 is obtained by the Moore-Penrose pseudo inverse. It can be computed from a singular value decomposition of the PSF matrix by decomposing it into a matrix product U ⋅ Σ ⋅ V * with U and V being unitary matrices. The diagonal matrix ∑ contains the singular values (decreasing series of real numbers; Figure  1C). Multiplying the pseudo-inverse, V ⋅ Σ −1 ⋅ U * , from the left to the image column vector ã yields a distortion-corrected and intensity-corrected representation of the object.
As deconvolutions are often ill-posed and ill-conditioned, 38 noise contributions are amplified, in particular, following inversion of small singular values (ie, due to voxels with low signal intensities increasing the sparsity of U and V) ( Figure 1C,D). Therefore, we applied a smoothing and interpolation procedure to the measured maps, 10 and Tikhonov regularization was integrated by introducing a parameter > 0 in the inversion of the singular values: This binds the condition number to (X) < ∑ 11 ∕ ∑ nn . 39

| Combination of images acquired with opposite phase blip polarity
While spatially correct coordinates can be straightforwardly reconstructed for a distorted image, this does not necessarily yield a correct intensity distribution. An accurate reconstruction is impossible if information is lost due to compression of signal contributions from different object positions. 20 However, directions of spatial signal shifts depend on the phase-blip polarity: Regions appearing compressed with bottom-up EPI will appear stretched with top-down EPI and vice versa. This provides additional information about the intensity distribution. In et al 27 presented a procedure for a weighted summation of pairs of PSF n -corrected images acquired with opposite blip polarities to resolve local loss of spatial information due to compression. It can be straightforwardly adapted to combine pairs of DecoDisCo-corrected GRE or SE-EPI acquisitions: The weights, ± m , are further scaled by an empirical exponent c. 27 They are obtained by normalizing the magnitude PSF matrices along their columns and calculating the row totals: A value of ± m > 1 reflects compression during image formation, whereas ± m < 1 reflects stretching. To give larger weight to information from originally stretched areas, c is chosen to be negative. The weights can also be computed from the inverse of the PSF matrix (ie, the matrix used to correct the images), leading to c > 0 as originally used by In et al. 27 For c = 0, the combination yields the mean of both contributions with an expected √ 2-fold improvement of the SNR. For c → −∞, a binary combination is obtained, reflecting the most accurate geometry, although without SNR gain. Other exponents lead to a trade-off between SNR and geometrical accuracy. Note that application of Equation 12 requires acquisition of the whole data set with both blip polarities (sometimes done in dw-MRI, seldom in fMRI).

| Image acquisition
Experiments were performed at 3 T (MAGNETOM Skyra Connectom and MAGNETOM Prisma fit ; Siemens Healthcare, Erlangen, Germany) using the body transmit coil and a 32-channel head receive coil. Pulse sequences were implemented under Siemens IDEA and tested with a spherical water phantom (1.25 g NiSO 4 ·6H 2 O per kg water). The deconvolutions were evaluated with a 3D-printed phantom (photopolymer) ( Figure 2D) filled with an MnCl 2 ·4H 2 O solution (61.9 mg per kilogram of water; T 1 ≈ 430 ms, T * 2 ≈ 27 ms measured in an unstructured container). 40 Rows of rods (nominal diameters increasing from 1.0 mm to 3.5 mm in steps of 0.25 mm) and two air-filled hollow cylinders for inducing distortions were integrated for evaluating the resolution. To measure PSF n , a PE table was added to a blipped GRE-EPI F I G U R E 2 Experimental and analytical PSF n in a resolution phantom. A, Measured PSF (with color-coded intensity a n ) in x-y EPI coordinates, where x is the (undistorted) frequency-encoding (FE) direction and y EPI is the distorted PE direction (the presentation of the PSF map and the notation of coordinates are consistent with sequence (TR = 2 seconds, TE = 25 ms, Δy = 1.5 mm, matrix = 128 × 128, f p = 5/8) as described by Zaitsev et al. 25 The PSF encoding was performed with full resolution (ie, no partial Fourier or parallel imaging were applied along this dimension). Experimental PSF n maps were extracted after inverse Fourier transform with an optional PSF deconvolution  For obtaining the PSF matrix in a geometrically undistorted space, a 2D-ME FLASH sequence 43 was used with 12 bipolar gradient echoes, an interecho time of 1.2 ms or 1.6 ms, and identical slice and resolution parameters as in the GE-EPI and DEPICTING scans. The value of Δ 0,n was extracted from a linear fit to the unwrapped phases of all echoes. The T * 2,n values were obtained by exponential fitting of the echo magnitudes. Maps of the goodness of fit of both parameters were stored for postprocessing. In selected investigations of T * 2 -related blurring, another ME-FLASH scan was recorded with 24 echoes and 4.5-ms spacing to improve the accuracy of the T * 2,n fits. Unless otherwise noted, GRE-EPI and DEPICTING acquisitions were corrected using the ME-FLASH data, whereas the dw-MRI scan was corrected using Δ 0,n estimates obtained with TOPUP from the gradientreversed acquisitions.

| Image processing
All postprocessing routines were implemented on a standard desktop computer using MATLAB (MathWorks, Natick, MA). Voxels of the Δ 0 maps with a goodness of fit smaller than the mean over the slice were labeled as uncertain, and their values were omitted by thresholding or, alternatively, interpolated using discrete cosine transforms. Finally, maps were smoothed using a small Gaussian kernel (FWHM = 2 pixels).
The PSF matrices were constructed from cyclically shifted PSF ± n y mn with measured input values Δ 0,n and (optionally) T * 2,n . After decomposing the PSF matrices using MATLAB's singular value decomposition function, the inversion was performed with variable . For comparison, the MFR approach 13 was also implemented in MATLAB, and TOPUP 20 was used as a stand-alone application.

| Simulations
For numerical investigations of the combination of pairs of PSF n -corrected images, synthetic fGE-EPI images were generated using measured image-domain data as input. A forward convolution was applied to the experimental images with a synthetic PSF matrix that was obtained using experimental Δ 0,n values rescaled to 1.5 ⋅ Δ 0,n to amplify geometric distortions and T * 2 → ∞ (ie, R * 2 = 0) to avoid additional T * 2 related blurring. Complex-valued Gaussian noise (average magnitude of 22% of the average signal magnitude of the masked brain region) was added, and a regularized PSF deconvolution ( = 0.01) was applied. The corrected image pairs with opposite blip polarities were combined using Equation 12 with different settings of c.

| RESULTS
Distortions were evaluated separately for fGE-EPI and pGE-EPI. Similarly, partial Fourier and T * 2 effects were evaluated separately. Unless otherwise noted, voxel-specific corrections of T * 2 effects were omitted by arbitrarily setting T * 2,n = 50 ms, a typical (uncorrected) value in cortical gray matter at 3 T. 44-46

| Comparison of experimental and analytical PSF
The PSF n mapping yielded a position-dependent shift of the experimental PSF due to local frequency offsets (Figure 2A). Application of DecoDisCo with a separately recorded Δ 0 map aligned the PSF at the correct position in undistorted PE direction ( Figure 2B). Note that T * 2 -related blurring was not corrected here. A remaining deviation from the undistorted position is visible at the edges, resulting from less accurate frequency offset estimates at the phantom border ( Figure 2E).
Profiles of the experimental and analytical PSF showed excellent agreement in the central part, whereas the weak oscillatory pattern due to the partial-Fourier reconstruction with zero filling was smoothed in the experimental PSF ( Figure 2C). This might be related to the scanner implementation of the parallel-imaging reconstruction. Figure 3 shows high-resolution fGE-EPI images and a corresponding Δ 0 map. The comparison with the MPRAGE image ( Figure 3A) demonstrates that distortions on the uncorrected EPI images ( Figure 3B,C,H,I) have expected opposite directions for inverted PE directions. Good agreement of the two corrected EPI images is obtained after PSF deconvolution ( Figure 3E,F,K,L). Both results agree well with the reference, indicating sufficient accuracy for correcting highresolution EPI images.

| Correction of distortions in full-Fourier and partial-Fourier GE-EPI
Applications to pGE-EPI are presented in Figure 4. Note that typical partial-Fourier reconstruction algorithms, such as the homodyne 47 or the projection onto convex set 48 techniques, use assumptions on the phase evolution extracted from the symmetrically acquired portion of k-space data. As their exact scanner implementations were not available, we used simpler k-space filling techniques during image reconstruction for a proof of concept: (1) zero filling (ie, missing k-space data filled with zeroes) (Equation 3) and (2) conjugate filling (Equation 6). No assumptions about the phase were made in Equation 6, and no additional convolutions were applied to avoid typical errors in partial-Fourier reconstructions. 49 Compared with the reference (Figure 4A), DecoDisCo achieved similar distortion corrections with zero and conjugate filling ( Figure 4C,F). However, sharper images (ie, reduced blurring) were obtained with conjugate filling ( Figure 4F). For further evaluation, we simulated exponential signal decays with T * 2 values of brain tissue and partial-Fourier sampling. Increased blurring exceeding T * 2 effects resulted after decreasing f p , due to the shortened effective acquisition window ( f p ⋅ N y ⋅ Δt e ) and correspondingly higher level of truncation with zero filling (data not shown).

| Correction of distortions in double-shot center-out EPI
Intensity errors occur in distorted images if multiple voxels are mapped onto the same location. As DEPICTING uses a segmented acquisition with opposite blip polarity, frequency offsets lead to opposite voxel shifts for the two tiles, and the PSF shows a double-peak pattern. 37 Therefore, DEPICTING is more prone to intensity errors than standard EPI. The very short TE achieved with its center-out trajectory yields high signal intensities with minimal phase evolution. Consequently, constructive interferences of shifted signal contributions produce characteristic hyperintensities. These images are particularly suited to evaluate the intensity distribution obtained after distortion correction. Figure 5 demonstrates that both the MFR and DecoDisCo restored the right brain shape. Spurious hyperintensities (indicated by arrows) remained after MFR ( Figure 5C,G), whereas DecoDisCo partly restored piled-up signal ( Figure 5D,H). These regions contain contributions of a corrected intensity distribution from previous voxel stretching plus a smeared average intensity from previous voxel compression. The improved intensity distribution might improve co-registration results.

-related blurring
Acquisitions with fGE-EPI, pGE-EPI (with zero filling), and DEPICTING in the resolution phantom are compared in Figure 6. Uncorrected images ( Figure 6A-C) showed substantial distortions (appearing as double contours with DEPICTING), which were corrected with MFR ( Figure 6E-G) or DecoDisCo and Δ 0 maps ( Figure 6H-J; no consideration of T * 2 ). Unrecovered signal voids due to signal compression appeared in regions affected by strong susceptibility gradients (regions of interest 1, 2, and 4 with fGE-EPI and pGE-EPI; regions of interest 1 and 4 with DEPICTING).
While the differently sized rods appeared sharp with fGE-EPI, substantial and moderate blurring (along PE direction) was found with pGE-EPI and DEPICTING, respectively ( Figure 6E-J). Deblurring by consideration of T * 2 maps had only minor effects for fGE-and pGE-EPI, but improved the sharpness of DEPICTING images ( Figure 6K-M; compare also profiles in Supporting Information Figure S1). Figure 7 demonstrates corrections of T * 2 -induced blurring in fGRE-EPI in vivo acquisitions. As a zeroth-order approximation for T * 2 deblurring, spatially uniform T * 2 maps (with different T * 2 values between 10 ms and 80 ms) were assumed in the in vivo example, yielding already a sufficient level of correction obtained with T * 2 ≥ 40 ms in most brain regions.

| Combination of EPI scans with opposite blip polarity
Despite good general agreement of the corrected bottomup and top-down images (see Figure 3E,F,K,L), the signal differed for different blip polarities, especially in areas of large Δ 0,n . Figure 8A shows one-dimensional profiles of simulated top-down and bottom-up fGRE-EPI acquisitions after PSF n correction (red and blue solid lines, respectively) in comparison to the undistorted input object (black dashed F I G U R E 5 Demonstration of the inherent intensity correction of the PSF deconvolution method. B,F, Uncorrected axial double-shot EPI with center-out trajectories and intrinsic navigation (DEPICTING) acquisitions (ie, double-shot center-out GRE-EPI trajectories sampled with very short TE of 4.6 ms) at different levels through the brain demonstrating both geometric distortions along the PE direction (anterior-posterior) in comparison to MP2RAGE reference images (A,E) as well as hyperintensities. C,G, Results obtained with a multifrequency reconstruction (MFR) that achieves distortion correction but fails to correct the hyperintensities. D,H, Results obtained with DecoDisCo with slightly superior distortion correction as well as an improved intensity profile (arrows). Both the MRF and DecoDisCo used the same Δ 0 map, which was obtained with a separate ME-FLASH acquisition. The value of T * 2,n was arbitrarily set to 50 ms without consideration of spatially specific variations line). Combined profiles obtained with Equation 12 are shown in Figure 8C. Three different c values were considered ( Figure 8B). Only minor differences were observed for the weighted (c = −4) and the either/or binary combination, whereas for the arithmetic mean (c = 0), a larger deviation from the input profile is evident, particularly in the left part.
Corresponding 2D images are shown in Figure 9. Red arrows in the uncombined images indicate areas that were not accurately reconstructed due to signal compression. These artifacts are also visible, with reduced intensity, in the combined image corresponding to the arithmetic mean ( Figure 9D). They are mitigated for c ≠ 0, which yields an improved signal

| Combination of DecoDisCo and TOPUP
Deconvolutions tend to create oversharpened images due to the inversion of noise, which requires additional regularization in comparison to a convolution-based voxel-shift used with TOPUP. A combination of the two methods may therefore result in sharper images without the need to measure a field map: If acquisitions with opposite blip polarities are available (as often performed in dw-MRI), an estimate of the underlying field map may be generated with TOPUP. Subsequently, distortion correction and pair-wise (weighted) combination of PSF n -corrected images may be performed with DecoDisCo. An example applied to SE-EPI with and without diffusion weighting is shown in Figure 10. Geometric distortions were corrected by achieving increased similarity of acquisitions with opposite gradient polarities (Figure 10, DecoDisCo). The combined images (here c = −2) were superior to the individually corrected images, suggesting that spatial accuracy is gained from the combination. A comparison to TOPUP corrections with the "least-squares restoration" mode yielded similar performance of both approaches, which were based on the same underlying field map. Visual inspection did not indicate different sharpness of the conjugatefilling reconstruction ( f p = 6∕8) over the TOPUP result.

| DISCUSSION
A versatile procedure is introduced to correct geometric distortions in EPI-type acquisitions. It is based on an analytical model of the PSF, a preparation scan to acquire maps of the local B 0 (eg, with multiple-PE or ME methods) and (optionally) T * 2 distributions as input parameters, and a regularized PSF deconvolution.
Strong gradients of Δ 0,n lead to insignificant rows in the PSF matrix, whereas short T * 2,n values result in insignificant columns. These effects decrease the conditioning of the system, which requires regularization of the inversion. Thikonov regularization with , depending on the acquisition scheme and noise level, improved the conditioning sufficiently to yield noncorrupted images.
Interestingly, the double-peak pattern of the DEPICTING PSF (Equation 8) also improved the conditioning of the associated inversion problem, as the corresponding PSF matrix is the sum of two (half-Fourier GE-EPI) PSFs of opposite blip polarity. This yields additional information, which may be exploited in future work for an improved reconstruction. 50 Due to multiple-to-one mappings, specific image regions may not be reconstructed correctly if only acquisitions with a single blip polarity are available. Combined GE-EPI or SE-EPI acquisitions with opposite polarity achieve improved corrections, as demonstrated with the adaptation of a previously published scheme for a weighted summation. 27 A potential advantage of DecoDisCo is that experimental mapping is not required, thus avoiding, for example, flow artifacts from CSF or ghosting, which may affect the measured PSF.
Evaluation of the deblurring efficiency requires careful consideration of the imaging readout. 51 With fGE-EPI, high spatial frequencies carrying information about fine structures are acquired early (ie, with minimal T * 2 -weighting) at one edge of k-space, but late at the other edge (ie, with maximal T * 2 -weighting). Despite the signal decay during the readout, information on high spatial frequencies (ie, information on "sharpness") is available from the initially acquired k-space lines, which efficiently mitigates related blurring. Consequently, further deblurring is relatively inefficient. This is not a limitation of the method, but more an intrinsic advantage of the well-posed conditions of fGE-EPI. With pGE-EPI (and, particularly, with small), high spatial frequencies are acquired late. Hence, no compensating information is available from the other k-space half, leading to degraded sharpness of partial-Fourier acquisitions (in particular, with zero-filling reconstructions). Of note, the full complex PSF is required to characterize this specific contribution to blurring, whereas a restriction to the magnitude as used in earlier discussions of pGE-EPI 35,37,52,53 is inadequate. 51 With DEPICTING, the entire k-space is acquired (albeit in two tiles), mitigating the specific loss of sharpness associated with pGE-EPI and zero filling. However, both k-space edges are sampled after significant relaxation. This specific T * 2 -related blurring can be corrected with DecoDisCo if a sufficiently accurate map is available.
Motion producing an inconsistent field-map estimate, field drifts, or higher-order field variations were not considered in the mathematical model and, hence, remain as potential sources of image artifacts. For example, respiration-induced dynamic field changes were shown to produce dynamic changes. 54,55 These are general limitations of static correction methods 37 and not specific to DecoDisCo. Such effects may be addressed by integrating prospective motion correction 56 or information F I G U R E 8 Results from combining synthetic PSF-based distortion-corrected fGE-EPI images acquired with opposite PE directions (onedimensional profiles through the slices shown in Figure 9). Red and blue solid lines indicate corrected profiles obtained with top-down ("EPI -") and bottom-up ("EPI + ") trajectories. obtained with a dynamic field camera. 57 Additionally, GE-EPI acquisitions with opposite blip polarities might be differently modulated by variations of the effective induced by echo shifting. 58,59 This might lower the accuracy of a pairwise combination of PSF n -corrected images. In dw-MRI, only a single pair of non-dw acquisitions with reversed PE is, in general, sufficient for obtaining a field map, which might be used to correct all dw data. A potential challenge is the shot-to-shot phase instability of acquisitions using strong diffusion-weighting gradients. Although previous work has shown that sophisticated phase correction is possible even under such conditions, 60 this strategy requires future investigations.
In our approach, T * 2 relaxation is modeled as a monoexponential decay. This is likely an oversimplification due to partial-volume effects or combined local and long-range contributions to the field distribution inside a voxel. More complex descriptions of the signal decay have been published (eg, Refs. 45,46,61), suggesting potential adaptation of the PSF. However, the results from Figure 7 demonstrate that T * 2 -related blurring can already be mitigated by assuming a uniform T * 2 according to the longest T * 2 (or longer) of the image volume. In our experiments, values of 40-50 ms served as a reasonable approximation for expected gray-matter values in the human brain at 3 T. 45,46 Comparing DecoDisCo with other distortion corrections, we note that the MFR technique also requires a separate field-map acquisition. In our experiments, the PSF-based corrections recovered more image details than an MFR with substantially reduced intensity artifacts ( Figure 5). On the other hand, the MFR is universally applicable to any k-space trajectory without further adaptation, whereas an analytical approximation of the PSF must be derived for DecoDisCo.

F I G U R E 9
Results from combining synthetic PSF-based distortion-corrected fGE-EPI images acquired with opposite PE directions (2D slices corresponding to the profiles shown in Figure 8). In addition to the input image (A) and the corrected single acquisitions with top-down (B) and bottom-up trajectories (C), results from weighted summations with different exponents (c = 0 [D], c = −4 [E], and c → −∞ [F]) are also presented. B,C, Areas that are not sufficiently corrected due to signal compression are indicated by red arrows. This information is (partly) recovered in weighted summations with c ≠ 0. Weighted averaging with c = −4 (empirically obtained) yielded a smaller MSE than with c = 0 (5.52 vs 6.17, respectively), whereas a larger MSE (7.36) resulted with c → −∞ due to the higher noise level in the binary combination TOPUP is a popular method for an offline field-map estimation. It relies on a more heuristic algorithm, aiming at mapping pairs of gradient-reversed images onto each other. This cannot correct multiple-to-one mappings. In principle, it can be beneficial to combine multiple-PE field-map estimation methods and PSF-based deconvolution to obtain a first-order estimate of the field map based on simplifying assumptions.
Whereas conventional EPI-based neuroimaging (in particular, routine fMRI) is often performed at relatively low resolution and involves spatial smoothing, EPI resolution is critical in specific applications including laminar fMRI, 6,8 fast quantitative susceptibility mapping, 62,63 or high-resolution dw-MRI. 60,64 Such areas of increasing current interest require high-quality distortion correction. Further experimental evaluation of potential benefits from additional deblurring is therefore warranted.

| CONCLUSIONS
We have shown that distortion correction of EPI acquisitions is feasible by deconvolution in image space following proper regularization. The analytic approach intrinsically corrects T * 2 -related blurring and intensity artifacts in addition to geometrical distortions. The method is applicable in combination with different field-mapping techniques, namely ME-FLASH or TOPUP, with overall improved image quality. F I G U R E 1 0 Application of DecoDisCo to diffusion-weighted spin-echo EPI. Images were acquired with opposite phase-blip polarities using anterior-posterior (AP; first column) and posterior-anterior (PA; second column) PE direction without (b = 0 s/mm 2 ; top panel) and with diffusion weighting (b = 1000 s/mm 2 ; bottom panel). The field map was estimated with TOPUP. For both b-values, successful correction of geometric distortions is evident from the increased similarity of images acquired with the different PE directions, although residual differences remain in heavily distorted regions (in particular, in frontal regions affected by signal compression). The weighted combination (c = −2) of both acquisitions (third column) suggests a superior image quality in these regions. Results obtained with TOPUP (fourth column) are also shown for comparison