Published in: Medical Imaging 1996: Physiology and Function from Multidimensional Images, Eric A. Hoffman, Editor, Proc. SPIE 2709, 180-196 (1996).

ASAP: Interactive quantification of 2D airway geometry

Neil D. D'Souza, Joseph M. Reinhardt, and Eric A. Hoffman

Division of Physiologic Imaging
Department of Radiology
University of Iowa College of Medicine
Iowa City, IA 52242

Table of Contents

Abstract

Introduction

Methods

Systems Description

Experimental Results

Summary

Acknowledgements

References

Abstract:

Keywords: cardio-pulmonary imaging, airways, medical imaging, interactive image segmentation, quantitative image analysis, multi-dimensional image processing, X-ray computed tomography, high-resolution CT.

1 INTRODUCTION

Peripheral bronchial airway geometry is a valuable tool for assessing pathology related to pulmonary function. Airway geometry measurements can be used to evaluate and track the progression of disease directed at both the airway and lung parenchyma such as cystic fibrosis, asthma, and emphysema. High-resolution X-Ray CT (HRCT) can provide detailed volumetric images of the lungs and bronchial tree. It is recognized that to evaluate the true geometry of the bronchial tree, one must measure the airway in planes perpendicular to the local long axis of the airway. The need to acquire three-dimensional image data sets to evaluate the heart and lungs was addressed through the development of the Dynamic Spatial Reconstructor (DSR) [1], which has remained a one-of-a-kind machine and is actually the only X-ray CT scanner capable of acquiring true dynamic volumetric images of the lung. Block and colleagues [2] used the DSR to demonstrate that because of partial volume effects, algorithms for the assessment of cross-sectional area of tubular structures (they used a glass branching tube model) are not applicable across all size tubes and that one must customize the measurement based upon an initial assessment of the area. Considerable efforts are currently underway to develop methods for acquiring volumetric image data sets of the lung using commercially available scanners with spiral (continuous rotation based upon slip ring technology) or electron beam CT (EBCT) [3] technology. Other work is underway to extract and visualize the three-dimensional airway tree [4,5,6,7]. While significant advances have been made in the research laboratory, the techniques currently available for image acquisition and image analysis are not yet ready for the rigors of a clinical environment. As a stop-gap measure, significant attention has been paid to the evaluation of intrathoracic airways using two-dimensional sections, while limiting the evaluation to those airway sections which are nearly round and thus have likely been sectioned, fortuitously, perpendicular to their long axis. Much of the analysis of the image data sets has been approached using non-digital, highly labor-intensive approaches. Some studies evaluating HRCT images have relied upon manual application of calipers placed on printed film. Other studies go through a process of filming the 11-bit digital CT data with a predetermined window and level (to either maintain a range of graylevels or to convert the image into a binary data set). The film is then either projected onto paper and manually traced, or digitized (note that the images were originally digital) and evaluated by manually identifying airway borders [8,9,10,11]. Other approaches have relied directly on display window and level to identify airway walls [11], an approach that is neither reliable nor repeatable.

Specialized software has been developed to automate region-of-interest (ROI) analysis on the original 11-bit CT data sets and to facilitate objective region measurement [12,13,14]. Fully digital image analysis, starting from the direct transfer of the CT sections to a stand-alone image analysis workstation, has, to date, relied upon the user providing at least an initial trace at the approximate airway border [15,16,14,17]. These results have still been subject to inter-observer and intra-observer variabilities, despite the fact that a final computer-controlled trace adjustment is made based upon the half-max principle.

Most recent analysis approaches have been general-purpose and not tailored specifically for airway analysis. The approach used by Amirav et al. [14] was actually developed to evaluate the upper airway [13]---a task that is considerably different than the evaluation of the intrathoracic airways. To give added objectivity to their results, Amirav et al. [14] averaged ten manual tracings of the airway borders. The method gave a good level of accuracy. However, the approach was still too cumbersome for routine use in a clinical setting.

This paper describes a semi-automatic, minimally labor-intensive, objective approach for airway segmentation and quantification using 2-D slices of CT images. In our approach, the user interacts with a GUI-based software system called ASAP (Airway Segmentation and Analysis Program). Using a pointing device, the user identifies the approximate locations of the centroids of the airways of interest on 2-D slices of the volume. ASAP then automatically detects airway walls and measures a number of parameters, including inner and outer airway diameters, lumen cross-sectional area, and airway wall thickness, along with statistical information related to variation of the measurements around the airway circumference. ASAP will assign unique labels to the segmented airways, which can be used to build a topological description of the 3D bronchial tree. In addition, ASAP has the ability to load multiple data sets simultaneously, so that airway geometry can be compared before and after intervention.

Particular attention has been paid to the user-interface design and extensibility of ASAP to allow for ease of use and to provide the important flexibility needed in a clinical research environment. Inexperienced users can measure airways simply by pointing and clicking with the mouse, while more experienced users can, if necessary, adjust measurement parameters and program layout through a series of configuration menus. ASAP can exchange data with a broad-based image analysis package, VIDA [12]. Using VIDA, the user can visualize the segmented airway tree in 3D, or apply additional specialized processing to the segmented data from ASAP.

This paper describes the overall analysis strategy and design of the ASAP program. We demonstrate the effectiveness of the semi-automatic approach using two specific test cases: (1) measurements made using ASAP on the CT scan of a known physical airway phantom; and (2) measurements made using ASAP on CT scan of a human subject. Our results show that the measurements obtained using our approach closely match the actual physical dimensions of the airway phantom, and the measurements made using ASAP have little inter-observer and intra-observer variability. Further, the semi-automatic method is much faster and easier than both manual analysis and our previous general-purpose approach.

2 METHODS

ASAP is designed to facilitate airway identification, airway registration between multiple data sets, and airway measurement in 2D data sets. ASAP can make geometric measurements such as inner and outer wall diameter, wall thickness, and lumen cross-sectional area. The methods used to make geometric measurements are tailored toward airways whose long axis is nearly perpendicular to the imaging plane, i.e., airways that appear circular (or nearly circular) on 2D slices of the data. If needed, measurements can be made to airways whose direction of travel is not perpendicular to the imaging plane, but an appropriate correction must be applied to account for the off-axis measurement.

Figure 1: Typical human airways viewed on 2D transverse slice of 3D volumetric CT data set. Adjacent to the airway at the center of the image is a blood vessel (bright, solid, circular blob).

As shown in Figure 1, airways traveling perpendicular to the imaging plane appear as ``ring-like'' structures on transverse slices of the CT data set. Blood vessels in the lung, which typically travel parallel to the bronchial passages, appear as bright, solid, circular objects adjacent to airways. The airway lumen is normally dark, although in small airways the partial-volume effect increases the average graylevel within the lumen. The lumen is surrounded by the solid airway wall.

In our approach, the user interacts with a graphical-user interface to identify a pixel within the lumen of an airway of interest. The system then automatically identifies the airway walls using one of two selectable approaches: (1) ray casting, which works best for circular and nearly-circular airways; and (2) threshold-based region growing, which can be applied to both circular and non-circular airways. Each of these methods is described individually in the next sections.

2.1 Ray casting

The ray-casting method can be used to estimate inner and outer diameter and lumen cross-sectional area. The ray-casting method sends rays outward from the centroid of the airway and determines the inner and outer airway walls by examining the individual graylevel profile of each ray. The method can be briefly summarized as follows: (1) user selects an airway by clicking at the approximate airway centroid; (2) a preliminary set of rays is cast from the user-specified centroid to get an approximate estimate of the inner and outer diameter; (3) the centroid is refined based on the estimated location of the inner wall; (4) wall detection thresholds are computed based on the estimate of the inner and outer walls; (5) a new set of rays is cast from the refined centroid; (6) final estimates of the inner and outer airway wall locations are determined by examining the graylevel profile along each ray. Next, the processing steps are described in detail.

Given an approximate user-specified airway centroid, the first step in the processing is to refine the centroid based on a rough estimate of inner airway wall locations. To refine the centroid location, N rays are cast from a user-specified centroid, with an angular spacing of radians between each ray. The graylevel profile along each ray is examined to estimate the location of the inner and outer wall.

Let (x_u1y_u) be the location of the user-specified centroid, which is assumed to lie within the airway lumen. Each ray is of length R units. The j^th ray, j = 1, . . ., N, is oriented at an angle with respect to the x axis and connects the points (x_u, y_u) and (see Figure 2).

  
Figure 2: Diagram illustrating how ray casting is used to determine the approximate inner and outer wall locations. Left-hand figure depicts a typical airway viewed on a 2-D slice. The local long axis of the airway is approximately perpendicular to the imaging plane --- as a result the airway lumen is approximately circular. Adjacent to the airway is a blood vessel (solid bulge at top left). N rays are cast from the approximate centroid of the airway (here N=8, with the angular ray spacing equal to radians). Right-hand figure shows the possible grayscale intensity profile along one ray. As discussed in the text, the maximum and minimum values along the ray are used to estimate the inner and outer wall locations.

Let f_j(r) be the graylevel along the j^th ray at a distance of r units from the centroid (at the point ). The minimum and maximum values of f_j(r), 0rR,are used to estimate the inner and outer wall locations. Let f_j(r_max) be the maximum intensity along the ray, f_j(r_min1) be the minimum intensity along the ray for 0r_min1r_max, and f_j(r_min2) be the minimum intensity along the ray for r_maxr_min1R. Thus, as illustrated in Figure 2, r = r_max is the location of the maximum along the entire ray, r_min1 is the location of the minimum between r=0 and r = r_max, and r_min2 is the location of the minimum between r = r_max and r=R.

Next, we compute the expected graylevel at the locations of the inner and outer wall. The expected graylevels at the wall are based on the maximum and minimum graylevels observed along the ray, and the shape of the ray profile. Let rⁱ_j be the inner wall radius and be the outer wall radius for the j^th ray, rⁱ_j and are computed as follows:

where t_i and are thresholds used to determine the graylevel at the inner and outer wall locations, and g(^.) = f^-1(^.), i.e., g(^.) is the mapping from graylevel to radial position along the ray. Per (1) and (2), the expected graylevel value is estimated to be t_i(f(r_max)-f(r_min1))+f(r_min1) at the inner wall and at the outer wall. During the centroid refinement step, we use . Thus, for this initial wall location estimate, we use the half-max criteria to select the edge. Note that g(^. is not explicitly computed, however, g(^.is determined by evaluating the graylevel profile along the ray. A linear search along the ray is used to detect the locations of rⁱ_j and . Interpolation is used to obtain a good estimate of the wall positions.

The next step is to refine the user-specified centroid based upon the initial wall location estimates. To refine the centroid estimate, we consider the polygon formed using the inner wall locations rⁱ_j, j = 1, . . ., N, as vertices. The refined centroid (x_c,y_c) is estimated as the centroid of the polygon:

As discussed earlier, the first estimates of the wall locations using (1) and (2) were computed using thresholds values of t_i and equal to 0.5, i.e., the wall was estimated using the half-max criteria. In practice, this approach can lead to poor estimates of wall position, especially for small airways. To address this problem, we have developed a method of adapting the thresholds t_i and based on the initial estimates of the inner and outer wall locations. For the results presented in this paper, the method uses the following empirically-derived thresholds (see Section 4.1 regarding phantom studies):

These thresholds were determined by scanning a phantom containing plexiglass tubes of known sizes (see Section 4.1 for details). The phantom scans were used to determine the threshold used to detect the inner wall, t_i, as a function of actual inner diameter, and the threshold used to detect the outer wall, , as a function of actual outer diameter. Since this empirical study was completed, we have performed a theoretical analysis showing that (4) is only a rough approximation to the ideal threshold selection strategy, however, (4) seems to gives good results over a wide range of common airway sizes.

After the refined centroid (x_c,y_c) and appropriate thresholds t_i and have been determined, a new set of N rays is cast outward from the refined centroid location. For each ray, the wall locations are computed using (1) and (2), however, using the new values of t_i and above. Again, we use a linear search along the ray to map grayscale intensity to radial position and to compute the locations of rⁱ_j and .

After the inner and outer wall locations are determined for each ray, three criteria are used to eliminate measurements that are likely to be in error. Erroneous measurements may be obtained if the ray crosses a wall that is especially weak, or passes through a neighboring structure (e.g., a blood vessel, as shown in Figure 3) leading to a bad estimate of the maximum along the ray. The measurements made using a ray j are accepted if they meet the following criteria:

1. wall thickness test: . Default values for t₁ and t₂ are 0.10 and 0.35. Thus, by default, ray measurements are accepted only if the wall thickness is between 10% and 35% of the inner radius.

2. median inner radius test: t₃rⁱ_medianrⁱ_jt₄rⁱ_median, where rⁱ_median is the median inner radius computed over all N rays. Default values for t₃ and t₄ are 0.85 and 1.15. Using these values, the measured radius must be within 15% of the median radius to be accepted.

3. mean inner radius test: rⁱ_mean - k rⁱ_j rⁱ_mean + k, where rⁱ_mean and are the mean inner radius and inner radius standard deviation computed over all N rays. The default value for k is 2. Thus, by default, the measured radius must be within 2 of the mean radius to be accepted.

Figure 3: Example illustrating the results of the ray-casting method. Figure shows rays used to detect inner and outer wall locations, plus splined wall boundary obtained using estimated wall locations as vertices. Light color rays have passed the ray acceptance criteria, the four dark rays have failed the ray acceptance tests and do not contribute to the measurements.

2.2 Threshold-based region growing

Threshold-based region growing uses a user-specified point inside the airway to identify the airway lumen and inner wall. This method can work well for both circular and non-circular airways, so it can be used to extract airways not traveling perpendicular to the imaging plane, or near airway bifurcation points. The average airway inner diameter is estimated from the lumen area.

Region-growing is used to identify the lumen by finding all pixels that are less than a specified graylevel threshold, and that are 8-connected to the seed point [18]. The inner wall of the airway is identified as the outer border of the lumen. Because the graylevel threshold is a critical parameter that directly affects the ability to obtain accurate and repeatable measurements, we have developed a method of automatically selecting an appropriate threshold. Our region-growing threshold is selected by using the ray-casting method to find the approximate location of the inner wall, and then using the minimum graylevel we find along that wall location. As described in 2.1, we cast rays from the user-defined centroid and estimate the inner wall locations based on the minimum and maximum graylevels observed along a ray. To select a threshold, we estimate the inner wall location by casting 4 rays from the centroid and apply (1) with (half-max criteria). We then examine the graylevel at each detected ray location and choose the minimum value for the region-growing threshold. If necessary, the region-growing threshold can also be directly specified by the user.

3 SYSTEM DESCRIPTION

ASAP is an event-driven GUI system designed to reduce the amount of user interaction required to accurately measure airway geometry. In the best case, the user is required to click only once in the lumen of each airway to be measured. As described above, the system uses two measurement methods: (1) ray casting; and (2) threshold-based region growing. The analysis process consists of the user loading a dataset from disk or shared memory, selecting a measurement method, then zooming and panning through the slices selecting airways of interest with the mouse, and finally generating a report of the measured parameters.

3.1 User interface design

ASAP is written in C using the OpenLook/Xview libraries and is intended for use in conjunction with VIDA [12], a broad-based system for quantitative image analysis. The program runs on multiple platforms, including Sun Workstations running Solaris 2.x, HP 700 series and Silicon Graphics machines.

  
Figure 4: ASAP main control window. This window allows the user to load and save data sets and reports, select measurement strategy, modify method parameters, and zoom and pan through the data set.

  
Figure 5: ASAP data slice display window. Figure shows the user applying the ray-casting method. This window allows the user to zoom, pan, and step through the data slice-by-slice.

  
Figure 6: ASAP measurement statistics window. Statistics shown in figure are slice number (SN), region number (RN), lumen cross-sectional area for threshold-based region growing (LA(T)) and ray-casting methods (LA(R)), inner diameter (ID), outer diameter (OD), and wall thickness (WT).

The User Interface Design of the ASAP software system consists of a main control window (Figure 4), a slice display window (Figure 5), a slice display layout editor, a statistics window (Figure 6), a ray profile window and windows for Input/Output of data.

3.1.1 Main control window

This window contains several subpanels for display of independent sets of variables used in the program, and a set of menus to control these subpanels. The subpanels are contained within the main control panel to provide better organization of windows, thereby reducing the complexity of the software resulting from several windows being displayed on the screen.

3.1.2 Slice display layout editor

The editor is used to customize the layout of images in the slice display window. Depending on the number of datasets being used (1 or 2), this editor window gives the capability of displaying up to six images (at most 3 images for each dataset). For each dataset, it is possible to display the current, previous and next slice. The option for viewing the previous and next slice is provided to enable the user to follow airways in a dataset. An additional feature of this editor is the ability to save different types of layouts depending on the type of study or analysis, into a resource file which the program loads at startup time.

3.1.3 Slice display window

The window displays up to six slices simultaneously (depending on the current layout selected in the slice display layout editor). Mouse interaction with the window, namely zooming (right mouse button) and panning (middle mouse button) of the images and selecting (left mouse button) and deleting (shift key plus mouse button) airways, occurs only in the canvas displaying the current slice of a dataset.

3.1.4 Statistics window

The window is used to report measurements and statistics, including slice number, region number, lumen area, mean wall thickness, mean inner diameter, mean outer diameter and centroid. The contents of this window can be sent to a printer or saved to a file.

3.1.5 Ray profile window

This window is only used when applying the ray-casting method. The window displays the intensity profile and wall locations for any selected ray. It is possible to accept or reject any ray in this window and override the ray rejection criteria described in Section 2.1. This window is shown in Figure 7.

  
Figure 7: ASAP Ray Profile Window. Here the user can examine the ray profile for a specific ray, or manually accept or reject rays.

3.1.6 Input/Output windows

These windows are used to read datasets from disk or shared memory, load and save regions created using the segmentation methods in the ROI (regions of interest) file format and, load and save analysis parameters for airways in an ASAP file format.

3.1.7 Zoom-Pan panel

This panel can be opened from the ``Properties'' menu in the main control window. It contains sliders for panning and zooming the images displayed in the slice display window.

3.1.8 Ray acceptance panel

This panel, which can be selected from the ``View'' menu, contains items associated with the ray acceptance criteria described in the ray casting method. Using this panel users can control the ray acceptance criteria described in Section 2.1.

3.1.9 Region labeling panel

The program by default uses sequential numeric labels for the airways within a slice. However, the airway labels can be modified using this panel, which can be opened from the ``Properties'' menu in the main control window.

3.1.10 Statistics options panel

This panel, which can be selected from the ``View'' menu in the main control window, contains a list of all the possible statistics that can be computed by the program.

3.1.11 Miscellaneous properties panel

This panel, which can be opened from the "Properties" menu in the main control window, contains a set of user selectable flags for modifying the segmentation process and customizing the display of the results

4 EXPERIMENTAL RESULTS

The methods and software system described earlier were tested by scanning a known plexiglass phantom and by analyzing data from two normal subjects and data from one patient presenting with mild cystic fibrosis (CF). ASAP was compared against the known physical dimensions of the phantom, and to concommitent measurements made using a program designed to facilitate general region-of-interest (ROI) analysis (the ROI program from VIDA [12]). The data from the human studies was used to test inter-observer and intra-observer variabilities.

4.1 Phantom tests

A plexiglass phantom was constructed to model the airways in the human lung. The phantom consists of seven plexiglass tubes, each with different inner and outer diameters. The tube dimensions are given in Table 1. The airspace around the tubes was filled with potato flakes to simulate the density of the lung held at functional residual capacity (approximately 60--70% air).

  
Table 1: Plexiglass phantom tube sizes (mm)

Figure 8: CT scan of plexiglass phantom used to validate airway geometry measurements. Phantom consists of seven plexiglass tubes, encased in a large plexiglass cylinder and surrounded by potato flakes. Tube dimensions are given in Table 1. Scanner parameters are given in the text.

The phantom was scanned at ten independent locations using an Imatron electron-beam CT (EBCT) scanner, with in-plane resolution 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. Five of the scans were used to calibrate the method utilized to select the wall detection thresholds t_i and (see discussion in Section 2.1). The remaining five scans were used to test the accuracy of the program. Figure 8 shows a single slice of the phantom scan.

The phantom data was analyzed using ASAP with the ray-casting method. The results obtained using ASAP were compared against the actual known physical dimensions of the plexiglass tubes in the phantom, and against measurements obtained manually using the Region of Interest (ROI) program in the VIDA image analysis system [12]. Figure 9 shows the comparison for the measurements of inner diameter, outer diameter, and lumen area. Each graph shows the mean of five measurements made on each tube in the phantom using ASAP and the ROI program. The mean measurement is plotted next to the actual physical dimensions of the phantom tube. The graphs use error bars to show the standard deviation of the five measurements made on each tube, although the standard deviation is so small compared to the actual tube dimensions that the error bars are difficult to see. Table 2 summarizes the mean measurement error for each of tubes in the phantom. The results show that, in general, using ASAP to analyze the data gives more accurate results than those obtained using manual tracing with the ROI program, although ROI gives measurements that are quite accurate relative to the size of the tubes being measured (the mean absolute measurement error is less using ASAP than using ROI, with p < 0.0025. p < 0.0005, p < 0.0005 for lumen area, inside diameter, and outside diameter). For ASAP, the mean absolute measurement errors (averaged across all measurements) for inner and outer diameter were 0.13 mm and 0.08 mm, or less than one-half of a pixel. For ROI, the mean absolute measurement errors (averaged across all measurements) for inner and outer diameter were 0.50 mm and 0.39 mm, or approximately one to one and one-half pixels. The graphs show that the variability of the measurements (as measured by the standard deviation of five independent measurements) is less using ASAP than when using ROI. ASAP requires significantly less analysis time since many of the operations and parameter selections are performed automatically.

  
Table 2: Plexiglass phantom: mean tube measurement error (mm)

It is interesting to note that the largest measurement errors for ASAP occur on phantom tube 3, which has an inside diameter of 6.5 mm and an outside diameter of 12.6 mm, corresponding to a wall thickness of about 3.1 mm. The method we use to select the wall thresholds t_i and was approximated to simplify the software implementation. Analysis shows that the approximation is least accurate for tube geometries like those in phantom tube 3. As a result, ASAP is selecting suboptimal wall detection thresholds for phantom tube 3, which results in slightly increased inner and outer diameter measurement errors. An improved implementation will use a more sophisticated look-up table to select thresholds so that the method can adapt to accurately measure a wider variety of airway geometries. If we eliminate tube 3 measurements from consideration, the maximum absolute error in measuring inner and outer diameter is reduced to about 0.2 mm (two-thirds of a pixel).

  
Figure 9: Mean error and standard deviation of inner diameter, outer diameter, and lumen cross-sectional area measurement obtained using ASAP and VIDA ROI [12], compared to actual physical dimensions of phantom.

4.2 Human data

Three human data sets were used to test the inter-observer and intra-observer variability of ASAP. The data consisted of two CT scans of normal human subjects and one CT scan of a subject with cystic fibrosis (CF). All subjects were scanned for other purposes---we simply utilized the data here for testing our analysis system. The normal subjects were scanned with the EBCT scanner, with in-plane resolution of 0.508 mm by 0.508 mm, and with a slice thickness of 3 mm. Scan aperture was 100 msec, and the reconstruction field-of-view was 26 cm. Reconstruction was done using a very sharp (high-spatial frequency) reconstruction kernel. The CF patient was also scanned on the EBCT scanner, with in-plane resolution of approximately 0.684 mm by 0.684 mm, and with a slice thickness of 3 mm. Scan aperture was 100 msec, and reconstruction field-of-view was 35 cm. Reconstruction used a very sharp reconstruction kernel.

4.2.1 Inter-observer variability

Two observers independently analyzed the two normal data sets and the CF data set using ASAP. The results show little inter-observer variability for any of the data sets. For each of the three measurements (inner diameter, outer diameter, and lumen cross-sectional area) the correlation coefficient, r, was greater than 0.98. The mean absolute difference (averaged over all observations) between the two observers was 0.18 mm and 0.25 mm for the inner and outer diameter measurements (less than one pixel). Figures 10 and 11 compare the measurements made by each observer.

  
Figure 10: Comparison of inner diameter, outer diameter, and lumen cross-sectional area measurements obtained using ASAP for Observer 1 and Observer 2; Normal data.

  
Figure 11: Comparison of inner diameter, outer diameter, and lumen cross-sectional area measurements obtained using ASAP for Observer 1 and Observer 2; CF data.

4.2.2 Intra-observer variability

A single observer analyzed the same normal data set using ASAP on two separate occasions. The results show little intra-observer variability. For each of the three measurements, the correlation coefficient, r, was greater than 0.98. The mean absolute difference (averaged over all observations) between the two analyses was 0.15 mm and 0.23 mm for the inner and outer diameter measurements (less than one pixel). Figure 12 compares the results of the these analyses.

  
Figure 12: Comparison of inner diameter, outer diameter, and lumen cross-sectional area measurements obtained using ASAP for Observer 1 for two distinct analyses; Normal data.

5 SUMMARY

We have described a system, ASAP, designed for interactive analysis of airway geometry for 2-D slices of a 3-D data set. While the most accurate airway measurements will certainly be based on 3-D data analysis, in many clinical cases it remains impractical to routinely gather and analyze 3-D data. Additionally, 3-D airway segmentation of the airway tree (including airways as small as 1--2 mm in diameter) is an area of ongoing research. ASAP is intended as a stop-gap until true 3-D techniques are fully developed.

In our method, we assume that the airway to be measured (such as before and after methacholine challenge, for instance) is scanned so that the long axis is approximately perpendicular to the imaging plane. Alternatively, ASAP can accept reformatted sections selected from a stack of slices using other modules of VIDA (sent to ASAP through the shared-memory facility in VIDA) [19].

Using a user-specified seed point in the airway lumen, ASAP will automatically detect and measure the inner and outer airway walls, and compute lumen cross-sectional area and wall thickness. ASAP uses two methods to identify and measure the airways: (1) a ray-casting approach with automatic selection of the thresholds used to determine wall locations; and (2) a threshold-based seeded region growing method, with automatic determination of the region-growing threshold. Both methods require minimal user interaction (in the best case, just one mouse click per airway), and provide accurate, repeatable measurements of airway geometry. The experimental results show the measurements made by the system are very accurate when applied to the CT scan of a known plexiglass phantom. Using ASAP, the mean absolute errors for the inner and outer diameter measurements of the phantom were less than one pixel (averaged across all seven tubes). The largest errors occurred for phantom tube 3, which has a geometry that causes the most inaccuracy in our approximation used to select wall detection thresholds. We have since developed a more sophisticated lookup-table method to eliminate much of this inaccuracy. This result indicates that the commonly-used half-maximum approach to edge localization must be adapted to reflect the approximate airway size and imaging point spread function. Further, the studies show the measurements are objective and repeatable: there is strong inter-observer (r=0.98) and intra-observer (r=0.98) agreement when applying the system to real human CT scans. For both the normal data and the CF data, the mean absolute difference (averaged over all observations) between the two observers using ASAP was less than one pixel for the inner and outer diameter measurements.

6 ACKNOWLEDGEMENTS

This project was supported in part by the National Library of Medicine (contract N01-LM-4-3511).

7 References

1

E. A. Hoffman and P. B. Heffernan, ``A computer graphics-aided 3-D analysis of heart/lung interaction reconstructed via DSR scanning,'' Proceed Comput Graphics, pp. 81--92, 1985.

2

M. Block, Y. H. Liu, D. Harris, R. A. Robb, and E. L. Ritman, ``Quantitative analysis of a vascular tree model with the dynamic spatial reconstructor,'' J. Comp. Assist. Tomogr., vol. 8, pp. 390--400, 1984.

3

D. P. Boyd and M. J. Lipton, ``Cardiac computed tomography,'' Proceedings of the IEEE, vol. 71, pp. 298--307, 1983.

4

M. Sonka, G. Sundaramoorthy, and E. A. Hoffman, ``Knowledge-based segmentation of intrathoracic airways from multidimensional high resolution CT images,'' in Proc. SPIE Conf. Medical Imaging, vol. 2168, pp. 73--85, 1994.

5

G. Sundaramoorthy, E. A. Hoffman, J. Hoford, T. Mitsa, J. Qian, and M. Sonka, ``Knowledge-based intrathoracic airway segmentation: Computer model and physical phantom based validation,'' in Proc. SPIE Conf. Medical Imaging, vol. 2168, pp. 86--97, 1994.

6

M. Sonka, W. Park, and E. A. Hoffman, ``Validation of an enhanced knowledge-based method for segmentation and quantitative analysis of intrathoracic airway trees from three-dimensional CT images,'' in Proc. SPIE Conf. Medical Imaging, vol. 2433, pp. 158--166, 1995.

7

S. A. Wood, E. A. Zerhouni, J. D. Hoford, E. A. Hoffman, and W. Mitzner, ``Measurement of three-dimensional lung tree structures by using computed tomography,'' J. Applied Physiology, vol. 79, no. 5, pp. 1687--1697, Nov. 1995.

8

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,'' Amer. Review Respir. Disease, vol. 144, pp. 208--212, 1991.

9

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.

10

A. E. McNamara, N. L. Muller, M. Okazawa, J. Arntorp, B. R. Wiggs, and P. D. Pare, ``Airway narrowing in excised canine lungs measured by high-resolution computed tomography,'' J. Appl. Physiol., vol. 73, pp. 307--316, 1992.

11

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 Journal, vol. 7, pp. 596--600, 1994.

12

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.

13

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, no. 4, pp. 1504--1514, 1993.

14

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, pp. 2239--2250, 1993.

15

R. H. Brown, C. J. Herold, C. A. Hirshman, E. A. Zerhouni, and W. Mitzner, ``Individual airway constrictor response heterogeneity to histamine assessed by high-resolution computed tomography,'' J. Appl. Physiol., vol. 74, pp. 2615--2620, 1993.

16

R. H. Brown, W. Mitzner, E. Zerhouni, and C. A. Hirshman, ``Direct in vivo visualization of bronchodilation induced by inhalational anesthesia using high-resolution computed tomography,'' Anesthesiology, vol. 78, pp. 295--300, 1993.

17

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, 1994.

18

R. C. Gonzalez and R. C. Woods, Digital Image Processing. Reading, MA: Addison-Wesley, 1992.

19

G. McLennan, S. Shamsolkottabi, and E. A. Hoffman, ``Assessment of major airway obstruction using image analysis of digital CT information,'' in Proc. SPIE Conf. Medical Imaging, vol. 2709, Newport Beach, CA, 10--15 Feb. 1996.