|Home | About | Journals | Submit | Contact Us | Français|
It is often desirable to separate voxels that contain signal from tissue along with measurement noise from those that contain purely measurement noise. Generally this separation called thresholding utilizes only the magnitude portion of the images. Recently methods have been developed that utilize both the magnitude and phase for thresholding voxels. This manuscript is an extension previous work and uses the bivariate normality of the real and imaginary values with phase coupled means. A likelihood ratio statistic is derived that simplifies to a more familiar form that is F-distributed in large samples. It is shown that in small samples, critical values from Monte Carlo simulation can be used to threshold this statistic with the proper Type I and Type II error rates. This method is applied to susceptibility weighted magnetic resonance images and shown to produce increased tissue contrast.
In magnetic resonance images, it is well known that noise manifests as independent and identically distributed normally distributed noise in the real and imaginary parts of the k-space measurements with a mean of zero and a constant variance [1,2,3]. It is also known that there is a linear relationship between complex-valued k-space measurement and complex-valued voxel measurements . From the normality of k-space measurements and their linear relationship with voxel measurements, the voxel measurements are also normally distributed [1,2,4]. Voxel measurements can be described as
where yR and yI are the measurements for the real and imaginary parts, εR and εI are the error terms for the real and imaginary parts, while ρ and θ are the population magnitude and phase. It is desirable to separate voxels that are pure noise from those that contain signal and noise . Assuming that the additive noise terms in Eq. (1) are normally distributed with a mean of zero and variance σ2, the joint probability distribution of the bivariate voxel observation (yR,yI) is
Upon conversion from Cartesian coordinates to polar coordinates in Eq. (2), the joint distribution of the observed magnitude and phase (m,) is
We would like to determine if the observed magnitude and phase in a voxel are signal or if they are noise. Given a collection of measurements (m1,1),…,(mn,n) from p(m,) in Eq. (3), the likelihood is
It is desirable to increase tissue contrast by thresholding noise voxels where the magnitude and phase are not statistically different from zero from those where the magnitude and phase are statistically different from zero. This separation of voxels that contain signal (within the object being imaged) has been described by other means that use an estimate of the variance from voxels that contain pure noise . We will show in the next section that a formal statistic can be derived from Eq. (4) and a statistical hypothesis test performed on the population magnitude and phase parameters.
In statistics the formal procedure to separate voxels that contain noise from those that also contain signal is to perform a hypothesis test on the magnitude and phase. In hypothesis testing there are four possible outcomes. One can imagine the usual 2×2 table where the first row is labeled “Reject H0,” the second row is labeled “Do Not Reject H0,” the first column is labeled “H0 True,” and the second column is labeled “H0 False.” In the top right and bottom left cells correct decisions are made. In the top left cell a Type I error is made in which the null hypothesis is rejected when it is true. The probability of a Type I error is the level of significance denoted by α. In the bottom right cell a Type II error is made in which the null hypothesis is not rejected when it is false. The probability of a Type II error is the level of significance denoted by β. These two error rates will be examined in more detail later.
Formally, this voxel separation procedure can be achieved by testing null and alternative hypotheses H0: ρ=0, θ=0 versus H1: ρ>0, θ≠0 where ρ is the population magnitude and θ is the population phase. In general, hypotheses can be evaluated with test statistics that are derived via a likelihood ratio statistic [6,7]. Under the constrained null hypothesis H0: ρ=0, θ=0, the maximum likelihood estimators (MLEs) for the magnitude and phase from Eq. (4) are
while under the unconstrained alternative hypothesis H1: ρ>0,θ≠0 the MLEs from Eq. (4) are
where is the mean of the real channel measurements and is the mean of the imaginary channel measurements. These estimates in Eq. (5) and Eq. (6) are then inserted back into the likelihoods and the ratio taken. Under some regularity conditions (that are satisfied in this case), − 2log(λ) is asymptotically χ2 distributed with degrees of freedom equal to the difference in the number of estimated parameters between H0 and H1, (two in this case).
Upon insertion of the MLEs in Eq. (5) and Eq. (6) back into the likelihood function in Eq. (4), the ratio of null hypothesis likelihood over alternative hypothesis taken. Some terms in the likelihood ratio can be cancelled to become
However, algebra can usually be performed to simplify this into a statistic with known distribution under the null hypothesis. This procedure can be applied to the general linear model to derive the usual t and F statistics. Applying this procedure here, the test statistic
can be arrived at as F = (1−λ1/n). The probability distribution of this statistic in Eq. (8) needs to be found. The steps through the logic for the derivation of the distribution for the numerator and the denominator of the F statistic are shown below in Table 1.
Since the numerator and denominator terms x1 and x are χ2 distributed with 2 and 2n degrees of freedom, dividing them by their degrees of freedom and taking the ratio should result in a statistic that under the null hypothesis has an F distribution with 2 and 2n degrees of freedom. However, this is not true in this case. These two χ2 statistics must be statistically independent for this to be true. The correlation between these two statistics can be derived. First, since these are χ2 distributed, their expectation is their degrees of freedom and their variances are twice their degrees of freedom E (x1) = 2, E (x2) = 2n, var (x1) = 4, var (x2) = 4n. The covariance can be found as cov (x1, x2) = E (x1 · x2) − E (x1) E (x2), where
The correlation between x1 and x2 is now , which is . This correlation tends to zero in large samples (large n) and the F statistic becomes F distributed with two numerator and 2n denominator degrees of freedom. In non-fMRI applications where there may be a very small number of repeated images if any, this asymptotic result does not hold. Thus, critical values from the F distribution do not apply to this F statistic. However, critical values for small n can be achieved by way of Monte Carlo simulation. For a given level of significance (Type I error rate α), we reject H0 (do not threshold voxel) if the test statistic F is larger than the critical value Fα(2,2n) and do not reject (threshold voxel) if F is smaller than the critical value Fα(2,2n).
To examine the convergence of the distribution of the F statistic to the F distribution, one million simulated data sets was created under the null hypothesis (ρ=0 and θ=0) for n=5, 9, 25, 50, 100, and 250. Normally distributed independent noise variates were generated for the real and imaginary parts with a variance of σ2=1. This corresponds to a signal-to-noise ratio (SNR) of zero. It was found (but not shown), that regardless of the sample size n, the two statistics x1 and x2 in Table 1 that make up the numerator and denominator are χ2 distributed. Additionally, the correlation between the numerator and denominator χ2 statistics was verified. For each sample size n, the F statistic was computed for each of the data sets. Cumulative distribution functions (CDFs) from the Monte Carlo empirical distributions (dashed curves) and corresponding CDF for the F distribution (solid curves) are presented in Fig. 1 for n=5 (red), 9 (black), 25 (yellow), 50 (green), 100 (blue), 250 (violet). Note the disparity between the Monte Carlo empirical distributions (dashed curves) and corresponding CDF for the F distribution (solid curves). One can see that it takes a relatively large sample size for this disparity to decrease. Therefore for small samples, Monte Carlo critical values need to be used.
For our application, we will be utilizing the test statistic denoted by F for a sample of size n=9. The Type I error rate was examined in Fig. 1 for several values of n. However, the Type II rate is also of interest. To examine the Type II error rate, one million simulated data sets for n=9 were created under the alternative hypothesis (ρ≠0 and θ≠0) with ρ=(0,1,2,3,5) and θ=0°. Normally distributed independent noise variates were generated for the real and imaginary parts with a mean of zero and a variance of σ2=1. Pairs of histograms (not shown) for the one million data sets when ρ=0 and each of the alternative hypothesis cases ρ=(1,2,3,5) were formed. A vertical line denoted by Fα is drawn with the false positive rate α being the fraction of variates in the ρ=0 histogram that are larger than the vertical line. The true positive rate 1−α is the fraction of variates in the ρ=0 histogram that are smaller than the vertical line at Fα. The fraction of variates in the alternative hypothesis case (ρ=1,2,3,5) larger than the vertical line at Fα is the true negative rate 1−β. The false negative rate β is the fraction of variates in the alternative hypothesis case (ρ=1,2,3,5) smaller than the vertical line at Fα. By sliding this vertical line left and right for varying false positive rate areas α to the right of it, we can determine the true positive rate 1−β. A plot of α on the horizontal axis and 1−β on the vertical axis is called a receiver operating characteristic (ROC) curve . This curve can be made for each combination of ρ and θ. In Fig. 2 are the ROC curves for ρ=(0,1,2,3,5) and θ=0°. When the θ values are not zero, minor variations from the curves for the corresponding ρ value are produced.
To determine more accurate critical values in the upper tail of the F statistic for n=9, 5×107 data sets were generated. For each set, the F statistic was computed. A histogram of these F statistics for n=9 is presented in Fig. 3. These F statistics were ordered and percentiles determined. For example, the 0.95×(5×107)th largest value was taken as the 95th percentile. Selected critical values are presented in Table 2 for n=5 and n=9. These critical values will be used for thresholding the magnitude and phase of voxels. Additional critical values can reliably be interpolated for significance levels between those that have been presented.
From a statistical point of view, we would like to have n repeated images. However, we are interested in thresholding high-resolution anatomical images where replicated images are rarely available. We will use the observed magnitude and phase values in each voxel and its surrounding voxels. Each voxel and its eight neighbors (n=9) will be used to estimate the magnitude and phase of every voxel with image wrap around. Alternatively, each voxel and its four neighbors (n=5) can be used. The estimated values for each voxel will be used to compute the previously described F statistic. The map of these F statistics is thresholded with the critical values in Table 2. A zero-one mask is produced from the thresholded F statistics then applied to the original magnitude and phase images.
Since only a single image is available, using a voxel and its neighboring voxel will result in very minor local correlation in the F statistics. A statistical critical value for a given α level might be slightly lower than the values we have used from the Monte Carlo simulation. To investigate this phenomenon, additional simulations were performed where images were created and the F statistics computed using a voxel and its four neighbors. Images of size 352×512 and 1000×1000 were created where each voxel had a magnitude of ρ=(0,1,2,3,5) and a phase of θ=0°. Normally distributed independent noise variates were generated and added to the real and imaginary parts with a mean of zero and a variance of σ2=1. ROC curves were produced for each magnitude (equivalently SNR) value. The fraction of true positives and false positives were determined. The false positive α and false negative rates β were consistent with previous values and yielded visually identical ROC curves as seen in Fig. 2. Simulations were also performed in which 512×512 images were produced with a circle of radius 128 pixels containing the same magnitudes of ρ=(0,1,2,3,5) and a phase of θ=0°. True and false positive rates were also consistent as were ROC curves. The very small local correlation yields the same critical values in Table 2 and does not affect the global image threshold as previously reported . If one were extremely concerned about a potential change in the image threshold due to local correlations (particularly in small images), one could perform larger simulations for more accuracy or use a thresholding method that incorporates a neighborhood about the voxel of interest .
The magnitude and phase statistical thresholding method described in this paper was applied to the human susceptibility weighted imaging (SWI) MRI data . A SWI brain volume was acquired on a 3T Siemens Trio with a matrix size of 352×512, FOV of 176 mm × 256 mm, and an in-plane resolution of 0.5×0.5 mm2, TR=26 ms, TE= 15 ms, flip angle (FA) = 11° . The results from this analysis are presented in Figs. 4, ,5,5, and and6.6. The scale for the magnitude images are in arbitrary units and the scale for the phase images are in radians. In Fig. 4A and Fig. 4B the original observed SWI magnitude and phase image data are respectively presented. Note that in Fig. 4B there is high noise in the phase image outside the head and in some internal areas while the magnitude noise in Fig. 4A is visually relatively low. The magnitude and phase model was applied to each voxel using its eight neighbors for samples of n=9. In Fig. 4C and Fig. 4D the estimated SWI magnitude and phase image data are respectively presented. Note that there is a slight reduction in the high frequency image content. In Fig. 4E the estimated voxel variance is presented. Note that some high frequency spatial content is present in Fig. 4E. In Fig. 4F the computed F statistic map that is used for thresholding is presented. Note that spatial anatomical structure is present in Fig. 4F where voxels with larger magnitude and or phase signal have larger F statistics.
In Fig. 5A and Fig. 5B histograms of the original observed SWI magnitude and phase image data are respectively presented. Note that in Fig. 5A, the observed magnitude contains two populations of voxels that visually are fairly easily distinguishable. The first population has high magnitude (SNR) values while the second has low magnitude (SNR) values. In Fig. 5B, the observed phase contains two overlapping populations of voxels that are much more difficult to visually distinguish. The first population of voxels has phase values that are distributed closely around zero. The second population of voxels has phase values that a more uniformly distributed in the −π to π interval. The marginal distribution of the phase is uniform when ρ=0. In Fig. 5C the two populations of voxels are still visually present but appear to have greater separation. In Fig. 5D, there appears to be slightly fewer voxels with phase values near ±π and more near zero. In Fig. 5E, it can be seen that most voxels have a small estimated variance but there are some with very large variances that span the entire horizontal axis. In Fig. 5F, it is obvious that there are F statistic values from two different distributions. The first distribution is for large values of the F statistic (corresponding to large magnitude and or phase voxels) that tapers for smaller F statistic values. The second distribution is on the smaller side for smaller F statistic values (corresponding to small magnitude and or phase voxels) that tapers for larger F statistic values. Note the similarity between this second distribution for small F statistic values and the distribution of F statistic values in Fig. 3 when no magnitude or phase signal is present (i.e. the null hypothesis is true).
The map of these F statistics in Fig. 4F is thresholded with the critical values in Table 2 for α = 0.05, 0.0001, and 0.05/512/352. A zero-one mask is produced from the thresholded F statistics then applied to the original magnitude and phase images in Fig. 4A and Fig. 4B respectively. Thresholded original observed images are presented in Fig. 6. The magnitude and phase of thresholded voxels are set to zero but for display the thresholded phase voxel values are set to -π. In Fig. 6A and Fig. 6B the images are thresholded at F=2.8102, in Fig. 6C and Fig. 6D the images are thresholded at F=6.1512, in Fig. 6C and Fig. 6D the images are thresholded at F=7.5869. A vertical line can be drawn in Fig. 3 (that contains the distribution of F statistics under the null hypothesis) and in Fig. 5F (that contains the observed F statistics with both null and alternative hypothesis statistics) for each of these threshold values. Note in Fig. 6 that as the false positive rate decreases from Fig. 6A and Fig 6B to Fig. 6E and Fig. 6F, the number of above threshold voxels outside of the head decreases but slightly more voxels within the head are also eliminated. This phenomenon is due to the relationship between Type I and Type II error rates. It is apparent that the magnitude image in Fig. 6E shows similar anatomical to the phase image in Fig. 6F indicating similar biological information.
One of the more recent approaches to suppress image noise also used a complex threshold method (CTM) . This method used not only complex thresholding but also connectivity to enhance suppression of noise or prevent the incorrect assignment of signal to noise in order to reduce Type I errors. We do not use any connectivity in this approach since we are looking at local variance on a pixel by pixel basis. Also, the previous CTM  suffers when the phase itself deviates from zero or any set zero point. This can lead to the suppression of areas where the phase is offset from flow effects or susceptibility effects (the later being acceptable in some cases when the goal is to also suppress veins). However, the phase has a rapid variation near the air/tissue interfaces, particularly at the front of the head in Fig. 4B for example. In that area the CTM method  would fail and that part of the brain will be suppressed. Not so with this approach where we look specifically at the local variance and keep the local value as the phase offset (i.e., there is not a global offset of zero in phase). This method can also easily be adapted for 3D images.
In summary, a magnitude and phase statistical thresholding procedure based upon a likelihood ratio test was presented. It was shown through Monte Carlo simulation that this method operates according to its theoretical statistical properties in terms of both false positives and false negatives. This statistical thresholding method was successfully applied to real human SWI data and shown to produce increased tissue contrast by eliminating false positives.
This work was supported in part by NIH EB000215 and EB007827.
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.