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.

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 [3], [4], [5], [6], [7], [8], [9], [10] 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 [12], [13], [14], [15], 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 [3], [4], [5], [6], [7], [8], [9]. 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 Fig. 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.

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 Fig. 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 Fig. 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 excellent methods are available 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.