1. Introduction
Magnetic Resonance Imaging (MRI) is a technique used in medical environment to produce high quality images of tissues inside human body. MRI is based on the principles of nuclear magnetic resonance (NMR), a spectroscopic technique used to obtain microscopic chemical and physical information about molecules. In MRI, a radiofrequency signal stimulates the hydrogen atoms of an object posed in a uniform magnetic field to emit complex-valued signals. These signals are coherently collected by MR coil, converted from analog to digital and recorded in the so called
k-space, providing MR raw data. Commonly known MR images are defined in the spatial domain, obtained by Fourier transforming the raw data [
1,
2].
MR images are characterized by the combination of three intrinsic parameters: the spin density of hydrogen atoms
ρ, the spin-lattice relaxation time T
1 and the spin-spin relaxation time T
2 [
2]. In particular, relaxation times can give useful information about local environment interaction and provide a major source of contrast in MR images, while, a quantitative analysis of spin-spin relaxation time T
2 can give useful information for cancer discrimination [
3].
Conventional relaxation parameter estimation techniques work on the amplitude of MR images. A commonly used approach consists in using Least Square (LS) estimator [
4]. This approach assumes additive Gaussian distributed noise on amplitude MR images. However, the noise probability density function (pdf) follows the Rice distribution, as amplitude MR images are computed from two independent Gauss distributed random processes (the real and imaginary parts) [
5–
7]. As far as we know, many papers present in literature consider the Rice distribution for the amplitude MR images [
8–
12], but only two papers, [
13] and [
14], have considered it for the T
2 estimation problem. In the first one [
13], although assuming Rice distribution for the model, the authors estimate spin-spin relaxation parameter still using LS. They justify the use of LS by stating that at high SNRs Rice distribution approaches a Gaussian pdf. This assumption, however, introduces a bias in the estimation, especially for low SNRs. In [
14], the authors propose a Maximum Likelihood Estimator (MLE) for retrieving relaxation parameters. Since the correct noise model is assumed, the proposed estimator is able to avoid the bias even at low SNRs in case of large data sets.
In this paper we propose an approach for spin-spin relaxation time estimation that works directly on the complex-valued MR images instead of amplitude. Working in complex domain allows us to implement LS estimation since the noise is Gaussian both on real and imaginary parts. Differently from [
13], we do not expect the estimator to be biased since we use the correct model. Compared to [
14], the proposed approach has two main advantages: first of all we exploit twice the available data, as the estimation is performed using the real and the imaginary parts. This allows us to use more information to reduce noise obtaining a more accurate estimation. Secondly, LS estimation ensures lower computational cost compared to the Rice based MLE proposed in [
14]. Note that the likelihood function in Rice case contains the Bessel function, which can be harder to be managed than an exponential function.
In Section 2 the statistical model is briefly addressed. In Section 3 the accuracy of the proposed model is evaluated exploiting Cramer Rao Lower Bounds (CRLB) and a comparison with other models present in literature is discussed. The performances of the proposed estimator are shown in Section 4. In Section 5, a fast version of the Non Linear LS estimator is proposed. Finally, we draw some conclusions about the presented technique.
2. Statistical Description of MR Images
In MRI the data are recorded in the
k-space, where they are corrupted by additive, zero mean and uncorrelated Gaussian noise samples [
15]. In order to obtain MR images in spatial domain, an inverse Fourier transform is required. Thanks to the linearity and orthogonality of this operation, complex-valued MR images are still corrupted by additive, zero mean and uncorrelated Gaussian noise in both real and imaginary parts. After the inverse Fourier transform, although the signal intensity and the noise variance could be modified, the ratio between the signal power and the noise power does not change.
Let us consider a pixel of a complex-valued MR image; in a noise free case, the signal expression is given by:
where
m represents the amplitude of the signal and
ϕ the phase. The amplitude
m depends on the used imaging sequence and on the unknown parameters
θ (the spin density of hydrogen atoms
ρ, the spin-lattice relaxation time T
1 and the spin-spin relaxation time T
2)
via:
while the phase
ϕ depends on the local intensity of the magnetic field (field map) [
16][
17].
Equation (1) can be alternatively written in terms of real and imaginary parts:
Let us now focus on the noisy case. As stated before, real (
yR) and imaginary (
yI) parts are corrupted by additive, zero mean and uncorrelated Gaussian noise
wR and
wI:
providing the following distributions:
where
σ2 represents noise variance. By multiplying the two equations
(5), we obtain the joint statistical distribution The proposed method is based on this statistical distribution of the measured data, called Gaussian Complex Model.
Other techniques work on the amplitude instead of the complex domain, assuming different distribution for the data [
14]. In particular, starting from
Equation (5), it is easy to show that the amplitude is corrupted by Rice distributed noise [
5–
7], leading to the following distribution:
where
I0(·) is the modified Bessel function of the first kind with order zero and
m̂ represents the measured noisy amplitude.
The term
f(
θ) is determined by the imaging sequence adopted during the MR scan. In this paper we consider a conventional spin echo imaging sequence which leads to [
2]:
where the instrumental variables T
E and T
R are the Echo Time and the Repetition Time of the sequence, respectively.
Since we are interested in T
2 estimation, we enclose T
1 and
ρ in a new parameter
k, called pseudodensity, as in [
13]. So
f(
θ) becomes:
which is commonly referred as monoexponential transversal magnetization decay model and is widely adopted [
13,
14]. In the next Section the achievable accuracy for T
2 estimation using the proposed model is addressed exploiting CRLB.
3. Cramer Rao Lower Bounds Evaluation
The Cramer Rao Lower Bounds provide the minimum variance for a given acquisition model that any unbiased (non polarized) estimator can reach [
18]. They provide a benchmark against which we can compare the performances of any unbiased estimator. Moreover, they alert us of the physical impossibility of finding an unbiased estimator with a variance smaller than the bounds. CRLB for amplitude estimation in MR Imaging have been presented in different papers [
8,
9]. In this section, a study on the CRLB for the specific T
2 estimation problem using the Gaussian Complex Model is conducted. To assess the performances of the proposed model, the CRLB are numerically evaluated using Monte Carlo method.
In the following simulation, a pixel with k = 80 and T2 = 100 msec is considered, corresponding to white matter tissue at 1.5 Tesla. The number of Monte Carlo iterations is fixed to 105.
First of all the dependency of CRLB on SNR is analyzed. SNR is defined as the ratio between the mean value of
f(
θ) and noise standard deviation
σ. A set of eight T
E equispaced in the interval [50 350] msec and a phase
ϕ = 60° are considered. The CRLB for this configuration are presented in the second column of
Table 1; for comparison, in the third column the CRLB for the same configuration but using the amplitude model (
Equation (6)) are reported. The CRLB behaviors for both acquisition models are also presented in
Figure 1.
Investigating the behavior of CLRB in
Table 1 and in
Figure 1, it can be noted that the Gaussian complex model
(5) is more accurate than the Rice amplitude
(6) one, in particular for low SNRs. As a matter of fact, for the considered SNR values range (SNR = [1–5]), the Minimum Variance Unbiased (MVU) estimator for the Rician model cannot reach the accuracy of the MVU estimator for the Gaussian case. For larger SNRs, the two models tend asymptotically to ensure the same accuracy.
The dependency of CRLB on the number of available acquisition (
i.e., the number of T
E) is then analyzed. The T
E values are equispaced in the interval [50 350] msec. The SNR is set to 3 and the phase
ϕ equal to 60°. The CRLB behaviors for both acquisition models are presented in
Figure 2.
Investigating the behavior of CLRB in
Figure 2, it can be noted that the number of acquisition has the same impact on both models. The behavior of the CRLB for the two considered models is the same, except for an offset. This offset is representative of the better accuracy achievable by the Gaussian model.
The CRLB evaluation can be useful not only to compare the achievable accuracy of different models, but also to investigate the best parameter configuration for a considered model.
In order to evaluate the impact of T
E values on the CRLB we perform two more analysis, using SNR = 4,
ϕ = 60° and eight acquisitions (
i.e., eight T
E values). The first one is related to the role of different spacing between T
E values. Such values are chosen in the interval [50 T
E_max] msec, varying T
E_max. The results of this simulation are shown in
Figure 3. As expected, increasing T
E_max (
i.e., increasing the spacing between T
E values) has a positive effect on the achievable accuracy. Anyway, saturation appears when T
E_max assumes high values.
Finally, the influence of phase
ϕ on the CRLB is investigated. The simulation has been conducted considering eight T
E values equispaced in the interval [50 350] msec, SNR = 4. The results, shown in
Figure 4, indicate that
ϕ does not substantially affect the behaviour of CRLB.
4. Non Linear Least Square Estimator: Results and Discussion
In this section, we present the estimator for the proposed model and we evaluate its performances in terms of mean and variance, compared to the CRLB, for different case studies.
Since we are working in the complex domain, as previously stated, signal statistical distributions for both real and imaginary parts are Gaussian. This assumption allows us to use a simple Least Square estimator. In particular we implement a Non Linear LS (NLLS) as the signal is characterized by an exponential decay
Equation (8). The NLLS is obtained by minimizing the square difference between the observed data and equation
(8) which contains the unknown parameter.
The presented results are obtained applying the method to realistically simulated data sets, using the same parameters of the previous section.
In the first case study the performances of the proposed estimator for different number of available complex-valued data sets
N are evaluated.
Figure 5(a) shows the trend of the mean value of the T
2 estimator for different number of T
E. As we can see the estimator results non-polarized since the bias is not present already with few data.
Figure 5(b) shows the behaviour of the variance of the T
2 estimator, compared to the CRLB. It can be noted that the estimator rapidly tends to be efficient, since its variance approaches the CRLB as the number of data grows.
For the second case study, we consider the performances of the T
2 estimator for different values of SNR. We fix the number of acquisitions to 8 (
N = 8). The performances of the estimator are satisfactory even at low SNRs. As it is shown in
Figures 6(a) and (b) the mean of the estimator tends to the true value and its variance approaches the CRLB already at low SNRs.
Finally the proposed algorithm has been tested on a 256 × 256 realistically simulated slice of the human head. The spin density and the relaxation times of the simulated tissues are reported in
Table 2, according to the Shepp-Logan phantom proposed in [
19]. The pseudodensity values
k have been obtained considering T
R = 5,000 msec.
Figure 7 shows the true spin density, the pseudodensity and the relaxation times maps of the considered tissues. Starting from these maps and fixing 8 T
E values equispaced in the interval [50 350] msec, 8 complex-valued images have been generated. Phase value has been set to 60° and noise has been added. The noise variance has been set equal to 3.5.
The results of the ML estimator applied to the Rician amplitude model and of the NLLS estimator with Gaussian complex model are reported in
Figures 8(a) and (b) and
Figures 8(c) and (d), respectively. In order to better appreciate the obtained results the reconstruction error map using both estimators have been reported in
Figure 8(e) and (f). All no signal areas have been masked after estimation. The simulation demonstrates the effectiveness of the proposed method and its better accuracy compared to the other considered method.
Concerning the computational cost of the algorithm, we have to underline that an optimization algorithm needs to iteratively compute a function to find the solution. Clearly, Gaussian Complex Model (
Equation (5)) shows a lower computational complexity compared to a Rician one (
Equation (6)), where the Bessel function has to be evaluated. This leads to a higher computational time for optimizing Rician distribution, which in some applications could be a problem. For the considered simulation, computational times were 660 sec for the LS estimator applied to the Gaussian model and 1,900 sec for the ML estimator applied to the Rician model using a SUN Ultra 40 workstation.
Figure 9 shows the reconstruction of the T2 map and the reconstruction error map for the same data set of
Figure 7, but with a noise variance of 1.5, showing the accuracy of the proposed method also with less noisy data; with this noise variance, the computation of the Bessel function in the Rician pdf becomes a not trivial task.
The better accuracy of the Gaussian Complex Model respect to the Rician one can be explained in terms of functions to be optimized. In the first case, the square error function between the observed data and
(7) (Least Square) has to be minimized. In the second one, the likelihood function, obtained by pdf
(6) after having observed the data, has to be maximized.
Figure 10 shows the normalized functions to be minimized in the Gaussian Complex Model square error function–
Figure 10(a), 10(b) and 10(c)] and in the Rician model [opposite likelihood function—
Figure 10(d)–(f)] simulating a pixel with k = 50 and T2 = 100, and considering eight TE values. As expected, changing the SNR from 1 [
Figure 10(a) and 10(d)] to 3 [
Figure 10(b) and 10(e)] and to 5 [
Figure 10(c) and 10(f)] has a positive effect on global minimum, allowing an easier minimization. Anyway in all cases the minimum related to the Rician model can be found with a lower accuracy respect to the Gaussian model case. This is due to the behavior of likelihood function
(6) which has a flat surface in the proximity of the minimum: even increasing the SNR from 1 to 5 the absolute minimum still remains located in a large and flat area.
5. Fast Non Linear Least Square Estimator
In Section 3 we presented the comparison between the Amplitude/Rice and Complex/Gaussian models, showing the best performances of the second one in both estimation accuracy and computational time. In this section we introduce a method to improve the speed of the second approach, making it quasi real time. This approach is based on Non Linear LS and is defined in [
18]:
where
y = [
y1R,
y1I,
y2R,
y2I,…
yNR, yNI]
T is the collection of the real and imaginary parts of measured data for the
N T
E values [T
E1, T
E2, …, TEN]
T, and:
where the components of
H(T
2) are reported twice in order to take into account both real and imaginary parts of
y.
The estimator of
Equation (9) allows us to perform a monodimensional minimization instead of a bidimensional one, greatly reducing the computational time. For example, spin-spin relaxation time estimation of 64 × 64 pixel image using eight acquisitions is performed in 71 seconds for classic Non Linear LS and in 10 seconds for Fast Non Linear LS estimator
Equation (9) with our workstation.
Figures 11(a) and (b) show the performances of this estimator varying the number of T
E values.
As we can see, in terms of performances, the fast estimator can be considered unbiased as the classic one while the accuracy is slightly worst compared to the classic one. This is due to the non linear relation between the observation matrix
H and T
2 in equation
(10). Anyway this loss of accuracy is balanced by the 7× speed improvements in computational time which can be very useful for quasi real time applications.
A final note about the possibility of estimating T
1 has to be marked. In the proposed approach we are interested in the estimation of T
2. Thus all our simulations and analysis are conducted considering only the spin-spin relaxation parameter. Anyway, it is easy to show that the proposed approach can be extended also for the estimation of T
1, expliciting
Equation (7) and
(8) in terms of spin-lattice relaxation parameter instead of spin-spin relaxation parameter.