Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2763057

Formats

Article sections

Authors

Related links

Magn Reson Imaging. Author manuscript; available in PMC 2010 November 1.

Published in final edited form as:

Published online 2009 June 23. doi: 10.1016/j.mri.2009.05.008

PMCID: PMC2763057

NIHMSID: NIHMS127510

The publisher's final edited version of this article is available at Magn Reson Imaging

See other articles in PMC that cite the published article.

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 [4]. 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

$$\begin{array}{cc}\hfill {y}_{R}& =\rho \phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\theta +{\epsilon}_{R}\hfill \\ \hfill {y}_{I}& =\rho \phantom{\rule{thickmathspace}{0ex}}\mathrm{sin}\theta +{\epsilon}_{I}\hfill \end{array}$$

(1)

where *y _{R}* and

$$p({y}_{R},{y}_{I}\rho ,\theta ,{\sigma}^{2})=\frac{1}{\sqrt{2\pi {\sigma}^{2}}}\mathrm{exp}\left[-\frac{({y}_{R}-\rho \phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\theta )}{2{\sigma}^{2}}\right]\frac{1}{\sqrt{2\pi {\sigma}^{2}}}\mathrm{exp}\left[-\frac{({y}_{I}-\rho \phantom{\rule{thickmathspace}{0ex}}\mathrm{sin}\theta )}{2{\sigma}^{2}}\right].$$

(2)

Upon conversion from Cartesian coordinates to polar coordinates in Eq. (2), the joint distribution of the observed magnitude and phase (*m*,) is

$$p(m,\mathrm{\rho ,\theta ,{\sigma}^{2})=\frac{m}{2\pi {\sigma}^{2}}\mathrm{exp}\{-\frac{1}{2{\sigma}^{2}}[{m}^{2}+{\rho}^{2}-2\rho m\phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}(\mathrm{-\theta )}]\}.}$$

(3)

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 (*m _{1}*,

$$L(\rho ,\theta ,{\sigma}^{2})={\left(2\pi {\sigma}^{2}\right)}^{-n}\left[{\Pi}_{i=1}^{n}{m}_{i}\right]\mathrm{exp}\{-\frac{1}{2{\sigma}^{2}}\sum _{i=1}^{n}[{m}_{i}^{2}+{\rho}^{2}-2\rho \phantom{\rule{thickmathspace}{0ex}}{m}_{i}\phantom{\rule{thickmathspace}{0ex}}\mathrm{cos}\left({i}_{-}\right]\}.$$

(4)

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 [5]. 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 H_{0},” the second row is labeled “Do Not Reject H_{0},” the first column is labeled “H_{0} True,” and the second column is labeled “H_{0} 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 H_{1}: ρ>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 H_{0}: ρ=0, θ=0, the maximum likelihood estimators (MLEs) for the magnitude and phase from Eq. (4) are

$$\stackrel{~}{\rho}=0,\stackrel{~}{\theta}=0,\phantom{\rule{thickmathspace}{0ex}}{\stackrel{~}{\sigma}}^{2}=\frac{1}{2n}\sum _{i=1}^{n}\left({y}_{Ri}^{2}+{y}_{Ii}^{2}\right)$$

(5)

while under the unconstrained alternative hypothesis H_{1}: ρ>0,θ≠0 the MLEs from Eq. (4) are

$$\widehat{\rho}={\left[{\left({\stackrel{\u2012}{y}}_{R}\right)}^{2}+{\left({\stackrel{\u2012}{y}}_{I}\right)}^{2}\right]}^{{\scriptstyle \frac{1}{2}}},\phantom{\rule{thickmathspace}{0ex}}\widehat{\theta}={\mathrm{tan}}^{-1}\left[\sum _{i=1}^{n}{y}_{Ii}/\sum _{i=1}^{n}{y}_{Ri}\right],\phantom{\rule{thickmathspace}{0ex}}{\widehat{\sigma}}^{2}={\scriptstyle \frac{1}{2}}\sum _{i=1}^{n}\left({y}_{Ri}^{2}+{y}_{Ii}^{2}\right)-{\scriptstyle \frac{1}{2}}\left[{\left({\stackrel{\u2012}{y}}_{R}\right)}^{2}+{\left({\stackrel{\u2012}{y}}_{I}\right)}^{2}\right]$$

(6)

where ${\stackrel{\u2012}{y}}_{R}$ is the mean of the real channel measurements and ${\stackrel{\u2012}{y}}_{I}$ 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 $\lambda =L(\stackrel{~}{\rho},\stackrel{~}{\theta},{\stackrel{~}{\sigma}}^{2})/L(\widehat{\rho},\widehat{\theta},{\widehat{\sigma}}^{2})$ 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 H_{0} and H_{1}, (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 $\lambda =L(\stackrel{~}{\rho},\stackrel{~}{\theta},{\stackrel{~}{\sigma}}^{2})/L(\widehat{\rho},\widehat{\theta},{\widehat{\sigma}}^{2})$ taken. Some terms in the likelihood ratio can be cancelled to become

$$\begin{array}{cc}\hfill \lambda & =\frac{{\left[\frac{1}{2n}\sum _{i=1}^{n}{m}_{i}\right]}^{-n}}{{\left[\frac{1}{2n}\sum _{i=1}^{n}{m}_{i}-\frac{1}{2}{\widehat{\rho}}^{2}\right]}^{-n}}\hfill \\ \hfill {\lambda}^{{\scriptstyle \frac{1}{n}}}& =\frac{\left[\frac{1}{2n}\sum _{i=1}^{n}({y}_{Ri}^{2}+{y}_{Ii}^{2})-\frac{1}{2}\left[{\left({\stackrel{\u2012}{y}}_{R}\right)}^{2}+{\left({\stackrel{\u2012}{y}}_{I}\right)}^{2}\right]\right]}{\left[\frac{1}{2}\sum _{i=1}^{n}({y}_{Ri}^{2}+{y}_{Ii}^{2})\right]}.\hfill \end{array}$$

(7)

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

$$F=\left(\frac{n\left[{\left({\stackrel{\u2012}{y}}_{R}\right)}^{2}+{\left({\stackrel{\u2012}{y}}_{I}\right)}^{2}\right]/{\sigma}^{2}}{2}\right)/\left(\frac{\left[\sum _{i=1}^{n}{y}_{Ri}^{2}+\sum _{i=1}^{n}{y}_{Ii}^{2}\right]/{\sigma}^{2}}{2n}\right)$$

(8)

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 *x _{1}* and

$$\begin{array}{c}\hfill E({x}_{1}{x}_{2})=E\{\frac{1}{n{\sigma}^{2}}\left[{\left(\sum _{i=1}^{n}{y}_{Ri}\right)}^{2}+{\left(\sum _{i=1}^{n}{y}_{Ii}\right)}^{2}\right]\frac{1}{{\sigma}^{2}}\left[\sum _{i=1}^{n}{y}_{Ri}^{2}+\sum _{i=1}^{n}{y}_{Ii}^{2}\right]\}\hfill & \hfill & =\frac{1}{n{\sigma}^{2}}E\left\{{\left(\sum _{i=1}^{n}{y}_{Ri}\right)}^{2}\sum _{i=1}^{n}{{y}_{Ri}}^{2}+{\left(\sum _{i=1}^{n}{y}_{Ri}\right)}^{2}\sum _{i=1}^{n}{{y}_{Ii}}^{2}+{\left(\sum _{i=1}^{n}{y}_{Ii}\right)}^{2}\sum _{i=1}^{n}{{y}_{Ri}}^{2}+{\left(\sum _{i=1}^{n}{y}_{Ii}\right)}^{2}\sum _{i=1}^{n}{{y}_{Ii}}^{2}\right\}\hfill \\ \hfill & =\frac{1}{n{\sigma}^{2}}E\left\{\left(\sum _{k=1}^{n}{y}_{Rk}\sum _{j=1}^{n}{y}_{Rj}\right)\sum _{i=1}^{n}{{y}_{Ri}}^{2}+\left(\sum _{k=1}^{n}{y}_{Rk}\sum _{j=1}^{n}{y}_{Rj}\right)\sum _{i=1}^{n}{{y}_{Ii}}^{2}+\left(\sum _{k=1}^{n}{y}_{Ik}\sum _{j=1}^{n}{y}_{Ij}\right)\sum _{i=1}^{n}{{y}_{Ri}}^{2}+\left(\sum _{k=1}^{n}{y}_{Ik}\sum _{j=1}^{n}{y}_{Ij}\right)\sum _{i=1}^{n}{{y}_{Ii}}^{2}\right\}\hfill \\ \hfill & =\{(n+2)+n+n+(n+2)\}.\hfill \end{array}$$

The correlation between *x _{1}* and

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 *x _{1}* and

CDF from Monte Carlo simulation (dashed lines) and F distribution (solid lines) for *n=5* (red), *9* (black), *25* (yellow), *50* (green), *100* (blue), *250* (violet).

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 [8]. 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×10^{7} 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×10^{7})^{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 [9]. 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 [10].

The magnitude and phase statistical thresholding method described in this paper was applied to the human susceptibility weighted imaging (SWI) MRI data [11]. 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 mm^{2}, TR=26 ms, TE= 15 ms, flip angle (FA) = 11° [8]. 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) [5]. 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 [5] 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 [5] 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.

[1] Henkelman RM. Measurement of signal intensities in the presence of noise in MR images. Med Phys. 1985;12:232–233. [PubMed]

[2] Bernstein MA, Thomasson DM, Perman WH. Improved detectability in low signal-to-noise ratio magnetic resonance images by means of phase corrected real reconstruction. Med Phys. 1989;16:813–817. [PubMed]

[3] Macovski A. Noise in MRI. Magn Reson Imaging. 1996;38:494–497. [PubMed]

[4] Rowe DB, Nencka AS, Hoffmann RG. Signal and noise of Fourier reconstructed fMRI data. J Neurosci Methods. 2007;159:361–369. [PubMed]

[5] Pandian D, Ciulla C, Haacke EM, Jiang J, Ayaz M. A complex threshold method for identifying pixels that contain predominantly noise in magnetic resonance images. J Magn Reson Imaging. 2008;28:727–735. [PubMed]

[6] Hogg RV, Craig AT. Introduction to Mathematical Statistics. Fourth Edition Macmillan Publishing Company; NY, NY: 1978.

[7] Bain LM, Engelhardt M. Introduction to Probability and Mathematical Statistics. Second Edition PSW-Kent Publishing Company; Boston, MA: 1992.

[8] Haacke EM, Brown RW, Thompson MR, Venkatesan R. MRI: Physical principles and sequence design. John Wiley and Sons; 1999.

[9] Logan BR, Rowe DB. An evaluation of thresholding techniques in fMRI analysis. NeuroImage. 2004;22:95–108. [PubMed]

[10] Logan BR, Geliazcova MP, Rowe DB. An evaluation of spatial thresholding technoques in fMRI analysis. Hum Brain Mapp. 2008;29:1379–1389. [PubMed]

[11] Haacke EM, Xu Y, Cheng YCN, Reichenbach J. Susceptibility weighted imaging (SWI) Magn Reson Med. 2004;52:612–618. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |