Abstract:

Measurements of the in vivo bronchial tree can be used to assess regional airway physiology. High-resolution CT (HRCT) provides detailed images of the lungs and has been used to evaluate bronchial airway geometry. Such measurements have been used to assess diseases affecting the airways, such as asthma and cystic fibrosis, to measure airway response to external stimuli, and to evaluate the mechanics of airway collapse in sleep apnea. To routinely use CT imaging in a clinical setting to evaluate the in vivo airway tree, there is a need for an objective, automatic technique for identifying the airway tree in the CT images and measuring airway geometry parameters. Manual or semi-automatic segmentation and measurement of the airway tree from a 3-D data set may require several man-hours of work, and the manual approaches suffer from inter-observer and intra-observer variabilities.

This paper describes a method for automatic airway tree analysis that combines accurate airway wall location estimation with a technique for optimal airway border smoothing. A fuzzy logic, rule-based system is used to identify the branches of the 3-D airway tree in thin-slice HRCT images. Raycasting is combined with a model-based parameter estimation technique to identify the approximate inner and outer airway wall borders in 2-D cross-sections through the image data set. Finally, a 2-D graph search is used optimize the estimated airway wall locations and obtain accurate airway borders. We demonstrate this technique using CT images of a plexiglass tube phantom.

Keywords: intra-thoracic airways, medical imaging, quantitative image analysis, X-ray computed tomography, high-resolution CT.

INTRODUCTION

Measurements of the in vivo bronchial tree can be used to assess regional airway physiology. High-resolution CT (HRCT) provides detailed images of the lungs and has been used to evaluate bronchial airway geometry. Such measurements have been used to assess diseases affecting the airways, such as asthma and cystic fibrosis, to measure airway response to external stimuli, and to evaluate the mechanics of airway collapse in sleep apnea. Airway measurement techniques have varied from manual and semi-automatic border tracing [1, 2, 3, 4, 5], to more quantitative semi-automatic and automatic approaches based on analyzing the image intensity near the airway border [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Specialized software packages have been designed to facilitate these types of airway geometry analyses [16, 17, 18].

To routinely use CT imaging in a clinical setting to evaluate the in vivo airway tree, there is a need for an objective, automatic technique for identifying the airway tree in the CT images and measuring airway geometry parameters. Manual or semi-automatic segmentation and measurement of the airway tree from a 3-D data set may require several man-hours of work, and the manual approaches suffer from inter-observer and intra-observer variabilities.

In this paper, we describe a method for automatic airway tree analysis that combines accurate airway wall location estimation [12, 13, 25] with a technique for optimal airway border smoothing [15]. First, a fuzzy logic, rule-based system is used to identify the branches of the 3-D airway tree in thin-slice HRCT images. Raycasting is combined with a model-based parameter estimation technique to estimate the approximate inner and outer airway borders in 2-D cross-sections through the image data set. Finally, a 2-D graph search is used optimize the estimated airway wall locations and obtain accurate airway borders. We demonstrate this technique using CT images of a plexiglass tube phantom.

DESCRIPTION

The proposed technique consists of three main components: (1) a rule-based system to identify the airways in the 3-D CT image data set; (2) a measurement approach that uses a model of an airway and the scanning process to estimate the inner and outer wall locations of the airway along the approximate airway tree centerline; and (3) a cost-based graph search that is used to globally smooth the airway wall location estimates to compute the airway wall borders. In this section, each of these components is described in detail.

Airway Tree Detection

The airway tree detection combines 3-D seeded region growing with a rule-based fuzzy logic classification system [15]. The region growing step identifies large diameter, 26-connected airways (the ``primary'' airway tree). The rule-based system classifies regions of homogeneous pixels into ``airway'' or ``non-airway'' in order to detect the small diameter airways. The entire method consists of the following steps:

1. Identify the primary tree using 3-D seeded region growing with automatically selected seed points.
2. Preprocess, enhance, and segment regions on the individual 2-D image slices.
3. Compute region features for candidate airways and candidate vessels.
4. Classify candidate regions into ``airway'' and ``non-airway'' using fuzzy if-then rules.
5. Extract the 26-connected 3-D airway tree.

Primary Airway Tree Detection

The primary airway tree is identified using 3-D 26-connected region growing [14, 17, 19, 20, 21, 22]. The seed points for region growing are automatically detected using thresholding and 2-D connectivity analysis to identify the largest, dark, round regions in the uppermost slices of the data set (these regions correspond to the trachea and/or the branches of the mainstem bronchi). A conservative, heuristically determined region growing threshold is used to find pixels that are 26-connected to the seed points. This set of pixels makes up the primary airway tree. The primary tree does not contain many small diameter airways due to the conservative region growing threshold.

Image Slice Preprocessing and Segmentation

The preprocessing and segmentation steps are used to identify the approximate borders of the lung field and to segment the lung parenchyma into homogeneous regions of pixels. These regions of pixels will be later classified into ``airway'' or ``non-airway'' categories.

Preprocessing starts with threshold-based lung segmentation. Non-lung structures are removed using 2-D connectivity and gray level information. A 13 X 13 hat-transform [23] is applied to the grayscale data to enhance the contrast around small diameter airways. The lung field (parenchyma) is then segmented into homogeneous regions of pixels using edge-based region growing [14]. The resulting segmentation will consist of dark regions and bright regions corresponding to possible airway and vessel locations. These regions are identified and labeled as candidate airways and candidate vessels.

Region Feature Calculation

• Brightness - the maximum gray level in the candidate airway region.
• Degree of wall existence - an airway is surrounded by a wall that appears as a bright ring in the individual slices. The degree of wall existence is expressed as a numerical feature and serves as a region property.
• Adjacency - the gray level of the brightest region that is adjacent to the candidate airway region.

Region Classification

Each region is classified using the region feature vector and the feature vectors of neighboring regions. The classification uses a fuzzy logic classifier with the membership functions and rules automatically trained using adaptive fuzzy partitioning [15]. The classification process labels each of the regions into either ``airway'' or ``non-airway'' categories by thresholding the fuzzy labeling confidence values.

3-D Tree Extraction and Airway Centroid Computation

The 3-D airway tree is constructed by finding the set of all labeled ``airway'' regions that are 26-connected to pixels in the primary airway tree. The 3-D tree is is thinned to a linked list of airway centroids, where the centroids are computed on each 2-D slice of the data set. This list of centroids captures the tree topology and will be used to identify the airway wall locations using the model-based parameter estimation technique. More sophisticated 3-D thinning approaches are being investigated [24].

Airway Wall Location Estimation

A popular technique for airway geometry measurement on 2-D image planes uses a ``ray casting'' technique. In ray casting, a set of rays oriented at different angles are cast outward from the approximate airway centroid. The gray-level profile along each ray is examined to estimate the radial location at which the ray crosses the airway wall. Each ray gives an independent estimate of the inner (and/or outer) wall location. If both inner and outer airway walls are detected, airway wall thickness can also be determined.

D'Souza describes a system, called ASAP, that uses ray casting for interactive airway geometry analysis [18]. ASAP requires an operator to indicate the approximate airway centroid with a mouse-driven video cursor. After the centroid is determined, rays are cast from the centroid and the inner and outer walls locations are estimated. Thus, in order to analyze an entire 3-D data set, the operator must manually inspect each slice and select the airways of interest.

Rather than relying on an operator to specify the approximate airway centroids, we use the airway tree centroid representation described in Section 2.2 to estimate the location of the airway centroids. For each airway centroid, a 2-D plane perpendicular to the airway centerline is computed. Image intensity values are computed (interpolated, if needed) within this plane. A set of rays lying in the 2-D plane is cast outward from the centroid through the airway wall.

For each ray that passes through the airway wall, the set of intensity values along the ray (the gray level profile) is used to estimate the inner and outer wall locations. The wall locations are estimated using a model-based parameter estimation technique [12, 13, 25]. We model the scanning process of an ideal airway and use the model to predict the shape of the gray level profiles of rays crossing the airway wall. The model parameters (inner and outer diameters, and parenchymal density) are adjusted to minimize the difference between the modeled profile and the actual profile observed in the data. The set of parameters that minimizes the difference between the model and the data yields an estimate of the airway geometry. For each of the estimated airway geometry parameters, an indication of measurement confidence (standard deviation) is computed. This confidence level is used in the next step to smooth the estimated airway wall locations.

Applying Wall Border Shape Constraints

Each ray that is cast from the the airway tree centroid provides an independent estimate of airway inner and outer wall location. Because each ray is assessed independently, there can be considerable ray-to-ray variation and measurement noise in the estimated wall locations--and there are no local shape continuity constraints applied during the model-matching process that will ensure the identification of a sensible airway border. To reduce ray-to-ray variation and enforce a smooth airway border, a contour smoothing step using a cost-based graph search [23] is applied to the raw wall location estimates.

For both the inner and outer walls, the model-based wall estimation method returns the estimated wall location, , and an indication of measurement accuracy, . A graph is constructed with the rows of the graph corresponding to the ray profiles measured at different angles. The node cost is computed using , and so the cost is proportional to the likelihood that a wall exists at that node location. A graph node (k, i) corresponds a radial location along ray k, where varies from 0 to R radial units (the ray is R units in length) and the ray is oriented at an angle of , where N is the total number of rays cast. The node cost function c(k,i) associated with each graph node (k,i) is defined as follows:

where and are the estimated wall location and measurement accuracy for the ray. The graph search is applied to find the minimum cost path through the graph from bottom (ray 0) to top (ray N-1). Figure 1 illustrates the use of the A-algorithm[26] to guide the graph search.

  
Figure 1: Example of a graph searching sequence using the A-algorithm[26]. Nodal costs are inversely proportional to the likelihood that the wall lies at the location of a given node. Individual rows of the graph correspond to individual ray profiles oriented at different angles in the image data. (a) Step 1; expansion of the start node; (b-e) Steps show the propagation of total path cost from one end of the graph to the other; (f) the optimal path is defined by backtracking. Figure adapted from Sonka, Hlavac, and Boyle [27].

After graph search has identified the minimum cost border, we apply a final local smoothing step to ensure that the border is closed (i.e., the estimated wall at the top of the graph meets the estimated wall at the bottom of the graph). The smoothed border points are then mapped back into the original image space to represent the airway wall boundary.

RESULTS

For our preliminary studies, the method has been quantitatively tested using CT scans of physical phantoms. The phantom consists of several ``off-the-shelf'' plexiglass tubes ranging in size from about 6.5 mm to 19.25 mm inside diameter, with tube wall thicknesses of approximately 3 mm. The measurement ``gold standard'' for the tube dimensions was obtained using micrometer measurements of the inner and outer tube diameters.

The phantom was scanned using an Imatron electron-beam CT (EBCT) scanner, with in-plane pixel size of 0.293 mm by 0.293 mm, and with a slice thickness of 3 mm. Scan aperture was 100 msec, and reconstruction field-of-view was 15 cm. Reconstruction used the normal kernel.

Our airway analysis method was applied to a stack of 10 2-D CT image slices. The method segmented the tubes (airways), identified approximate tube centerlines, and found the inner airway wall border using raycasting with the model-based wall location estimation followed by the cost-driven graph search. For this application, we did not attempt to detect the location of the outer tube wall. Table 1 compares the automatic method with semi-automatic measurements made using ASAP [18]. For the semi-automatic method the operator interactively specified airway centroids on each slice of the data and used ray-casting with the model-based wall location estimation to find the inner airway wall. This semi-automatic technique was previously validated on phantoms [18], and recently on actual lung tissue [25]. dD..-2

  

Title Number	True Inner Radius (mm)	Seni-Auto. Inner Radius Error (mm)	Seni-Auto. Error Std. Dev. (mm)	Automatic Inner Radius Error (mm)	Automatic Error Std. Dev. (mm)
1	9.63	-0.08	0.71	-0.15	0.26
2	4.75	-0.01	0.85	-0.18	0.14
3	3.25	0.02	0.64	-0.01	0.16

Table 1: Comparison of wall estimation accuracy for the semi-automatic approach and automatic method. Semi-automatic method uses operator selected centroid and model-based wall estimation. Automatic method uses automatic segmentation, centerline/centroid detection, model-based wall estimation, and cost-based graph search. Results are averaged over 10 slices of CT data set. A negative estimation error indicates the inner radius is underestimated with respect to the gold standard.

  
Figure 2: Results of applying automatic airway border detection to image of plexiglass tube phantom. Figure shows estimated airway border location from the automatic technique. Tube one is upper left, tube two is middle right, and tube three is bottom left.

DISCUSSION

The results in Table 1 show that the mean estimation error is slightly increased using the automatic technique compared to the semi-automatic analysis using ASAP. However, the increase in estimation error is small--on the order of 0.1 mm (about 1/3 of a pixel)--and the automatic method has a greatly reduced measurement error standard deviation compared to the semi-automatic technique. This is because the automatic method employs a graph search step to enforce local shape constraints that act to smooth the airway boundary. For the three tube sizes discussed in Table 1, the error standard deviation using the automatic method is three to four times smaller than the semi-automatic approach. As a result, the airway boundaries identified using the automatic approach are smoother and less noisy than the borders identified using the model-based technique alone.

One explanation for the increased measurement error using the automatic method is that while the model-based technique provides wall location estimates with sub-pixel accuracy, the graph searching step rounds all values to the nearest pixel when computing the optimal boundary. Thus, no matter how accurately the model-based approach can estimate the wall location on an individual ray, the final airway border will be forced to pass through the nearest image pixel. This problem can be addressed by oversampling the rays when constructing the graph. In the oversampled graph, a ray of length R units will be mapped to a graph row of length M R nodes, where M determines the amount of oversampling(M>=1 ). The graph search can be applied to oversampled graph, and the resulting optimal path around the airway border can be mapped back to the image data with sub-pixel accuracy.

SUMMARY

Quantitative analysis of the in vivo bronchial tree can be used to non-invasively assess regional airway physiology. Automatic measurement of airway geometry parameters from HRCT images can provide fast, accurate indices of lung pathology, and is free from the inter-observer and intra-observer variabilities present with manual analyses. We have described a technique for automatic 3-D airway tree analysis. This technique is the integration of several well-tested algorithms for image analysis: (1) fuzzy rule-based airway tree segmentation; (2) airway centroid/centerline detection; (3) airway wall location estimation; (4) border smoothing using a cost-driven graph search.

We have presented results demonstrating our method on CT scans of a plexiglass tube phantom. We compare the new automatic technique against a previously validated semi-automatic method [18, 12, 13, 25]. The results show that the automatic technique has a slightly increased mean measurement error (about 1/3 pixel), but greatly reduced measurement noise (three to four times reduction in measurement error standard deviation). The reduced measurement noise translates to smoother, more realistic airway border estimates. Future work will focus on implementing the oversampled graph to retain sub-pixel precision during the airway border estimation step. Final validation must be done using a realistic branching airway phantom.

ACKNOWLEDGEMENTS

This project was supported in part by the National Library of Medicine (contract N01-LM-4-3511) and the National Science Foundation (IRI 96-16747).

References

1

W. R. Webb, G. Gamsu, C. E. Cann, and E. Proctor, ``CT of a bronchial phantom: Factors affecting appearance and size measurements,'' Invest. Radiol., vol. 19, pp. 394-398, Sept.-Oct. 1984.

2

R. H. Brown, C. J. Herold, C. A. Hirshman, E. A. Zerhouni, and W. Mitzner, ``In vivo measurements of airway reactivity using high-resolution computed tomography,'' Am. Rev. Resp. Dis., vol. 144, no. 1, pp. 208-212, 1991.

3

C. J. Herold, R. H. Brown, W. Mitzner, J. M. Links, C. A. Hirshman, and E. A. Zerhouni, ``Assessment of pulmonary airway reactivity with high-resolution CT,'' Radiology, vol. 181, pp. 369-374, Nov 1991.

4

A. E. McNamara, N. L. Müller, M. Okazawa, J. Arntorp, B. R. Wiggs, and P. D. Paré, ``Airway narrowing in excised canine lungs measured by high-resolution computed tomography,'' J. Applied Physiology, vol. 73, no. 1, pp. 307-316, 1992.

5

E. Senéterre, F. Paganin, J. M. Bruel, F. B. Michel, and J. Bousqet, ``Measurement of the internal size of bronchi using high resolution computed tomography (HRCT),'' European Respiratory J., vol. 7, pp. 596-600, 1994.

6

B. S. Baxter and J. A. Sorenson, ``Factors affecting the measurement of size and CT number in computed tomography,'' Invest. Radiol., vol. 4, pp. 337-341, 1981.

7

J. H. C. Reiber, C. J. Kooijman, C. J. Slager, J. J. Gerbrands, J. H. C. Schuurbiers, A. den Boer, W. Wijns, P. W. Serruys, and P. G. Hugenholtz, ``Coronary artery dimensions from cineangiograms--methodology and validation of a computer-assisted analysis procedure,'' IEEE Trans. Medical Imaging, vol. 3, pp. 131-141, Sept. 1984.

8

S. R. Fleagle, M. R. Johnson, C. J. Wilbricht, D. J. Skorton, R. F. Wilson, C. W. White, M. L. Marcus, and S. M. Collins, ``Automated analysis of coronary arterial morphology in cineangiograms: Geometric and physiologic validation in humans,'' IEEE Trans. Medical Imaging, vol. 8, pp. 387-400, Dec. 1989.

9

R. J. Schwab, W. B. Gefter, L. R. Kline, A. I. Pack, and E. A. Hoffman, ``Dynamic imaging of the upper airway during respiration in normal subjects,'' J. Applied Physiology, vol. 74, pp. 1504-1514, 1993.

10

I. Amirav, S. S. Kramer, M. M. Grunstein, and E. A. Hoffman, ``Assessment of methacholine-induced airway constriction by ultrafast high-resolution computed tomography,'' J. Applied Physiology, vol. 75, no. 5, pp. 2239-2250, 1993.

11

S. S. Kramer, E. A. Hoffman, and I. Amirav, ``High resolution CT assessment of pediatric airways: Structure and function,'' in Proc. SPIE Conf. Medical Imaging, vol. 2168, (Newport Beach, CA), 13-18 Feb. 1994.

12

J. M. Reinhardt, N. D. D'Souza, and E. A. Hoffman, ``Improved quantitation of airway geometry,'' Ann. Biomed. Eng., vol. 24, pp. S-65, Oct. 1996.

13

J. M. Reinhardt, N. D. D'Souza, and E. A. Hoffman, ``Accurate measurement of intra-thoracic airways.'' Submitted to IEEE Trans. Medical Imaging.

14

M. Sonka, W. Park, and E. A. Hoffman, ``Rule-based detection of intrathoracic airway trees,'' IEEE Trans. Medical Imaging, vol. 15, no. 3, pp. 314-326, 1996.

15

W. Park, E. A. Hoffman, and M. Sonka, ``Segmentation of intrathoracic airway trees: A fuzzy logic approach.'' Submitted to IEEE Trans. Medical Imaging.

16

R. A. Robb, ``A software system for interactive and quantitative analysis of biomedical images,'' in 3D Imaging in Medicine (K.-H. Höhne, H. Fuchs, and S. M. Pizer, eds.), vol. F60 of NATO ASI Series, pp. 333-361, New York: Springer-Verlag, 1990.

17

E. A. Hoffman, D. Gnanaprakasam, K. B. Gupta, J. D. Hoford, S. D. Kugelmass, and R. S. Kulawiec, ``VIDA: An environment for multidimensional image display and analysis,'' in SPIE Conf. Biomedical Image Processing and Three-Dimensional Microscopy, vol. 1660, pp. 1-18, 1992.

18

N. D. D'Souza, J. M. Reinhardt, and E. A. Hoffman, ``ASAP: Interactive quantification of 2D airway geometry,'' in Proc. SPIE Conf. Medical Imaging, vol. 2709, (Newport Beach, CA), pp. 180-196, 10-15 Feb. 1996.

19

J. K. Udupa, ``Interactive segmentation and boundary surface formation for 3-D digital images,'' Comp. Graphics and Image Proc., vol. 18, pp. 213-235, 1982.

20

E. A. Hoffman, L. J. Sinak, R. A. Robb, and E. L. Ritman, ``Noninvasive quantitative imaging of shape and volume of lungs,'' J. Applied Physiology, vol. 54, no. 5, pp. 1414-1421, 1983.

21

S. Wood, J. Hoford, E. Zerhouni, E. Hoffman, and W. Mitzner, ``Quantitative 3-D reconstruction of airway and pulmonary vascular trees using HRCT,'' in Biomedical Image Processing and Biomedical Visualization, vol. 1905, (Bellingham, WA), pp. 316-323, SPIE, 1993.

22

R. A. Chiplunkar, J. M. Reinhardt, and E. A. Hoffman, ``Segmentation and quantification of the primary human airway tree from 3-D X-ray CT,'' in Proc. SPIE Conf. Medical Imaging, vol. 3033, (Newport Beach, CA), 23-28 Feb. 1997.

23

M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis and Machine Vision. London: Chapman & Hall, 1  ed., 1993.

24

C. Min Ma and M. Sonka, ``A fully parallel 3D thinning algorithm and its applications,'' Comp. Vis. and Image Understanding, vol. 64, pp. 420-433, 1996.

25

J. M. Reinhardt, S. A. Raab, N. D. D'Souza, and E. A. Hoffman, ``Intra-thoracic airway measurement: Ex vivo validation,'' in Proc. SPIE Conf. Medical Imaging, vol. 3033, (Newport Beach, CA), 23-28 Feb. 1997.

26

N. J. Nilsson, Problem-Solving Methods in Artificial Intelligence. New York: McGraw-Hill, 1971.

27

M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis and Machine Vision. London: Chapman & Hall, 2  ed., 1998. (In press.).