PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 April 19.
Published in final edited form as:
PMCID: PMC2856353
NIHMSID: NIHMS160353

Regularized field map estimation in MRI

Amanda K. Funai, Student Member, IEEE, Jeffrey A. Fessler, Fellow, IEEE, Desmond T.B. Yeo, Student Member, IEEE, Valur T. Olafsson, Student Member, IEEE, and Douglas C. Noll, Senior Member, IEEE

Abstract

In fast MR imaging with long readout times, such as echo-planar imaging (EPI) and spiral scans, it is important to correct for the effects of field inhomogeneity to reduce image distortion and blurring. Such corrections require an accurate field map, a map of the off-resonance frequency at each voxel. Standard field map estimation methods yield noisy field maps, particularly in image regions with low spin density. This paper, describes regularized methods for field map estimation from two or more MR scans having different echo times. These methods exploit the fact that field maps are generally smooth functions. The methods use algorithms that decrease monotonically a regularized least-squares cost function, even though the problem is highly nonlinear. Results show that the proposed regularized methods significantly improve the quality of field map estimates relative to conventional unregularized methods.

I. Introduction

MR imaging techniques with long readout times, such as echo-planar imaging (EPI) and spiral scans, suffer from the effects of field inhomogeneity that cause blur and image distortion. To reduce these effects via field-corrected MR image reconstruction, e.g., [1]–[5], one must have available an accurate estimate of the field map. A common approach to measuring field maps is to acquire two scans with different echo times, and then to reconstruct the images (without field correction) from those two scans. The conventional method is then to compute their phase difference and divide by the echo time difference Δ1. This model makes no account for noise and creates field maps that are very noisy in voxels with low spin density. Section II first introduces this model and then reviews standard approaches for this problem. A limitation of the standard two-scan approach to field mapping is that selecting the echo-time-difference Δ1 involves a trade off: if Δ1 is too large, then undesirable phase wrapping will occur, but if Δ1 is too small, the variance of the field map is large. One way to reduce the variance while also avoiding phase unwrapping procedures is to acquire more than two scans, e.g., one pair with a small echo time difference and a third scan with a larger echo time difference. By using multiple echo readouts, the scan times may remain reasonable, at least for the modest spatial resolutions needed in fMRI. Therefore, we present a general model that accommodates more than two scans and describe a regularized least-squares field map estimation method using those scans. Section III shows the improvements both in the estimated field maps and the reconstructed images using multiple scans. This is shown first with simulated results in Section III-A and then using real MR data in Section III-B.

II. Multiple Scan Fieldmap Estimation - Theory

A. Reconstructed Image Model

The usual approach to measuring field maps in MRI is to acquire two scans of the object with slightly different echo times, and then to reconstruct images y0 and y1 (without field correction) from those two scans, e.g., [6]–[8]. We assume the following model for those undistorted reconstructed images is

yj0=fj+εj0yj1=fjeωjΔ1+εj1,
(1)

where Δ1 denotes the echo-time difference, fj denotes the underlying complex transverse magnetization in the jth voxel which is a function of the spin density, and εj denotes (complex) noise. The goal in field mapping is to estimate an (undistorted) field map ω = (ω1, …,ωN) from y0 and y1, whereas f = (f1, …, fN) is a nuisance parameter vector. This section reviews the standard approach for this problem, other approaches in the literature, and then describes a new and improved method.

B. Conventional Field Map Estimator

Based on (1), the usual field map estimator [omega with circumflex]j uses the phase difference of the two images, computed as follows [1], [9]:

ω^j=(yj0*yj1)/Δ1.
(2)

This expression is a method-of-moments approach that would work perfectly in the absence of noise and phase wrapping, within voxels where |fj| > 0. However, (2) can be very sensitive to noise in voxels where the image magnitude |fj| is small relative to the noise deviations. Furthermore, that estimate ignores our a priori knowledge that field maps tend to be smooth or piecewise smooth. Although one could try to smooth the above estimate using a low pass filter, usually many of the [omega with circumflex]j values are severely corrupted so smoothing would further propagate such errors (see Fig. 2 top right). Instead, we propose below to integrate the smoothing into the estimation of ω in the first place rather than trying to “fix” the noise in [omega with circumflex] by post processing.

Fig. 2
Top row: magnitude image |yj|, conventional field map estimate(2). Middle row: (field map estimates): penalized-likelihood estimate using (13) with β = 2−6 (left) and β = 2−3 (right). Bottom row: Maps of the spatial resolution ...

C. Other Field Map Estimators

Although the conventional estimate (2) is most common, other methods for estimating field maps have appeared in the literature.

Different techniques have been proposed that incorporate field map acquisition with image acquisition ([7] for projection reconstruction and [10] for spiral scans). Chen et al. in [11] used multiple echos during EPI acquisition and used these distorted scans to create a final corrected undistorted image. Priest et al. in [12] used a two-shot EPI technique to obtain a field map for each image; this could prevent changes in the field map due to subject motion from being propagated through an entire fMRI time series.

Stand alone field map acquisition techniques have also been proposed. Windischberger et al. [13] used three echos and corrected for phase wrapping by classifying the degree of phase wrapping into seven categories. They then used linear regression to create a field map that is followed by median and Gaussian filtering. Reber et al. [14] used ten separate echo times and acquired distorted EPI images. They used a standard phase unwrapping technique of adding multiples of 2π and then spatially smoothed the image with a Gaussian filter. While these techniques both seek to use more echos to increase the accuracy of the field map, they have several disadvantages. Neither are based on a statistical model and, thus, do not consider any noise in developing their estimator. The filtering suggested by both techniques also adds additional blur. Aksit et al. [15] used three scans, the first two with a small echo time and no phase unwrapping and the third with a larger echo time. Two techniques were tried: 1) phase unwrapping by using the first two sets of data and 2) taking a Fourier transform to determine the EPI shift experienced. In phantom studies, using three scans yielded half to a third of the error of two scans. Because the estimator uses a linear fit, there is still error in voxels near phase discontinuities and along areas of large susceptibility differences.

An additional technique used to improve the conventional estimate is local (non-linear) fitting, e.g. [16], [17]. While this can improve the conventional estimate, we desire a more statistical approach.

Our technique is unique in that it uses a statistical model using multiple scans and operates without the constraint of linearity. By using a penalized-likelihood cost function, we can easily adjust the regularization parameter to control the amount of smoothing without any additional filtering step. By using a field map derived from the first two echos as the initialization for the iterative method (assuming the two echos are close enough together), no phase unwrapping is required. Our model also takes into account R2* decay which was ignored in previous multiple echo techniques.

D. Multiple Scan Model

We now generalize the conventional model (1) to the case of multiple scans, i.e., with more than one echo time difference. The reconstructed images are denoted here by y0, …, yL, where L is the number of echo time differences. Because we are using multiple echo time differences, R2* decay may no longer be negligible and should be included in our model. Our model for these images is:

yjl=fjeωjΔleRjΔl+εjl,
(3)

for l = 0, …, L, where Δl denotes the echo time difference of the lth scan relative to the original scan i.e., (Δ0 = 0), where j denotes the voxel number and where Rj denotes the value of R2* for the jth voxel. As in most field map estimation methods, this model assumes implicitly there is no motion between the scans. As in (1), fj denotes the complex transverse magnetization and εjl denotes the (complex) noise. If we choose the Δl values carefully, this data model allows for a scan that is free or largely free of phase wraps but which gives a phase difference lower in SNR, as well as scan(s) with wrapped phase but higher in SNR. Including the scan(s) with a larger echo time difference should help reduce noise in [omega with circumflex]j whereas the wrap-free scan helps avoid the need for phase unwrapping tools.

E. Maximum-Likelihood Field Map Estimation

The conventional estimate (2) appears to disregard noise effects, so a natural alternative approach is to estimate ω using a maximum likelihood (ML) method based on a statistical model for the measurements y. In MR, the k-space measurements have zero-mean white gaussian complex noise [18], and we furthermore assume here that the additive noise values in y in (3) have independent gaussian distributions1 with the same variance σ2. Under these assumptions, the joint log likelihood for f and ω given y = (y0, …, yL) is

logp(y;f,ω)=l=0Llogp(yl;f,ω)12σ2j=1Nl=0L|yjlfjeωjΔleRjΔl|2
(4)

where “[equivalent]” denotes equality to within constants independent of f and ω. If the Rj values were known, the joint ML estimate of f and ω could be solved by the following minimization problem:

argminωN,fNj=1N[yj0yj1yjL][1eωjΔ1eRjΔ1eωjΔLeRjΔL]fj2.
(5)

This problem is quadratic in fj; minimizing over fj yields the following ML estimate:

f^j=l=0LyjleωjΔleRjΔll=0Le2RjΔl.
(6)

Substituting this estimate back into the cost function (5) and simplifying considerably yields the following cost function used for ML estimation of ω:

ΨML(ω)j=1Nm=0Ln=0L|yjmyjn|·wjm,n·[1cos(yjnyjmωj(ΔnΔm))],
(7)

where wjm,n is a weighting factor that depends on R2* as follows:

wjm,n=eRj(Δm+Δn)l=0Le2RjΔl.
(8)

Similar weighting appeared in the weighted phase estimate proposed in [19] for angiography. The ML cost function ΨML(ω) is periodic, similar to cost functions used in phase unwrapping problems, e.g., [20]. The cost function (7) appears to require either knowledge of or a good estimate of R2*. However, we note that:

|E[yjm]|=|fj||eRjΔm|;

therefore, hereafter, we approximate wjm,n, as follows:

wjm,n|yjm||yjn|l=0L|yjl|2.
(9)

This approximation does not require knowledge of R2* values.

There is no analytical solution for the minimizer, ω in (7), except in the L = 1 case. Thus, iterative minimization methods are required, even for the ML estimator.

F. Special Case: L = 1 (Conventional Two Scans)

In the case where L = 1 usually Δ1 is chosen small enough that we can ignore R2* decay (i.e., let R2*=0) and the ML cost function in (7) simplifies to

ΨML(ω)j=1N|yj0yj1|[1cos(yj1yj0ωjΔ1)].
(10)

The ML estimate is not unique here due to the possibility of phase wrapping. But ignoring that issue, the ML estimate of ω is ω^j=(yj1yj0)/Δ1, because 1 − cos(t) has a minimum at zero. This ML estimate is simply the usual estimate (2) once again to within multiples of 2π. Thus the usual field mapping method (for L = 1) is in fact an ML estimator under the white gaussian noise model assuming R2*=0. The more general cost function (7) for the field map ML estimator for L > 1 is new to our knowledge.

G. Penalized-Likelihood Field Map Estimation

The ML estimator ignores our a priori knowledge that field maps tend to be spatially smooth functions due to the physical nature of main field inhomogeneity and susceptibility effects2. (We note that this assumption does not address the presence of signal from fat). A natural approach to incorporating this characteristic is to add a regularizing roughness penalty to the cost function. Here we regularize only the phase map ω and not the magnetization map f; we expect f to be far less smooth because it contains anatomical details. Such regularization is equivalent to replacing ML estimation with the following penalized-likelihood estimator:

(ω^,f^)=argmaxω,fl=0Llogp(yl;f)βR(ω),

where R(ω) is a spatial roughness penalty (or log prior in a Bayesian MAP philosophy). Based on (6) and (7), after solving for f and substituting it back in, the resulting regularized cost function has the form

ΨPL(ω)ΨML(ω)+βR(ω),
(11)

where we use the approximation (9) for ΨML(ω). This cost function automatically gives low weight to any voxels where the magnitude |yjmyjn| is small. For such voxels, the regularization term will have the effect of smoothing or extrapolating the neighboring values. Thus, this approach avoids the phase “outlier” problem that plagues the usual estimate (2) in voxels with low signal magnitude. If ω, corresponds to a N1 × N2 field map ωn,m, then a typical regularizing roughness penalty uses the second-order finite differences between horizontal and vertical neighboring voxel values as follows:

R(ω)=n=1N11m=0N21ψ(2ωn,mωn1,mωn+1,m)+n=0N11m=1N21ψ(2ωn,mωn,m1ωn,m+1)
(12)

where ψ is a convex “potential function.” Here, we use the quadratic potential function, ψ(t) = t2/2. In this paper, we used second-order differences for all results; we found that second-order finite differences are preferable to first-order differences because the resulting PSF tails decrease more rapidly even when the FWHM values are identical. A quadratic potential function has the advantage of being differentiable and easy to analyze, especially with Gaussian noise. Although quadratic regularization blurs edges, we assume the field map is smooth, so a more complicated potential function, such as Huber, is not considered here.

Usually ψ is differentiable, so we can minimize the cost function Ψ(ω) either by conventional gradient descent methods or by optimization transfer methods [21]–[23]. In particular, in the usual case where [psi](t)/t is bounded by unity, then the following iteration is guaranteed to decrease Ψ(ω) monotonically:

ω(n+1)=ω(n)diag{1dj+β·c}Ψ(ω(n)),
(13)

where [nabla] is the gradient of the cost function,

c{4,regularizationwith1st-orderdifferences16,regularizationwith2nd-orderdifferences
(14)

and

djm=0Ln=0L|yjmyjn|·wjm,n·(ΔnΔm)2,
(15)

using the approximation for wj shown in (9). For examples in this paper, we used a similar minimization algorithm described in Appendix A because of its faster convergence properties.

To initialize ω(0), we used the regularized ML estimate (11) based on the first two sets of data y0 and y1. We choose the echo times to avoid phase wrapping between these sets of data (this same idea is used in [15] in their three-point method). Therefore, there is no need to apply phase unwrapping algorithms - the algorithm will converge to a local minimizer in the “basin” of the initial estimate [21].

In [24], we considered approximating the 1 − cos term in (11) with its second-order Taylor series to create a penalized weighted least squares (PWLS) cost function. A simplified PWLS approach where the weights were thresholded was also considered. Those models ignore any phase wrap that may occur when evaluating (2). They also have increased error with little computational benefit. Therefore, those simplified methods are not explored further in this paper.

H. Spatial Resolution Analysis of Field Map Estimation

To use the regularized method (11) the user must select the regularization parameter β, which could seem tedious if one used trial-and-error methods. Fortunately, it is particularly simple to analyze the spatial resolution properties for this problem, using the methods in [25] for example. We make the second-order Taylor series approximation for this analysis. The local frequency response of the estimator using second-order finite differences at the jth voxel can be shown to be:

H(ΩX,ΩY)11+βdj(ΩX2+ΩY2)p,
(16)

where ΩX and ΩY are the Discrete Space Fourier Transform (DSFT) frequency variables. and where p = 1 for regularization based on first-order differences and p = 2 for second-order finite differences as in (12). (See [26] for related analyses.) From (16) we see that the spatial resolution at each voxel depends on the data through dj. In areas with small signal magnitudes, there will be more smoothing, as desired. The spatial resolution (16) also depends on the Δl values being used. Data from scans with larger Δl values have lower [omega with dot above]j variance (see (17) below), and will be smoothed less. However, data from these scans will also be affected by R2* decay through wjm,n if the data is not scaled to compensate for this factor. To simplify selecting β, we normalize the data by the median of the square root of (15) using the approximation (9) for ωj. Normalizing by this factor allows us to create a standard β to FWHM table or graph (e.g., Fig. 1). If this normalization were not applied, a similar figure would need to be calculated with each new data set (or at least with each new set of Δl values) or β would need to be chosen empirically. Normalizing based on the analytical result (16) enables us to use the same β for all scans.

Fig. 1
Angularly averaged FWHM of PSF for field map estimation as a function of log2β for dj = 1 in (16).

We used the inverse 2D DSFT of (16) to compute the PSF h[n, m] and tabulate its FWHM as a function of β, assuming the previous corrections were made and that the pixel j has dj = 1. Fig. 1 shows this FWHM as a function of log2(β), for both p = 1 and p = 2. The FWHM increases monotonically with β, as expected, although the “knees” in the curve are curious. Nevertheless, one can use this graph to select the appropriate β given the desired spatial resolution in the estimated field map. The resulting spatial resolution will be inherently nonuniform, with more smoothing in the regions with low magnitudes and vice versa. One could explore modified regularization methods [25] to make the resolution uniform, but in this application nonuniform resolution seems appropriate since the goals include “interpolating” across signal voids.

I. Qualitative Example: L = 1

Fig. 2 shows an example of the data magnitude |yj0| and the usual phase estimate based on L = 1 (2) which is very noisy. This is real data taken from a 3T MR scanner with Δ1 = 2 ms. The maximum value of |ωj · Δ1| is 1.61 radians in nonzero voxels, making the scan free of any phase wraps. Fig. 2 also shows the penalized-likelihood estimate based on (13) using two different values for β and using 150 iterations. Here, we can see the improvement from using a regularized estimator versus the conventional ML estimator. The effect of β on the smoothness of the estimate is also seen. The improvement seen is analyzed quantitatively in Section III. Fig. 2 also shows the effective FWHM (in pixels) of the regularized algorithm based on (16) for both values of β. Most of the image has a FWHM corresponding to the chosen β based on Fig. 1. Areas of low magnitude have a much higher FWHM (such as the sinuses) and areas of high magnitude have the lowest FWHM.

J. Theoretical Improvements Over 2 Data Sets

Using more than two sets of data requires a longer data acquisition and also involves choosing the Δl values. Analyzing the theoretical improvements that may be attained by using multiple data sets can help determine when the increased acquisition time is warranted and can guide our choice of the Δl values. Therefore, we calculated the Cramér-Rao bound (CRB) for the model (3). This bound expresses the lowest achievable variance possible for an unbiased estimator based on a given model. Although a biased estimator (the penalized-likelihood estimator) is used in our implementation, the bound quantifies the maximal improvement possible based on the model and allows for a comparison on how close our implementation is to the ideal, unbiased case.

Because there are multiple unknown parameters in these models θ = (ωj, |fj|, [for all]fj), the multiple parameter CRB must be used. In that case, the matrix CRB is

Covθ{θ^}F1(θ)

where F(θ)=E[θ2lnp(Y;θ)] is the Fisher information. Because fj is a nuisance parameter, we focus on the CRB for the variance of ωj, although the effect of fj will be felt through the inversion of the Fisher matrix. For simplicity, we initially set R2* to 0 in the CRB derivations shown below.

Applying the CRB to the L echo-time difference model (3) yields, after considerable simplification, the expression:

VarL{ω^j}σ2(L+1)Δ12|fj|2λL,
(17)

where, defining αl = Δl / Δ1:

λL(1L+1l=0Lαl2)(1L+1l=0Lαl)2.

The variance reduces, in general, as L is increased. The expression for λL is the “variance” of {α0, α1, (...) αL}, measuring the variance between the echo time differences. Increasing the variance (spread) of the Δl values will decrease the overall variance of the field map estimate.

For the L = 1 (2 sets of data) model, λ1 = 1/4 and (17) simplifies to:

CRB12σ2Δ12|fj|2.

As expected, the field map variance decreases when the signal strength |fj|, or echo time difference Δ1, increase. For an unbiased estimator based on the model (3) with L = 2 (3 sets of data) one can show:

CRB2CRB14/3(α22α2+1).
(18)

Interestingly, simply using three scans, but using Δ2 = Δ1 (or α2 = 1), would reduce the variance by only 4/3.

From (18), increasing α2 should decrease the variance for an unbiased estimator. Making α2 arbitrarily large, however, is not advisable for many reasons. A larger α2 creates more phase-wrapping. Eventually, the wrapping will lead to intra-voxel aliasing and the desired improvement would be unattainable. Another problem with large values of αl is the effect on the MR pulse sequence length. A large α2 also causes much more R2* decay in the signal as shown in (7). Choosing optimal Δl values requires some knowledge of R2* decay. This can be seen more clearly in the CRB bounds for the model (3) with R2* decay included. For the L = 1 model, one can show:

Var1{ω^j}CRB1·1+e2RjΔ12.
(19)

For the L = 2 (3 sets of data) model:

Var2{ω^j}σ2Δ12|fj|21+e2RjΔ1+e2RjΔ1α2b,
(20)

where

be2RjΔ1+α22e2RjΔ1α2+(1+α222α2)e2RjΔ1(1+α2).

Using these expressions, we can optimize the Δl values, which will be inversely proportional to the value of R2*. In fact, for L = 1, one can show that the optimal choice is Δ1opt=1.11/Rj. Therefore, small values of αl based on the amount of R2* decay expected should be used.

III. Experiments

A. Simulation: Comparison of L = 1 and L = 2 Methods

We compared the L = 1 and L = 2 methods with two examples. First, we used a simulated Gaussian true field map (Fig. 3) with a magnitude map equal to unity at all points. Second, we simulated a brain example. For the magnitude, we used a simulated normal T1-weighted brain image [27], [28]. We generated a simple field map consisting of a 4.8 cm diameter sphere of air (centered around the nasal cavity) embedded in water using simple geometrical equations [29], [30], using a slice slightly above the sphere. Fig. 4 shows the field map and magnitude image |fj|. We added complex Gaussian noise at many levels of SNR to the images. For this paper, we used the following definition of SNR:

SNR=20logfy0f.
(21)

The SNR remains consistent even when varying R2*, L, or α.

Fig. 3
Top: “True” field map for Gaussian example in Hz; Noisy (SNR = 10dB) wrapped phase yj2 with α2 = 3, Noisy (SNR = 10dB) wrapped phase with α2 = 7. Bottom: Conventional estimate for L = 1, PL estimates for L = 1, ...
Fig. 4
Top: True field map and magnitude for brain example and mask, (SNR = 8.5dB) wrapped phase for α2 = 3 and α2 = 5 images. Center and Bottom: Conventional, Conventional convolved with a Gaussian filter, PL with 2 sets (L = 1), and PL with ...

We used Δ1 = 2 msec for both cases. For the L = 2 case we also varied α2 to produce several Δ2 values. We used a uniform value of R2* = 20 sec−1 in generating our simulations.

The field map was reconstructed using the penalized-likelihood method (11) using normalization as described in Section II-H for both L = 1 and L = 2. The algorithm (13) was run at each SNR level for the L = 1 case and for the L = 2 case of data with varying values of α2 using 5 realizations. We ran 300 iterations of the algorithm, using β = 2−3.

We also applied the conventional estimator to our data. To reduce the noise, we convolved the conventional estimate with Gaussian filters of varying widths (σ = 0.0625, 0.1250, …, 3.125). We chose the “optimal” σ based on the minimum masked RMSE. Choosing the optimal σ using the true field map gives the conventional estimate an advantage in this example unavailable in practice.

The RMS error (in Hz) was computed between the “true” field map and the field map reconstructed using the PL method (11) and the conventional estimate. This RMSE was calculated in a masked region (pixels with magnitudes at least 20% of the maximum true magnitude).

Fig. 3 shows an example of the PL with L = 1 estimate compared to the PL with L = 2 estimate at α2 = 3 and α2 = 7 at an SNR of 10dB. Qualitatively, we can see improvements with increases in both L and α2. Fig. 4 shows similar results for the brain example.

The largest errors in these field maps occur where the magnitude is smallest. The RMSE is much higher using only the conventional method. We also calculated the RMSE in the sinus region of the brain (the ROI is shown in Fig. 4). We chose this ROI because the low magnitude makes the field map difficult to estimate here although the field inhomogeneity is also greatest here. The RMSE in this ROI was 61.1 Hz for the conventional estimate, 11.6 Hz for the Gaussian filtered estimate, 3.4 Hz for the L = 1 regularized estimate, and 1.9 Hz for the L = 2 α2 = 3 regularized estimate and 1.7 Hz for the L = 2 α2 = 5 regularized estimate. Overall, the filtered conventional estimate performed similar to the PL method with L = 1 over the masked region, but had higher error in the ROI. The PL method with L = 2 showed a decreased error in both the masked region and the ROI. We would expect even higher improvement over any practical Gaussian filtered estimate because a suboptimal σ would be used. The proposed regularized estimators are more accurate in pixels with low magnitude. Adding additional scans (L > 1) makes the PL estimate even more accurate.

Fig. 5 shows the improvement (defined as the RMS error for PL estimate with L = 1 divided by the error for PL estimate with L = 2) gained by using an additional set of data for the Gaussian example. For comparison, we also plotted the predicted improvement, given by the square root of the ratio of the expressions (19) and (20). The experimental gains are actually higher than the improvements anticipated as shown by the dotted lines (the predicted improvement) for some SNR values. Because this is a ratio of RMSEs and the amount of bias can vary between L = 1 and L = 2, the unbiased CRB provides a benchmark of expected ratios rather than an exact upper limit. Also, recall that (19) and (20) considered R2* to be a known value when, in fact, R2* is unknown and approximated through (9). The RMSE is low (in voxels with large magnitudes) at high SNRs using either L = 1 or L = 2. At lower SNRs, however, including in voxels with low magnitudes, using L = 2 and higher values of α2 greatly reduces RMS error. We repeated these simulations with R2* = 0 (results not shown) and the empirical improvement almost exactly matched (18).

Fig. 5
Improvement in the RMSE for the Gaussian example by using 3 data sets rather than 2 sets. Expected improvements shown by dotted lines.

Fig. 6 shows the improvement gained by using an additional set of data for the brain image. For a low SNR (for example 10 dB), the improvements are close to expected. The brain image has some areas where the magnitude is very low, making estimation using any method quite challenging. In addition, the field map phase itself is less smooth than in the Gaussian case, making the estimation more difficult. For a higher SNR (for example 20 dB), the 3-set case still outperforms the 2-set case substantially but by less than predicted by (18).

Fig. 6
Improvement in the RMSE for the brain example by using 3 data sets rather than 2 sets. Expected improvements shown by dotted lines.

The RMSE has components of both bias error and variance in it, as shown below:

RMSE(X)=Var{X}+bias2(X).

Therefore, we analyzed the bias and the standard deviation at a single representative SNR = 20 dB and at α2 = 1, 2, … 7 using 500 iterations and 100 realizations for each factor. Fig. 7 compares the standard deviation for each α2 relative to that at α2 =1 and the empirical improvements were compared to those predicted by the CRB (20) for the Gaussian example. As expected, the improvements in variance are very close to predicted. Here, the bias is also very low at all levels of SNR - explaining the improvement seen in RMSE in Fig. 5.

Fig. 7
Bias and RMSE improvement for Gaussian example. Top: Space-averaged σ and absolute bias for several α2 values; Bottom: RMSE improvement, empirical and expected, over α2 = 1 for several α2 values.

Fig. 8 shows the bias and standard deviation for a signle SNR = 20 for the brain example. The empirical variances were close to those expected. The bias, however, introduced in part by the regularization, was nearly constant (independent of α). So for large values of α2, the bias begins to dominate the variance in RMSE calculations, explaining Fig. 6.

Fig. 8
Bias and RMSE improvement for brain example. Top: Space-averaged σ and absolute bias for several α2 values; Bottom: RMSE improvement, empirical and expected, over α2 = 1 for several α2 values.

Overall, the variance reductions in both examples due to using three echo times were close to the results predicted by the CRB. For low values of α2 (i.e., five or less), the expected benefit using L > 1 holds even with a moderate value of R2*. The RMSE reductions are largest at lower SNRs. For phase estimation, the local SNR depends on the spin density of each voxel as seen in (17). Voxels with lower spin density effectively have lower SNR. It is precisely in these voxels where using 3 or more scans has the greatest benefit.

B. MR phantom data: Application to Spiral Trajectories

To illustrate how improved field map estimation leads to improved reconstructed images, we used field maps produced by the conventional method (2) and produced by the PL method with three scans (11) to correct real spiral MR data for field inhomogeneities. We imaged a phantom with large field inhomogeneity. We used a spiral-out trajectory with a TE of 30 ms, TR of 2 sec, and a flip angle of 90 degrees. We took six slices spaced 5 cm apart over the 15 cm field of view. First, we collected data to create the field maps (using eight interleaves to minimize the effect of the field inhomogeneity) at the original 30 ms, as well as at 32 ms (Δ1 = 2 ms) and at 34 ms (Δ2 = 4 ms) and at 40 ms (Δ2 = 10 ms). We took ten realizations for each echo difference. We reconstructed iteratively the resulting 64 × 64 pixel images in a masked region using [31]. Then, we used these images to create (for each slice) a conventional field map (2), a conventional field map blurred with a Gaussian filter, a PL field map with L = 1, a PL field map with L = 2 and α2 = 2, a PL field map with L = 2 and α2 = 5, and a PL field map with L = 3, (11). We used β = 2−6 for the regularized iterative algorithm and σ = .5 for the Gaussian filter approach, approximately matching the FWHM of the two approaches. Finally, we collected one-shot spiral out data with TE = 30 ms. This scan is thus much more affected by final inhomogeneity. We collected two realizations and then averaged them in k-space. We first reconstructed this data iteratively without a field map as in [31]. Uncorrelated field inhomogeneity causes a blurred image for spiral trajectories. Finally, we iteratively reconstructed this one-shot data with each of the field maps previously created as in [2].

Fig. 9 shows one representative slice. The regularized field maps are less noisy than the conventional one, especially in areas of low magnitude and along the edges. Fig. 9 illustrates the blur and distortion in the one-shot image reconstructed without a field map. The images reconstructed with a field map do not have this blur. Nevertheless, a noisy field map can cause error in the reconstructed image. For example, in Fig. 9, the image reconstructed with the conventional field map shows more artifacts than the eight-shot data or either of the images reconstructed with regularized field maps. Using the eight-shot data as “truth”, we computed the NRMSE of each image and Table I shows the mean and variance over the ten realizations. We include data from two representative slices to show a range of values, although slice three is not shown. In addition, we calculated the NRMSE in the one-shot reconstructed images in pixels where the magnitude is less than .2 times the maximum pixel value of the eight-shot reconstructed image to see if the regularized field maps reduce errors in areas of the image with low magnitude. This is also reported in Table I. We use the norm of the eight-shot 30 ms image for normalization. The regularized iterative PL methods have a lower RMSE and much less variability than the other methods. Therefore, these regularized methods (especially with more than one echo time) give a very reliable estimate of the field map with little variability.

Fig. 9
First Slice - Top: Reconstructed 8-shot image, Conventional field map, Gaussian filtered field map, regularized field map L=1, regularized field map L=2; α2 = 2, regularized field map L=2 α2 = 5, regularized field map L=3. The field maps ...
TABLE I
Phantom NRMSE for two representative slices

IV. Discussion

We described a regularized method for field map estimation using two or more scans: the penalized-likelihood method (11). This method yields field maps that interpolate smoothly over regions with low spin density, thereby avoiding phase outliers that plague the conventional estimate (2). The method has been used with L = 1 (without full description) in [3], [32], [33].

Our analysis also shows that the conventional estimate (2) is in fact the ML estimate, a property that has previously gone unnoticed to our knowledge.

We also analyzed the spatial resolution properties of this method, leading to a practical procedure for choosing the regularization parameter to achieve a given desired spatial resolution.

We studied the CRB on the variance of the estimate for this method and found that our empirical simulation results for the PL method compared favorably, showing a reduction in the RMSE in comparison to using only two scans.

We collected real MR phantom data and created conventional and PL estimates of the field map which were used to reconstruct final images. The PL estimate reduces image artifacts caused by the field inhomogeneity and has a reduced RMSE, especially in areas of very low magnitude where the conventional estimate has many errors. Omitting or using a poor field map estimate for image reconstruction can dramatically affect the final image quality.

As noted in Section II-D, our cost function assumes, as do most other field map estimation problems, that there is no motion between scans. While our analysis indicated that a larger L is better in terms of variance, motion could be a problem during the larger time required for L echo time differences. Practically, L = 1 or L = 2 are the most likely choices for L and here motion is less likely to be an issue. If a larger number of echo differences are desired, then the cost function could be further generalized to include a joint estimation of the field map and rigid motion parameters.

We have focused here on the case of a single receive coil. It is straightforward to generalize the method for phased array coils, cf. [34].

Although we did not estimate R2* we used a simple weighting (9) in our algorithm to partially account for R2* decay; the improvements seen over estimation with two scans are still large, especially when using a small value of α2.

While this method assumed the first two echo time differences were close enough to prevent phase wrapping, this method could, with proper initialization, extend to data with larger echo time differences and some phase wrapping. This is especially interesting at higher field strengths where wrapping still exists at low echo time differences.

Overall, this method has potential to be a reliable estimator for MR field maps, able to utilize many scans to produce a good estimate. The general penalized-likelihood approach in this work is also applicable to estimating other parametric maps in MRI, such as relaxation maps [35] and sensitivity maps [36]. It may also be useful for phase unwrapping problems with noisy data. In some cases, it may be preferable to use edge-preserving regularization in (12), such as the Huber potential function [37].

Ultimately, this method is a tool that may help answer the main question of field mapping: how to best allocate scan time to achieve the most accurate field map. The preliminary CRB analysis guides choice of echo times given a set number of scans. In future work, we wish to further explore the relationship between number of echoes, signal to noise ratio, and spatial resolution.

Acknowledgments

Supported in part by NIH grants EB002683 and DA15410 and NSF-GRFP.

APPENDIX A MINIMIZATION ALGORITHMS

To minimize the cost function (11) developed in this paper, we need a method that will decrease it monotonically. The simple minimization algorithm shown in (13) is guaranteed to decrease Ψ(ω) monotonically; the proof that ensures mono-tonicity uses the fact that the second derivative of 1 – cos t is bounded above by unity. This algorithm will converge to a local minimizer of Ψ(ω) within the “basin” that contains the initial estimate [38].

However, this simple minimization algorithm shown in (13) is only one possible option to minimize the cost function given in (11). In our implementation, we used an optimization transfer approach to refine the iterative algorithm [22], [38]. First express (11) as shown below:

Ψ(ω)j=1Nm=0Ll=0Lφjml(ωj)+βR(ω),
(22)

where we define

φjml(ω)|yjmyjl|wjm,lφ(ω(ΔlΔm)+yjmyjl)

with

φ(t)1cos(t).

To minimize this cost function, we adopt an optimization transfer approach, for which we need a surrogate function for [var phi](s). In particular, we use the following parabola surrogate for [var phi]:

φ(t)q(t;s)=Δφ(s)+φ˙(s)(ts)+12κφ({s}2π)(ts)2

where {s} denotes the principle value of s. Huber stated that parabola surrogate functions (which he called a comparison function) exist for [var phi] that satisfy Huber’s conditions [39, p.184-5]; the functions must be differentiable, symmetric, and have curvatures (κ[var phi](s)) that are bounded and monotone non-increasing for s > 0. For |s| ≤ π, [var phi](s) shown above satisfies Huber’s conditions. We note

φ˙(s)=sin(s)

and

κφ(s)φ˙(s)s=sin(s)s.

Substituting this curvature κ[var phi](s) into the expression for [var phi]jml(ω) gives us the following curvature for the parabola surrogate

κφ,jml(s)=Δφ˙jml(s)s=|yjmyjl|wjm,l(ΔlΔm)2sin(s)s,

which is bounded as s → 0 and decreasing as |s| increases. For values of |s| > π, we exploit the periodicity of [var phi] and find an integer k such that |s − k2π| ≤ π, i.e., the principal value of the phase s. Fig. 10 shows [var phi] and parabola surrogates for several values of s. When s is an even multiple of π, the curvature κ[var phi] is the maximum curvature of [var phi]. When s is an odd multiple of π, the curvature κ[var phi] is zero, and [phi] is also zero, so the surrogate function is a constant.

Fig. 10
Illustration of [var phi](t) and quadratic surrogates for several values of s.

Aggregating such surrogates leads to the following surrogate function for the cost function Ψ(ω):

ϕ(n)(ω)j=1Nm=0Ll=0Lqjml(n)(ωj)+βR(ω)

where

qjml(n)(ω)φjml(ωj(n))+φ˙jml(ωj(n))(ωωj(n))+12κφ,jml(sr(n))(ωωj(n))2,

and where

sr(n)(ωj(n)|ΔlΔm|+yjmyjl)modπ[π,π].

If the roughness penalty R(ω) is a quadratic function, which is the natural choice for smooth phase maps, then the surrogate ϕ(n) above is a quadratic function that can be minimized easily by any classical method such as the conjugate gradient algorithm.

In our implementation, we used a separable quadratic surrogate algorithm to minimize this cost function [40]. Then, the following iteration, similar to that of (13), is guaranteed to decrease Ψ(ω) monotonically:

ω(n+1)=ω(n)diag{1d˜j(n)+β·c}Ψ(ω(n)),
(23)

where c was defined in (14) and where

d˜j(n)=m=0Ll=0Lκφ,jml(sr(n)).

The advantage of (23) over (13) is that d˜j(n)dj (15), so (23) will converge faster [41].

Footnotes

1Independence in image space is an approximation. In k-space, data points are independently distributed, but reconstruction may produce correlations, especially in scans with non-Cartesian k-space imaging.

2There may be discontinuities at air/water boundaries. Even in this case, sharp boundaries can be problematic if there is motion between scans, further motivating the use of regularization.

REFERENCES

1. Sekihara K, Matsui S, Kohno H. NMR imaging for magnets with large nonuniformities. IEEE Trans. Med. Imag. 1985 Dec;vol. 4(no. 4):193–199. [PubMed]
2. Sutton BP, Noll DC, Fessler JA. Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities. IEEE Trans. Med. Imag. 2003 Feb;vol. 22(no. 2):178–188. [PubMed]
3. Noll DC, Fessler JA, Sutton BP. Conjugate phase MRI reconstruction with spatially variant sample density correction. IEEE Trans. Med. Imag. 2005 Mar;vol. 24(no. 3):325–336. [PubMed]
4. Twieg DB. Parsing local signal evolution directly from a single-shot MRI signal: A new approach for fMRI. Mag. Res. Med. 2003 Nov;vol. 50(no. 5):1043–1052. [PubMed]
5. Barnet C, Tsao J, Pruessmann KP. Efficient iterative reconstruction for MRI in strongly inhomogeneous B0. Proc. Intl. Soc. Mag. Res. Med. 2004:347.
6. Jezzard P, Balaban RS. Correction for geometric distortion in echo planar images from B0 field variations. Mag. Res. Med. 1995 Jul;vol. 34(no. 1):65–73. [PubMed]
7. Nayak KS, Nishimura DG. Automatic field map generation and off-resonance correction for projection reconstruction imaging. Mag. Res. Med. 2000 Jan;vol. 43(no. 1):151–154. [PubMed]
8. Cusack R, Brett M, Osswald K. An evaluation of the use of magnetic field-maps to undistort echo-planar images. NeuroImage. 2003 Jan;vol. 18(no. 1):127–142. [PubMed]
9. Glover G. Multipoint Dixon technique for water and fat proton and susceptibility imaging. J. Mag. Res. Im. 1991 Sep;vol. 1(no. 5):521–530. [PubMed]
10. Nayak KS, Tsai C-M, Meyer CH, Nishimura DG. Efficient off-resonance correction for spiral imaging. Mag. Res. Med. 2001 Mar;vol. 45(no. 3):521–524. [PubMed]
11. Chen N, Wyrwicz AM. Correction for EPI distortions using multi-echo gradient-echo imaging. Mag. Res. Med. 1999 Jun;vol. 41(no. 6):1206–1213. [PubMed]
12. Priest AN, Vita ED, Thomas DL, Ordidge RJ. EPI distortion correction from a simultaneously acquired distortion map using TRAIL. J. Mag. Res. Im. 2006 Apr;vol. 23(no. 4):597–603. [PubMed]
13. Windischberger C, Robinson S, Rauscher A, Barth M, Moser E. Robust field map generation using a triple-echo acquisition. J. Mag. Res. Im. 2004 Oct;vol. 20(no. 4):730–734. [PubMed]
14. Reber PJ, Wong EC, Buxton RB, Frank LR. Correction of off resonance-related distortion in echo-planar imaging using EPI-based field maps. Mag. Res. Med. 1998 Feb;vol. 39(no. 2):328–330. [PubMed]
15. Aksit P, Derbyshire JA, Prince JL. Three-point method for fast and robust field mapping for EPI geometric distortion correction. Proc. IEEE Intl. Symp. Biomed. Imag. 2007:141–144.
16. Schneider E, Glover G. Rapid in vivo proton shimming. Mag. Res. Med. 1991 Apr;vol. 18(no. 2):335–347. [PubMed]
17. Irarrazabal P, Meyer CH, Nishimura DG, Macovski A. Inhomogeneity correction using an estimated linear field map. Mag. Res. Med. 1996 Feb;vol. 35(no. 2):278–282. [PubMed]
18. McVeigh ER, Henkelman RM, Bronskill MJ. Noise and filtration in magnetic resonance imaging. Med. Phys. 1985 Sep;vol. 12(no. 5):586–591. [PMC free article] [PubMed]
19. Bernstein MA, Grgic M, Brosnan TJ, Pelc NJ. Reconstructions of phase contrast, phased array multicoil data. Mag. Res. Med. 1994 Sep;vol. 32(no. 3):330–334. [PubMed]
20. Leitao JMN, Figueiredo MAT. Absolute phase image reconstruction: a stochastic nonlinear filtering approach. IEEE Trans. Im. Proc. 1998 Jun;vol. 7(no. 6):868–882. [PubMed]
21. Jacobson MW, Fessler JA. An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms. IEEE Trans. Im. Proc. 2007 Oct;vol. 16(no. 10):2411–2422. [PMC free article] [PubMed]
22. Böhning D, Lindsay BG. Monotonicity of quadratic approximation algorithms. Ann. Inst. Stat. Math. 1988 Dec;vol. 40(no. 4):641–663.
23. Lange K, Hunter DR, Yang I. Optimization transfer using surrogate objective functions. J. Computational and Graphical Stat. 2000 Mar;vol. 9(no. 1):1–20.
24. Fessler JA, Yeo D, Noll DC. Regularized fieldmap estimation in MRI. Proc. IEEE Intl. Symp. Biomed. Imag. 2006:706–709.
25. Fessler JA, Rogers WL. Spatial resolution properties of penalized-likelihood image reconstruction methods: Space-invariant tomographs. IEEE Trans. Im. Proc. 1996 Sep;vol. 5(no. 9):1346–1358. [PubMed]
26. Unser M, Aldroubi A, Eden M. Recursive regularization filters: design, properties, and applications. IEEE Trans. Patt. Anal. Mach. Int. 1991 Mar;vol. 13(no. 3):272–277.
27. Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, Evans AC. Design and construction of a realistic digital brain phantom. IEEE Trans. Med. Imag. 1998 Jun;vol. 17(no. 3):463–468. [PubMed]
28. Kwan RK-S, Evans AC, Pike GB. MRI simulation-based evaluation of image-processing and classification methods. IEEE Trans. Med. Imag. 1999 Nov;vol. 18(no. 11):1085–1097. [PubMed]
29. Haacke EM, Brown RW, Thompson MR, Venkatesan R. Magnetic resonance imaging: Physical principles and sequence design. New York: Wiley; 1999.
30. Schenck JF. The role of magnetic susceptibility in magnetic resonance imaging: MRI magnetic compatibility of the first and second kinds. Med. Phys. 1996 Jun;vol. 23(no. 6):815–850. [PubMed]
31. Fessler JA, Sutton BP. Nonuniform fast Fourier transforms using min-max interpolation. IEEE Trans. Sig. Proc. 2003 Feb;vol. 51(no. 2):560–574.
32. Yip C, Fessler JA, Noll DC. Iterative RF pulse design for multidimensional, small-tip-angle selective excitation. Mag. Res. Med. 2005 Oct;vol. 54(no. 4):908–917. [PubMed]
33. Sutton BP, Noll DC, Fessler JA. Dynamic field map estimation using a spiral-in / spiral-out acquisition. Mag. Res. Med. 2004 Jun;vol. 51(no. 6):1194–1204. [PubMed]
34. Lu K, Liu TT, Bydder M. Optimal phase difference reconstruction: comparison of two methods. Mag. Res. Im. 2007 [PubMed]
35. Gamier SJ, Bilbro GL, Gault JW, Snyder WE. Magnetic resonance image restoration. J. Math. Im. Vision. 1995 Feb;vol. 5(no. 1):7–19.
36. Ying L, Sheng J, Liu B. Joint estimation of image and coil sensitivities in parallel MRI. Proc. IEEE Intl. Symp. Biomed. Imag. 2006:17–20.
37. Yu DF, Fessler JA. Edge-preserving tomographic reconstruction with nonlocal regularization. IEEE Trans. Med. Imag. 2002 Feb;vol. 21(no. 2):159–173. [PubMed]
38. Jacobson MW, Fessler JA. Comm. and Sign. Proc. Lab. Ann Arbor, MI: Dept. of EECS, Univ. of Michigan; 2004. Nov, Properties of MM algorithms on convex feasible sets: extended version. 48109-2122, Tech. Rep. 353.
39. Huber PJ. Robust statistics. New York: Wiley; 1981.
40. Erdoğan H, Fessler JA. Ordered subsets algorithms for transmission tomography. Phys. Med. Biol. 1999 Nov;vol. 44(no. 11):2835–2851. [PubMed]
41. Fessler JA, Clinthorne NH, Rogers WL. On complete data spaces for PET reconstruction algorithms. IEEE Trans. Nuc. Sci. 1993 Aug;vol. 40(no. 4):1055–1061.