Journal of Magnetic Resonance 249 (2014) 100–107

Contents lists available at ScienceDirect

Journal of Magnetic Resonance journal homepage: www.elsevier.com/locate/jmr

New approach to phase correction in multi-echo T 2 relaxometry q Marcus Björk ⇑, Petre Stoica Department of Information Technology, Uppsala University, Uppsala, Sweden

a r t i c l e

i n f o

Article history: Received 9 July 2014 Revised 27 September 2014 Available online 16 October 2014 Keywords: Phase correction Parameter estimation algorithm Multi-echo T2 relaxometry Multi-component T2 In-vivo brain data

a b s t r a c t Estimation of the transverse relaxation time, T 2 , from multi-echo spin-echo images is usually performed using the magnitude of the noisy data, and a least squares (LS) approach. The noise in these magnitude images is Rice distributed, which can lead to a considerable bias in the LS-based T 2 estimates. One way to avoid this bias problem is to estimate a real-valued and Gaussian distributed dataset from the complex data, rather than using the magnitude. In this paper, we propose two algorithms for phase correction which can be used to generate real-valued data suitable for LS-based parameter estimation approaches. The first is a Weighted Linear Phase Estimation algorithm, abbreviated as WELPE. This method provides an improvement over a previously published algorithm, while simplifying the estimation procedure and extending it to support multi-coil input. The algorithm fits a linearly parameterized function to the multi-echo phase-data in each voxel and, based on this estimated phase, projects the data onto the real axis. The second method is a maximum likelihood estimator of the true decaying signal magnitude, which can be efficiently implemented when the phase variation is linear in time. The performance of the algorithms is demonstrated via Monte Carlo simulations, by comparing the accuracy of the estimates. Furthermore, it is shown that using one of the proposed algorithms enables more accurate T 2 estimates; in particular, phase corrected data significantly reduces the estimation bias in multi-component T 2 relaxometry example, compared to when using magnitude data. WELPE is also applied to a 32-echo in vivo brain dataset, to show its practical feasibility. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Relaxometry provides quantitative tissue characterization and is useful in many branches of MRI, cf. [1,2]. A standard approach to estimating T 2 is based on acquiring several images at different echo times using a spin echo sequence, and fitting one or several decaying exponentials to the magnitude data in each voxel individually [3–6]. A single exponential decay is often assumed [7–9], but in practice, the data usually consists of several decaying components [10]. Multi-component T 2 relaxometry has application in, for example, myelin water imaging [11,12]; but it is also of general importance, as the single-component T 2 estimates can suffer from significant bias when the assumed model is wrong [10]. The magnitude images commonly used for T 2 relaxometry are Rice distributed [13], which means that the least squares (LS) parameter estimates will be suboptimal. A few suggestions on how to solve this sub-optimality problem are available in the q This work is partly supported by the European Research Council (ERC) Advanced Grant 247035. ⇑ Corresponding author at: Division of Systems and Control, Department of Information Technology, Uppsala University, Box 337, SE-751 05 Uppsala, Sweden. E-mail address: [email protected] (M. Björk).

http://dx.doi.org/10.1016/j.jmr.2014.09.025 1090-7807/Ó 2014 Elsevier Inc. All rights reserved.

literature [14,15], where the authors apply maximum likelihood (ML) methods, taking the Rice distribution into account; however, these approaches are nonlinear and rather involved, which complicates the implementation and can lead to convergence problems. To make the problem more tractable, many algorithms for T 2 estimation are based on minimizing the LS criterion; however, it has been suggested that LS-based approaches can lead to tissue mischaracterization caused by the Rician noise [16]. Another approach, that was applied in [17], estimates the signal decay and the noise properties from the magnitude images, and uses this information to transform the magnitude data into a Gaussian signal that can be used for T 2 estimation. This method is shown to improve the T 2 estimates when only the magnitude images are available. In this paper, we present two methods that, based on complexvalued data, compute Gaussian-distributed real-valued images, which can be used in LS-based T 2 estimation without inflicting any bias. By splitting the problem into two parts: one of phase correction to obtain an estimate of the underlying magnitude decay, and another of estimating the relaxation components, the highly nonlinear Rician estimation problem is avoided. Moreover, this enables increased accuracy when using the non-negative least squares (NNLS) algorithm [11,18], or the method called

101

M. Björk, P. Stoica / Journal of Magnetic Resonance 249 (2014) 100–107

Exponential Analysis via System Identification based on Steiglitz– McBride (EASI-SM), for multi-component exponential estimation [19]. It should be noted that the phase correction methods presented herein can be used for other problems than T 2 relaxometry, as only the phase of the data is modeled; however, this possibility will not be considered here. A method for temporal phase correction (TPC) has been proposed in [16], where the phase of the signal in an arbitrary voxel is fitted by two fourth-order polynomials using weighted LS. The estimated phase is then removed from the data, effectively projecting the signal onto the real axis. However, the TPC method does not make the best use of the data, and can be both simplified and generalized. Moreover, no simulations were performed in [16] to show the statistical performance of TPC. In the following, our aim is to statistically improve upon the method in [16], generalize it to any linearly parameterized phase function, and extend it to multi-coil data. Furthermore, we present a ML estimator of the true magnitude decay. This results in a nonlinear problem; however, when the phase changes linearly with time, the ML estimator can be efficiently implemented using the Fast Fourier Transform (FFT). This case is of special interest, as linear phase variations were observed in the in vivo data used both in this paper and in [16]. The main drawbacks of TPC from a signal processing perspective are: (1) the data is split into two parts to counter any even/odd-echo oscillations, which effectively reduces the number of samples used for estimation; (2) a suboptimal weighting based on the echo number is applied when solving the LS problem, which can be problematic, particularly for nonuniform echo spacing; and (3) estimating the parameters of high-order polynomials, such as the two fourth-order polynomials used in TPC, is rather sensitive to noise when only a few samples are available, which can lead to spurious oscillations in the phase. Typically, a simpler description of the phase is sufficient. The method presented herein extends to multi-coil acquisitions, and for each voxel, makes use of all the available samples to fit a single function to the phase of the data by weighted linear LS. The weights are chosen based on the magnitude of the measured signal, which directly relates to the variance of the phase noise. Using these weights in the LS problem corresponds to an approximation of the best linear unbiased estimator (BLUE) [20]. The proposed method also enables the use of any linearly parameterized function, rather than merely the set of polynomials, to fit the phase of the data. Moreover, the number of parameters can be chosen based on the observed variations in the data, or other prior information, avoiding the potential over fitting associated with always using a fourth-order polynomial. The resulting method is called Weighted Linear Phase Estimation (WELPE). The alternating phase phenomenon described in [16], which causes the phase to change sign from echo to echo, can be absorbed into the known parts of the equations. Doing so simplifies the algorithm, and makes it more robust compared to splitting the data into even and odd echoes and performing separate phase estimation for each of these datasets. In this work, we compare the phase correction performance of WELPE, ML, and TPC using simulated data; and illustrate the difference in accuracy when estimating multiple T 2 components using phase corrected data and magnitude data. Furthermore, WELPE is applied to a multi-echo spin-echo dataset comprising 32 in vivo brain images, to evaluate the practical feasibility of the algorithm. As the focus of this paper is to compare phase correction methods, and not approaches to estimate the model parameters, only EASI-SM will be used for parameter estimation. This also simplifies the bias comparison, as EASI-SM estimates discrete components rather than a continuous spectrum of T 2 values as, for example, NNLS does.

The paper is structured as follows: In the next section the signal and noise models are presented. Section 3 contains details on the two algorithms proposed, while Section 4 describes the simulation and data acquisition procedure. The estimation results for both simulated and in vivo data are shown in Section 5, and the discussion of the results follows in Section 6. Finally, some concluding remarks are found in Section 7. 2. Theory 2.1. Signal model For a multi-echo spin-echo acquisition, the received time domain signal from coil i, in an arbitrary voxel, can be modeled as a sum of complex exponentials

~si ðtn Þ ¼ ki ejPðtn Þ

M X

cm etn =T 2m þ ~i ðtn Þ ¼ g~i ðt n Þ þ ~i ðt n Þ;

n ¼ 1; . . . ; N;

m¼1

ð1Þ where ki 2 C is the coil sensitivity, M is the number of exponential components, cm 2 Rþ and T 2m 2 Rþ are thepamplitude and relaxaffiffiffiffiffiffiffi tion time of component m, respectively, j ¼ 1; Pðt n Þ is a function describing the phase change over time, ~ i ðtn Þ is zero mean complexvalued Gaussian white noise with variance 2r2 (i.e. r2 in both the real- and imaginary part), t n is the sampling time relative to excitation, and N is the number of samples (echoes). It is important to note that the factor ejPðtn Þ is common to all damped exponentials in a voxel, and all coils, as this phase is mainly attributed to field inhomogeneity and drift [13]. It is also possible to model the signal using a continuous distribution of T 2 values [21], rather than the discrete sum of components in (1). The phase-correction algorithms presented herein can be used with either of these models; however, to keep the description concise, the continuous model will not be considered here. As mentioned, the magnitude data commonly used in T 2 relaxometry is Rice distributed, and can be described by the following probability distribution function [8,9]:

pR ðj~si ðtn ÞjÞ ¼

j~si ðtn Þj

r2

exp 

!   j~si ðt n Þj2 þ g 2i ðt n Þ j~si ðt n Þjg i ðtn Þ I ; 0 r2 2r 2 ð2Þ

where g i ðt n Þ ¼ jg~i ðtn Þj is the true signal magnitude, and I0 is the modified Bessel function of the first kind and order zero. However, by estimating the phase function PðtÞ in (1), including the constant phase contribution from the coil, and removing it from the data, the real part of the resulting signal can be modeled as

si ðtn Þ ¼ jki j

M X cm etn =T 2m þ i ðt n Þ ¼ g i ðt n Þ þ i ðtn Þ;

ð3Þ

m¼1

where now, i ðt n Þ is a Gaussian distributed noise with variance r2 . Unlike the magnitude data, si ðt n Þ is suitable for use in any LS-based method, like NNLS or EASI-SM, without causing bias problems. In the next section two methods for estimating g i ðt n Þ from complexvalued data are suggested. These estimated samples will be Gaussian distributed, and can hence be described by the model in (3). In this paper Pðt n Þ is assumed to be linearly parameterized, that is

Pðtn Þ ¼ p0 þ

Q X pq wq ðt n Þ;

ð4Þ

q¼1

where fwq ðt n ÞgQq¼1 are a set of predefined basis functions, Q is the order of the basis expansion, and p0 captures any constant phase offset in the data. A wide class of phase variations are described by this parameterization, which should enable general use. For

102

M. Björk, P. Stoica / Journal of Magnetic Resonance 249 (2014) 100–107

example, wq ðt n Þ ¼ t qn gives the set of polynomials of order Q, as used in [16] with Q ¼ 4. In most cases, using a polynomial basis is sufficient, for example when no other assumption than smoothness of the phase can be made. Therefore, this is the general recommendation; however, given some prior information on the expected phase variation, other basis functions could be useful. The signal-to-noise ratio (SNR) is defined as

PN SNR ¼

n¼1 gðt n Þ

Nr

ð5Þ

;

that is, the mean intensity of the signal over the noise standard deviation. This is the common SNR definition for MR images; however, here it is used across the echoes. 3. Methods 3.1. WELPE

ð6Þ

where argfg is the phase (argument) of a complex number, /i is the phase of the coil sensitivity ki , and v i ðt n Þ is uncorrelated noise with time dependent variance. Note that /i can be absorbed into p0 , giving the coil dependent parameter p0i ¼ p0 þ /i . The phase noise, v i ðtn Þ, in (6) is random and cannot be unwrapped; therefore, as the wrapped phase is confined to the interval ½p; p; v i ðtn Þ will be wrapped-normal distributed rather than normal distributed [22]. Constructing the BLUE estimator of fpq gQq¼0 from (6) only requires knowledge of the variance of v i ðt n Þ for each tn , independent of the distribution [20]; however, BLUE assumes that the noise is zero mean, which does not hold for the wrapped-normal distribution at low SNR. It is possible to make use of circular statistics to avoid this problem, cf. [22], but for practical purposes it will be sufficient to use the variance of v i ðt n Þ to obtain an approximation of the BLUE. The model in (6) can be written as

3 p0i 6 7  6 p 1 7 7 ai ðtn Þ ¼ 1 w1 ðt n Þ . . . wQ ðtn Þ 6 6 .. 7 þ v i ðtn Þ: 4 . 5 2

ð7Þ

pQ Since only one parameter depends on the coil index i, there is a benefit in using the data from all coils to simultaneously estimate the model parameters. Using Nc coils, this reduces the number of unknowns to ðNc þ Q Þ, compared to ðN c þ QNc Þ when solving the problem for each coil separately; a significant decrease, particularly when Q > 1. By stacking the time samples and the basis functions into vectors, we can define ai ¼ ½ai ðt 1 Þ; . . . ; ai ðt N ÞT ; v i ¼ ½v i ðt1 Þ; . . . ; v i ðtN ÞT , and Wq ¼ ½wq ðt1 Þ; . . . ; wq ðtN ÞT , to obtain the multi-coil matrix model

0

   0 W1 .. . .. 1 . .. . .. .. .. . . . 0 0    0 1 W1

2p 3 01 36 . 7    WQ 6 .. 7 7 6 .. 7 7 6 p0Nc 7 . 7 76 7 þ v ¼ Rh þ v ; 6 7 .. 76 p1 7 7 . 56 6 . 7 6 . 7    WQ 4 . 5

ð8Þ

pQ T

ð9Þ

where b is a uniformly distributed random phase, and the dependence on i and t n has been dropped for notational convenience. As the phase P þ / is to be estimated, the disturbing phase is given by argfg þ jjejb g. It can be shown that, for high SNR, this phase disturbance can be approximated as v  Refg=g, where Refg indicates the real part, and hence, the variance becomes varðv Þ  r2 =g 2 [23]. Therefore, the BLUE weights at high SNR can be approximated as

ð10Þ

where the constant r has been omitted as it has no effect on the BLUE. However, as gðt n Þ is not known, we propose to use the magnitude of the noisy data from each coil i to approximate the quantity in (10), that is 2

ai ðt n Þ ¼ argf~si ðt n Þg ¼ Pðtn Þ þ /i þ v i ðt n Þ;

1 6 60 6 a¼6. 6. 4.

  ~s ¼ ejðPþ/Þ g þ jjejb ;

w ðtn Þ ¼ g 2 ðtn Þ;

The unwrapped angle of the model in (1) is given by

2

problem is linear in h, and by introducing weights proportional to the inverse of the variance of each sample, we can obtain the BLUE estimate of h [20]. A full derivation of the variance of the samples, and hence the weights, is mathematically rather involved and not particularly useful in practice; however, to provide some guidance for how to choose the weights, we rewrite (1) as

T

where a ¼ ½aT1 ; . . . ; aTNc  ; v ¼ ½v T1 ; . . . ; v TNc  ; R 2 RNNc ðNc þQ Þ , and 0; 1 2 RN1 are column vectors with all elements equal to zero or one, respectively. By construction, the corresponding estimation

wi ðtn Þ ¼ j~si ðt n Þj2 :

ð11Þ

For details on the accuracy of this approximation, particularly at lower SNR, see Section 5. The weighted LS solution can now be written as

^h ¼ ðRT WRÞ1 RT Wa;

ð12Þ

where the hat indicates an estimate, and the weight matrix W is given by

  W ¼ diag ½wT1 ; . . . ; wTNc  ;

ð13Þ T

where wi ¼ ½wi ðt 1 Þ; . . . ; wi ðt N Þ , and diagðxÞ is a diagonal matrix with the elements of the vector x along the diagonal. If the data is collected in a manner that yields alternating phases, as was the case in [16], this is easily detected and compensated for by conjugating every other sample of the data prior to running the estimation method. This procedure does not alter the noise properties, and is in fact equivalent to including the alternating sign in the model. Next, the estimated phase function is removed from the measured data, giving an estimate of the true signal magnitude:

n o ^i ¼ Re diagð~si ÞejRi ^hi ; g

i ¼ 1; . . . ; Nc ;

ð14Þ

hi is the part of (8) that corresponds to the estimated where Ri ^ T phase function for coil i; ~si ¼ ½~si ðt 1 Þ; . . . ; ~si ðt N Þ , and gi ¼ ½g i ðt 1 Þ; T . . . ; g i ðt N Þ . In the noise-free case, phase unwrapping over time is one dimensional, and therefore straight-forward, assuming that the phase varies smoothly. In practice, however, the phases of the low-SNR samples cannot be successfully unwrapped as these are practically noise. This unwrapping problem occurs when the noise intensity is of the same order of magnitude as the signal. In WELPE, this corresponds to a very small weight in the criterion, meaning that the impact on estimation is typically negligible. Moreover, unwrapping the phase noise does not affect the noise properties, and therefore, phase unwrapping can generally be applied to the full dataset without effecting the outcome of the estimation procedure. The order of the basis expansion, Q, and the basis functions, fwq gQq¼1 , can be chosen based on the data, and as mentioned, the set of polynomials of first or second order is often sufficient. It is

103

M. Björk, P. Stoica / Journal of Magnetic Resonance 249 (2014) 100–107

also possible to choose Q automatically on a voxel basis by utilizing an order selection method, for example, the Bayesian Information Criterion [24]. This is, however, beyond the scope of this paper. The fitting performance of WELPE depends on the SNR, the model order, and the choice of basis functions. For low-order polynomials the phase function can be accurately estimated at rather low SNR values. In particular, WELPE typically works at a much lower SNR than what is required for the application of multi-component T 2 estimation [6,25]. The estimation and phase correction procedure is summarized in Algorithm 1. Algorithm 1. Summary of WELPE for phase correction

e L¼2

N X 

2 Re ejPðtn Þ~sðtn Þ :

ð18Þ

n¼1

Using the fact that for any complex number z

½Refzg2 ¼

1 2 ½jzj þ Re z2 ; 2

ð19Þ

we can rewrite (18) as

e L¼

N h X

i j~sðt n Þj2 þ Re e2jPðtn Þ~s2 ðt n Þ

n¼1 

 ¼ const: þ hðp1 ; . . . ; pQ Þ cos arg hðp1 ; . . . ; pQ Þ  2p0 ;

ð20Þ

where we have defined

1: Inputs: Complex-valued data: ~sðtn Þ Sampling times: tn model order: Q P 1

hðp1 ; . . . ; pQ Þ ¼

! Q N X X ~s2 ðtn Þ exp 2j pq Wq ðt n Þ :

Basis functions: fwq ðtn ÞgQq¼1

ð21Þ

q¼1

n¼1

It now follows directly that the maximum with respect to p0 is given by

2: Construct the matrix in (8) 3: for each voxel do 4: Compute the unwrapped angle of the data (6) 5: Construct the matrix in (13) 6: Compute the WELPE parameter estimates (12) 7: Project data to the real axis (14) 8: end for

^0 ¼ p

1 ^1 ; . . . ; p ^Q Þg; argfhðp 2

ð22Þ

where the estimates of fpq gQq¼1 are

^q gQq¼1 ¼ argmax jhðp1 ; . . . ; pQ Þj: fp

ð23Þ

fpq gQ q¼1

Q

3.2. Maximum likelihood estimator In the case when only one coil is used, it is possible to derive a special form of the ML estimator of gðt n Þ. This method can also be used when multiple coil images have been merged into one image using a method that preserves the data model. For example, taking the mean of the complex-valued images across the coils, or combining the data corresponding to the maximum intensity voxels across the coils into one image, preserves the model in (1), while the sum-of-squares approach does not have this property. To derive the ML estimator, the true signal magnitude, gðt n Þ 2 Rþ is seen as a time varying amplitude, that is, no model of gðtn Þ is assumed. By rewriting (1) as

~sðt n Þ ¼ gðt n ÞejPðtn Þ þ ~ðtn Þ;

ð15Þ

where the coil index i has been dropped, and the coil sensitivity parameters have been absorbed into gðt n Þ and PðtÞ, the nonlinear least squares (NLS) estimates of gðt n Þ and Pðt n Þ can be obtained by minimizing the criterion function



N X ~sðt n Þ  gðtn ÞejPðtn Þ 2 n¼1

N n X 

2  jPðtn Þ 2 o ~sðtn Þ ¼ j~sðt n Þj2 þ gðt n Þ  Re ejPðtn Þ~sðt n Þ  Re e ; n¼1

ð16Þ where the second expression follows by expanding and completing the square. Estimates obtained by globally minimizing (16) are ML ^q gQq¼0 and the corunder the assumption of Gaussian noise. Given fp b n Þ, the estimate of gðt n Þ is immediresponding function estimate Pðt ately found to be

n o g^ðtn Þ ¼ Re ej bPðtn Þ~sðtn Þ :

ð17Þ

By substituting (17) into (16), it can be seen that the NLS estimate of fpq gQq¼0 can be obtained by maximizing the following function

^q gq¼1 from (23), and back-substituting into (22) to By computing fp obtain p0 , (17) can be used to obtain the ML estimate of gðtn Þ. Note that this approach does not require phase unwrapping, nor is there a need to use a potentially sub-optimal weighting in the criterion. For Q ¼ 1, finding a solution to (23) is a one-dimensional problem, and the maximization can be performed by brute force search on a grid of potential p1 values. Moreover, when the phase function Pðtn Þ can be modeled as a linear function of t n , that is

Pðtn Þ ¼ p0 þ p1 t n ;

ð24Þ

the expression in (21) is similar to a discrete Fourier transform, and for uniform echo spacing T s , (23) can be efficiently implemented using the FFT:

^1 ¼ p

1 argmax FFTK f~s2 ðt n Þg ; 2T s f 2½0; 2p

ð25Þ

where f corresponds to the normalized frequency in the transform, and K is the number of evaluated frequencies. Typically, K  10N (rounded to the closest power of two) provides sufficient accuracy in practical applications. For nonuniform echo times, the corresponding nonuniform FFT can be used, cf. [26]; alternatively, it is possible to precompute a nonuniform Fourier matrix based on grid of p1 values, and use this matrix at each voxel. For higher orders of the phase function, that is Q P 2, the estimates in (23) can be obtained by nonlinear maximization, given proper initialization; however, in such a case we can only expect convergence to a local maximum of jhðp1 ; . . . ; pQ Þj, and hence, the obtained estimates may be suboptimal. MATLAB implementations of both WELPE and the ML algorithm are available at: https://github.com/AAAArcus/PhaseCorrection. 4. Simulation and data acquisition As the phase variation observed in the in vivo data was approximately linear over time (see Section 5), the simulations will focus on this scenario. Inspired by [16], Monte Carlo (MC) simulations were performed at SNR ¼ 70, using datasets with 48 samples and linear alternating phase, where the first 32 echoes were spaced

104

M. Björk, P. Stoica / Journal of Magnetic Resonance 249 (2014) 100–107

by 10 ms, and the following 16 by 50 ms. However, to enable reliable estimation of more than two T 2 components at SNR ¼ 70, 128 echoes, uniformly spaced by 10 ms, were used in the T 2 estimation example. The data was generated using the model in (1), and complex-valued Gaussian white noise of variance 2r2 was added to achieve the appropriate SNR, as defined by (5). To enable direct comparison with TPC, only a single coil was simulated, and the coil sensitivity k was set to unity (without loss of generality). Furthermore, to obtain a high accuracy of the sample means and variances, 10,000 MC simulations were performed in all examples. The chosen parameter set used to generate the data in the simulation examples was

T 2 ¼ ½20; 80; 200 ms;

c ¼ ½0:4; 1; 0:1:

ð26Þ

Single slice, multi-echo spin-echo in vivo brain data was collected at Uppsala University Hospital using a 1.5 T Philips scanner equipped with a 8 element SENSE headcoil. For an echo spacing of 10 ms and TR ¼ 2200, 32 images with echo times ranging from 10 ms to 320 ms, were acquired. No averaging was used. An imaging matrix of 240  188 voxels was used to obtain a resolution of 1  1  5 mm. The total scan time was 7 min and 11 s. 5. Results 5.1. Simulation To illustrate the accuracy of the proposed criterion weights, given by (11), the variance of the phase noise was empirically approximated by the means of MC simulation. Noisy datasets were generated using the model parameters in (26), and the variance of the phase at each sampling instance was computed. The corresponding weights are shown in Fig. 1, along with the high SNR approximation, g 2 ðt n Þ, and the proposed square magnitude of a noisy realization given by (11). The 1=n weights used in TPC [16] were also included for comparison. As can be seen, the empirical weights are relatively well approximated by the proposed weights in this example (note the logarithmic scale). Moreover, (11) provides a better approximation to the empirical weights than g 2 ðtn Þ at later samples, that is, at lower SNR. TPC, on the other hand, gives a relatively low weight to the initial, most informative, samples, and too high weight to the noisy samples at the end of the decay. An example of a simulated dataset with linear alternating phase is shown in Fig. 2, in terms of the magnitude and phase of the noisy data. The generated data closely resembles the measurements from a white matter voxel shown in [16]. As can be seen, the phase

information is very noisy when the signal magnitude is low, even though the true phase is linear and alternating. To show the statistical improvement associated with using WELPE or ML, compared to TPC, when estimating the true magnitude decay, MC simulations were performed. The average of the estimated phase-corrected signals are shown in Fig. 3. As can be seen, both WELPE and ML are on the average close to the true noise-free magnitude decay, while TPC gives a signal that is essentially more similar to the magnitude of the noisy data. To illustrate the multi-component T 2 estimation performance based on WELPE, TPC, and magnitude data, MC simulations were performed using the parameters in (26). The root mean square errors (rMSE) of the T 2 estimates, obtained using EASI-SM on the different datasets, are listed in Table 1. As can be seen, the WELPE phase-corrected data provides the highest parameter estimation accuracy, followed by TPC, and the magnitude data. The estimates from the first 200 MC simulations, plotted in the fc; Tg-plane, are shown in Fig. 4, together with the true parameter values. As can be seen, using the magnitude data causes a bias in the estimates, particularly for the slowest T 2 component. TPC reduces this bias, but the estimates are slightly more spread out compared to when using the data corrected by WELPE. Moreover, TPC leads to a few outliers in the parameter estimates, shown as squares clearly separated from the main cluster of estimates. ML provided similar results to WELPE, but was omitted from Fig. 4 for clarity. 5.2. In-vivo To reduce the computation time, the in vivo images were initially masked by excluding voxels below 20% of the maximum intensity voxel. This masked region corresponds to the low signal background of the images, mainly located outside of the skull. The magnitude and phase of the 80 ms echo from the collected

(a)

(b)

Fig. 1. An example of the LS weights given by the inverse of the variance of the phase (approximate BLUE weights), the proposed weights obtained directly from the data, and the TPC weights given by the inverse of the echo number.

Fig. 2. (a) Magnitude and (b) phase (rad) over time, for a simulated noisy dataset where the true phase is linear and alternating.

M. Björk, P. Stoica / Journal of Magnetic Resonance 249 (2014) 100–107

105

are approximately Gaussian distributed with zero mean and a relatively small standard deviation. 5.3. Computational details

Fig. 3. The average estimated magnitude decay g^ðt n Þ from 10,000 MC simulations for WELPE, ML, TPC, and the magnitude data, compared to the true decay gðtn Þ.

Table 1 Root mean square errors (rMSE) of the fT 2m g3m¼1 estimates using magnitude, WELPE and TPC data, computed from 10,000 MC simulations. The simulated data consisted of 128 echoes with uniform spacing of 10 ms, and SNR ¼ 70. Methods

Magnitude WELPE TPC

rMSE (ms) T 21

T 22

T 23

1.1 0.5 1.2

4.9 1.6 3.6

87.0 13.3 16.4

Fig. 4. Estimates of c and T 2 from 200 MC simulations, obtained by applying EASISM to phase-corrected data generated by WELPE and TPC, as well as to the magnitude of the noisy data. The true parameter values are indicated by the stars. The estimates correspond to the first 200 realizations of the simulations in Table 1.

32 echo dataset is shown in Fig. 5, together with the estimated magnitude and phase using WELPE, and the corresponding differb As can be seen, the magnitude ence images (j~sj  g^ and argf~sg  P). image and the WELPE corrected data are close to one another. Moreover, the phase is accurately modeled as no major structure remains in the difference image, and a clear denoising effect is observed in the estimated phase image. An example of the WELPE fitting of the phase for a single voxel is shown in Fig. 6. The phase initially varies approximately linearly, while as the signal magnitude decays, the phase noise grows larger, leading to a more random behavior. Because of the weighting in the criterion, WELPE is able to fit to the initial linear trend, while effectively disregarding the heavily distorted samples from the later echoes. In Fig. 7, a histogram of the imaginary part of the phase corrected image from Fig. 5 is shown. As can be seen, the samples

The computational burden of the considered methods varies depending on the settings and the implementation; however, we can still give an indication of the run times to be expected. The computation times for all methods when performing 10,000 MC simulations on a single processor thread using an Intel i7 860 at 2.93 GHz, are shown in Table 2. The simulated data consisted of 128 echoes with uniform spacing of 10 ms, similar to the example in Table 1. A uniform sampling was used to be able to show the computational performance of the ML estimator based on the FFT implementation. WELPE is the fastest method, running approximately twice as fast as TPC. ML has a slightly longer runtime than WELPE for K ¼ 1024, but is still more computationally efficient than TPC in this case. 6. Discussion 6.1. Simulation The proposed weights used in the WELPE criterion are data adaptive, provide a good approximation of the empirical BLUE weights, are easily computed, and lead to superior statistical performance compared to the 1=n weights used in TPC. The simulated data with linear and alternating phase displayed in Fig. 2 showed that a linear phase model, rather than a fourthorder polynomial, would have been more suitable to model the phase of the data in [16]. This motivates both WELPE, which has the option to choose the number of parameters used to model the phase based on the data, and the ML method, which provides a simple FFT-based solution for the case of linear phase variations. As was seen in Fig. 3, the TPC data is on the average rather similar to the magnitude data, in this example with non-uniformly spaced echoes. This is partially due to the suboptimal weights used, which puts too high weight to the noisy samples at the end of the decay. The data adaptive weights used by WELPE, on the other hand, provide an accurate estimate of the true signal magnitude. The jump discontinuity shown in the last two samples of the TPC estimate in Fig. 3 is due to robustness issues in the phase fitting. The two fourth-order polynomials used in TPC usually fits the initial linear phase, where the SNR is high, quite well. However, as a consequence of the high order model, several extreme points of the polynomials are placed in the low SNR region, essentially fitting to the noise. This can cause rapid variations in the estimated phase function, and in turn, can lead to an increase in the bias. As was shown in Fig. 4, using the magnitude of the noisy data and an LS-based approach can lead to a significant bias in the T 2 estimates, and therefore, mischaracterization of the tissue. This can be expected, as the Rician noise will effectively raise the tail of the decay curve, which is interpreted as a slower decay. By applying WELPE to correct the phase, however, the data becomes Gaussian distributed, which eliminates the bias problem. The TPC approach also improves the T 2 estimation significantly, but results in a higher rMSE than WELPE, and also causes outliers in the estimates. The outliers are partially due to the robustness issues in the polynomial fitting mentioned above, which for some noise realizations can lead to poor estimates of the true phase function. 6.2. In-vivo For the collected 32 echo in vivo dataset, WELPE is able to accurately model the phase in the whole image, and provides an

106

M. Björk, P. Stoica / Journal of Magnetic Resonance 249 (2014) 100–107

Fig. 5. (a) Magnitude, and (b) phase (rad) at echo time 80 ms of the single slice in vivo dataset; and the estimated (c) magnitude, and (d) phase, provided by WELPE, together with the error in (e) magnitude (amplified by a factor of 1000), and (f) phase (amplified by a factor of 100).

Fig. 6. Phase (rad) over time for a single voxel of the in vivo dataset together with the linear fit provided by WELPE.

Fig. 7. Distribution of imaginary part of the phase corrected image at echo time 80 ms.

M. Björk, P. Stoica / Journal of Magnetic Resonance 249 (2014) 100–107 Table 2 MATLAB run times in seconds for phase correcting 10,000 datasets using WELPE, ML, and TPC, at different SNRs. Methods

WELPE ML (K ¼ 1024) TPC

SNR 50

70

150

1.2 1.9 2.8

1.2 1.9 2.7

1.2 1.8 2.7

estimate of the true magnitude decay that is close to the magnitude of the noisy data for high SNR, as was shown in Fig. 5. Moreover, Fig. 6 showed that the linear fit to the phase variation in time was relatively good, and that it would be hard to motivate the fitting of a higher order phase function in this example. The histogram in Fig. 7, showed that the imaginary part of the data after phase correction is small and approximately Gaussian distributed. This would be expected in the case of successful phase correction, and indicates that WELPE is able to project a large proportion of the signal energy to the real axis, leaving mainly noise in the imaginary part. 6.3. Computational details Since TPC requires two systems of equations to be solved and the results recombined, the computation times for WELPE are typically lower. The FFT implementation of the ML estimator is also fast for K  10N, which provides reasonable accuracy. Note that it is possible to trade off computation time and accuracy by setting the number of evaluation points, K, in the FFT. 7. Conclusion Two methods for phase correction have been presented, and through simulations, the algorithms have been shown to be useful for avoiding bias in multi-component T 2 estimation by accurately estimating the true magnitude decay. WELPE is statistically sound and is easy to implement; moreover it works with multi-coil data, general sampling schemes, and a wide range of phase functions. The ML estimator is optimal, and does not require phase unwrapping or weighting of the criterion. Moreover, it can be rather efficiently implemented for linear phase variations in time; however, in the general case, finding the ML estimates is computationally more intensive than using WELPE. Acknowledgments Dr. Joel Kullberg and Anders Lundberg are gratefully acknowledged for their assistance during the image acquisition stage. References [1] P. Tofts, Quantitative MRI of the Brain: Measuring Changes Caused by Disease, Wiley, Chichester, UK, 2003. [2] H.-L. Margaret Cheng, N. Stikov, N.R. Ghugre, G.A. Wright, Practical medical applications of quantitative MR relaxometry, J. Magn. Reson. Imaging 36 (4) (2012) 805–824, http://dx.doi.org/10.1002/jmri.23718.

107

[3] E.L. Hahn, Spin echoes, Phys. Rev. 80 (1950) 580–594, http://dx.doi.org/ 10.1103/PhysRev.80.580. [4] H.Y. Carr, E.M. Purcell, Effects of diffusion on free precession in nuclear magnetic resonance experiments, Phys. Rev. 94 (1954) 630–638, http:// dx.doi.org/10.1103/PhysRev.94.630. [5] E.M. Haacke, R.W. Brown, M.R. Thompson, R. Venkatesan, Magnetic Resonance Imaging: Physical Principles and Sequence Design, Wiley, New York, NY, USA, 1999. [6] E. Alonso-Ortiz, I.R. Levesque, G.B. Pike, Mri-based myelin water imaging: a technical review, Magn. Reson. Med. (2014), http://dx.doi.org/10.1002/ mrm.25198. [7] M. Björk, P. Stoica, Fast denoising techniques for transverse relaxation time estimation in MRI, in: Proc. 21st European Signal Processing Conference (EUSIPCO-2013), No. 1569743663, Marrakech, Morocco, 2013. . [8] J. Sijbers, A.J. den Dekker, M. Verhoye, E. Raman, D. Van Dyck, Optimal estimation of T2 maps from magnitude MR images, Proc. SPIE Med. Imaging 3338 (1998) 384–390. [9] J.-M. Bonny, M. Zanca, J.-Y. Boire, A. Veyre, T2 maximum likelihood estimation from multiple spin-echo magnitude images, Magn. Reson. Med. 36 (2) (1996) 287–293, http://dx.doi.org/10.1002/mrm.1910360216. [10] K.P. Whittall, A.L. MacKay, D.K. Li, Are mono-exponential fits to a few echoes sufficient to determine T2 relaxation for in vivo human brain?, Magn Reson. Med. 41 (6) (1999) 1255–1257, http://dx.doi.org/10.1002/(SICI)15222594(199906)41:63.0.CO;2-I. [11] A. Mackay, K. Whittall, J. Adler, D. Li, D. Paty, D. Graeb, In vivo visualization of myelin water in brain by magnetic resonance, Magn. Reson. Med. 31 (6) (1994) 673–677, http://dx.doi.org/10.1002/mrm.1910310614. [12] I.R. Levesque, P.S. Giacomini, S. Narayanan, L.T. Ribeiro, J.G. Sled, D.L. Arnold, G.B. Pike, Quantitative magnetization transfer and myelin water imaging of the evolution of acute multiple sclerosis lesions, Magn. Reson. Med. 63 (3) (2010) 633–640, http://dx.doi.org/10.1002/mrm.22244. [13] R.M. Henkelman, Measurement of signal intensities in the presence of noise in mr images, Med. Phys. 12 (2) (1985) 232–233, http://dx.doi.org/10.1118/ 1.595711. [14] J. Sijbers, A.J. den Dekker, E. Raman, D. Van Dyck, Parameter estimation from magnitude MR images, Int. J. Imaging Syst. Technol. 10 (2) (1999) 109–114. [15] S. Walker-Samuel, M. Orton, L.D. McPhail, J.K.R. Boult, G. Box, S.A. Eccles, S.P. Robinson, Bayesian estimation of changes in transverse relaxation rates, Magn. Reson. Med. 64 (3) (2010) 914–921, http://dx.doi.org/10.1002/mrm.22478. [16] T.A. Bjarnason, C. Laule, J. Bluman, P. Kozlowski, Temporal phase correction of multiple echo T2 magnetic resonance images, J. Magn. Reson. 231 (2013) 22– 31, http://dx.doi.org/10.1016/j.jmr.2013.02.019. [17] R. Bai, C.G. Koay, E. Hutchinson, P.J. Basser, A framework for accurate determination of the T2 distribution from multiple echo magnitude MRI images, J. Magn. Reson. 244 (2014) 53–63. http://dx.doi.org/10.1016/ j.jmr.2014.04.016. [18] C.L. Lawson, R.J. Hanson, Solving Least Squares Problems, Prentice Hall, Englewood Cliffs, NJ, USA, 1974. [19] P. Stoica, P. Babu, Parameter estimation of exponential signals: A system identification approach, Digital Signal Process. 23 (5) (2013) 1565–1577, http://dx.doi.org/10.1016/j.dsp.2013.05.003. [20] T. Söderström, P. Stoica, System Identification, Prentice Hall, Cambridge, UK, 1989. [21] R.M. Kroeker, R.M. Henkelman, Analysis of biological NMR relaxation data with continuous distributions of relaxation times, J. Magn. Reson. 69 (2) (1986) 218–235, http://dx.doi.org/10.1016/0022-2364(86)90074-0. [22] K.V. Mardia, P.E. Jupp, Directional Statistics, Vol. 494, John Wiley & Sons, 2009. [23] S. Tretter, Estimating the frequency of a noisy sinusoid by linear regression, IEEE Trans. Inf. Theory 31 (6) (1985) 832–835, http://dx.doi.org/10.1109/ TIT.1985.1057115. [24] P. Stoica, Y. Selen, Model-order selection: a review of information criterion rules, IEEE Signal Process. Mag. 21 (4) (2004) 36–47, http://dx.doi.org/10.1109/ MSP.2004.1311138. [25] S.J. Graham, P.L. Stanchev, M.J. Bronskill, Criteria for analysis of multicomponent tissue T2 relaxation data, Magn. Reson. Med. 35 (3) (1996) 370–378, http://dx.doi.org/10.1002/mrm.1910350315. [26] J.A. Fessler, B.P. Sutton, Nonuniform fast fourier transforms using min–max interpolation, IEEE Trans. Signal Process. 51 (2) (2003) 560–574, http:// dx.doi.org/10.1109/TSP.2002.807005.

New approach to phase correction in multi-echo T2 relaxometry.

Estimation of the transverse relaxation time, T2, from multi-echo spin-echo images is usually performed using the magnitude of the noisy data, and a l...
1MB Sizes 2 Downloads 4 Views