Using synthetic MR images for distortion correction

2.1. Creating an undistorted synthetic EPI image using Synth

Our approach to correcting EPI distortion is to create a synthetic EPI image based on information combined from a participant’s undistorted T1w and T2w images and then use this synthetic image as a reference for nonlinearly aligning a participants’ real EPI image. Here, we refer to any real image that Synth attempts to synthesize as a target image. The ideal synthetic EPI image will have three properties: (1) it will match the real EPI target image in terms of overall signal intensity; (2) it will exhibit similar contrast between fluids and tissue types that closely correspond to what is observed in the EPI target image and; (3) It will account for the difference in spatial resolution between typically high-resolution anatomical images and typically low-resolution EPI images. Thus, the synthetic image can be represented by the following model:

The left hand side of Eq. (1), $\mathbf{s}$, represents a final synthetic image. Working from inside to outside, the right hand side of Eq. (1) represents the following: the matrix, $\mathbf{F}$, is the decomposition of pre-aligned T1w and T2w anatomical images into a set of basis vectors (e.g., radial basis functions — as in Fig. 1a, or b-splines, etc.) in order to model the continuous relationship between the voxel intensities of the anatomical images and EPI target image; $\mathit{\theta }$ comprises the weights for each column of $\mathbf{F}$. $\mathbf{E}$ is a blurring operator modeling the effective difference in spatial resolution between the anatomical images and EPI target image. In our implementation of Synth, the blurring operator, $\mathbf{E}$, is effected with an Epanachnikov smoothing kernel, chosen for its optimal noise-resolution properties (Gureyev et al., 2020). Finally, $\mathbf{B}$ is a global image contrast operator modeled as a cumulative beta distribution, parameterized by scalars $\alpha$ and $\beta$. When $\mathit{\theta }$, $\alpha$, and $\beta$ are at the optimum solution, $\mathbf{s}$ represents a synthetic image with the geometry of the undistorted source images and the contrast properties of the EPI target image. Solving for the optimal $\mathit{\theta }$, $\alpha$, and $\beta$ requires solving a joint optimization problem, which is discussed in the next section.

The columns of the matrix, $\mathbf{F}$, are composed of a radial basis function (RBF) decomposition of the source T1w and T2w anatomical images. This decomposition is visualized in Fig. 1a, where the T1w and T2w input images have been divided into smooth voxel intensity ‘bins’, each of which corresponds to an RBF component. Each of these components are weighted by the parameter, $\mathit{\theta }$, such that the weighted sum of these components replicates the contrast properties of the EPI target image as visualized in Fig. 1b. Synth allows for a user-defined number of intensity ‘bins’ into which the T1w/T2w images are decomposed. The choice of number of intensity bins is dictated by a tradeoff between computational complexity and the accuracy of the synthetic image. As the number of RBF components in $\mathbf{F}$ increases, so too do the memory requirements. For the synthetic images in the presented results, a 24 component RBF decomposition of both T1w/T2w images (12 components of T1w; 12 components of T2w) at $1\mathrm{mm}$ isotropic resolution consumed $\sim 20$ GB of RAM.

2.2. Estimation of synthetic image and distortion-correction warp

To relate the synthetic EPI image to the real EPI target image, we model the distortion of the synthetic image as a non-linear warp, $f$, constrained to operate in the phase-encoding direction. This model is described by:

Here, $\mathbf{y}$ represents the target image (e.g., a distorted EPI image), $\mathit{\varphi }$ represents the underlying parameterization of the nonlinear transformation, $f$, and $\mathit{\eta }$ represents additive gaussian noise. By combining Eqs. (1), (2), and solving for the parameters that maximize the correlation coefficient between the real and distorted synthetic image, we obtain the following joint optimization problem:

To solve Eq. (3), we employ an alternating minimization approach, where one parameter is optimized while the others are held fixed over each iteration. $\mathit{\varphi }$ is initialized with a rigid-body transform so that both synthetic and real EPI images are closely registered. Once an updated value for $\mathit{\varphi }$ is found through non-linear registration, the inverse of the estimated nonlinear warp, $f^{-1}$, is then computed. $f^{-1}$ is analogous to a traditional field map correction and maps voxels in a distorted EPI image, $\mathbf{y}$, back to their correct locations. Using the corrected target image, $f^{-1}\left(\mathbf{y}\right)$, an updated value for $\mathit{\theta }$ can be found by solving the linear system $f^{-1}\left(\mathbf{y}\right)=\mathbf{E}\mathbf{F}\mathit{\theta }$. This updated $\mathit{\theta }$ is used to generate an updated intermediary synthetic image. Finally, the parameters, $\alpha$ and $\beta$, that control the global contrast for the intermediary synthetic image, are computed by minimizing the least squares difference between the synthetic image and distortion corrected EPI image. The resulting synthetic image is then used for another iteration of non-linear registration.

In regions of the brain with significant distortion and signal dropout or pile-up affecting the EPI image, there exists no valid mapping between a participant’s anatomical and EPI images. Therefore Synth also allows for the inclusion of a weight volume to reduce the contributions of these areas when estimating the parameters, $\mathit{\theta }$, $\alpha$ and $\beta$. For the presented results, we down-weighted the contributions of voxels in high-distortion areas (Supplemental Figure 3).

With the application of $f^{-1}$ to the real EPI data during each iteration, geometric correspondences between anatomical and functional images are improved. Because the quality and accuracy of the intermediary synthetic EPI images, depends on how closely registered the anatomical and true EPI images are, modest improvements to the synthetic image can be achieved after each iteration of the Synth algorithm. By default, Synth repeats this process for 3 iterations. The full procedure for solving Eq. (3) is outlined in Algorithm 1 (see Section 6.8 in Supplemental Methods) and Fig. 1c.

In principle, many of the pre-existing nonlinear registration utilities are suitable for estimating $f$ and $f^{-1}$, owing either to their methods of diffeomorphic warp construction which guarantee the existence of an invertible warp —as is the case with AFNI’s 3dQwarp (Zhark the Grotesquely Warped (but still strangely handsome), 2021) and ANTs SyN (Avants et al., 2008)— or their ability to project a potentially non-invertible warp onto a “nearest invertible” space— as is the case with FSL’s FNIRT (Paul, 2021). For the results presented here and in the reference preprocessing scripts associated with this manuscript, we used the ANTs SyN algorithm with a local cross correlation metric for estimating all non-linear warps (see Appendix A.7 in Supplemental Methods). ANTs SyN was chosen for its reliable high performance in non-linear warp estimation (Avants et al., 2011, Klein et al., 2009, Wang et al., 2017).

2.3. Description of MRI datasets and processing

We assessed the quality of distortion corrections produced by direct mapping to Synth-generated synthetic images in two primary contexts representing important use cases for fMRI data. First, we examined a subset of 100 participants selected randomly from the Adolescent Brain Cognitive Development (ABCD) study dataset (median age: 9.85 years; min: 9 years; max: 11 years). Our random sample used multi-band data acquired using GE, Philips, and Siemens scanners (see Table 1 in Supplemental Material). Second, we evaluated Synth’s performance on the Midnight Scan Club (MSC) dataset, which consists of resting state fMRI scans acquired from 10 participants on 10 separate occasions (300 min of resting state fMRI data/participant). The MSC precision functional mapping (PFM) (Braga and Buckner, 2017, Gordon et al., 2017, Gordon et al., 2018, Gordon et al., 2020, Gratton et al., 2018, Greene et al., 2020, Laumann et al., 2015, Lynch et al., 2020, Marek et al., 2018, Newbold et al., 2020b, Newbold et al., 2020a, Sylvester et al., 2020, Zheng et al., 2020) dataset allowed us to assess session-to-session reliability of EPI distortion correction schemes across multiple sessions for the same participant. Importantly, image acquisition parameters for ABCD and MSC datasets differ significantly allowing us to assess the performance of Synth distortion correction on images with a range of contrast and levels of detail. Image acquisition parameters have been reported in detail elsewhere (Casey et al., 2018, Gordon et al., 2017).

In the ABCD and MSC datasets, we compared the effectiveness of this approach against five widely used distortion correction algorithms. Three of these approaches rely on separately acquired field map data (FSL fugue; FSL topup; and AFNI’s 3dQwarp) while the fourth and fifth approaches implement alternative field map-less distortion correction methods, SyN-SDC (Esteban et al., 2019, Esteban et al., 2021) and SynB0-DisCo. FSL field map corrections were constructed using either fugue (used for double echo field maps acquired in the MSC dataset) or topup (used for opposite direction phase encoded field maps in the ABCD dataset). In AFNI (Cox, 1996), the distortion correction is estimated from opposite phase encoded field maps using a “meet in the middle” nonlinear warp estimation implemented in 3dQwarp (Zhark the Grotesquely Warped (but still strangely handsome), 2021). This estimated warp is used to correct the image and combined with a separately estimated rigid body transform to align the functional to the anatomical image. Because 3dQwarp does not operate on double echo field map data, we could not assess its performance on the MSC data set. We therefore only compared Synth’s performance to ${\text{FSL}}_{fugue}$ and SyN SDC, SynB0-DisCo.

For field map-less approaches, we assessed the SyN-SDC field map-less method that corrects distortion by non-linearly aligning the participant’s EPI and anatomical images while constraining allowable warps to regions known to be strongly affected by distortion (Huntenburg, 2014, Treiber et al., 2016, Wang et al., 2017). Additionally, we evaluated the SynB0-DisCo method, which uses deep learning methods to synthesize a pseudo infinite bandwidth EPI image as an input to FSL’s topup to estimate the distortion correcting warp (Schilling et al., 2020).

Performance of the different distortion correction approaches was assessed in matched datasets, thus all statistical comparisons were performed within a multi-level modeling framework (fitlme, MATLAB). For these analyses, the metrics produced by the Synth-based registration pipeline were modeled as the baseline and differences in performance metrics associated with other approaches were modeled with main effect factors. To account for participant-specific variability, independent of the registration approach, participant identity was modeled as a random effect.

Studies in which the phase encoding direction is chosen to align with the anterior–posterior axis represent the most common use case. However, researchers have elected in some instances to choose other phase encoding directions for a variety of theoretical and practical reasons, for instance in an attempt to minimize total distortion or to preserve signal in orbitofrontal regions. The Human Connectome Project (HCP), in which the phase encoding direction was chosen to align with the left–right axis, represents a notable example of such a study (Van Essen et al., 2013, Glasser et al., 2016). While this approach is decidedly less common, it does represent an important use case. We therefore applied Synth to a 10 subject data set chosen at random from the publicly available HCP data set. We compared corrections produced by Synth against those produced by FSL’s topup. Instances in which HCP results differ in character from those observed in the ABCD and MSC data set are noted in the relevant portion of this manuscript and are otherwise included in the supplement.

2.4. Evaluation metrics

No consensus exists as to what global metric best quantifies alignment quality. Many metrics used to measure alignment quality tend to bear some relation to the cost functions that are used during registration and therefore directly introduce a risk of confounding circularity. Thus, we computed several metrics comparing the similarity of different features of the field map corrected EPI images to their associated anatomical images.

3.1. Synth images consistently match the contrast of EPI images

We reasoned that anatomically-based reference images with greater contrast similarity to EPI images would serve as a more reliable reference for distortion correction than images with lower contrast similarity. To assess differences in contrast similarity between EPI and T1w, T2w and synthetic images, we calculated the linear correlation coefficient between participants’ EPI images and each of their T1w, T2w and synthetic images within a whole brain mask. By default, Synth uses a model that includes a 12-component radial basis function (RBF) decomposition of the T1w/T2w source images along with pair-wise T1w/T2w interaction terms to create synthetic EPI images (Fig. 1). This model was sufficient to produce synthetic images with comparable EPI contrast similarity in both ABCD and MSC datasets (Fig. 2a,b). In both datasets, synthetic images exhibited significantly increased contrast similarity to EPI images compared to T1w and T2w images (Fig. 2c). Contrast similarity between the EPI and Synth images was consistently highest. Mixed-effects model comparisons in the ABCD data set revealed that the differences between Synth contrast similarity and T1w/T2w contrast similarity were both significant: (T1w-Synth$=-0.417$; $p<0.001$; $t=-46.18$; $df=297$) and (T2w-Synth$=-0.18$; $p<0.001$; $t=-19.62$; $df=297$). The same pattern held in the MSC dataset: (T1w-Synth$=-0.374$; $p<0.001$; $t=-49.64$; $df=297$) and (T2w-Synth$=-0.274$; $p<0.001$; $t=-36.4$; $df=297$). Synthetic EPI images derived from the HCP data set exhibited greater contrast similarity to their respective EPI reference images than their associated T1w/T2w images (Supplemental Figure 4)

3.2. Synth outperforms existing field map-less methods on global measures of image alignment

We assessed Synth’s distortion correction performance (Fig. 3) against existing field map-based (AFNI 3dQWarp; FSL topup; FSL fugue) and field map-less (SyN-SDC; SynB0-DisCo) distortion correction methods using both established and novel metrics that quantify global image similarity between the EPI image and the associated T1w/T2w anatomical images (Fig. 4). These metrics included normalized mutual information; the correlation coefficient between the image gradients; and segmentation alignment. Here, we outline Synth registration performance for each metric described above as compared to the highest performing competitor. Synth registration tended to produce EPI images with the highest registration quality metrics for both ABCD and MSC data. Comparisons are reported in the form of statistical contrasts corresponding to the mean difference between registration metrics (e.g., Method A-Method B). Full statistical tables comparing all methods are included in section 5.4 and 5.5 of the Supplemental Materials. Illustrative comparisons between FSL fugue/topup and Synth are included in Supplemental Figures 1 and 2.

First, we examined a well-established image similarity metric, normalized mutual information (NMI). In the ABCD dataset, NMI between EPI and T1w images was not significantly different (${\text{FSL}}_{topup}$-Synth$=9.25\text{e-6}$; $p=0.98$; $t=0.025$; $df=495$). In contrast, NMI for T2w images was significantly higher for Synth images than those aligned by ${\text{FSL}}_{topup}$ (${\text{FSL}}_{topup}$-Synth$=-1.49\text{e-3}$; $p<0.001$; $t=-2.85$; $df=495$). Synth alignment produced EPI images with the highest NMI shared between T1w/T2w images in the MSC dataset as well (${\text{FSL}}_{fugue}$-Synth$=-4.02\text{e-3}$; $p<0.001$; $t=-10.74$; $df=396$) and (${\text{FSL}}_{fugue}$-Synth$=-3.22\text{e-3}$; $p<0.001$; $t=-8.53$; $df=396$) (Fig. 4a). We observed no significant differences in these alignment metrics in the HCP data set (Supplemental Figure 5, Supplemental Section 5.6).

Next, we examined the quality of alignments emphasizing registration of high contrast boundaries. EPI-T1w edge alignments produced by Synth were outperformed by those produced by SyN SDC (SyN SDC-Synth$=3.19\text{e-2}$; $p<0.001$; $t=4.93$; $df=495$) in the ABCD dataset; for T2w gradient magnitude images, there was no significant difference between Synth and the highest performing competitor, ${\text{FSL}}_{topup}$ (${\text{FSL}}_{topup}$-Synth$=-3.12\text{e-3}$; $p<0.001$; $t=-0.54$; $df=495$). In the MSC dataset, its closest competitor varied by modality. For EPI-T1w edge alignment, ${\text{FSL}}_{fugue}$ was the highest performing competitor, but still produced edge alignments that were significantly lower than those produced by Synth distortion correction (${\text{FSL}}_{fugue}$-Synth$=-1.70\text{e-2}$; $p<0.001$; $t=-3.22$; $df=396$). For T2w images in the MSC dataset, ${\text{FSL}}_{fugue}$ was once again the closest competitor but did not perform significantly better than Synth, (${\text{FSL}}_{fugue}$-Synth$=4.88\text{e-3}$; $p=0.39$; $t=0.87$; $df=396$) (Fig. 4b). In the HCP data set, these metrics revealed no significant differences between Synth and ${\text{FSL}}_{topup}$.

In addition, we examined the quality of the alignment of EPI images to anatomical images segmented by tissue type ($gw$: gray/white matter; $vw$: ventricles/white matter; $ie$: interior/exterior of brain) based on the participants’ freesurfer segmentation. These metrics quantify the separability of EPI image voxel intensities based on anatomically defined segmentation (higher AUC [area under the receiver operating characteristic curve] is better). ${\text{AUC}}_{gw}$ image segmentation alignment metrics for Synth were significantly lower than those produced by ${\text{FSL}}_{topup}$ in the ABCD dataset (${\text{FSL}}_{topup}$-Synth$=1.62\text{e-2}$; $p<0.001$; $t=5.86$; $df=495$) and MSC (${\text{FSL}}_{topup}$-Synth$=1.31\text{e-2}$; $p<0.001$; $t=6.03$; $df=396$) datasets. Synth-aligned images exhibited greater ${\text{AUC}}_{ie}$ in both ABCD (${\text{FSL}}_{topup}$-Synth$=-1.60\text{e-2}$; $p<0.001$; $t=-5.60$; $df=495$) and MSC (${\text{FSL}}_{fugue}$-Synth$=-2.16\text{e-2}$; $p<0.001$; $t=-8.82$; $df=396$) datasets. We also found that the SyN SDC pipeline produced the greatest ${\text{AUC}}_{vw}$ values (SyN SDC-Synth$=2.14\text{e-2}$; $p<0.001$; $t=3.77$; $df=495$) for the ABCD dataset. However, Synth performed better than the highest performing competitor, ${\text{FSL}}_{fugue}$ (${\text{FSL}}_{fugue}$-Synth$=-1.02\text{e-2}$; $p=0.03$; $t=-2.21$; $df=396$) in the MSC dataset (Fig. 4c). In the HCP data set, ${\text{FSL}}_{topup}$ produced significantly greater ${\text{AUC}}_{gw}$ values (${\text{FSL}}_{topup}$-Synth=1$.63\text{e-2}$; $p<0.001$; $t=-4.52$; $df=18$) whereas corrections produced by Synth exhibited greater ${\text{AUC}}_{vw}$ values (${\text{FSL}}_{topup}$-Synth=$-1.29\text{e-2}$; $p=0.001$; $t=-2.82$; $df=18$)

3.3. Synth improves local measures of anatomical and functional image alignment

Existing experimental evidence indicates that field map correction quality varies regionally and depends on the method used to estimate the B0 inhomogeneity (Schallmo et al., 2021). Consequently, to identify whether Synth was better at correcting certain brain regions, we carried out a winner-take-all (WTA) analysis across all voxels by determining which method produced the greatest average local similarity between T1w and T2w anatomical images ($R^2$ computed within a 7x7x7 voxel spotlight). Representative spotlight $R^2$ images underlying the WTA analysis are provided in Supplemental Figures 6 and 7. Each distortion correction method exhibited regions of the brain for which they performed consistently better than other distortion correction methods. Although the structure of the WTA maps varied between datasets, they shared some consistent features (Fig. 5).

In ABCD and MSC datasets, Synth produced either the largest, or second largest regions of high quality local alignment, regardless of anatomical image modality (i.e., T1w, T2w). In the MSC data set, Synth generally had the best regional correction performance across the entire brain. For the ABCD dataset, the best regional distortion correction performance was split between ${\text{FSL}}_{topup}$ and Synth, where ${\text{FSL}}_{topup}$ provided better local alignment to T1w images over a large swathe of the superior frontal cortex while Synth’s particular strengths were in improving local registration of cerebellar, brainstem, orbitofrontal and occipital regions with respect to a participant’s T2w image. Overall, ${\text{FSL}}_{topup}$ won the greatest percentage of voxels in the T1w and T2w comparisons for the ABCD dataset, (T1w: ${\text{FSL}}_{topup}$ $=43\text{\%}$; Synth$=33\text{\%}$; T2w: ${\text{FSL}}_{topup}$ $=51\text{\%}$; Synth$=46\text{\%}$). While Synth won the greatest percentage of voxels in the T1w and T2w comparisons over ${\text{FSL}}_{fugue}$ (T1w: ${\text{FSL}}_{fugue}$ $=21\text{\%}$; Synth$=65\text{\%}$; T2w: ${\text{FSL}}_{fugue}$ $=31\text{\%}$; Synth$=54\text{\%}$) for the MSC dataset. In the HCP data set, ${\text{FSL}}_{topup}$ produced better measures of local alignment over a modestly larger fraction of the image (T1w: ${\text{FSL}}_{topup}=55\text{\%}$; Synth=45%; T2w: ${\text{FSL}}_{topup}=57\text{\%}$; Synth=43%) (See Supplemental Figures 8, 9, and 10). These results indicate that Synth can perform at parity with state-of-the-art field map-based approaches.

The analysis of functional images depends to a greater extent on the quality of gray matter alignment and less so on the alignment quality of other voxels representing other tissue types. Accordingly, we summarized local gray matter alignment performance of each distortion correction method by computing the average spotlight $R^2$ across only the gray matter voxels. We observed that Synth and ${\text{FSL}}_{topup}$ registration pipelines produced similar quality of local registration to the participants’ T1w image (${\text{FSL}}_{topup}$-Synth$=-2.51\text{e-3};p=0.17;t=-1.38;df=495$), in the ABCD dataset. Synth performance did not differ significantly from the highest performing competitor, ${\text{FSL}}_{topup}$, for average local similarity between EPI and T2w images (${\text{FSL}}_{topup}$-Synth$=-4.14\text{e-3};p=0.21;t=-1.27;df=495$). Synth registration produced EPI images with the greatest average local similarity for associated T1w images (${\text{FSL}}_{fugue}$-Synth$=-1.93\text{e-2};$p$<$ 0.001$;t=-8.64;df=396$) and did not differ significantly from the highest performing competitor, ${\text{FSL}}_{fugue}$ for T2w images (${\text{FSL}}_{fugue}$-Synth$=-3.27\text{e-3};p=0.17;t=-1.37;df=396$) in the MSC dataset (Fig. 6).

3.4. Synth improves intra-participant consistency across sessions

The primary goal of applying distortion correction to EPI images is to improve data quality by reducing variability caused by poor alignment. It is therefore important to know the extent to which a particular approach to distortion correction introduces additional measurement variability. By applying different distortion correction strategies to the repeated-measures in the MSC dataset, we assessed how each method contributed to within-participant structural and functional variability.

First, we examined each voxel’s average intensity across all sessions normalized by its standard deviation across sessions. We term this stability metric, session-to-session SNR (s-SNR). We observed that for Synth-aligned images, s-SNR was not significantly different than that produced by the highest performing competitor, ${\text{FSL}}_{fugue}$ (${\text{FSL}}_{fugue}$-Synth$=0.80$; $p=0.55$; $t=0.60$; $df=36$) (Fig. 7a).

Next, we quantified global similarity of time averaged EPI images for each pair of sessions contributed by a participant. We extracted the region of the aligned time-averaged EPI volumes containing just brain matter (as classified by the FreeSurfer segmentation of the participant’s anatomical images) from each session. Then, for each participant, we computed the average linear correlation coefficient across all pairs of sessions. Owing to the general stability of a given participant’s brain geometry on the time scales over which the MSC dataset was collected, inter-session linear correlation coefficients between scan sessions tended to be quite high across all methods ($0.94-0.99$). While Synth produced EPI images with the greatest session-to-session similarity, it did not produce images that were significantly more stable than those produced by ${\text{FSL}}_{fugue}$ (${\text{FSL}}_{fugue}$-Synth$=-5.56\text{e-3};p=0.23;t=-1.20;df=1796$) (Fig. 7b).

Finally, we examined how the choice of distortion correction approach affects the stability of resting state functional connectivity (RSFC) data. To do this, we created voxel-level RSFC matrices for each session of data and computed RSFC similarity between pairs of sessions as the correlation between upper triangular portions of the RSFC matrices. Mixed effects analysis of similarity between unique session pairs revealed RSFC correlation matrices from Synth aligned datasets were not significantly less stable than those produced by the highest performing competitor, ${\text{FSL}}_{fugue}$ (${\text{FSL}}_{fugue}$-Synth$=5.47\text{e-3}$; $p=0.12$; $t=1.55$; $df=1796$) (Fig. 7c).

4.1. Synth achieves parity with or exceeds the distortion correction performance of alternative methods on most measures of global alignment

We began our exploration of distortion correction performance by examining the commonly used global image similarity metric, normalized mutual information. We found that Synth performed comparably to or better than state of the art field map corrections provided by FSL’s topup and fugue in both ABCD and MSC data sets. One of the challenges to assessing the quality of multi-modal image registration and distortion correction is determining which global metric serves as the optimal summary of final alignment quality. To address this concern, we assessed the performance of each distortion correction approach with several additional global image similarity metrics, each intended to assess alignment with respect to different image features.

The first alternative global similarity metric we examined was the quality of high-contrast edge alignment. We reasoned that these features would be less influenced by residual low spatial frequency bias fields and the differences in tissue contrast properties between images acquired with different modalities (T1w, T2w, EPI). Here, Synth’s performance on cross-modal edge alignment metrics was generally high, either achieving the highest performance or placing in the top two. The single exception to this outcome was observed in the T1w-EPI edge alignment metrics in the ABCD data set. There, FSL’s topup and SyN SDC produced the highest alignment metrics. This result may arise from the fact that the Synth pipeline initially registers the EPI image into the anatomical image space using a T2w reference image, while the FSL and SyN SDC pipelines both register to the anatomical image space using the T1w image as a reference. It may be the case that the choice of image modality (i.e. T1w/T2w) used for the initial EPI-anatomical alignment biases the edge alignment metric for those two image modalities. Significantly, we observed that Synth’s high performance generalized well to the HCP data set in which, due to the left–right orientation of the phase encoding direction, EPI distortion exhibits very different structure.

A second notable exception to Synth’s typically high performance was the outcome of the segmentation alignment metric focusing on the alignment accuracy of the ventricles and white matter (${\text{AUC}}_{vw}$) in the ABCD dataset. In that instance, the field map-less SyN SDC method provided the best ventricles and white matter alignment metrics (${\text{AUC}}_{vw}$) in the ABCD dataset followed by FSL’s topup. One potential reason for this outcome is that both FSL’s topup and SyN SDC pipelines used brain-based registration for their EPI to anatomical alignments (Greve and Fischl, 2009). This registration technique incorporates anatomical segmentation information when guiding anatomical alignments and may improve ventricle/white matter correspondence between the EPI and anatomical images. This interpretation is complicated somewhat by the fact that ventricular alignment metrics were greatest for Synth-based alignments in the HCP data set, while FSL’s topup tended to produce better gray/white matter alignment.

4.2. Synth exhibits high performance in regional measures of image alignment

In addition to the challenge of selecting and justifying a particular global metric for measuring multi-modal image alignment quality, prior work has shown it does not always provide a sufficient condition for demonstrating optimal alignment (Greve and Fischl, 2009). This is because global metrics, while flexible, often cannot fully capture the complicated mapping of voxel intensity values between images acquired with different modalities. While achieving high measures of global image similarity is a necessary condition for quality image alignment, it is important to consider them alongside non-global measures to fully assess alignment (Saad et al., 2009). For this reason, we examined small regions of the image (i.e. spotlight), which contain a limited range of tissue types and minimal residual bias field. In this way, we created an alignment metric that assesses image similarity while reducing the influence of non-linearities in the relationship between voxel intensity values of images acquired with different modalities.

Our winner-take-all analyses of local image alignment revealed Synth to be among the top performing distortion correction approaches, producing the highest local image quality metrics in large swathes of the cortex, but also performing particularly well in correcting distortions affecting the cerebellum and brain stem. That different distortion correction approaches excelled in correcting different areas of the brain is in line with existing evidence suggesting that the accuracy of standard field map distortion correction methods varies across the brain (Schallmo et al., 2021, Graham et al., 2017).

One possible reason for Synth’s superior performance in the MSC dataset overall, compared to the ABCD dataset, may be due to the dual echo field maps used in the MSC dataset. Prior work has shown that distortion correction with reverse-phase encoding field maps outperform corrections produced by dual echo field maps (Graham et al., 2017). Our findings are consistent with these observations, and suggest that Synth may provide superior correction to dual-echo field map corrections in datasets in which image acquisition parameters share sufficient similarity to those used for the MSC study (Gordon et al., 2017).

4.3. Synth does not introduce session-to-session variability like existing field map-less approaches

We observed that as a general rule, field map-less approaches tended to introduce more session-to-session structural variability than traditional field map-based approaches. Synth stood out as a clear exception to this trend by introducing no more session-to-session variability than a simple rigid-body alignment procedure or field map correction. Importantly, because Synth distortion correction did not measurably increase the session-to-session variability of RSFC matrices we can conclude that the improvement in EPI-to-anatomical alignment provided by Synth does not come at a hidden cost of increasing the variability of the underlying functional data.

4.4. Accounting for differing performance of existing field map-less approaches

Synth tended to produce higher quality and more reliable corrections than existing field map-less approaches (SyN SDC and SynB0-DisCo) with which it shares many conceptual similarities. This prompts us to consider the potential reasons why Synth performs with greater reliability and accuracy.

In the case of SyN SDC, one reason for decreased reliability may be the use of a group average field map template in order to constrain its solutions. Such a template may not allow sufficient flexibility to align high spatial resolution participant-specific features. SyN SDC also relies on an intermediate atlas alignment stage, which may be a source of additional variability. Because participant motion can influence the structure of the EPI distortion, the use of a field map template that does not account for this effect may reduce overall efficacy.

Both Synth and SyN SDC employ a direct mapping approach to correct EPI distortion. The fundamental difference in their approaches is that SyN SDC attempts to directly align to a participants’ T1w image –whose contrast properties are very different from the EPI image– while Synth uses an intermediate synthetic image whose contrast properties are more similar to the EPI. That Synth proved more effective and reliable than the similar approach implemented in SyN SDC would seem to indicate the importance of matching the contrast properties of target and reference images when attempting to use direct mapping approaches for distortion correction.

With SynB0-DisCo, a synthetic pseudo-infinite bandwidth intermediate EPI image facilitates distortion correction using FSL’s topup. The corrections, though, were less accurate and less reliable than those produced by Synth and SyN SDC. On the surface, it would seem that the use of a contrast matched intermediate image along with FSL’s state of the art field map estimation software should produce the same reliable distortion correction observed with Synth. Two factors may potentially contribute to the variable performance of this approach. The tissue contrast properties of EPI data can vary significantly depending on image acquisition parameters. Therefore, one possibility is that SynB0-DisCo’s deep learning model may generalize poorly to unseen datasets whose resolution and contrast are quite different from its training data. A second possibility is that synthetic images produced by deep learning models are prone to introducing spurious artifacts into their outputs (Antun et al., 2020, Bhadra et al., 2021). The presence of such artifacts may be a significant source of registration variability. Both of these effects may have contributed to reducing SynB0-DisCo’s performance in the ABCD and MSC datasets. Synth’s ability to generate effective synthetic images based solely on the participant’s anatomical data without relying on a large training dataset avoids both of these potential issues by flexibly adapting to the contrast properties of a specific acquisition and minimizing the possibility of spurious artifacts.

While our results indicate that Synth is able to correct distortion comparably well to current field map based techniques, it is possible that advances in theory or software may improve field map-based corrections even further so that they reliably exceed Synth’s performance. We therefore do not advocate abandoning the collection of field map data entirely. Rather we propose that Synth distortion correction is a reliable substitute for field maps at the present time and given the current state of the art. We can unreservedly recommend using Synth as an effective tool for correcting distorted fMRI data when no field map data is available.

4.5. General guidelines for deploying Synth

Our results indicate that synthetic images with contrast similarity (linear correlation coefficient) of 0.4 with the target EPI image are sufficient for producing high quality corrections. In order to accommodate datasets with different contrast, resolution, and fields-of-view, Synth provides a variety of command line options: $source\text{\_}images$, is a list of pre-aligned undistorted anatomical images from which the target image is modeled. Typically this is the participant’s T1w and T2w images, but any undistorted images can be used in principle; $interaction\text{\_}model$ is a string that specifies both main-effect terms and interaction terms between a user-specified number of source images, each decomposable into a user-specified number of RBF components. The number of RBF components applied to each source image enables more complex models that will allow Synth to capture more complicated relationships between the voxel intensities of the source images and the EPI images. This flexibility comes at the cost of increased memory requirements and computational time, but can produce more accurate synthetic images that may ultimately improve registration quality. In our results, the RBF model parameters were chosen to produce a synthetic image of qualitatively sufficient visual similarity to their associated EPI images. By default, Synth uses the RBF model parameters optimized for the ABCD and MSC datasets (i.e., a 12 component decomposition on each of the participants’ T1w/T2w images, along with pairwise interaction terms between T1w/T2w components). In general, Synth attempts to map the contrast of high resolution source images to a low resolution EPI image, which requires modeling the difference in point spread function between the two. Modifying the bandwidth parameter enables Synth to be applied to datasets acquired at different resolutions. Users acquiring high resolution EPI images will decrease this parameter and users acquiring low resolution EPI images will increase it. Lastly, the $weight\text{\_}mask$ parameter is a binary mask that indicates the region of the EPI image in which a valid mapping between source image (T1w/T2w) intensities and target image (EPI) intensity values exists. This mask will exclude areas of extreme signal dropout. Only voxels within this mask are included when estimating the model that produces the synthetic image.

We provide a reference pipeline implementing Synth which only requires the user to organize their data into the BIDS format standard (Gorgolewski et al., 2016). Using this pipeline, all of Synth’s parameters will be automatically determined. Users seeking to integrate Synth into their own pipelines, will need to construct their own weight masks and determine the appropriate Epanachnikov smoothing kernel bandwidth independently. In general, we recommend that inputs to Synth should be deobliqued and bias field corrected, but otherwise minimally preprocessed. For walkthroughs detailing the Synth’s installation and use as a standalone utility or as part of our BIDS compliant reference pipeline, users are invited to explore https://gitlab.com/vanandrew/omni.

4.6. Future directions for field map-less distortion correction

Although the distortion corrections produced by Synth are high quality, several open questions remain relating to how this general approach can be implemented most effectively. Chief among these is whether Synth’s performance can be augmented further by including existing field map information derived using traditional approaches. The focus of the present report is assessing the quality of the field map corrections produced by Synth in a situation in which there is no initial estimate of the EPI image distortion provided to the underlying nonlinear warping software. Including an existing field map as an initial estimate may further increase Synth’s performance. Future iterations of Synth will include this ability and may produce even greater correction fidelity between functional and anatomical images than presented here.

Even though Synth’s underlying RBF model provides a great degree of flexibility in fitting the functions that map T1w/T2w voxel intensities to EPI image voxel intensities, other data collection parameters and preprocessing steps may improve the quality of the synthetic images or reduce the needed complexity of the RBF model. For example, adjusting echo times during data acquisition may produce images that can be more accurately modeled with fewer RBF components while still retaining sufficient sensitivity to BOLD contrast. Alternatively or in addition, performing contrast enhancing preprocessing of the anatomical and EPI images so that the resulting function relating their voxel intensities can be modeled more efficiently by the chosen RBF model may allow for simpler RBF models to perform equivalently. In our own preliminary tests and for the results reported here, we observed that synthetic images were more accurate when a bias field was estimated and removed from the anatomical images and image contrast was increased through histogram normalization (see Appendix A.2.2 in Supplemental Methods).

4.7. Conclusion

The results reported here have demonstrated that it is possible to achieve high-quality EPI distortion correction without the need to collect separate field map data. We have shown that field map-less approaches, such as Synth, can perform comparably to existing gold-standard distortion correction approaches, and in some cases, may surpass them on measures of global and local image alignment quality. Importantly, Synth produced high quality alignments even within the HCP dataset, indicating that its high performance generalizes not only across differences in EPI resolution and contrast properties, but is robust to differences in phase encoding direction. Removing the reliance on field maps to correct EPI distortion will allow researchers to recover samples with missing or corrupted field maps while maintaining high quality alignment between their anatomical and EPI images. Field map-less distortion correction may prove to be a particular asset to researchers studying high motion cohorts, such as pediatric or neuropsychiatric populations (Greene et al., 2016), where acquiring high quality field map data is a significant challenge. Reliable field map-less distortion correction shows great promise for overcoming the limitations arising from the acquisition, processing, and quality checking of field map data and has the potential to greatly simplify data processing in the neuroimaging field.

A.1. Description of common pipeline

Distortion correction methods were compared using a common registration pipeline, which differed only in the distortion correction method and functional-to-T1w alignment. This includes alignment of anatomical images to atlas, functional framewise alignment and resampling of the final aligned functional data to $3\mathrm{mm}$. Details of each common step are explained in the following sections.

A.2. HCP dataset

Ten randomly selected participants from the Human Connectome Project data set were used for the supplemental analyses detailed in this manuscript. EPI data was drawn from the “Emotional Processing” portion of the dataset, along with each participants’ T1w/T2w images. A thorough account of the acquisition details for task and anatomical imaging data are located in Barch et al. (2013) and Van Essen et al. (2013). Briefly, data was acquired using a modified SIEMENS Skyra using a TR of 720 ms, TE 3.31 ms, flip angle of 52 degrees and multiband acceleration factor of 8. EPI images were acquired with 2.0 mm isotropic voxels.

A.3. Description of radial basis function model used to estimate the synthetic EPI image

Synth allows the user to specify a regression model that maps source image intensities to target image intensities. The radial basis components were constructed using the following procedure: Let ${\mathbf{I}}_m$ refer to an individual source image, e.g.,T1w or T2w; and $\left\{c_1,c_2,\dots ,c_J\in C\right\}$ refers to a set of $J$ evenly spaced values that lie between the maximum and minimum values of ${\mathbf{I}}_m$. Then each radial basis component, ${\mathbf{r}}_j$, of an image, ${\mathbf{I}}_m$, is constructed according to the following equation:

A visual representation of the images produced by applying Eq. (4) to T1w and T2w images is depicted in Fig. 1a. The final design matrix fit by Synth consists of each radial basis component concatenated into a single matrix $\mathbf{F}$, along with intercept terms, ${\mathbf{e}}_m$, (a vector of all ones) for each source image, such that:

For both datasets, $\mathbf{F}$ was constructed from a 12 component decomposition on each of the participants’ T1w/T2w images, along with all T1w/T2w pairwise interaction terms. As in traditional linear models, these interaction terms are simply the pairwise product of the columns of $\mathbf{F}$.

A.4. The blur operator

While Synth estimates the field map corrections at the $1\mathrm{mm}$ isotropic voxel resolution of the T1w/T2w and synthetic EPI images, the true EPI images are acquired at much lower resolutions. In the case of the MSC dataset, this corresponds to a $4\mathrm{mm}$ native voxel resolution resampled to $1\mathrm{mm}$ isotropic resolution; and for the ABCD dataset this corresponds to $2.6\mathrm{mm}$ native resolution resampled to $1\mathrm{mm}$. When lower resolution EPI images are upsampled to the higher resolution grid, their lower resolution manifests as apparent blur. Synth models this blur ($\mathbf{E}$ in Eq. (1)) as the convolution of a source image with an Epanachnikov kernel, whose bandwidth is related to the resolution difference between the EPI and anatomical images. This kernel is given by the following expression:

where $\mathbf{x}\in {\mathbb{R}}^3$ is a location on the kernel, and $h$ is the bandwidth parameter. For the results in this paper, the bandwidths of the Epanechnikov kernels were set to $h=6$ and $h=10$ for the ABCD and MSC dataset respectively.

Synth models the difference in effective spatial resolution between T1w/T2w and EPI images by blurring the radial basis function (RBF) components of the anatomical images with an Epanachnikov smoothing kernel prior to fitting the model to the target EPI image. In principle, the optical blurring should occur only after the field map deformation is applied, rather than prior to it. However, implementing the model in this way would require repeatedly blurring the high spatial resolution radial basis components of the T1w/T2w source images after every iterative update of the estimated field map parameters. A naive implementation of this approach would be computationally prohibitive. Prior research has shown that in practice the sort of approximation implemented in Synth in many cases does not introduce a significant amount of error (Zhang et al., 2012). However, while the quality of the field map corrections produced by Synth appears to be quite high, some improvement may yet be gleaned by using the “correct model”, assuming the technical challenges to its implementation can be surmounted.

A.5. Contrast correction for error induced image bias

When source and target images are acquired at substantially different spatial resolutions, or in the early phases of alignment when source and target images are far out of register, the synthetic images produced by Synth exhibit significantly “flattened” contrast, even if the rank order voxel intensities comprising the synthetic image largely correspond to those observed in the target image. This is due to the fact that errors in registration or large differences between source and target image resolution bias estimates of the radial basis function parameters, $\mathit{\theta }$, toward the target image mean. In an extreme, but illustrative, case where image registration is so far off that there is no meaningful relationship between corresponding source and target image voxels, the synthetic image will simply be a volume in which all voxels have the same intensity as the target image mean. That is, it will have perfectly ‘flat’ contrast. As registration improves some, the synthetic image will begin to exhibit more accurately the contrast properties of the target image. In general, however, the contrast of naively estimated synthetic images will tend to be somewhat flattened due to error. To address this, Synth includes a tone curve adjustment option that increases contrast of the final synthetic image. This is accomplished by passing the synthetic image voxel intensities through a monotonic non-linear curve, a process similar to tone curve adjustment in traditional image processing. The curve is constrained to reside within the set of beta function cumulative distributions –distributions well suited for modeling tone curves that tend to be sigmoidal– and its parameters, $\alpha$ and $\beta$, are optimized by minimizing the negative linear correlation coefficient between the target and synthetic images.

A.6. Estimating field maps with an undistorted synthetic EPI image

The warp parameters that map the undistorted synthetic EPI image to the reference EPI image were estimated using a three step iterative procedure. The final iteration of this procedure produced the estimated field map corrections used for our analyses. This procedure was conducted as follows:

• 1.

  Using Synth, construct an initial synthetic EPI image based on the current best affine alignment between T1w, T2w and a reference EPI volume. An EPI signal-to-noise mask was used to reduce the contribution of low signal-to-noise areas when constructing the synthetic EPI (see Appendix A.7).

• 2.

  Estimate a warp that maps the current synthetic EPI to the reference EPI image. With the exception of the first iteration, each warp estimate is initialized from the warp estimate of the previous iteration.

• 3.

  Apply the inverse of this warp, which approximates a field map correction, to the reference EPI image.

• 4.

  Repeat steps 1–3, at each iteration, using Synth to re-estimate the synthetic EPI image using the most recently corrected EPI image (see Eq. (1)).

For the results presented here, we estimated the field maps using the ANTs SyN algorithm. Each subsequent iteration used a smaller gradient step size, 3, 1, 0.1 for each iteration respectively. SyN parameters were set to use a neighborhood cross-correlation cost function (radius = 2), update and total field variance parameters of 0, winsorize limits of [0.005, 0.995], histogram matching, BSpline interpolation of order 5 (Lehmann et al., 1999), convergence parameters of 500x0, 1e-6 and 15, smoothing parameters of 0x0, and shrink factor parameters 2x1. In order to leverage the detailed spatial information available in the high resolution synthetic EPI images during the estimation of the field map, the warps at each iteration were constructed in the 1 mm $x$ 1 mm $x$ 1 mm native resolution of the T1w/T2w source images.

A.7. Creating the EPI signal-to-noise mask

To construct a binary mask that would serve to label EPI voxels suffering from a high degree of noise or signal dropout we implemented a linear discriminant analysis procedure. We began by using the brain mask derived from the skullstripping stage of the anatomical volume and aligned it to the functional dataset using the affine transformation generated from Synth. This mask served as an initial guess, or prior estimate, for which voxels were likely to represent meaningful fMRI signals in the resting state EPI dataset, and those which represent noise and areas of signal dropout Supplemental Figure 3. Using this initial labeling, we implemented a between-class linear discriminant analysis (LDA) based, first, on the raw un-centered and un-scaled voxel time series and then computed the projection of each voxel’s time series onto the resulting principal discriminant vector. We repeated this LDA procedure on the centered and normalized time series from each voxel, again projecting the times series back onto the principal discriminant vector. Thus, each voxel was associated with two projection values. We then created a third value associated with each voxel by computing the product of the two principal discriminant vector projections. A second level of LDA was done on these three values. Lastly, Otsu’s method (Otsu, 1979) was used to label each voxel as ‘signal’ or ‘noise’, depending upon its projection onto this final 3-element principal discriminant vector. We refined this labeling by repeating the entire LDA procedure over 20 iterations. During each subsequent iteration the labels produced by the previous iteration were used during the LDA procedure. Sufficiently iterated, this approach converges to a stable binary signal-or-noise label for each voxel. Lastly, to construct a mask containing only within-brain noise voxels –voxels most likely to be affected by signal dropout– we masked the final result of this iterative process to remove non-brain voxels (defined by the skullstrip mask) and morphologically open and dilate the image by 1 voxel. This mask was used as a weighting volume to further refine the estimate of the synthetic functional EPI by reducing the contributions of brain regions where no correspondence could exist between EPI images and T1w/T2w images due to signal drop out.

A.8. Distortion correction using alternating descent optimization with synth

Distortion correction with Synth is accomplished through an alternating descent optimization approach, which solves for each parameter individually while keeping the others fixed. At the end of each iteration of alternating descent, updates to the distortion correction estimation improves the accuracy of the synthetic image constructed during the subsequent iteration. These sequential improvements allow for iteratively more accurate distortion correction. The complete procedure is described in Algorithm 1.

Prior to the alternating descent procedure, initial warp parameters, ${\mathit{\varphi }}_0$, are initialized by a rigid-body transform that aligns the T2w anatomical image and EPI image using a mutual information cost function. The rigid-body parameters from this step are then converted into a displacement field and applied to the EPI image, aligning it to the anatomical images that underlie the matrix, $\mathbf{F}$. Algorithm 1 begins by constructing the radial basis function matrix, $\mathbf{F}$, from the T1w/T2w anatomical images. Then the main alternating minimization loop is run for each iteration $n$. First, ${\mathit{\theta }}_n$ is computed by solving the linear system, $f^{-1}\left(\mathbf{y};{\mathit{\varphi }}_{n-1}\right)=\mathbf{E}\mathbf{F}{\mathit{\theta }}_n$. Here, $f^{-1}$ denotes the inverse displacement field that non-linearly transforms the real EPI image to the space of the synthetic EPI image. ${\mathit{\theta }}_n$ is then used to construct an intermediate synthetic image. Next, global contrast parameters, ${\alpha }_n$ and ${\beta }_n$ are estimated to improve overall similarity between the EPI image and the intermediate synthetic image. Finally, a new estimate for the displacement field parameters, ${\mathit{\varphi }}_n$, are estimated using the ANTs SyN non-linear warping algorithm. Each iteration updates ${\mathit{\varphi }}_n$, ${\mathit{\theta }}_n$, ${\alpha }_n$, and ${\beta }_n$ until a specified number of iterations, $N$, is reached.