PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Magn Reson. Author manuscript; available in PMC 2010 April 1.
Published in final edited form as:
PMCID: PMC2765718
NIHMSID: NIHMS83011

A Signal-Transformational Framework for Breaking the Noise Floor and Its Applications in MRI

Abstract

A long-standing problem in Magnetic Resonance Imaging (MRI) is the noise-induced bias in the magnitude signals. This problem is particularly pressing in diffusion MRI at high diffusion-weighting. In this paper, we present a three-stage scheme to solve this problem by transforming noisy nonCentral Chi signals to noisy Gaussian signals. A special case of nonCentral Chi distribution is the Rician distribution. In general, the Gaussian-distributed signals are of interest rather than the Gaussian-derived (e.g., Rayleigh, Rician, and nonCentral Chi) signals because the Gaussian-distributed signals are generally more amenable to statistical treatment through the principle of least squares. Monte Carlo simulations were used to validate the statistical properties of the proposed framework. This scheme opens up the possibility of investigating the low signal regime (or high diffusion-weighting regime in the case of diffusion MRI) that contains potentially important information about biophysical processes and structures of the brain.

Keywords: MR noise floor, signal transformation, Rician, nonCentral Chi, Gaussian, MRI

I. INTRODUCTION

Magnetic resonance imaging (MRI) (1) is a rapidly expanding field and a widely used medical imaging modality—possessing many noninvasive techniques capable of probing functional activities (2) and anatomical structures (310) of the brain in vivo. In quantitative MRI, important parameters of biophysical relevance are typically estimated from a collection of MR signals that are related to one another through a function of one or more experimentally controlled variables. As ever higher sensitivity and specificity to biophysical processes are achieved in MRI through improved spatial or temporal resolution, the adverse effect of noise on the overall accuracy of MRI-based quantitative findings also increases.

MR signals are complex numbers where the real and imaginary components are independently Gaussian distributed (11). The phase of the complex MRI signal is highly sensitive to many experimental factors, e.g., see (11,12), and as such, the magnitude of the complex MR signal is used instead in most quantitative studies. Although several techniques have been proposed to correct the phase error (1215), the magnitude of the complex MR signal (hereafter, magnitude MR signal) remains the most commonly used measure in MRI. While the magnitude MR signal is not affected by the phase error, it is not an optimal estimate of the underlying signal intensity when the signal-to-noise ratio is low (11) because it follows a nonCentral Chi distribution (16,17) rather than a Gaussian distribution. We should note that the Rician distribution (18,19) is a special case of the nonCentral Chi distribution. It is also well known that a Rician distribution (20) reduces to a Rayleigh distribution when the underlying signal intensity is zero, and the first moment of a Rayleigh distribution is usually known as the “noise floor” (21).

It is increasingly apparent that a resolution of the noise-induced bias in the magnitude MR signals could make it possible to gain further insights into the low signal regime that contains potentially important information about intrinsic functional activity (22) and tissue microstructure (39). Although several correction methods have been proposed (11,16,19,23,24) to address this problem, these methods do not produce corrected data that are Gaussian distributed.

A simple means of assessing Gaussianity in the corrected data when the noisy magnitude signals are drawn from the same distribution, e.g., see Figure 1, is to check if the corrected data follow a Gaussian distribution. In practice, this type of data is rare. Rather, we usually have MRI data that are drawn from a family of distributions all of which are characterized by different location parameters (e.g., the location parameter of a Gaussian distribution is the first moment and the location parameter of a nonCentral Chi distribution will be pointed out later). For example, each of the noisy magnitude signals of interest may be acquired under a slightly different experimentally controlled condition so that each noisy magnitude signal is actually drawn from a slightly different distribution. The proposed scheme is the first method capable of obtaining corrected data that are distributed evenly in both the positive and negative axes when the signal-to-noise ratio is very close to zero, which is a very important but simple criterion for testing the accuracy or lack thereof of a correction scheme. We should point out that none of the previously published methods (11,16,19,23,24) satisfies this criterion because these methods cannot produce corrected data that have negative values.

Fig. 1
(A) A schematic diagram of the proposed scheme. (B) A schematic diagram of the two possible approaches that can be used to map nonCentral Chi signals to Gaussian signals. The approach using the samples that are drawn from the same distribution (the bottom-left ...

In this work, we present a framework for making the magnitude signals Gaussian-distributed. A simple example illustrates the idea behind the proposed framework: suppose the noisy magnitude signals are drawn from a family of nonCentral Chi distributions all of which are characterized by different location parameters but with the same scale parameter. The proposed framework attempts to transform the noisy magnitude signals such that each noisy transformed signal may be thought of as if it were drawn from a Gaussian distribution with a different mean but the same standard deviation. Note that the location and scale parameters that characterize a nonCentral Chi distribution are exactly the mean and the standard deviation of the Gaussian distribution that characterize the transformed signal.

Three important considerations will have to be taken into account in order to construct such a framework. First, we need a method that can find an estimate of the first moment of a nonCentral Chi distribution from which the datum is drawn. Second, we need a method that can find an estimate of the first moment of the Gaussian distribution if an estimate of the first moment of a nonCentral Chi distribution is provided. Third, we need a method that can find a noisy Gaussian-distributed signal for each of the magnitude signals if the first moment of the nonCentral Chi distribution, the first moment and the standard deviation of the Gaussian distribution are provided. Each consideration above constitutes a separate procedure or stage.

Therefore, it is necessary to have a procedure in the first stage that can find an “average value” for each datum. In other words, the first moment of a nonCentral Chi distribution from which the datum is drawn is estimated in the first stage. Once an estimate of the first moment of a nonCentral Chi distribution is known, a procedure in the second stage must be able to produce the “average value” of the underlying signal intensity, which is an estimate of the first moment of a Gaussian distribution. A procedure in the third stage must be able to use each original noisy datum, which is nonCentral Chi-distributed, to find the corresponding transformed noisy signal that is Gaussian-distributed. The schematic representation of the three stages of the proposed framework is shown in Figure 1A.

Specifically, in the first stage, a data smoothing or fitting method may be used to obtain the average values of the noisy magnitude signals. The data may be fitted with some parametric functions (single exponentially or bi-exponentially decaying functions) or smoothed with a variety of smoothing methods. Although a comparison of various fitting or smoothing methods is of interest, such a comparison, if thoroughly investigated, would take us too far afield. Here, we use a penalized or smoothing spline model (25,26), to obtain the “average values”. The penalized spline model is chosen for its ease of implementation and use. The degree of smoothness is selected based on the method of generalized cross-validation (GCV) (26,27). Again, other methods may be used to select the degree of smoothness, see e.g., (28).

In the second stage, we propose an iterative method that takes in an “average value” of a noisy magnitude signal as an input and returns an “average value” of the underlying signal intensity as an output. This iterative method is closely related to but different from our previously proposed fixed point formula of the signal-to-noise ratio (SNR) because it is a fixed point formula of the underlying signal intensity, see Figure 1B. Specifically, the present iterative method treats the estimations of the underlying signal intensity and of the Gaussian noise standard deviation (SD) separately rather than simultaneously. The key advantage of such an approach is that there exists excellent methods for estimating the Gaussian noise SD from a much larger sample (29,30). Consequently, a more precise estimate of the Gaussian noise SD will result in a more precise estimate of the underlying signal intensity.

In the third stage, the corresponding noisy Gaussian signal of each of the noisy magnitude signals is found through a composition of the inverse cumulative probability function of a Gaussian random variable and the cumulative probability function of a nonCentral Chi random variable. Both the inverse cumulative probability function of a Gaussian random variable and the cumulative probability function of a nonCentral Chi random variable depend on the “average value” of the underlying signal intensity and the Gaussian noise SD. The third stage is exactly a Gaussian random number generator if the input data are Rician-distributed.

The statistical properties of the proposed framework is investigated using Monte Carlo simulations. Experimental data is also used to illustrate the proposed framework.

II. METHODS

Since the first stage of the proposed scheme is readily available (25,26), our focus in this paper will be on the latter stages. For completeness and notational consistency, we have included a brief discussion of one-dimensional penalized splines in Appendix A, and of spherical harmonics splines in Appendix B. These spline models share the same matrix structure, and therefore, the computation of this matrix structure is briefly touched on in Appendix C.

A. Theoretical preliminary

The probability density function (PDF) and the cumulative distribution function (CDF) of a nonCentral Chi random variable, m, are needed respectively in the second and third stages of the proposed scheme. It is known that magnitude MR signals obtained from an N-receiver-coil MRI system follow a nonCentral Chi, [chi], distribution of 2N degrees of freedom and the corresponding PDF can be expressed as (16,17):

pχ(mη,σg,N)dm=mNσg2ηN1exp(m2+η22σg2)IN1(mησg2)dm,m0
(1)

where the PDF is zero when m < 0, η is the underlying (combined) signal intensity (also known as the location parameter of the nonCentral Chi distribution), σg is the Gaussian noise standard deviation, and Ik is the kth-order modified Bessel function.

The corresponding (CDF) can be expressed as:

Pχ(αη,σg,N)=0αpχ(mη,σg,N)dm.
(2)

In practice, it is more convenient to compute Eq. (2) in terms of series representations of the generalized Marcum-Q function (31), QN. It can be shown that Eq. (2) can be simplified to:

Pχ(αη,σg,N)=1αpχ(mη,σg,N)dm=1QN(η/σg,α/σg),
(3)

where the definition of the generalized Marcum-Q function is:

QN(λ,γ)=1λN1γsNexp(λ2+s22)IN1(λs)ds.
(4)

When the underlying signal is zero, i.e., η = 0, the PDF and the CDF are given by (30):

pχ(m0,σg,N)dm=m2N12N1σg2N(N1)!exp(m22σg2)dm,
(5)

and

Pχ(α0,σg,N)dm=11(N1)!Γ(N,α2/(2σg2)),
(6)

where the incomplete Gamma function is defined as Γ(N,x)=xtN1exp(t)dt. The complete Gamma function is Γ(N,0) and is typically written simply as Γ(N).

B. Fixed point formula of the underlying signal intensity

The derivation of the fixed point formula of the underlying signal intensity, η, which is needed in this work, is closely related to that of the fixed point formula of the signal-to-noise ratio, θ [equivalent] η/σg, shown in our previous work (16). The main difference is in the separation of the underlying signal intensity and the Gaussian noise SD, σg. This separation is conceptually very important because the Gaussian noise SD, σg, is held as a fixed constant during the iterative process of successively estimating η. The key advantage of this approach is that the precision in the estimate of η is higher than its counterpart through our previous approach because the estimate of σg computed from a much large sample is less variable.

Here, we present the derivation of the fixed point formula of η. We begin with the first two moments of a nonCentral Chi distribution, Eq. (1), and they are given by:

m=σgβNF11(1/2,N,η2/(2σg2)),
(7)

and

m2=η2+2Nσg2,
(8)

respectively, where βN=π/2(2N1)!!2N1(N1)!, the double factorial is defined as: n!!= n(n −2)(n −4) × …, and 1F1 is the confluent hypergeometric function.

The variance of a nonCentral Chi random variable is defined as:

σχ2m2m2=ξ(ησg,N)σg2,
(9)

where the scaling factor, ξ, is given by:

ξ(ησg,N)=2N+η2σg2[βNF11(1/2,N,η2/(2σg2))]2.
(10)

The fixed point formula of the underlying signal intensity can be obtained by substituting the expression in Eq. (8) into Eq. (9). This leads to the following expressions:

η=g(ηm,σg,N)m2+[ξ(ησg,N)2N]σg2.
(11)

Note that the implementation of the fixed point formula of η, which is based on Newton’s method of root finding and is described in Appendix D, has important differences compared to that of the fixed point formula of θ [equivalent] η/σg (16).

To find the fixed point estimate, denoted by [eta w/ hat], in Eq. (11), left angle bracketmright angle bracket and σg are replaced by their corresponding estimates, denoted by m and [sigma with hat]g, respectively. In general, m may be taken to be the smoothed estimate obtained from the smoothing spline and [sigma with hat]g may be taken to be the estimated Gaussian noise SD obtained through various techniques mentioned above (11,29,30,32,33).

In short, the fixed point formula maps m to [eta w/ hat]. Fixed point formulae are powerful methods of successive approximation because their convergence can be tested under a very simple and general assumption (34). Specifically, let [eta w/ hat] be the fixed point that satisfies Eq. (11), i.e., g([eta w/ hat] | m, [sigma with hat]g, N) = [eta w/ hat], and η0 be the initial approximation, if both [eta w/ hat]0 and [eta w/ hat] belong to an interval in which dg(ηm^,σ^g,N)dη<1 for all η in that interval then η0 will always converge to [eta w/ hat]. In the context of Eq. (11), the interval does not have to be specified because [eta w/ hat] is always less than m, and therefore, the iteration with η0 = m as the initial approximation will always be convergent if m is exactly at or above the level of the noise floor, i.e., m ≥ βN[sigma with hat]g. Due to statistical variations, however, m may occasionally be below the noise floor. If this is the case, then [eta w/ hat] is set to −[eta w/ hat]δ where [eta w/ hat]δ is the fixed point estimate derived from a new estimate of left angle bracketmright angle bracket, which is defined by mδ = βN[sigma with hat]g + δ with δ = βN[sigma with hat]gm. This particular choice is needed to ensure the symmetry of the resultant distribution of [eta w/ hat] at a zero signal-to-noise ratio. Finally, an implementation of the fixed point formula is provided in Appendix D.

C. Mapping nonCentral Chi to Gaussian signals

Mapping a nonCentral Chi random variable, m, to a Gaussian random variable, x, can be achieved by a composition of the inverse cumulative distribution function of a Gaussian random variable and the cumulative probability function of a nonCentral Chi random variable, i.e.,

x=PG1(Pχ(mη,σg,N)η,σg),
(12)

where the inverse cumulative distribution function of a Gaussian random variable is given by

PG1(yη,σg)=η+σg2erf1(2y1).
(13)

Note that erf −1 is the inverse of the error function. We should mention that, in practice, an outlier-rejection step is recommended in Eq. (12). Specifically, we shall identify x in Eq. (12) as an outlier if the following inequalities do not hold: (α/2) ≤ P[chi] (m | η,σg, N) <1−(α/2) where α may be between 0 and 1 inclusively but it is usually set to a user-specified value of 0.005, 0.001 or 0.0005.

The method of mapping an arbitrary distribution to a Gaussian distribution is well known, e.g., (35,36). In general, however, this type of mapping is of limited value without a priori knowledge of both η and σg, except for those that map from a Gaussian-derived distribution to a Gaussian distribution in which η and σg may be estimated through the proposed technique. Therefore, in the context of the present work, the parameters, η and σ, in Eq. 12) are replaced by their corresponding estimates, [eta w/ hat]g and [sigma with hat]g, which can be obtained through the techniques discussed in the previous section and in (30), respectively.

III. RESULTS

The validity of the proposed scheme is analyzed with several simulation tests.

A. NonCentral Chi random samples drawing from the same distribution

We will begin with the simplest case—that is, the mapping of noisy nonCentral Chi signals, which are drawn from the same distribution characterized by constant η and σg, to noisy Gaussian signals. Without loss of generality, we take N to be unity.

This type of data in which samples are drawn from the same distribution is rare in practice but is useful for illustrating the basic idea of the mapping between nonCentral Chi and Gaussian distributions. Note that this type of data is not an ordered sequence, and therefore, does not require a smoothing spline to estimate the “average value”—the sample mean of the data is sufficient in this case.

We should also note that the Gaussian noise SD cannot be estimated from this type of data using the noise variance estimation techniques discussed in (11,29,30,32,33) because there is no “background” in this type of data to estimate noise variance. Fortunately, other approaches can estimate both the underlying signal intensity and the Gaussian noise SD. Here, we note two approaches—our previously proposed analytically exact scheme (16), and the maximum likelihood approach as discussed in (37). One of the notable differences between these two approaches is that the former is a 1-D optimization procedure while the latter is a 2-D optimization procedure.

In this example, we will use the analytically exact scheme (16) to estimate both the underlying signal intensity and the Gaussian noise SD. Figure 2A shows the histogram of 20000 random samples that were drawn from a Rician distribution with η = 25 and σg = 50 (or η/σg = 0.5). The sample mean and the standard deviation of these random samples were 66.727 and 34.729, respectively. The magnitude SNR was 66.727/34.729 = 1.921, which corresponded to an estimated SNR, [eta w/ hat]/[sigma with hat]g, of 0.515. The estimated underlying signal intensity and the estimated Gaussian noise SD were found to be [eta w/ hat] = 25.74 and [sigma with hat]g = 49.98, respectively.

Fig. 2
(A) Histogram of 20000 random signals generated from a Rician distribution. (B) Histogram of the transformed signals.

Based on the estimated values of η and σg, the noisy Rician samples were then transformed to noisy Gaussian samples through the third stage of the proposed scheme. The histogram of the transformed signals is shown in Figure 2B. The sample mean and the standard deviation of these random transformed samples were 25.72 and 50.01, respectively.

B. Medium samples generated from a 1-D exponentially decaying model with a large number of repeated measurements

In this and the next examples, we investigate the statistical properties of the proposed scheme with data generated from a simple exponentially decaying model of the following form, s0ebD, taken from diffusion-weighted MRI. The b -value is an experimentally controlled variable that determines the level of diffusion weighting, which affects the level of attenuation of the non-diffusion-weighted signal, and D is the unknown diffusion coefficient.

Data generated from an exponentially decaying model are particularly useful for testing the proposed scheme because each measurement obtained at a different b-value is in fact drawn from a different distribution. Since there is only one measurement at each b-value, using the sample mean as the “average value” at each b-value would be too variable. Therefore, the “average value” at each b-value has to be estimated from a smoothing method such as the penalized spline where a collection of measurements at different b-values is treated as a whole to estimate the “average values” at all b-values.

Here, we generated 50000 sets of 30 measurements (Rician signals) from the following expression (s0exp(bD)+ε1)2+ε22 with s0 = 1000, D = 2.1×10−3 mm2/s, and ε’s are the Gaussian random variables with mean zero and standard deviation of 100.

The 30 measurements are sampled uniformly from b-value of 50 s/mm2 to 1993 s/mm2 with a gap of 67 s/mm2. Figure 3A shows the sample mean and the sample standard deviation of the 50000 measurements at each b-value. The error bar denotes one standard deviation away in both directions from the sample mean. The blue curves in Figures 3A and 3B are the expected value computed from the first moment of the Rician random variables.

Fig 3
(A) The expected value of the magnitude signal evaluated with a known Gaussian noise SD of 100 unit is shown as a blue curve, and the gray box and the error bar at each b-value represent the sample mean and the sample standard deviation that are obtained ...

Each set of 30 measurements is analyzed through the proposed scheme using the penalized spline with truncated polynomial basis of degree 4 and with 3 knots at {452, 988., 1457} s/mm2. The results of these 50000 sets for each stage of the proposed scheme are shown in Figures 3B, 3C and 3D. Figures 3B and 3C show the sample mean and the sample standard deviation of the spline estimates and of the fixed point estimates, respectively. The red curves in Figures 3C and 3D are the ground truth, i.e., s0 exp(−bD). Figures 3D and 3E show the sample mean and the sample standard deviation of the transformed signals obtained through the proposed framework and the method of Gudbjartsson and Patz (19), respectively.

In Figure 3D, it is clear that the sample mean at each b-value is close to the ground truth value but the variance (or SD) increases as the SNR decreases. The increase in SD is mainly due to a lack of sufficient samples because the ideal or expected behavior is that the variance should be constant (Figure 1B). As an example, we compare the result from the above simulation to that of another simulation in which the number of sampling points on the b-value axis was increased to 98, see Figure 4. It is clear from Figure 4 that the Gaussian noise SD estimates of the 98-point fit are collectively much closer to the ground truth value of 100 (arbitrary unit) that those of the 30-point fit.

Fig 4
The estimated Gaussian noise SD as a function of b-value with two Monte Carlo runs under the same simulation conditions but with two different sample sizes—30 and 98.

C. Large samples generated from a 1-D exponentially decaying model without repeated measurements

The same exponentially decaying model in diffusion-weighted MRI and the same set of parameters, D = 2.1×10−3 mm2/s and the Gaussian noise SD of 100, are used in this example. Here, we have only one set of 2476 measurements sampled from 50 s/mm2 to 5000 s/mm2 with a gap of 2 s/mm2. The penalized spline with a truncated polynomial basis of degree 4 and with 5 knots at {872, 1698, 2524, 3348, 4174} s/mm2 was used in this example.

The goal of this example is to show the qualitative features of the noisy Rician signals and of the transformed signals obtained through the proposed framework and the method of Gudbjartsson and Patz (19). We also compare and contrast the results from the parametric fits (mono-exponential and bi-exponential fits) to both the noisy signals and the transformed signals.

Figure 5A shows the noisy Rician signals. Figures 5B and 5C show the transformed signals obtained through the proposed framework and the method of Gudbjartsson and Patz (19), respectively. The results of both a mono-exponential fit and a bi-exponential fit to the noisy Rician signals are shown in Figure 5D. It is interesting to note that a bi-exponential model fits the noisy Rician signals rather well—the bi-exponential model is almost superimposed upon the expected curve. Figure 5E shows the result of a mono-exponential fit to the transformed signals obtained through the proposed framework; the resultant curve is close to the ground truth. The results of both a mono-exponential fit and a bi-exponential fit to the transformed signals obtained through the method of Gudbjartsson and Patz (19) are shown in Figure 5F. The estimates of the parameters, (s0, D), obtained through a mono-exponential fit of the noisy Rician signals, the transformed signals based on the proposed framework, and the corrected signals based on the method of Gudbjartsson and Patz (19) were found to be (597.4, 7.3 × 10−4 mm2/s), (966.4, 2.0 × 10−3 mm2/s), and (774.3, 1.3 × 10−3 mm2/s), respectively. In the bi-exponential fit of the noisy Rician signals and of the corrected signals based on the method of Gudbjartsson and Patz (19), we found (ŝ0 = 1036.5, D1 = −1.8×10−5 mm2/s, D2 = 2.7×10−3 mm2/s, 0.11) and (ŝ0 = 1040.9, D1 = −3.0×10−5 mm2/s, D2 = 2.7×10−3 mm2/s, 0.087), respectively. Note that the last item in each of the lists above is the (volume) fraction associated with D1.

Fig 5
(A) A collection of 2476 noisy magnitude signals sampled from 50 s/mm2 to 5000 s/mm2 with a gap of 2 s/mm2. These noisy Rician signals are generated with a known Gaussian noise SD of 100. (B) The transformed signals obtained through the proposed method. ...

D. Medium samples generated from a 3-D exponentially decaying model with a large number of repeated measurements

In this example, we will illustrate the proposed scheme with data sampled on a unit sphere. The spherical harmonic spline model will be used to transform the nonCentral Chi signals to Gaussian signals. A brief introduction to the spherical spline is provided in Appendix B.

For simplicity, the noisy Rician signals will be generated from a single tensor model according to the following expression, (s0exp(bgTDg)+ε1)2+ε22 where s0 = 1000, D is the diffusion tensor, g is a unit gradient vector, T denotes matrix or vector transposition, and ε’s are the Gaussian random variables with mean zero and SD of 100. Further, the synthetic tensor is given by:

D=(9.51.11.61.16.70.51.60.54.8)×104mm2/s.

For visualization purposes, we first parametrize the unit gradient vector in terms of spherical coordinates, i.e., g = [sin(θ)cos([var phi]), sin(θ)sin([var phi]), cos(θ)]T. With this parametrization, we can plot the underlying signal intensity and the expected value of the Rician random variables as functions of the spherical coordinates. Figure 6A shows the underlying signal intensity as a function of the spherical coordinates at a b-value of 3000 s/mm2, and Figure 6B is the corresponding expected value in magnitude obtained from the known Gaussian noise SD of 100.

Fig 6
(A) The underlying signal intensity from a single tensor model as a function of spherical coordinates evaluated with a constant b-value of 3000 s/mm2. (B) The expected value of the Rician signals (with the known Gaussian noise SD of 100) as a function ...

Similar to the one-dimensional case, we chose 30 unit gradient vectors that are uniformly distributed on the sphere (based on the electrostatic repulsion scheme (38)), and the spherical coordinates of each of the gradient vectors are color-coded in Figure 6C. Figure 6D shows the color-coded underlying signal intensity in ascending order and their respective expected values (the first moment of the Rician random variables) with a Gaussian noise SD of 100. There are 50000 sets of 30 measurements and each measurement in the set is a sample on the unit sphere obtained through one of the gradient vectors. The sample mean and the sample SD of the noisy Rician signals of all the spherical coordinates are shown in Figure 6E. Finally, each set of measurements is analyzed through both the proposed scheme using the spherical spline with spherical harmonics of even degree up to l = 6 and the method of Gudbjartsson and Patz. The results are shown in Figures 6F and 6G, respectively.

It is clear from the results shown in Figure 6F that the sample means are close to the ground truth values but the variance increases slightly as the SNR decreases. The increase in variance is to be expected since only 30 gradient directions are used. More importantly, we can expect the variance to get closer to a constant value that is independent of the SNR level as the size of the samples becomes larger.

E. Illustration using experimental data

We illustrate the performance of our approach on an excised rat hippocampus data set. The data set contains a series of diffusion-weighted images obtained by varying the diffusion gradient strength. The rat was perfusion-fixed with 4% paraformaldehyde in phosphate buffered-saline (PBS), the hippocampus was dissected and kept in fixative for more than 8 days. Prior to imaging, the sample was washed overnight in PBS. The imaging was performed using a 14.1T narrow-bore spectrometer where a pulsed gradient stimulated echo pulse sequence was employed. The imaging parameters were: TE=12.6ms, TR=1000ms, resolution=(78×78×500)μm3, matrix size=(64×64×3), number of repetitions=4, diffusion gradient pulse duration (δ)=2ms, and diffusion gradient separation (Δ)=24.54ms. The data set contains a total of 33 images with different diffusion gradient strengths increasing from 0 to 2935mT/m in steps of 91.75mT/m. One diffusion weighted image is shown in Figure 7A.

Fig 7
Experimental data. (A) A diffusion-weighted image of a hippocampus with a red square indicating the four different pixel locations where the noisy magnitude signals of each pixel (with different b-values) are analyzed using the proposed method. The results ...

Four neighboring pixels indicated with a red square were selected for further analyses. The noisy magnitude signals and the noisy transformed signals of each of the pixels as a function of b-value are shown in Figures 7B–7E as blue and red dots, respectively. The blue curve in each of the panels is obtained through a least squares fit of a bi-exponential function to the noisy magnitude signals. The red curve in each of the panels is obtained through a least square fit of a bi-exponential function to the noisy transformed signals produced by the proposed framework. Note that the penalized spline with a truncated polynomial basis of degree 4 and with 4 knots was used in this example. The estimated Gaussian noise standard deviation was 0.88. Further, the estimated parameters obtained from a least squares fit of a bi-exponential function to both the noisy magnitude signals and noisy transformed signals are shown below:

Bi-exponential fit to the noisy magnitude signalsŝ0 (a.u.)D1 (×10−5 mm2/s)D2 (×10−4 mm2/s)Volume fraction associated with D1
Fig. 7B62.480.825.30.027
Fig. 7C63.102.06.20.037
Fig. 7D64.280.816.00.026
Fig. 7E64.361.45.50.027
Bi-exponential fit to the noisy transformed signals
Fig. 7B62.69.05.50.060
Fig. 7C63.310.96.60.077
Fig. 7D64.411.36.20.056
Fig. 7E64.49.95.70.048

If both the estimated Gaussian noise SD and each of the red curves are assumed to be the ground truth values then the expected value (or the first moment) of a Rician distribution as a function of b-value can be computed and is shown in dark gray; these expected values are in good agreement with the blue curve, which is an indication that the red curve is a good approximation of the underlying signal intensities.

IV. DISCUSSION

In this work, our main objective is to demonstrate that nonCentral Chi signals can be transformed into Gaussian signals and present as clearly as possible the basic ideas as well as the nuts and bolts of the proposed scheme.

This paper can be thought of as a sequel to but independent of our recent paper on the probabilistic and self-consistent approach to the identification and estimation of noise (PIESNO) (30) because the noise estimate on which the proposed framework depends can estimated through other techniques. The fixed point formula of the underlying signal intensity and the technique proposed in (30) represent our major attempt to decouple the fixed point formula of SNR (16) into two self-consistent approaches for estimating the underlying signal and the Gaussian noise SD.

The advantage of this decoupling is substantial because the estimation of the Gaussian noise SD can be obtained from a much larger collection of samples (30). As a consequence, the precision of the Gaussian noise SD estimate will be significantly increased, and in turn, the precision of the underlying signal intensity estimate will also be increased. As discussed above, the decoupling is more useful and practical than the fixed point formula of SNR because we do not usually have many data that are drawn from the same distribution, Fig 1B. It is interesting to note that the way in which the present scheme is realized is due in part to this practical constraint.

The combination of these stages presented here is, to the best of our knowledge, unique and novel. Moreover, the formulation of the second stage is conceptually very different from our previous approach (16), see Figure 1B.

The first and third stages of the proposed scheme are well known but these stages, alone or together, are not sufficient for mapping nonCentral Chi signals to Gaussian signals without the second stage. The three stages used in the proposed scheme in a sense form an irreducible set of steps that is necessary to map noisy nonCentral Chi signals to noisy Gaussian signals. While different fitting or smoothing methods may be used in the first stage, the last two stages are strictly mathematical and fixed even though there exists several different iteration scheme for finding the fixed point of the underlying signal intensity, see Appendix D.

In the second stage, we should point out that the suggested modification to the fixed point formula of the underlying signal intensity for the special case in which the “average value” of the magnitude signals is below the noise floor can be further improved. Although we have provided a theoretical justification for this modification, we believe further studies are needed to investigate other approaches to find the fixed point estimate for this particular situation.

The examples illustrated above clearly show the feasibility and effectiveness of the proposed scheme in mapping noisy magnitude signals to noisy signal intensities. The proposed scheme can be extended to transforming any Gaussian-derived noisy signals, e.g., Rayleigh, Rician, nonCentral Chi, and nonCentral Chi-squared distributed signals, to noisy Gaussian signals by finding the specific fixed point formula used in the second stage.

The basic idea of our approach is general and can be easily adapted to many MRI and non-MRI applications, e.g., the Laser Interferometric Gravitational Wave Observatory (LIGO) (39,40) and communication systems (31), by selecting an appropriate data smoothing method that is optimal for the application-specific sampling space. For example, the penalized spline model or the wavelet smoothing spline may be useful in the analysis of functional MRI data while spherical splines are particularly useful to diffusion tensor imaging and high angular diffusion imaging techniques (39,21). We should also point out that some algorithms of least squares estimation may need to be modified in order to handle negative values in the transformed data. For example, the nonnegative least squares approach, e.g., (41), may be needed to analyze the transformed signals.

Spline models are known for their flexibility in capturing unknown trends in the data but they come at the cost of slightly higher susceptibility to noise such as spurious oscillatory trends in the spline estimates. Therefore, optimal performance cannot be expected of any spline or regression models when the number of samples is very small, and simulation studies may be needed to get an initial assessment of the number of samples needed for a particular experimental design. In this work, the GCV function was used as a smoothing criterion because it has several desirable properties, the most notable of which is that as the number of samples increases, the spline estimate obtained via the GCV becomes closer to the estimate that is obtained by minimizing the mean square error between the estimate and the unknown ground truth (27). Finally, the spurious trends in the spline estimates mentioned above can be partially removed if the transformed signals are fitted with some parametric functions based on a priori physical or mathematical model that is less flexible than the smoothing spline, e.g. mono-, bi- or tri-exponential functions for the one-dimensional diffusion data or the diffusion tensor model for the three dimensional diffusion data.

In quantitative MRI, anatomically or physiologically relevant parameters are usually estimated from a least squares model. As noted in the introduction, the Gaussian-distributed noisy signals are of interest here rather than the Gaussian-derived random signals because the Gaussian-distributed noisy signals are generally more amenable to statistical treatment based on the principle of least squares, e.g., (4244). It is important to point out that one of the basic assumptions in a least squares model is that random errors follow a Gaussian distribution. The principle of least squares is very powerful because of its mathematical tractability not only in parameter estimation but also in hypothesis testing and confidence interval estimation. Further, the least squares and maximum likelihood estimators are equivalent under the assumption of normality of random errors (45).

In this work, we have presented a novel approach for transforming noisy nonCentral Chi signals to noisy Gaussian signals, thus making least squares approaches uniformly applicable for analyzing MRI data. The present approach is a major advance in facilitating and improving all subsequent data analysis and processing steps in a quantitative MRI pipeline.

Table 1
The algorithm for finding the fixed point estimate of the underlying signal intensity.

Acknowledgments

We are grateful to Liz Salak for reviewing the paper. C. G. Koay was supported in part by the Eunice Kennedy Shriver National Institute of Child Health and Human Development, the National Institute on Drug Abuse, the National Institute of Mental Health, and the National Institute of Neurological Disorders and Stroke as part of the NIH MRI Study of Normal Pediatric Brain Development with supplemental funding from the NIH Neuroscience Blueprint. We would like to thank Drs. Timothy M. Shepherd and Stephen J. Blackband for the MRI data set.

APPENDIX A

PENALIZED SPLINE

A penalized spline function with a truncated polynomial basis (25) of degree p and K knots at {κ1, …, κK} is given by:

f(x)=β0+i=1pβixi+j=1Kβpj(xκj)+p,
(A1)

where the operation, (x)+, returns x if x > 0, and zero otherwise, (xκj)+p are the spline basis functions and κj are the knots.

If there are n observations, {y1, …, yn}, sampled at {x1, …, xn}, then Eq. (A1) can be expressed in matrix notation as follows:

y=Xβ,
(A2)

where the design matrix, X, is given by:

(1x1x1p(x1κ1)+p(x1κK)+p1xnxnp(xnκ1)+p(xnκK)+p).
(A3)

In practice, we usually normalize the coordinates,{x1, …, xn}, by the maximum of the absolute value of the elements in {x1, …, xn} to avoid numerical instability associated with the QR decomposition of X. Therefore, the construction of X is based on the normalized coordinates, and so are the knots. Therefore, the knots reported in the Results section have to be scaled accordingly.

In the ordinary least squares estimation, the goal is to find β that minimizes ||y||2 while, in the penalized spline estimation, the goal is to find β that minimizes

yXβ2+λβTDβ,
(A4)

where T denotes matrix or vector transposition, D is a diagonal matrix whose first p +1 diagonal elements are zero and the rest of the diagonal elements are unity and λ is the penalty parameter (or the smoothing parameter). The smoothed observation vector, ŷλ, estimated from the penalized spline can be expressed as follows:

y^λ=Sλy,
(A5)

where

Sλ=X(XTX+λD)1XT
(A6)

is known as the smoother matrix.

The procedure presented thus far does not provide a means to find an optimal λ. Here, we use the GCV function (27) to select an optimal λ, which will be denoted by λGCV; note that λGCV is a minimizer of the GCV function and the GCV function is given by:

GCV(λ)=RSS(λ)/(1tr(Sλ)/n)2,
(A7)

where RSS(λ) = ||yŷλ||2 is the residual sum of squares, tr denotes the matrix trace operation, and n is the number of observations. For a numerically stable implementation of the penalized spline estimation, see Appendix C.

APPENDIX B

SPHERICAL SPLINE

According to the expansion theorem of the spherical harmonics (46), any continuous function, f (θ,[var phi]), on the unit sphere together with continuous derivatives up to second order can be expanded in terms of the Laplace series of the spherical harmonics:

f(θ,φ)=l=0m=llβlmYlm(θ,φ)
(B1)

where Ylm(θ,φ) is the spherical harmonic of lth degree and of mth order. The spherical harmonic can be expressed as a real rather than complex function, and this is given by (9,46):

Ylm(θ,φ)={2l+12π(l+m)!(lm)!sin(mφ)Plm(cos(θ)):lm12l+14πPlm(cos(θ)):m=02l+12π(lm)!(l+m)!cos(mφ)Plm(cos(θ)):1ml.

Note that Plm is the associated Legendre polynomial of mth order, and the arguments of the spherical harmonic function are defined within these intervals: 0 ≤ θ < π and 0 ≤ [var phi] < 2π.

The smoothing spherical spline (9,26) is built on the Laplace series with finite number of terms as well as on the following linear matrix structure:

y=Xβ,
(B2)

where y is an array of measurements sampled at {(θ1,[var phi]1), …,(θn, [var phi]n)}, and the design matrix, X, is made up of the spherical harmonics up to some maximum degree, lmax, and in a specific order that corresponds exactly to how the coefficients of the Laplace series are ordered in β, e.g.,

β=[β00,β11,β10,β11,β22,,β22,,βlmaxlmax]T.
(B3)

The goal in the smoothing spherical spline estimation (9,26) is to find β that minimizes

yXβ2+λβTDβ,
(B4)

where D is a diagonal matrix with each diagonal element takes on the value of l2 (l +1)2 where l is the degree associated with the corresponding element, βlm, in β.

The solution of the above estimation has the same matrix structure as that of the penalized spline estimation in Appendix A. Note that in diffusion MRI, only spherical harmonics of even degree are of interest because of the assumption that the diffusion process has antipodal symmetry. Therefore, Eq. (B3) has to be modified accordingly.

APPENDIX C

A COMMON COMPUTATIONAL METHOD FOR THE PENALIZED AND SPHERICAL SPLINES

The key computational problem in penalized spline estimation is to find an efficient matrix decomposition of the smoother matrix:

Sλ=X(XTX+λD)1XT.
(C1)

Our approach in computing the smoother matrix is slightly different from that of (25) in that we use the QR decomposition to factor X rather than the Cholesky decomposition to factor XT X.

Let the QR decomposition of X be Q R where Q is an orthogonal matrix, i.e., QT Q = I, and R is an upper triangular matrix. Note that I is the identity matrix. Substituting Q R into Eq. (C1), we have:

Sλ=QR(RTR+λD)1RTQT=Q(I+λRTDR1)1QT.
(C2)

At this stage, the singular value decomposition (SVD) of Σ [equivalent] RTDR−1 is needed, which will be denoted by Σ [equivalent] UΔVT. Note that Δ is a diagonal matrix and its diagonal elements are the singular values of Σ. Further note that U = V because Σ is a symmetric matrix. Finally, the smoother matrix is given by:

Sλ=QU(I+λΔ)1(QU)T,=MWMT.
(C3)

where M = Q U is an orthogonal matrix and W is a diagonal matrix and its diagonal elements are defined by Wii=11+λΔii. Since M is an orthogonal matrix, tr(Sλ) is simply tr(W)=i11+λΔii. In practice, the factor M may be precomputed and only the diagonal matrix W needs to be updated during the optimization search for λGCV.

APPENDIX D

AN IMPLEMENTATION OF THE FIXED POINT FORMULA OF THE UNDERLYING SIGNAL INTENSITY

In this appendix, we provide an implementation of the fixed point formula of the underlying signal intensity, which is based on Newton’s method of root finding. It begins with an iteration scheme of the following form:

ηk+1K(ηkm^,σ^g,N)=ηkf(ηkm^,σ^g,N)f(ηkm^,σ^g,N),
(D1)

where f (η| m, [sigma with hat]g, N) = g(η| m, [sigma with hat]g, N) − η, f′ denotes the first order derivative of f with respect to η, and g is defined in Eq. (11). The function, K(η| m, σ, N), can be further simplified to the following expression:

ηg(ηm,σ,N)(g(ηm,σ,N)η)η(1βN22NF11(12,N,η22σ2)F11(12,N+1,η22σ2))g(ηm,σ,N).
(D2)

The basic algorithm of the above iteration is given in Table 1. It is clear from Eq. (D2) that the expression is different from that of (16).

We should note that there are other iteration schemes for finding the fixed point of the underlying signal. Here, we provide another iteration scheme also based on the Newton’s method:

η+2Nσ(mβNσF11(12,N,η22σ2))βNηF11(12,N+1,η22σ2).
(D3)

The expression above is derived directly from Eq. (7). Note that Table 1 can be easily adapted for Eq. (D3) instead of Eq. (D2).

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

1. Lauterbur PC. Image formation by induced local interactions: Examples employing nuclear magnetic resonance. Nature. 1973;242(5394):190–191. [PubMed]
2. Ogawa S, Lee TM, Kay AR, Tank DW. Brain Magnetic Resonance Imaging with Contrast Dependent on Blood Oxygenation. Proceedings of the National Academy of Sciences. 1990;87(24):9868–9872. [PubMed]
3. Basser PJ, Mattiello J, Le Bihan D. MR diffusion tensor spectroscopy and imaging. Biophys J. 1994;66(1):259–267. [PubMed]
4. Tuch DS, Reese TG, Wiegell MR, Makris N, Belliveau JW, Wedeen VJ. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magnetic Resonance in Medicine. 2002;48(4):577–582. [PubMed]
5. Frank LR. Characterization of anisotropy in high angular resolution diffusion-weighted MRI. Magnetic Resonance in Medicine. 2002;47(6):1083–1099. [PubMed]
6. Anderson AW. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magnetic Resonance in Medicine. 2005;54(5):1194–1206. [PubMed]
7. Hess CP, Mukherjee P, Han ET, Xu D, Vigneron DB. Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis. Magnetic Resonance in Medicine. 2006;56(1):104–117. [PubMed]
8. Özarslan E, Shepherd TM, Vemuri BC, Blackband SJ, Mareci TH. Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT) NeuroImage. 2006;31(3):1086–1103. [PubMed]
9. Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Regularized, fast, and robust analytical Q-ball imaging. Magnetic Resonance in Medicine. 2007;58(3):497–510. [PubMed]
10. Wu YC, Alexander AL. Hybrid diffusion imaging. Neuroimage. 2007;36(3):617–629. [PMC free article] [PubMed]
11. Henkelman RM. Measurement of signal intensities in the presence of noise in MR images. Med Phys. 1985;12(2):232–233. [PubMed]
12. Liu J, Koenig JL. An automatic phase correction method in nuclear magnetic resonance imaging. Journal of Magnetic Resonance. 1990;86:593–604.
13. Chen L, Weng Z, Goh L, Garland M. An efficient algorithm for automatic phase correction of NMR spectra based on entropy minimization. Journal of Magnetic Resonance. 2002;158:164–168.
14. Bretthorst GL. Automatic phasing of MR images. Part I: Linearly varying phase. Journal of Magnetic Resonance. 2008;191:184–192. [PubMed]
15. Bretthorst GL. Automatic phasing of MR images. Part II: Voxel-wise phase estimation. Journal of Magnetic Resonance. 2008;191:193–201. [PubMed]
16. Koay CG, Basser PJ. Analytically exact correction scheme for signal extraction from noisy magnitude MR signals. Journal of Magnetic Resonance. 2006;179(2):317–322. [PubMed]
17. Constantinides CD, Atalar E, McVeigh ER. Signal-to-noise measurements in magnitude images from NMR phased arrays. Magnetic Resonance in Medicine. 1997;38(5):852–857. [PMC free article] [PubMed]
18. Bernstein MA, Thomasson DM, Perman WH. Improved detectability in low signal-to-noise ratio magnetic resonance images by means of a phase-corrected real reconstruction. Med Phys. 1989;16(5):813–817. [PubMed]
19. Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magnetic Resonance in Medicine. 1995;34(6):910–914. [PMC free article] [PubMed]
20. Rice SO. Mathematical analysis of random noise. Bell System Technical Journal. 1944;23 and 24
21. Jones DK, Basser PJ. Squashing peanuts and smashing pumpkins: How noise distorts diffusion-weighted MR data. Magnetic Resonance in Medicine. 2004;52(5):979–993. [PubMed]
22. Vincent JL, Patel GH, Fox MD, Snyder AZ, Baker JT, Van Essen DC, Zempel JM, Snyder LH, Corbetta M, Raichle ME. Intrinsic functional architecture in the anaesthetized monkey brain. Nature. 2007;447(7140):83–86. [PubMed]
23. McGibney G, Smith MR. An unbiased signal-to-noise ratio measure for magnetic resonance images. Med Phys. 1993;20(4):1077–1078. [PubMed]
24. Miller AJ, Joseph PM. The use of power images to perform quantitative analysis on low SNR MR images. Magnetic Resonance Imaging. 1993;11(7):1051–1056. [PubMed]
25. Ruppert D, Wand MP, Carroll RJ. Semiparametric regression. Cambridge University Press; 2003.
26. Wahba G. Spline models for observational data. SIAM; 1990.
27. Craven P, Wahba G. Smoothing noisy data with spline functions. Numerische Mathematik. 1978;31(4):377–403.
28. Bertero M, Boccacci P. Introduction to inverse problems in imaging. Philadelphia: Institute of Physics Publishing; 1998.
29. Sijbers J, Poot D, den Dekker AJ, Pintjens W. Automatic estimation of the noise variance from the histogram of a magnetic resonance image. Physics in Medicine and Biology. 2007;52(5):1335–1348. [PubMed]
30. Koay CG, Özarslan E, Pierpaoli C. Probabilistic Identification and Estimation Noise (PIESNO): A self-consistent approach via the median method and its applications in MRI. Journal of Magnetic Resonance Submitted [PMC free article] [PubMed]
31. Proakis J. Digital communications. New York: McGraw-Hill; 2001.
32. Edelstein WA, Bottomley PA, Pfeifer LM. A signal-to-noise calibration procedure for NMR imaging systems. Med Phys. 1984;11(2):180–185. [PubMed]
33. Chang L-C, Rohde GK, Pierpaoli C. An automatic method for estimating noise-induced signal variance in magnitude-reconstructed magnetic resonance images. SPIE Medical Imaging: Image processing. 2005;5747:1136–1142.
34. Courant R, John F. Introduction to Calculus and Analyis I. New York: Springer; 1989.
35. Liu P, Der Kiureghian A. Multivariate distribution models with prescribed marginals and covariances. Prob Eng Mech. 1986;1:105–112.
36. van Albada S, Robinson P. Transformation of arbitrary distributions to the normal distribution with application to EEG test-retest reliability. J Neu Meth. 2007;161:205–211. [PubMed]
37. Sijbers J, Den Dekker AJ. Maximum likelihood estimation of signal amplitude and noise variance from MR data. Magn Reson Med. 2004;51:586–594. [PubMed]
38. Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magnetic Resonance in Medicine. 1999;42(3):515–525. [PubMed]
39. Abramovici A, Althouse WE, Drever RWP, Gursel Y, Kawamura S, Raab FJ, Shoemaker D, Sievers L, Spero RE, Thorne KS, Vogt RE, Weiss R, Whitcomb SE, Zucker ME. LIGO: The Laser Interferometer Gravitational-Wave Observatory. Science. 1992;256(5055):325–333. [PubMed]
40. Cutler C, Flanagan EE. Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral waveform? Phys Rev D. 1994;49(6):2658–2697. [PubMed]
41. Lawson CL, Hanson RJ. Solving least squares problems. Philadelphia: SIAM; 1974.
42. Koay CG, Chang L-C, Carew JD, Pierpaoli C, Basser PJ. A unifying theoretical and algorithmic framework for least squares methods of estimation in diffusion tensor imaging. Journal of Magnetic Resonance. 2006;182(1):115–125. [PubMed]
43. Koay CG, Chang LC, Pierpaoli C, Basser PJ. Error propagation framework for diffusion tensor imaging via diffusion tensor representations. IEEE Transactions on Medical Imaging. 2007;26(8):1017–1034. [PubMed]
44. Koay CG, Nevo U, Chang L-C, Pierpaoli C, Basser PJ. The elliptical cone of uncertainty and its normalized measures in diffusion tensor imaging. IEEE Transactions on Medical Imaging Accepted [PMC free article] [PubMed]
45. Seber GAF, Lee AJ. Linear Regression Analysis. New York: Wiley; 2003.
46. Courant R, Hilbert D. Methods of Mathematical Physics. New York: Wiley; 1989.