[0]Address correspondence to eric-hoffman@uiowa.edu
2Dept. of Radiology
The University of Iowa, Iowa City, IA-52242
[0]Address correspondence to eric-hoffman@uiowa.edu
There has been an increased interest in automatic segmentation of volumetric medical image data. One of the reasons is that automated segmentation takes away the variability which exists when data is segmented manually. It also reduces processing time significantly. However, because of the stochastic nature of biological structures and the fact that no two data sets and scanner models are alike, it is very important to develop automated methods which process images in an adaptive manner and use a priori information to simplify the process. The method which we present here adaptively determines thresholds in order to segment out the primary human airway tree and uses some a priori information about the manner in which branching occurs, specifically the order in which the upward and downward branches arise from the right and left bronchi. We present preliminary results from this method, which automatically segments out the first four generations of the airway tree reliably, in data sets from both normal and airway compromised subjects and present comparisons with the current `gold standard' of manual segmentation.
Keywords: image segmentation, volumetric pulmonary imaging, region growing, virtual bronchoscopy
Reliable and timely assessment of lung anatomy and physiology is crucial, both in the clinical and research arenas. The high spatial resolution of volumetric HRCT, in conjunction with instrumentation to objectively control patient lung volume, has made it possible to resolve extremely small anatomic structures without the annoying motion artifacts. An area which is receiving increased attention and which forms the scope of this work is the automatic segmentation of pulmonary tree structures, specifically the primary airway tree.
Accurate segmentation of the airway tree is an essential step in so called Virtual Bronchoscopy[1, 2, 3, 4, 5, 6]: a newly popularized technique which volume renders airways in real time, allowing the clinician to visualize the inside of the airways in a bronchoscopic like paradigm. We make a strong distinction between visualization and detection. It is our belief that visualization without detection is of limited clinical use. Visualization with detection leads the clinician towards development of novel interventions. Geometric measurements of the airways are frequently made in order to study and plan surgical interventions for pathological conditions like tracheal and subglottic stenosis and include laser ablation of tumors, balloon dilatation procedures and stent design. Manual methods are available in the form of image processing tools like VIDA [7], to allow segmentation of a 3-D data set by performing 2-D analysis. Since manual extraction requires the user to interactively manipulate parameters such as region grow thresholds, it allows for introduction of operator specific variability. Furthermore, manual analysis is time consuming, with segmentation time being largely dependent on the skill level of the technician and data set size. Any automated method used to aid clinicians and researchers alike, needs to address two important criteria: accuracy and time taken for the procedure.
Conventional 26-connected region grow methods use a conservative fixed threshold and hence do not allow for accumulation of the entire airway lumen and also suffer from the inability to handle loss in airway connectivity, as illustrated in figure 1: a problem which may occur in real world situations and cannot be overcome by 26-connected region growing. The threshold value is empirical and does not take into account the fact that radiopacity for a given anatomic structure is not only machine dependent, but varies spatially. Partial volume effects and incompletely defined in-plane branch airway walls exacerbate the problem.
Figure 1: This figure indicates a case commonly encountered wherein airway connectivity is lost. (a) indicates an extremely narrowed trachea at the center. (b), (c), (d) and (e) demonstrates a complete loss of the trachea due to swallow. (f) shows a reappearance of the trachea. From this case it is apparent that conventional 26 connected region grow would fail.
Clearly, there is a crucial need to develop a robust method to allow accurate segmentation of the primary airway tree in patient data: a method which circumvents problems of the conventional hand tracing approach and does it in a manner to make it of clinical significance. We present a method which automatically extracts the primary airway tree, reliably, with minimum operator interaction and present results from clinical cases and comparisons with currently employed manual methods.
Any image processing technique is not only highly modality dependent, but also depends on the specific protocol followed during image acquisition. Pulmonary image data can be corrupted by cardiogenic artifacts, mis-registration between breath-holds and minimum patient motion (respiration during scanning yields unacceptable data sets). Furthermore, partial volume effects and finite slice spacing cause a reduction in resolvability of the modality, this problem is easily seen in the blurring of smaller airways, through incompletely defined airway walls when the branching is in-plane. Following a well defined protocol during image acquisition such as cardiac gating, lung volume control[8], yields data sets better suited for automated segmentation.
The role of pulmonary HRCT has been previously highlighted in studies [9] which have indicated that HRCT can compensate for short comings in conventional flexible bronchoscopy as well as offer a more objective view about pathology when compared to pulmonary function tests, which at best can provide a subjective global picture about the presence or absence of pathology. They have established the invaluable role which CT, in conjunction with image processing can play in surgical planning and pathology specific tasks like laser ablation and stent design. A high degree of correlation has been established between the diagnostic results of the two modalities ( Bronchoscopy and Helical CT) [10]. Results published by McLennan et al. [9] were generated by skilled operators, who extracted 2-D airways using established manual medical image processing tool sets [11]. Although this manual method, coupled with computer based adjustments of the borders, provides useful information, it does suffer from operator variability and long analysis times. A skilled technician requires one hour or more to segment airways of clinical significance from a 3 mm thick, 40 slice deep data set.
Previous work in the area of automatic airway segmentation includes simple 26-connected region growing to extract geometric information of branching in excised canine lungs[12]. However, the method exhibits shortcomings, since the threshold used needs to be as small as possible, governed primarily by the gray scale value of the weakest wall encountered, which would correspond to the gray scale value of the wall in in-plane airway branches.
Using a priori knowledge about the pulmonary anatomy and the anatomic and functional relationships between the pulmonary vasculature and the airway tree, peripheral airways[13] can be marked with reliable accuracy. However, this relationship may be undermined by disease and pathology and increased accuracy demands reduction of flexibility or development of more restrictive rules. D'Souza et al. [14] have demonstrated a 2-D based method which allows a user to segment out the airways on a slice by slice basis, by providing a seed point in the lumen. The results indicated a high degree of correlation in studies carried out with phantoms. However, increased interactiveness allows some inter-operator variability and has a considerable time cost. This served as the motivation to develop a robust method to segment the primary airway tree, with minimal operator interaction and in a reasonable amount of time to allow this method to be viable in the clinical setting and aid in the areas of patient management as cited by McLennan et al.[9]
The primary human airway tree can be modeled in most cases as a bipodal branching structure, as seen in figure 2. The main anatomic structures present within the lungs are the branches of the pulmonary circulation, lung parenchyma and the airways. The success of CT in imaging the thorax stems from the fact, that the resulting image provides excellent contrast between lung parenchyma, the blood filled vessels and the airways, because of the greatly different attenuation constants. The airways can be modeled as dark areas (air) surrounded by a bright wall. This statement is generally true, except in cases when the airway is oriented parallel to the imaging plane. In order to overcome this problem, it is essential to provide some means of adapting the region grow threshold based on local information and gray scale relationships between the air in the airway, the airway wall and the surrounding lung parenchyma.
Figure 2: schematic of primary airway tree, as used by pulmonologists at the University of Iowa Hospitals and Clinics. This figure of the branching structure of the normal primary airway tree indicates the order in which the LU and RU bronchi branch from the Left and Right main bronchus.
This method addresses these problems by using both adaptive methods and a priori knowledge regarding the anatomy and by trying to achieve a balance so as to be robust without being conservative in it's approach; allowing accumulations of as many conducting airways as possible and avoiding the pitfalls of marking non airway regions. The trachea is automatically detected using a predetermined threshold . The threshold selected is close to the gray level of air. Using the image geometry, connectivity analysis carried out near the center of the first few slices of the image space, by searching for a large dark region allows detection of the trachea. The centroid of these points is computed by averaging the X and Y coordinates of the region marked as the trachea. This centroid is used as the center of two rays which are cast out at right angles to each other. The gray level of the pixels in the original image which lie on this ray are examined to determine the maximum intensity. This maximum value indicates the gray level of the airway wall. However, in order to avoid spurious noise within the airway, from being detected as the airway wall, the intensity gradient profile is examined. A true airway wall exhibits a profile as shown in the figure 3 (a) and (b). In the case of the airway completely enclosed within the airway wall, four maximum and hence four half maximum values can be determined. For a given slice these four values should be close to each other. The average half maximum value and the standard deviation (SD) is computed. Half maximums which lies within 2 SD of the mean are considered legal values. If the value is beyond +2 SD, it is neglected in order to avoid region grow failure by whereby the grown region extends outside the desired branch.
Figure 3: Ray profile through airway (a) and (b). Note the distinct edges and the relatively flat intensity profile within the airway. Ray profile through parenchyma (c) and (d). Note variations and moderately increased intensity values. This results in a higher coefficient of variation along the profile, allowing discrimination of parenchyma.
The mean of the remaining values is used as the threshold for region growing. 2-D 8-connected region grow is used to accumulate all points within the airway. Biasing the threshold to the lowest half maximum value may result in a valid pixels not being marked as airway pixels. However, methods are available [15], to subsequently detect the correct edge. The mean gray scale of all these pixels is calculated, and each is projected to the next slice. Pixels in the next slice are considered as valid luminal pixels if their individual gray level is within +2 SD of the mean value of the accumulated pixels in the previous slice. Connectivity analysis is used to determine the number of contiguous objects present. These objects represent possible luminal pixels, except in cases where the resulting subset was due to pixels in the lung parenchyma having gray scale values within 2 SD of the average gray scale intensity in the previous slice. Considering such a region as a the lumen, will result in marking of the lung parenchyma as the airway. To avoid this pitfall, the resulting overlap regions are preprocessed to ensure that the pixels are those within the lumen. The method employed to do this involves, determination of the centroid of the overlap region.
This centroid is used as the center to cast two rays at right angles to each other. The texture of the pixels which lie on these rays is investigated. If the pixels belong to the lung parenchyma, the ray will demonstrate a continuously varying profile, because of the bright blood vessels which it crosses. If the subset of pixels lie within a valid airway, a smooth texture will be encountered with a gradual increase in intensity values near the airway wall. Thus overlaps which result between the previous airway pixels and those in the present slice can be discriminated as being within a valid airway or in the lung parenchyma. Cases may be encountered where no overlap exists because of a stenotic airway or a physiological activity such as swallowing. This situation is overcome, by skipping to the next slice. It is interesting to note that in this case the radiologist examining the data set, slice by slice, felt this to be the airway narrowing and missed the true stenosis. Also, to account for changes in the direction of the airway, the search area is increased with increasing slice numbers. This particular condition has been illustrated in figure 1. This process of threshold detection, region growing and branch point detection continues until all downward going branches are detected.
The detection of the upward going branches proceeds in the same manner. The only difference is that a priori knowledge regarding the branching structure of the primary airway tree is used to determine the slices from which to began the upward search. It is evident from the schematic of the airway tree in figure 2, that a search for the LU and RU bronchi should begin after the branching of the trachea into the right and left main bronchi. From the geometry of the airway tree, it is known that the second branch results in the Right and Left Upper bronchi. Also, a branch point can be characterized by a sudden increase in area and an increase in pixel variability within the airway. We use the largest area between the first and second branch point, detected during the first pass as the starting point to begin the upward search for airways. At each step statistics related to airway geometry are logged.
Development of this method was directed towards segmentation of the primary airway tree from a volumetric scan whose slice thickness was of the order of 3 mm, with no slice overlap and a scan aperture of 0.1 sec. This ensured sufficiently good image quality, with a minimum radiation dose. Previous attempts have involved volume (stacked thin sections) scans with slice thickness ranging from 2 to 3 mm, with an overlap equivalent to half the slice thickness, resulting in a significant increase in patient exposure. Data acquisition was ECG gated to reduce mis-registration due to cardiogenic motion. Additionally, scan aperture was limited to 100 msec to reduce cardiogenic blurring.
The image processing method described above was tested on both normal subjects and subjects exhibiting airway pathology. The pathology exhibited by the diseased subjects, included benign and malignant causes of bronchial and tracheal stenoses as well as bronchiectasis and airway wall thickening due to cystic fibrosis. We studied two patients diagnosed with bronchial and tracheal stenoses and one subject diagnosed with cystic fibrosis. The one data set of normal airway geometry came from a previous study, carried out to verify the relationship between airway reactivity and airway receptor density. The `normal subjects' exhibited no form of anatomically based airway pathology and were not in a state of airway constriction. The normal data set used in the study has a slice thickness of 3mm and the slices were contiguous, from just above the carina. Data were acquired using an Electron Beam CT, with cardiac gating and lung volume control, lung volume being held at 65 % vital capacity. The subjects exhibiting forms of airway pathology, were scanned using the Electron Beam CT, at total lung capacity, with cardiac gated image acquisition. All subjects were either scanned for clinical reasons or as part of other experimental protocols approved by the hospital IRB. The four test cases presented here were selected to demonstrate:
The accuracy of our method in segmenting out the airways can be measured at two levels
Table 1: Comparision of airways obtained by the manual and automatic method. bp1, bp2 are abbreviations for branch points. Results were obtained using a normal,a subject suffering from cystic fibrosis and 2 subjects suffering from tracheal stenosis.A 
represents area determined by our automatic method and A 
represents area determined from the manually traced airways ways. Pixel counting was employed to count number of lumen pixels. Conversion to metric units was achieved using pixel dimensions. 
and 
represent contours of manually and automatically segmented airways on each slice.
| Generation | Percentage Match |
|---|---|
| 1 | 100 |
| 2 | 100 |
| 3 | 100 |
| 4 | 80 |
| Airway Generation |
Euclidcan Distance mm |
|---|---|
| 1 | 0.865 |
| bp1 | 1.047 |
| 2 | 0.783 |
| bp2 | 1.236 |
| 3 | 0.461 |
| 4 | 0.583 |
To verify how closely our airway contours (C 
) matched the manually extracted contours (C 
), the mask images were thinned. Each point in C was compared with the points in C , to determine the minimum distance between each point C with contour C . The maximum was computed for each these minimas. Also, the mean distance between the two curves was computed. The results were subdivided based on the airway generation. Further, branch points were classified separately, to highlight the fact that the highest degree of mis-match occurs at a branch point. Simple pixel counting was used to compute the area of the airways in each slice for all subjects. Because airway dimensions vary from person to person based on age and other vital statistics a spread in airway areas is expected. However, the fact that areas with very small dimensions matched well with those traced with the manual method, serves as a clear indicator of the ability of our method to work well with airway compromised data sets. Table 1 is a compilation of these results.
Measurement of an accurate region growing threshold is a moot point, which determines the accuracy with which the method's results come close to the original airway area and shape. Our method used the half maximum criteria in cases where the airway was not oriented parallel with the imaging plane. However, in most cases the half maximum criteria is not an accurate measure of the presence of an airway wall. Reinhardt et al [15] have shown that accurate measurement of the airway wall is also a function of airway dimension, airway generation and characteristics of the scanner itself. It is thus possible to further process the resulting segmented airway to come as close as possible to the original airway tree.
Although the distance between contours shows degraded values at the branch points, table 3 clearly demonstrates the close match between centroids of the manually traced and automatically segmented airways. This further strengthens are belief that a threshold selected using additional information as cited in Reinhardt et al. [15], will show a dramatic improvement in the degree of match between the contours. Figure 4, is a surface rendered view of the airway tree segmented using our method. Close inspection of the images, will clearly reveal the pathology in cases where it exists. The maximum number of generations which our method detects is greatly dependent on slice thickness, as with other methods. The primary reason for this limitation is the fact that connectivity is lost after generation four or five. Our method will produce better results if slice thickness is reduced or if slice thickness is kept at 3mm with table increments of 1.5mm. In the present method only two perpendicular rays were used to determine the threshold and to detect the presence of a weak airway wall. Increasing the number of rays would have a positive impact on the method in terms of it's deterministic ability with reference to threshold detection. The down side of this would be increased processing time. The manner in which the method proceeds to segment out the various branches of the primary airway tree, makes it an excellent candidate for implementation of a parallel processing method.
Figure:4 (a) Airway tree from a normal subject segmented by this method, (c) Airway tree exhibiting stenosis and loss in tracheal connectivity as indicated in figure 1, segmented automatically by our method,(d) Airway tree of a subject suffering from cystic fibrosis, our method segmented out the first four and in some cases the fifth generations.
Figure 5: (a) White area represents automatically segmented airway, (b) demonstrates the overlap which exists with (c), the next slice. Overlap between (a) and (c) is computed by including all pixels in (c) which have gray scale statistics which match those of the segmented airway lumen in (a)
We have developed an automated method which segments the primary airway tree. The method has been successfully tested on both normal subjects and subjects exhibiting stenotic airway segments and cystic fibrosis. Although additional edge finding methods need to be incorporating, we believe that the segmented airway tree using our method can be used by pulmonologists in surgical planning and treatment of subjects suffering from various forms of airway stenosis.
Probably one of the main hurdles which needs to be overcome when testing any medical image processing methodology is access to patient data sets. Dr. Geoff McLennan who has a keen interest in exploring the feasibility of Volumetric CT as a complement to bronchoscopy allowed us a sufficient amount of data to test our method. Special thanks to Suzanne Shamsholkotabbi, Shereen Chang and Janice Cook-Granroth, for the manual segmentation of the data sets which have served as a gold standard. This work was supported in part by the National Library of Medicine (contract N01-LM-4-3511).
Papers |
DPI Homepage |
VIDA |
NLM |
Contact Us |
Search