|Home | About | Journals | Submit | Contact Us | Français|
To develop a robust algorithm for tissue-air segmentation in MRI using the statistics of phase and magnitude of the images.
A multivariate measure based on the statistics of phase and magnitude was constructed for tissue-air volume segmentation. The standard deviation of first-order phase difference and the standard deviation of magnitude were calculated in a 3×3×3 kernel in the image domain. To improve differentiation accuracy, the uniformity of phase distribution in the kernel was also calculated and linear background phase introduced by field inhomogeneity was corrected. The effectiveness of the proposed volume segmentation technique was compared to a conventional approach which uses the magnitude data alone.
The proposed algorithm was shown to be more effective and robust in volume segmentation in both synthetic phantom and susceptibility weighted images of human brain. Using our proposed volume segmentation method, veins in the peripheral regions of the brain were well depicted in the minimum-intensity projection of the susceptibility weighted images.
Using the additional statistics of phase, tissue-air volume segmentation can be substantially improved compared to that using the statistics of magnitude data alone.
Tissue-air volume segmentation is a necessary pre-processing step for quantitative volume analysis and display of brain images in many applications. Since manual volume segmentation is laborious and can be operator biased, automatic and robust volume segmentation is highly desirable. One of the crucial tasks in volume segmentation of the brain is the reliable segmentation of the brain tissue from air and the skull. MR signal in the boney skull has a very low intensity in conventional MRI. This allows for reliable segmentation of brain tissue and skull. Subsequent removal of scalp tissue is thus made much easier because of the spatial separation between scalp and brain by the skull.
Tissue-air segmentation is commonly conducted by using the intensity in the magnitude images alone, due in large part to the lack of access to phase images on clinical scanners. Despite their simplicity, these magnitude-based approaches are, however, prone to several pitfalls. The intensity of a magnitude image can have a high spatial variation caused by tissue contrast, spatial variation of coil sensitivity, B0 inhomogeneity, and dielectric resonance effect. In susceptibility weighted imaging (SWI) (1–4), for example, veins and regions with rich iron deposits can have considerably reduced magnitude in the data. On the other hand, for images with a reasonably high signal-to-noise ratio (SNR), the distribution of phase is a weak function of intensity. In low SNR images, it was shown in (5) that the standard deviation (STD) of phase is inversely proportional to the SNR of the magnitude data (5). The phase distribution is expected to be uniform between −π to π in air and substantially narrower in tissue even at a low SNR. Inclusion of the statistics of phase is, therefore, expected to substantially improve the differentiation between tissue and air in volume segmentation. Nevertheless, several technical challenges arise from the use of phase for volume segmentations. In particular, the following two factors introduce substantial background phase variation in images: a) local field gradients in gradient-echo imaging affect regions with severe field inhomogeneities (6) and b) phase aliasing which can occur anywhere in the images.
Pandian et al. have recently proposed a complex threshold method to identify voxels that contain high noise levels in MRI data by applying thresholds to both the magnitude and phase images (7). The combination of magnitude and phase thresholds, along with the use of connectivity criteria, provides substantially improved differentiation between tissue and air. In this study, we present a novel multivariate approach which utilizes the statistics of both phase and magnitude data in the images for volume segmentation. The first-order phase difference (FPD), instead of the original phase in the complex images, is used to address technical challenges a) and b) outlined above. Moreover, a correction scheme is developed to remove the effect of background phase variation introduced by local field gradient. The effectiveness of the proposed algorithm is demonstrated by applying the algorithm to a synthetic phantom and 3D SWI dataset of the brain.
Conventional tissue-air segmentation algorithms use the difference of intensity between tissue and air in magnitude images. In this study, the statistics of magnitude and phase was calculated in a cubic kernel. A kernel of 3×3×3 voxels is illustrated in Fig. 1. A parameter of standard-deviation-to-mean ratio (SMR) normalized to the SMR in air was defined as:
where Si,n is the magnitude intensity of the nth voxel and is the averaged signal in the kernel centered at location i, N is the number of voxels in the kernel, meanair and STDair are the mean and standard-deviation (STD) of the magnitude intensity in air calculated in the kernel, respectively. It is noted that the ratio of meanair and STDair presents the SMR in air. For image data acquired with a single receiving coil, the ratio meanair/STDair is expected to be 1.912 (8, 9). The parameter SMR is expected to have a value near unity in air and a small value near zero in tissue. In contrast to the direct use of magnitude intensity of the images, such normalization enables us to define a threshold independent of signal intensity for volume segmentation. A kernel size of 3×3×3 voxels with N=27 was used in the rest of this study.
It has been observed that the phase of voxels has distinct statistical properties in tissue and air (7). The phase in air is expected to have a uniform distribution between −π and π. Conversely, the STD of phase in tissue, in the unit of radians, is given by:
where SNRmag is the SNR of the magnitude (5, 7). The phase in tissue is expected to be tightly distributed in a narrow range in regions with a reasonably high SNR. For example, the estimated σphase in tissue is less than 10% of that in air at SNRmag≥6.0.
As mentioned earlier, there are factors which contribute to additional background phase variation in images, thus resulting in higher STD values. One of the factors is phase aliasing, which commonly occurs in phase images acquired with gradient-echo techniques. Phase aliasing in tissue can possibly cause the phase values to get clustered near both ends of the (−π, π) range. Another factor is local field gradient (LFG) which can introduce a linear background phase and a bias in the STD calculation in gradient-echo images. We propose first-order phase difference (FPD) to circumvent both of the issues above.
An FPD vector can be presented as an arrow, with complex division of the voxel on the tip of the arrow by the voxel on the tail of the arrow. In the example shown in Figure 1, the FPD between the (j−1)th and jth voxels along the xy direction is given by:
where Ai,j−1 and Ai,j are the complex signal at the (j−1)th and jth voxel in the kernel, respectively; arg indicates an operation to take the polar angle of a complex number.
FPD is expected to be uniformly distributed between −π and π in air and is expected to be more tightly distributed near zero in tissue. The STD of FPD (stdFPD) at a voxel i, which is the center of the kernel, normalized by the stdFPD in air is given by:
where m is the index of the FPD pair, is the averaged FPD, and M is the number of FPD values calculated in the cubic kernel. Given a uniform distribution of FPD in a range of (−π, π), the STD of FPD in air is estimated to be . Therefore, stdFPDi is expected to have a value near 1/(0.577π·SNRmag) in tissue and near unity in air.
Although stdFPD is a good measure of the broadness of the phase distribution, it is not necessary a good measure of the uniformity of phase distribution. For example, a Gaussian distribution and rectangular distribution can have the same stdFPD value despite of their gross difference in the shape of distribution. In this study, a second parameter was employed to characterize the uniformity of phase distribution:
where Cout,i is the counts of FPD in the outer subdivisions (−π, −π/2) and (π/2, π) at the voxel i. Since the parameter θFPD is zero for uniformly distributed FPD and unity for tightly distributed FPD, θFPD is expected to have a value near zero in tissue and near unity in air.
In our approach, FPD values were calculated for two adjacent voxels in the kernel. Different approaches can be used in selecting the pairs of voxels in the kernel for the calculations of stdFPD and θFPD. A simple approach would calculate FPD between the central voxel and the other 26 voxels in a 3×3×3 kernel, resulting 26 FPD values. In this study, we selected a more thorough approach which used a total of 158 pairs of any two voxels among the 27 voxels: 18 pairs along x, 18 pairs along y, 18 pairs along z, 12 pairs along xy, 12 pairs along yz, 12 pairs along zx, 12 pairs along (−x)y, 12 pairs along (−y)z, 12 pairs along (−z)x, 8 pairs along xyz, 8 pairs along x(−y)z, 8 pairs along xy(−z), and 8 pairs along (−x)yz. Using such a large number of FPD values would reduce the variation of stdFPD and θFPD in the kernel.
The simple division in FPD directly resolves the issue with phase aliasing by subtracting the corresponding large values of phase. The stdFPD would thus have a lower value than the STD of the phase values.
As mentioned earlier, LFG can introduce a linear background phase in the kernel and a bias in the FPD calculation in gradient-echo images. The phase dispersion in a kernel along a direction u is given by
where γ is the gyromagnetic ratio with γ/2π=42.58MHz/T, LFGu is the LFG along the direction u, Δu is the dimension of the kernel, and TE is the echo time of data acquisition. LFGu can be up to 0.745mT/m at the orbitofrontal region (10) at 3T. With a TE=20ms commonly used in SWI at 3T, LFGz=0.5mT/m and a kernel size of 1.15×1.15×1.6 mm3 (twice of voxel size in a 3×3×3 kernel), the phase dispersion in the kernel along the z direction is 1.4π. Unless the effect of LFG is corrected, voxels in brain regions with severe field inhomogeneity could be mis-categorized as air. The correction of LFG in the calculation of FPD is described next.
The FPD along the LFG direction is presented by FPDi,LFG = arg(Ai,1/Ai,0). The FPD with the correction of LFG along this direction is calculated using complex division:
The corrected FPD′ is expected to have a narrow distribution near zero in tissue and a uniform distribution between −π and π in air. After correcting LFG along 13 directions, a total number of 145 FPD′ were obtained for the calculation of stdFPD and θFPD. It is noteworthy that the use of complex division in the calculation of FPD′ was necessary to circumvent phase aliasing and ensure that the range of FPD′ distribution remains to be (−π, π). If a simple subtraction of (FPDi,j − FPDi,LFG) was used, the distribution of FPD′ would be shifted by an amount of FPDi,LFG, resulting in erroneous elevation of the stdFPD and θFPD values.
A multivariate parametric Ω was constructed using the multiplication of stdFPD, θFPD, and SMR for air-tissue segmentation:
Ωi is expected to have a value near 1.0 in air and a small value near zero in tissue. The 1D histogram Ω presents a projection of the 3D histogram onto a line that connects between the peak of tissue and the peak of air if the tissue peak can be considered at the origin. Voxels in air and tissue are expected to be more distinguishable using the multivariate Ω than using any of the single parameter alone. A binary mask was generated by applying a threshold to the Ω map for tissue-air segmentation. In this study, the global image threshold was determined to minimize the intra-class variance of the black and white voxels (11).
The proposed algorithm was applied to a synthetic phantom that simulates the SWI data of the brain with a size of 288×384, as shown in Figure 2. This phantom contains a 2-pixel layer of CSF surrounding the brain (A), a small size vein with a width of one pixel (B), a median size vein with a width of two pixels (C), a major vein with a width of four pixels (E), and a pial vein on the surface of the brain with a width of 2 pixels (F). The veins have reduced intensity because of the increased T2* decay of venous blood compared to the brain tissue in SWI. The veins have negative phase shifts compared to the background phase of the brain which was set to zero. The smallest vein has the smallest shift and the major vein has the largest phase shift. A layer surrounds the upper-left quarter of the phantom has a reduced intensity without phase-shift, mimicking the intracranial cerebrospinal fluid (CSF) which has a reduced intensity because of the T1 saturation commonly found in SWI. The circular region (D) near the center of the brain represents the off-resonance effect at the orbitofrontal region of the brain in SWI (10). The amplitude reduction and phase shift both have a quadratic shape, as shown in (c) and (d). Seven levels of random Gaussian noise were added in the real (e) and imaginary (f) part of the phantom with the resulting SNR of 3.5, 4.0, 4.5, 5.0, 5.5, and 6.0, respectively, in the phantom. Three slices from the same phantom but with different instances of independently added random noise were generated for the calculation of the SMR, stdFPD, θFPD maps, since they require a 3×3×3 kernel for the calculations. Linear background phase was corrected in the calculation of stdFPD, θFPD maps. In this study, the type I and II errors were defined as follows (7):
The dependency of the type I and II errors on SNR was investigated.
A healthy volunteer was scanned with a 3D spoiled gradient-recalled echo (SPGR) pulse sequence on a GE 3T scanner (Milwaukee, WI) with a standard head volume coil with written informed consent approved by the local Institutional Review Board. The following imaging parameters were used in the scans: field of view (FOV)=22×17.8cm2, 1.0 mm slice thickness, 384×288 matrix, full echo acquisition, readout bandwidth=±5.6 kHz, and TE/TR/α=16ms/32ms/15°s. Sixty-four slices were acquired. The imaging time was 9 minutes and 53 seconds. Flow compensation was applied to the readout and phase-encoding directions to reduce the phase variation caused by flowing blood. Complex 3D image data were reconstruction off-line using MATLAB (MathWorks, Natick, MA). Six additional levels of random Gaussian noise were added to the in vivo data. The dependency of the type I and II errors on SNR was assessed.
Figure 3 shows the results of applying the proposed algorithm to the synthetic phantom data: the SMR map (a), the stdFPD map with the correction of linear background phase (b), the θFPD with the correction of linear background phase (c), and the Ω map (d). Figure 4 shows the histograms of the original magnitude (a), the SMR map (b), the stdFPD map (c), the θFPD map (d), and the Ω map (e), as shown in Figure 3. The air peaks in (b–e) were all near 1.0, due to the normalization of these statistic values to that in air. Figure 4f shows the details of the air and tissue distributions by removing the large tissue peak in (e). It was observed that the tissue and air peaks were much better separated in (f) than in (a) and (b).
A threshold of 0.58, near the minimal between the tissue and air peaks obtained by an algorithm described in ref. 11, was applied to the original magnitude (SNR=5.0) for volume segmentation. The resultant binary mask is shown in Fig. 5a, with the type I error of 0.0248 and the type II error of 0.0142. The type I error pixels were densely distributed at the regions of the veins and the region with severe off-resonance effect. A threshold of 0.46, near the minimal between the tissue and air peaks obtained by an algorithm described in ref. 11, was applied to the Ω map for volume segmentation. The resultant binary mask is shown in Fig. 5b, with the type I error of 0.0003 and the type II error of 0.0007. Most of the type I error pixels in (f) located on the major veins. The type I and II errors at other SNR levels are shown in Table 1. The total error with the use of the proposed algorithm was only a small fraction of that with the use of the conventional algorithm. The advantage of the proposed algorithm compared to applying a threshold to the original magnitude images was clearly demonstrated in this phantom study.
The proposed algorithm was applied to the entire 3D SWI dataset with 58 slices of the subject. Figure 6 shows the results in a slice located slightly superior to the Circle of Willis. The magnitude image is shown in Fig. 6a. SNR was 15.5 in a selected region of the brain. The arteries, veins, and CSF in the brain have higher intensity in the SMR map. The off-resonance effect caused elevated intensity at the orbitofrontal region in the stdFPD and θFPD maps, as indicated by the arrows in Figs. 6c and 6d, and can potentially introduce type I errors in volume segmentation. After the correction of the linear background phase, the off-resonance effect at the orbitofrontal region was effectively suppressed in the stdFPD and θFPD maps, as shown in Figs. 7a and 7b, respectively. The Ω map, as shown in Fig. 7c, has substantially improved tissue-air differentiation than any other maps shown in Fig. 6 and Fig 7. The computation time was 21 minutes for a matrix of 384×288×64.
The binary mask obtained by applying a threshold directly to the magnitude image contains many dark holes inside the brain due to the low intensity in veins and CSF, as shown in Fig. 8a. The binary mask obtained by applying a threshold to the Ω map, as shown in Fig. 8b, shows substantially reduced dark holes in the brain compared to that in Fig. 8a. In each of these two cases, the threshold was selected to be near the lowest point between the tissue and air peaks on the histogram. After removing the scalp tissue using a region growing algorithm, the resultant binary masks are shown in Figs. 8c and 8d, respectively. Figures 8e, 8g, and 8i show the binary masks obtained by applying a threshold directly to the magnitude images followed by a region-growing algorithm at reduced SNR level of 12.4, 9.3, and 6.2, respectively. Figures 8f, 8h, and 8j show binary masks obtained by applying a threshold to the Ω map followed by a region-growing algorithm at the same reduced SNR levels. The binary masks in the bottom row were constantly superior to that in the top row at all SNR levels.
A 3D standard binary mask with 58 slices was obtained with the original in vivo data using the proposed algorithm followed by a region growing procedure to remove the scalp and an imfill procedure provided in the MATLAB to fill in the holes introduced by major vessels and CSF. This standard binary mask was used to assess the type I and II errors in the 3D binary masks obtained without using the imfill procedure at different SNR levels. The type I error voxels include the ones shown as holes on major vessels and CSF. The region of scalp was excluded in the counting of the type I and II voxels. As shown in Table II, the use of the proposed algorithm resulted in substantially reduced type I and II errors compared to the conventional algorithm.
The mIP of the same 3D SWI dataset without tissue-air volume segmentation is shown in Fig. 9a. While the venous vasculature was well depicted in the middle and parietal regions of the brain in the projection, the veins in the frontal region and the peripheral region of the brain was largely obscured due to the passing of air region of the projection path (10). Tissue-air volume segmentation of the 3D complex SWI dataset was obtained by multiplying the dataset with the 3D binary mask. A region growing algorithm seeded in the brain tissue was applied to strip the scalp tissue and the skull. The residual holes in the 3D binary mask introduced by major vessels and CSF were removed by an imfill procedure provided in the MATLAB prior to volume segmentation. The mIP of the 3D volume segmented SWI dataset is shown in Fig. 9b. The lost signal at the peripheral region of the brain was effectively recovered using the proposed volume segmentation algorithm along with the scalp-striping and hole-filling procedures. The veins in the frontal and peripheral regions of the brain were well depicted in the mIP display of the 3D data.
The phase has different distribution patterns in air and tissue in magnetic resonance images. The use of phase information provides an additional measure for more accurate differentiation of tissue and air. The joint use of both magnitude and phase images has demonstrated to be helpful to improve the removal of noise pixels in magnetic resonance images (7). The concept of joint use of magnitude and phase images was extended in this study. In stead of using voxel-wise magnitude and phase values, we have used the statistical properties of magnitude and phase in a cubical kernel to further improve the differentiation of tissue and air. The calculations of SMR, stdFPD, and θFPD were normalized to the theoretically predicted values in air without additional measurement of SNR. The multivariate parameter Ω was, therefore, normalized and had an expected value of unity in air. With such normalization, it becomes feasible to use similar threshold values of Ω for tissue-air segmentation of images with a wide rage of SNR and contrasts.
The distribution of phase can be severely distorted in regions with severe field inhomogeneity in gradient-echo images (12,13). The signal intensity can be substantially reduced due to phase dispersion in these regions (6). These brain regions are prone to mis-categorization as air during volume segmentation because the voxels in these regions have similar statistics of phase and magnitude distributions as that in air. The linear background phase introduced by the LFG was effectively removed in the calculation of stdFPD and θFPD, resulting in more accurate tissue-air segmentation in regions with severe field inhomogeneity.
Different approaches can be used for the selection of voxel pairs in a 3×3×3 kernel for the calculations of stdFPD and θFPD. For example, the statistics of phase could be calculated using the 26 FPDs along 13 directions in the kernel. Such calculation, however, will make the correction of background linear phase difficult. In this study, we opted for the use of all the possible pairs of voxels for the calculation of FPD. In addition to the anticipation of improved accurate of statistics due to a larger number of samples, the linear background phase can be readily corrected by subtracting one FPD from the other FPDs along the same direction.
The proposed algorithm has superior performance compared to the conventional algorithm at all SNR levels in both phantom and in vivo data, as shown in Table 1–Table 2 and Fig. 8. The type I error was increased rapidly at reduced SNRs, suggesting that a hole-filling procedure is more necessary in the volume segmentation of data with a low SNR. Conversely, the effect of SNR on the type II error was much less than that on the type I error. It is noteworthy that a decrease in threshold would reduce the type I error and increase the type II error. The ratio of type I error and type II error can change considerably depending on the selection of the threshold. Table 1 and Table 2, however, show a clear advantage of the proposed algorithm in reducing the overall errors compared to the conventional algorithm.
In conventional mIP implementation, the minimum intensity in voxels along a projection path is selected for the display. Using mIP for the display of veins that have negative contrast introduce signal loss in regions where the projection encounters air or bone. As a result, the peripheral regions of the brain can be obscured in through-plane mIP display. Robust tissue-air segmentation is desirable to mitigate the signal loss in peripheral regions of the brain in the mIP display of 3D SWI data (10). The signal intensity in SWI, however, is highly heterogeneous due to the T2* signal drop in veins and regions with high iron concentration. The signal difference among white matter, gray matter, and CSF can also be substantial due to the heterogeneity of T2*. Manual volume segmentation can be a tedious task subject to errors due to the fact that many of the pial veins lay on the surface of the brain and the veins appear dark in the magnitude images. The joint uses of both magnitude and phase data in this study have shown to be effective in tissue-air segmentation of SWI data. With volume segmentation, signal loss in peripheral regions of the brain during in the mIP display of the SWI data was substantially mitigated. Since the pial veins are distributed on the surface of brain, the remnant edge errors in volume segmentation can introduce artifacts in the mIP display of 3D SWI data. The inclusion of air or bone in the mask can generate vein-like structures in the mIP display, whereas the exclusion of voxels on the surface of the brain in the mask can obscure some of the pial veins in the mIP display.
In summary, this study demonstrates that robust volume segmentation of 3D MR images can be obtained by the joint uses of the local statistics of both phase and magnitude distributions. The use of first-order phase difference in a small volume of 3×3×3 voxels circumvents the adverse effects of phase aliasing and local field gradient to the calculations of phase distribution. Robust air-tissue volume segmentation was achieved by applying the multivariate approach to 3D SWI data. In addition to volume segmentation of brain tissue, potential applications of the proposed algorithm include tissue-air segmentation of the extremities and the delineation of the spinal core and spinal nerves from their surrounding vertebral column.
We thank Yuzheng Hu at Zhejiang University and Dr. David Miller at the University of Colorado for their assistance in data analysis during the initial stage of this study. We thank Dr. Joseph Dagher for his assistance in editing the manuscript. ZJ is grateful to Prof. Ling Xia at Zhejiang University for his support.
Grant support: NIH grants MH68582, MH070037, MH47476, DA009842, MH079485, DA011015 (YPD).