Abstract
A segmentation algorithm to isolate areas of ventilation from hyperpolarized helium-3 magnetic resonance imaging (HP 3He MRI) is described. The algorithm was tested with HP 3He MRI data from four healthy and six asthmatic subjects. Ventilated lung volume (VLV) measured using our semiautomated technique was compared to that obtained from manual outlining of ventilated lung regions and to standard spirometric measurements. VLVs from both approaches were highly correlated (; ) with a mean difference of 3.8 mL and 95% agreement indices of −30.8 mL and 38.4 mL. There was no significant difference between the VLVs obtained through the semiautomatic approach and the manual approach. A Dice coefficient which quantified the intersection of the two datasets was calculated and ranged from 0.95 to 0.97 with a mean of 0.96 ± 0.01 (mean ± SD). VLVs obtained through the semiautomatic algorithm were also highly correlated with measurements of forced expiratory volume in one second (FEV1) (; ) and forced vital capacity (FVC) (; ). The technique may open new pathways toward advancing more quantitative characterization of ventilation for routine clinical assessment for asthma severity as well as a number of other respiratory diseases.
1. Introduction
Recent advancements in hyperpolarized helium-3 magnetic resonance imaging (HP 3He MRI) enable direct visualization of ventilation in the lung [1, 2]. While normally ventilated lungs have been found to exhibit a homogeneous distribution of gas signal, obstructed lungs such as in asthma show areas of signal depletion, often referred to as ventilation defects [1–5]. It is increasingly accepted that quantifying spatial patterns in the ventilation distribution can provide rich insight on the severity of asthma and how well a specific patient responds to a prescribed therapy [1, 2]. Additionally, such information may provide novel perspectives in the fundamental nature of asthma with regard to whether it is a localized airway pathology or a global lung disease.
Traditional analysis of HP 3He MRI has primarily been qualitative in nature, largely restricted to a scoring system that required a radiologist to visually estimate the number of ventilation defects [1–4]. These approaches were subjective and were likely inconsistent and time intensive. A number of quantitative methods have emerged for the segmentation of ventilated airspaces. Initial attempts by Kauczor et al. [6] relied on a thresholding scheme which assumed a Gaussian distribution of noise. However, such an assumption leads to an approximately 60% underestimation of the true noise power [7]. Later efforts by Tzeng et al. [5] and Woodhouse et al. [8] applied a threshold value that relied on a signal-to-noise threshold but still required rigorous manually outlined lung boundaries. More recent work using class-based algorithms with lung partitioning using a Gaussian mixture model [9] and methods that employ fuzzy C-means and K-means clustering [10–12] have also been introduced. These methods were automated but required additional manual removal of the trachea. By default, the trachea and associated large airways comprise a majority of the anatomic dead space which contains the largest percentage of HP 3He gas [5]. Since our goal was to target gas exchange regions, removal of the trachea and associated large airways would result in a more accurate assessment of ventilation.
In this study, we introduced a robust, semiautomatic algorithm for rapid segmentation of HP 3He MRI into distinct regions based on ventilation. The ventilated lung volume (VLV) quantified using our method was compared to that measured using a conventional manual analysis by a trained technician to determine the accuracy of our segmentation. As spirometry still remains as the gold standard for measurement of airway obstruction, we compared measurements of lung volume from HP 3He MRI using our method to forced expiratory volume in one second (FEV1) and forced vital capacity (FVC). The scope of this paper is to introduce the methodology and a preliminary study with data from four healthy and six asthmatic subjects. The intent is to provide proof-of-principle, in a fashion that indicates the capability of this approach in analyzing spatial distributions for ventilation [5, 12, 13] and future modeling studies [14, 15] for asthma [1–5] as well as the potential to be streamlined to other respiratory diseases such as chronic obstructive pulmonary disease [8, 11, 12] and cystic fibrosis [11, 16, 17].
2. Materials and Methods
2.1. Subject Enrollment
The Health Insurance Portability and Accountability Act-Compliant research protocol in this study was approved by both Boston University and Brigham and Women’s Hospital Institutional Review Boards. Written informed consent was obtained from all recruits, which consisted of four healthy subjects (two men and two women: age range 21–23 years; mean age 22 years) and six asthmatic subjects (one man and five women: age range 19–23 years; mean age 22 years). Before the first study visit, each subject participated in a screening day visit during which a methacholine challenge was administered to determine a PC20 dose that elicited a 20% drop in baseline FEV1. This index was used to separate healthy from asthmatic subjects. For our protocol, healthy subjects were nonsmokers with no history of respiratory diseases and exhibited PC20 values of >25 mg/mL. Asthmatic subjects consisted of those with a history of asthma who exhibited PC20 values of <8 mg/mL. The demographics are detailed in Table 1.
2.2. Image Acquisition Protocol
Standard spirometry measurements were recorded with the subject in supine position. Each subject was instructed to inhale a ~1 liter mixture of ~33% HP 3He-67% N2 from functional residual capacity (FRC). Images were acquired on a General Electric Signa LX 1.5 MRI scanner equipped with a heterodyne system which included frequency mixers to image at the 3He NMR frequency of 48.65 Hz. The system interfaced with a flexible quadrature lung coil (Clinical MR Solutions, Brookfield, WI) tuned to the same frequency. Hyperpolarization of the 3He gas was initiated through a collision spin exchange with vaporized rubidium optically pumped using a custom-built polarizer. The scans employed a Fast Gradient Echo pulse sequence that compiled coronal multislice images with a field of view (FOV) of 46 cm, matrix dimensions (zero-padded to ), 13 mm slice thickness, 0 mm gap between slices, 1.8 mm in-slice resolution, 31.25 kHz bandwidth, 14–18° flip angle, TE/TR 1.228 ms/50–75 ms, and interleaved data acquisition. Typically, 8–14 slices were obtained for each subject, depending on the anterior to posterior depth of the lung.
2.3. MR Image Processing
A detailed schematic of our semiautomatic segmentation method is illustrated in Figure 1. Our methods will refer to various panels in Figure 1. There are three steps to our semiautomatic segmentation method. (1) A preprocessing routine is applied involving statistical noise subtraction. (2) The image pixels are correspondingly clustered into ventilation classes to refine our initial segmentation. (3) The trachea and major airways are removed to obtain a final binary image representative of ventilated airspaces.
2.3.1. Statistical Noise Subtraction
HP 3He MR images were first preprocessed through a denoising scheme by determining an optimal threshold from a sampled background noise distribution located outside of the lung field. This space comprised an automated pixel box in the bottom center of each image slice (Figures 1(a) and 1(b)). The distribution is fitted through a nonlinear regression with an adjusted Rayleigh curve where is the intensity of background noise, with parameters and . In contrast to a similar technique previously applied in brain tissue segmentation [18, 19], our approach employed an additional shifting parameter, , which accounted for horizontal shifts in curve-fitting and provided a much stronger fit to the sampled data. An optimal threshold, , was subsequently derived from the minimization of an error term where the function, , constituted the subtracted distribution calculated by the difference between the best-fit adjusted Rayleigh curve, , and the pixel intensity distribution of the sampled background noise, . Consider
The purpose of the preprocessing was to automatically remove discernible sites of noise artifacts to construct an initial binary mask (see Figure 1(a)).
2.3.2. Segmentation Refinement
Pixel intensities across the entire image space were correspondingly partitioned through a clustering scheme. Here, we describe the clustering using fuzzy C-means (FCM) clustering [10, 11]. However, this step can also be replaced by a K-means clustering algorithm [12, 13] as both these algorithms use the same cost function. Briefly, the algorithm initializes four random cluster centers in which a corresponding membership function, , is calculated. The membership function is based on a distance measure which describes the degree of similarity between each data point and each cluster center given by where is the number of distinct clusters, and is a weighing parameter used to control the level of fuzziness in the classification scheme, typically initialized to 2 [20]. The variables, and , constitute the distance between point to the cluster center of clusters and , respectively. From the resultant calculation of the membership, , a new cluster center for each class, , is calculated across all data points, given by the following relationship [20]:
Using these new cluster centers, , the membership is updated, and the process is iteratively repeated, based on minimization of the following objective function [20]:
A predefined criterion, , between 0 and 1, is set such that when reached, the algorithm is terminated. Following previous publications [10, 11], we split the ventilated lung region into four clusters that corresponded to negligible ventilation, low ventilation, intermediate ventilation, and high ventilation (Figures 1(c) and 1(d)). Since our eventual areas of interest comprised ventilated regions within the lung, pixel intensities designated to the negligible ventilation class were treated as part of the background.
2.3.3. Semiautomatic Trachea Removal
Typically, we acquired 8–14 image slices anterior to posterior for each subject. Some of these images contain the trachea and the main stem bronchi, particularly in the middle slices, which, by default, hold the largest percentage of HP 3He gas. Since the trachea and associated large airways are not directly involved in gas exchange, it became crucial to remove them for an accurate assessment of ventilation (Figure 1(e)). Therefore, we employed a slice-by-slice seeded region-growing algorithm [21, 22] which detected edges based on the intensity levels of connected pixels.
In its simplest form, the technique requires an initiation point, known as a seed, which is often manually selected by the user [21, 22]. Each connected component of the seed is then flagged, and a difference measure is calculated by a predefined criterion at each iteration. The goal of the algorithm is to enable a final segmentation of regions as homogeneous as possible while constrained by each pixel’s connectivity to the initial seed point. A basic model was defined by Adams and Bischof, using a running mean calculated at each iteration starting at a designated seed point [21]. Given as the set of all as-yet unallocated pixels which border at least one of the regions,
The difference measure, , bound by the running mean is described by the following expression, where is the gray value of the image point . Consider
From the set of unallocated pixels, , which border at least one of the regions connected to the seed point, a minimum distance, , was set as the segmented space [21]. Consider
A detailed schematic is illustrated in Figure 2. The technique was knowledge based and required two inputs: a user-defined bounding box to limit the processing space to the trachea and an initial seed point composed of a single pixel manually selected inside the trachea (Figures 2(a)–2(c)). Following interrogation of each pixel within the isolated bounding box, an outlined space was obtained (Figure 2(d)). Each element within the image space was correspondingly labeled based on connectivity (Figure 2(e)), and an area filter was applied to isolate the trachea (Figure 2(f)). The area filter was based on pixel connectivity in which connected areas of fewer than 50 pixels were selectively removed. A simple binary subtraction between the input (with the trachea) and output images (without the trachea) yielded a binary image with the trachea selectively removed (Figure 2(g)).
We discarded residual artifacts from the crude binary image subtraction through morphological operations. A binary erosion (Figure 2(h)) was first used. Then, each connected element was labeled, and another area filter was applied to selectively remove connected areas of fewer than 50 pixels (Figure 2(i)). Finally, a binary dilation was applied (Figure 2(j)). The rationale was to target weakly connected areas usually comprising the larger associated airways extending from the main stem bronchi. The details on the operators are outlined by Serra [23]. Briefly, given the mask as a discrete Euclidian image, , dilation of by a structural element, , is expressed as follows:
The erosion of by is given as
For our processing scheme, we used a disk structural element for both binary erosion and binary dilation. We applied these operations to a fixed template based on the maximum width and height of the segmented lung slice as illustrated in Figure 3. This was done to maximize removal of the larger attached airways while minimizing morphological distortion particularly along the concave lung base.
2.4. Statistical Data Analysis
To assess the accuracy of our method with the manual analysis currently employed to assess ventilation heterogeneity from HP 3He MRI, all images were processed by a trained lab technician (5 years experience with HP 3He MRI) using a MATLAB-coded software (MathWorks, Natick, MA). In the manual analysis, lung contours and ventilation defects were outlined manually, and the VLV was calculated from the number of the pixels in the regions identified as being ventilated. A paired t-test was used to compare the VLVs from the manual and semiautomated methods. An unpaired t-test was used to compare the VLVs between the healthy and the asthmatic subjects. A Dice coefficient was also calculated to measure the agreement or similarity between the VLV using our approach, , and the manual approach, [24]. The Dice coefficient ranges from 0 to 1, with 1 indicating perfect agreement. Consider
We performed a linear regression analysis across all subjects and calculated a correlation coefficient and the slope between semiautomatic and manual methods. Bland Altman analysis [25] was used to determine the 95% limits of agreement calculated from the mean and standard deviation of the volume difference between the two methods of segmentation. VLVs through our algorithm were compared to PFTs, specifically, to FEV1 and FVC, functional measures that vary with the size and level of lung obstruction. Finally, an unpaired t-test was also used to compare FEV1 and FVC between the healthy and asthmatic subjects.
3. Results
A typical segmentation of the ventilated regions in the HP 3He MRI into three distinct classes of ventilation is shown in Figure 4. In our limited subset of six asthmatics and four healthy subjects, the asthmatics showed increased predominance in pockets of hypointense areas indicating low ventilation in the asthmatics as compared to the healthy subjects. We compiled a total of 109 coronal slices for ten subjects. For healthy subjects, segmentation with the semiautomatic approach yielded a mean VLV of L (mean ± SD), compared to the manual approach, which gave a mean VLV of L. For asthmatics, the mean VLVs of the semiautomatic approach were L and L with the manual approach. There was not a statistically significant difference in FEV1 between the semiautomatic and manual approaches (). The resulting dice coefficients for each subject are illustrated in Table 2. The coefficients ranged from 0.95 to 0.97 with a mean of .
(a)
(b)
(c)
(d)
Across each coronal slice of the lung, VLV measurements obtained through both methods were highly correlated (; ; ) (Figure 5(a)). From the Bland-Altman analysis, the mean VLV difference was mL. The lower and upper 95% limits of agreement were −30.8 mL and 38.4 mL, respectively (Figure 5(b)). Comparison to spirometry yielded a high correlation to measurements of FVC (; ; ) (Figure 6(a)) and FEV1 (; ; ) (Figure 6(b)). For healthy subjects, mean FEV1 was L and mean FVC was L; for the asthmatic subjects, mean FEV1 was 3.23 ± 0.69 L and mean FVC was L. There was no statistically significant difference in FEV1 () and FVC () between the healthy and asthmatic subjects.
(a)
(b)
(a)
(b)
4. Discussion
To this day, HP 3He MRI has confirmed and advanced a number of new perspectives in asthma. For one, when exposed to airway smooth muscle provocation, the lungs will constrict heterogeneously with the number and size of ventilation defects directly correlating to the level of clinical severity [1–4]. In cases of very severe asthma, heterogeneously distributed ventilation defects may be even present at baseline [1, 2]. More recently, there is even some evidence that the size and location of many of these ventilation defects in asthmatic lungs tend to not change with time or repeated bronchoconstriction [1, 3]. These notions are primarily qualitative as they are based on visual inspection of ventilation images. However, together they raise intriguing clinical and structure-function questions regarding whether one could apply a quantitatively robust method for diagnosing the severity of baseline asthma and for evaluating the efficacy of treatment.
Up until now, much effort has been devoted to extracting detailed structural information from HP 3He MRI. While qualitative methods of analysis [1–4] have raised the concerns about consistency, quantitative methods [5, 6, 8–13] have paved new insights in the characterization of ventilation. However, many of these segmentation approaches did not include an extraction and removal of the trachea and mainstem bronchi with associated large airways [5, 6, 8, 10–13]. To our knowledge, our approach is the first to enable both a segmentation of ventilated airspaces and a direct selective removal of these components that constitute the anatomic dead space. When we applied our method to just a small pilot-study number of healthy and asthmatic subjects, we did not find a statistically significant difference in the VLV between our semiautomatic method and our manual tracings that served as the ground truth for our analytical comparisons. The results obtained from our method were able to produce high correlations to those obtained through manual processing and showed high degree of similarity and agreement through the Dice coefficients and Bland Altman analysis, respectively.
However, there are some limitations to our technique. Static scans do not represent real-time ventilation but instead represent snapshots in real time. True ventilation would require a multibreath technique [26] necessitating a greater amount of 3He gas. Hence, our technique was only capable of calculating a total VLV based on the gas distributive patterns at breath-hold. Another limitation is the dependence in user input in the semiautomatic trachea removal. Varying the size of the bounding box can certainly impact the semiautomatic trachea removal algorithm. A crucial element of the approach is application of an area filter that isolates the trachea. If a large region outside the trachea was chosen, then it would indeed be more difficult to adjust the area filter to discard small regions of connected pixels and large regions of connected pixels as opposed to the status quo of just removal of small regions. The seed point, we believe, should not impact the trachea because, for the most part, it is nearly homogenous in signal intensity. However, sensitivity studies in varying both the seed point and the bounding box can certainly be done in the future. A final limitation was in the thickness of the slices. Because these images were acquired at breath-hold, thick slices of 13.13 mm were compiled in order to cover the entire extent of the lung while trying to minimize discomfort suffered by the subject. To compare lung volumes between each subject, we recommend the use of 1H proton MRI scans to determine the volume of the thoracic cavity to normalize for lung size.
5. Conclusion
In conclusion, our work outlines a novel statistically and quantitatively driven imaging analysis that may provide a powerful and valuable additional tool for the clinical assessment of asthma severity. With the emergence of modeling approaches to combine imaging modalities to construct patient-specific models [14, 15], segmentation of lung ventilation becomes more important than ever. These methods may provide new perspectives in structure-function relations and hold the potential to be extrapolated to other respiratory diseases.
Acknowledgments
The authors would like to thank W. C. Karl, J. Konrad, Y-S. Tzeng, and J. P. Roche for their helpful advice and knowledge in image processing. The authors have been supported by the National Institutes of Health.