## Abstract

To objectively quantify airway geometry from three-dimensional computed tomographic (CT) images, an idealized (circular cross section) airway model is parameterized by airway luminal caliber, wall thickness, and tilt angle. Using a two-dimensional CT slice, an initial guess for the airway center, and the full-width-half-maximum principle, we form an estimate of inner and outer airway wall locations. We then fit ellipses to the inner and outer airway walls via a direct least squares fit and use the major and minor axes of the ellipses to estimate the tilt and in-plane rotation angles. Convolving the airway model, initialized with these estimates, with the three-dimensional scanner point-spread function forms the predicted image. The difference between predicted and actual images is minimized by refining the model parameter estimates via a multidimensional, unconstrained, nonlinear minimization routine. When optimization converges, airway model parameters estimate the airway inner and outer radii and tilt angle. Results using a Plexiglas phantom show that tilt angle is estimated to within ±4° and both inner and outer radii to within one-half pixel when a “standard” CT reconstruction kernel is used. By opening up the ability to measure airways that are not oriented perpendicular to the scanning plane, this method allows evaluation of a greater sampling of airways in a two-dimensional CT slice than previously possible. In addition, by combining the tilt-angle compensation with the deconvolution method, we provide significant improvement over the previous full-width-half-maximum method for assessing location of the luminal edge but not the outer edge of the airway wall.

- airway area
- bronchial tree
- lung imaging
- quantitative computed tomography
- deconvolution

multidetector row computerized tomography (MDCT) can produce high-resolution images that allow direct measurements of the bronchial tree geometry. MDCT imaging provides localized information and direct measurement of lung physiology and structure compared with global information and indirect measurements obtained by pulmonary function tests. In addition, MDCT allows in vivo measurements repeated across lung volumes and across short and long time periods, limited only by acceptable radiation exposures. Such methodology now allows for measures previously limited to ex vivo analysis of airway casts (12, 31, 35) or in vitro muscle reactivity studies.

Numerous methods have been proposed and used for the measurement of airway size and geometry from X-ray computed tomography (CT) images. Simpler, less objective methods for the assessment of airway dimensions include manual tracing and measurement either on film, computer screen, or projected image with the adjustment of window and level for better visualization (2, 4, 21, 29, 34). More sophisticated automatic and semiautomated methods have also been proposed, such as the full-width-half-maximum (FWHM) method (1, 3, 9, 37). FWHM has been used in numerous X-ray CT-based imaging studies (1, 5). The FWHM method assumes that the true edge position is halfway between the maximum and minimum values of a line profile across the edge. Region of interest-based methods using thresholding and region growing (13, 15, 18, 27, 36) select thresholds to determine edge location. Model-based approaches (24) use models for the airways and the scanner and try to match the model results to the actual measurement. Heuristic analysis and morphological filtering (7, 15, 23) use a priori information to form an intelligent guess about the geometry of airways.

Studies to date have largely limited their assessment of airway dimensions to those airways that happen to be imaged perpendicular to their local long axis so as to avoid the overestimation of wall thicknesses or luminal diameters associated with oblique sectioning of a tubular structure (1, 2, 4, 6, 17). Efforts have been underway for a number of years to devise methods of automatically extracting lung structure, including airway tree geometry, from volumetric CT images (11, 22, 24, 37). To further advance our understanding of how to best take advantage of the unique in vivo airway information offered by advanced CT imaging technology, in the present study, we have developed and tested a new method for the estimation of airway geometry from two-dimensional (2D) CT slices and have compared the results to the less computationally intense tilt-adjusted FWHM criterion. FWHM uses a simple ray-casting technique, and the wall locations are estimated on the basis of the gray-scale profile along the rays radiating from a central point outward. Our new method incorporates a three-dimensional (3D) airway model and the 3D point spread function (PSF) of the scanner, following the steps of King et al. (14) and Wang and colleagues (28, 32, 33), similar to the method used in a 2D form by Reinhardt et al. (24). Our new method is used to estimate airway inner and outer radii from 2D CT slices as well as to estimate the airway long axis tilt angle relative to the scan plane. The method allows for the quantitation of airways obliquely oriented to the scan plane. Assuming an ideal circular airway model arbitrarily oriented in space, we use mathematical deconvolution to extract the model parameters from the actual image knowing the scanner (PSF) characteristics. This is done in two phases: *1*) the initialization phase; and *2*) the optimization phase that generates an estimate for the inner and outer wall locations and the airway axis tilt angle.

## METHODS AND MATERIALS

The model-based and the FWHM tilt-compensated methods were evaluated by using a Plexiglas phantom in which tubes of known geometry were embedded in a pulmonary parenchymal-equivalent medium (potato flakes). Tubes were selected with varying inner and outer diameters. Results of the two quantitative methods for estimating the inner and outer diameters of the tubes are compared. In addition to presenting the two methods, we will describe the calibration steps needed to estimate the scanner model parameters. Reinhardt and colleagues (24) have previously validated the 2D model-based method alone by comparison with microscopy-based measures in pathology specimens (26). The two methods share some initial steps that will be described next (*phase 1*) followed by a description of the specific steps of each method (*phase 2*). See Fig. 1 for a complete overview of the methods.

### Phase 1

*Center determination*. The center of the airway is determined manually by using a point-and-click device on a computer screen. Although, for the purposes of this study, we did not concentrate on a “user friendly” integration of methods, the center point could be read from a file containing center coordinates generated automatically by a centerline-finding software (30, 36, 37) or by other automated means of identifying the airway locations automatically (15, 22). From this initial center, equally spaced (2π/*n, n* = 40) rays independent of airway size are cast out radially with enough length to cross both the inner and outer walls of the airway. The length of each ray is set to fixed length unless no outer wall is detected, which would increase the preset length until the outer wall is found.

*FWHM*. The finite bandwidth of the CT imaging system causes edges to appear blurred. Thus the location of an edge cannot be determined precisely from an image, because the sharp intensity transition at the edge is replaced by a more gradual gray-level transition. The true location is somewhere between the maximum and minimum gray levels along a ray across the edge. In the case of airways, there are two edges across the wall, one at the luminal surface and one at the boundary between the airway and other pulmonary structures that include parenchyma and associated pulmonary arteries. A line profile across the airway wall (exclusive of regions where the airway borders with pulmonary arteries) would look similar to the profile shown in Fig. 2. Let the line profile be a function *Y*(·) of the distance *r* from center of airway, *Y*(*r*). Define the location of the inner edge to be at a level *y*_{i} between *Y*(*r*_{max}) and *Y*(*r*_{min1}) and similarly that for the outer edge to be *y*_{o} between *Y*(*r*_{max}) and *Y*(*r*_{min2}). The FWHM method defines the levels *y*_{i} and *y*_{o} to be halfway between the maximum and minimum observed gray-level values (1, 3, 9, 24, 36) 1a 1b These correspond to the locations of the true edges being at a distance *r*_{i} and *r*_{o} from the center. Mathematically we can express this relation as *r*_{i} = *Y*^{-}^{1}(*y*_{i}) [or *r*_{o} = Y^{-}^{1}(*y*_{o})], where Y^{-1}(*r*) is the inverse of the function *Y*(*r*). Because the function *Y*(*r*) is not defined mathematically for real image data, the values *r*_{i} and *r*_{o} are found by searching along the ray. Final results might require interpolation for better *r*_{i} and *r*_{o} estimates.

The locations of the inner and outer airway walls are estimated by using the FWHM method along each ray of the set of rays generating two sets of points: one for the inner wall and the second representing the outer wall. To date, this method has been limited to measuring airways that perpendicularly cross the image plane. The following procedures correct for tilt angle to allow estimation of tube radii that are not perpendicular to the image plane.

*Ellipse fitting*. In CT images of the lung, airway appearance ranges from nearly circular to elliptical. The ellipse is the general geometric shape of the cross section of an ideal circular airway intersecting with the image plane. In the special case of the ideal airway being perpendicular to the image plane, the circular cross section can be viewed as a special-case ellipse.

To mathematically describe airways, an ellipse model can be fit to the edge points determined from the 40 rays and the FWHM estimation of edges. Ellipses can be defined as general conic sections (10) by using an implicit second-order polynomial *F*(**a**, **x**) 2 where **a** = [*a b c d e f*] and **x** = [*x*^{2} *xy y*^{2} *x y 1*]^{T}. Fitting a general conic to a set of points can be accomplished by minimizing the sum of squared algebraic distances between the conic *F*(**a**, **x**) = 0 and the set of *n* airway edge points (*x*_{i}, *y*_{i}), with some constraints applied to vector **a**. This method is ellipse specific, invariant to affine transformation, robust to noise, and numerically very efficient (10).

The fitted ellipses provide estimates that are used to calculate a refined approximation of the inner (*R*_{i}) and outer (*R*_{o}) radii of the airway, in-plane rotation angle (θ), and the long-axis tilt angle relative to the scan plane (ϕ) (see Fig. 3). In addition, it provides an adjustment to the initial center in both the *x* and *y* directions. The in-plane rotation and tilt angles are estimated by using the following formulae 3a 3b where ϕ is the tilt angle, *r*_{major} and *r*_{minor} are the semimajor and semiminor axes of the ellipse, θ is the in-plane rotation angle, and *a, b*, and *c* are defined in a general quadratic curve equation (*Eq. 2*). In this study, we compensate for the tilt angle (ϕ) by including the estimated values in the airway model. Because the model assumes the major and minor axes of the ellipse to be along the horizontal and vertical axes, we compensate for the in-plane rotation angle (θ) by rotating the original image by (-θ) and applying the model on the rotation-adjusted image.

### Phase 2

*FWHM method*. Because the initial FWHM measurements were obtained from the actual 2D cross sections, we adjust those measurements by the tilt angle estimated from the ellipse-fitting step. To reduce the noise-induced bias in the final measurements, we assume that the airway is symmetrical around the center and along the horizontal and vertical axes and average the two horizontal estimates and the two vertical estimates of radii. Airway circularity has never been conclusively proven; however, several researchers have assumed the airways to be circular or near circular either locally or globally and have applied measures of circularity to quantitatively describe airways (12, 13, 16, 24, 27, 35). This assumption is not valid in some severe pathological states. However, we seek to devise methods for the assessment of early lung disease before dramatic disease-based disruption of normal airway geometry.

*Model-based method*. In the model-based method, in addition to modeling the airway cross-sectional geometry as a circular tube that is tilted as it intersects the scan plane, we also model the scanner and the scanning process. We model the scanner by its 3D PSF described by a 3D Gaussian with standard deviations in the *x, y*, and *z* directions. The scanning process is modeled as a linear convolution of the airway and the scanner models (see *Eq. 5*). However, as discussed earlier, the cross-sectional appearance of the airway is, in general, elliptical in shape rotated by an angle θ (see Fig. 3). To determine *r*_{i} and *r*_{o} of the embedded tubes using the model-based method, we first generate a 3D data set using the tilt angle parameters from ellipse fitting (see *Ellipse fitting, Eqs. 3a* and *3b*), using 0 Hounsfield units (HU) for airway wall, -1,000 HU for lumen, and estimated values for lung parenchyma (α_{i}) outside the lumen. This sharp image is then convolved with the 3D PSF of the scanner to produce a predicted image. The profiles along the major and minor axes of the ellipse thus produced are compared with the actual profiles from the original image, and a sum of squared deviations is computed. The final values for *r*_{i} and *r*_{o} are then determined by adjusting these values to minimize the sum of squared deviations between model and actual image line profiles (see Fig. 1) by using the multidimensional unconstrained nonlinear minimization simplex method (Nelder-Mead simplex method as implemented in Matlab, The Mathworks, Natick, MA, and described in Ref. 19).

psf of the scanner. The PSF represents the resolution-limiting factors that cause images not to be exact replicas of the real objects (13). An ideal PSF is an impulse, and any deviation from this ideal function causes the PSF to widen. The increase in PSF width increases the blur that images incur as a result of being generated by this nonideal system. In short, the PSF or its frequency domain version, the modulation transfer function, specifies the resolution of the scanner (32). The model Reinhardt et al. (24) used considered only perpendicular airways, and therefore the scanner can be modeled by a 2D PSF. Because of the added degree of complexity in our airway model, the use of a 2D PSF does not fully represent the image-formation process. In the case of obliquely oriented airways, the effect of slice thickness influences the 2D reconstructed image because of partial volume effect. Using a 3D PSF better captures this phenomenon. Therefore, we modeled the 3D PSF as a spatially invariant Gaussian function (24, 28, 32, 33) of the form shown in *Eq. 4* 4 where *psf*(*x, y, z*) is the 3D PSF, and σ_{x}, σ_{y}, and σ_{z} are the standard deviations in the *x, y*, and *z* directions, respectively. The PSF parameters σ_{x}, σ_{y}, and σ_{z} must be estimated by a calibration procedure (24, 32).

The scanning process is modeled as the 3D linear convolution of the 3D ideal airway model with the 3D ideal scanner model (24, 33) as shown in *Eq. 5* 5 where *g*(*x, y, z*) is the 3D ideal airway model, *psf*(*x, y, z*) is the 3D PSF of the scanner, and *f*(*x, y, z*) is the resultant 3D blurred image.

psf estimation. System impulse response is the output of the system when stimulated by an impulse function. In 3D, we use a small spherical bead of a known size to simulate a 3D impulse function, and the resulting image provides the information needed to estimate the PSF of the scanner. Three different types of beads were used: 1.4 mm metal; 1.0 mm zirconium; and 0.794 mm sapphire (Small Parts, Miami Lakes, FL). The sapphire beads were used to estimate the in-plane standard deviation (σ_{x} and σ_{y}), and the larger beads were used to estimate the standard deviation in the *z* direction (σ_{z}). A Gaussian function (*Eq. 6*) is used to model the scanner characteristic function (PSF) (14, 24, 28) as shown below 6 where σ_{x}, σ_{y}, and σ_{z} are the standard deviations and *x*_{0}, *y*_{0}, and *z*_{0} are the bead center offsets in the *x, y*, and *z* directions, respectively; *C* is a constant offset representing the image background (in this case to be -1,000 HU representing air); and *A* is a scaling factor.

To solve for the PSF parameters, we treated the system as a nonseparable system and fitted the data obtained in 3D to the model of *Eq. 6*. This method has the advantage of providing a more realistic model of the PSF than the separable model because the slice thickness for the Imatron electron beam CT (EBCT) scanner is dependent on in- and out-of-plane factors that render the image voxels acquired to be nonisotropic. One of these factors is the X-ray beam shape and the relative location of the targets, detectors, and object imaged relative to the gantry (8).

Because of the limited dynamic range of the scanner and the high density of the beads, the images of the metal and the zirconium beads suffered from saturation. To model this saturation, we clipped the output *f*(*x,y,z*) when the simulated image value exceeded the upper limit of the scanner (3,095 HU). The clipping occurs at the peak of the Gaussian PSF while the transitional parts are preserved. Considering the other sources of error introduced because of the geometry and the nature of the Imatron scanner, this loss of information is not significant, and the PSF of the scanner remained recoverable by using the approach mentioned above.

The beads were scanned in the Imatron C-150 scanner (Imatron, South San Francisco, CA) using the smallest pixel size of 0.1758 mm corresponding to a field of view (FOV) of 9.0 cm. The slice thickness was 1.5 mm and an overlap of 1 mm, making the effective slice thickness 0.5 mm. The scan aperture was 100 ms, and normal, sharp, and very sharp reconstruction kernels were used.

### Airway Phantom

The airway measurement techniques were validated by using a Plexiglas phantom. The phantom consists of seven Plexiglas tubes resembling airways with various inner diameters and wall thickness values. The tubes are contained in a Plexiglas cylinder and surrounded by potato flakes to simulate the lung parenchyma because potato flakes were found to have a similar HU value (24) to that of the lung parenchyma at functional residual capacity (approximately -650 HU). Figure 4 shows representative CT images of the Plexiglas phantom at two different tilt angles (0 and 45°) with respect to the scan plane.

The five smallest tubes were used for validation. The tube dimensions are shown in Table 1. The images were acquired using the Imatron C-150 EBCT scanner. Three FOVs were reconstructed for each data set. The FOVs were 15, 26, and 35 cm, corresponding to in-plane pixel size of 0.29 × 0.29, 0.51 × 0.51, and 0.68 × 0.68 mm, respectively. Slice thickness was 3 mm for all cases, consistent with the protocol used at the time for human scans. The scan aperture was 100 ms, and normal, sharp, and very sharp reconstruction kernels were used at the 15-cm FOV and with no tilt. In addition, the phantom was scanned at three different tilt angles (0, 22.5, and 45° relative to the scan plane) while reconstructed with the normal reconstruction kernel and a 15-cm FOV. Only those at 0° tilt were reconstructed with the three different reconstruction kernels and with the three different FOVs (15, 26, and 35 cm). The FOV of the other tilt angles (22.5 and 45°) was 15 cm, and they were reconstructed with the normal kernel. Ten different slices were analyzed for each tube.

To compare the performance of the new model-based and the tilt-adjusted FWHM methods, the absolute error of the two methods for each combination of imaging parameters was analyzed by using the standard paired *t*-test and its nonparametric counterpart the Wilcoxon signed-rank test. Nonparametric tests are more robust for small samples (in this case *n* = 10) and do not rely on assumptions of normality for the underlying distributions. The null hypothesis, that the absolute errors were equal, was tested against the alternate hypothesis that the absolute errors from the model-based approach were less than the errors from the tilt-compensated FWHM method (a one-sided test). The results of both statistical tests were evaluated against an alternative hypothesis with α = 0.05 (probability of a type I error, false positive, equal to 5%).

## RESULTS

### Tilt Angle Estimation

The error in the estimation of the tilt angle is important because of the utilization of tilt angle information in the model to estimate the airway size. In numerical simulations, we were able to estimate the tilt angle in the range of 0-90°. However, actual phantom data analyzed was for 0, 22.5 and 45° of tilt. The error in tilt angle estimation is shown in Fig. 5.

Figure 5*A* shows errors in range of tilt angles, 15-cm FOV, and normal reconstruction kernel. The figure shows the error in estimating the tilt angle to be within 4.0° over the range of tilt angles considered (0-45°) for all the tubes except the first tube (the smallest), especially at the largest tilt angle (45°) where the lumen of the tubes almost disappears (Fig. 4*B*). The results show that the method tends to overestimate the tilt angle.

Similarly, Fig. 5*B* shows the effect of reconstruction kernels. The error is estimated to be <4° for all the tubes and slightly higher for the smallest tube. In this case, we notice that the tilt angle is overestimated for thin-walled tubes and slightly underestimated in the case of the thick-walled tube. In addition, it appears that the mean error differences for the same tube using all the kernels are small compared with the difference in case of different FOVs or tilt angles.

Figure 5*C* shows the effect of altering the FOV. The tilt angle estimation error is within 2.3° for all FOVs except for the first and the fourth tubes. It is worthy to note that at these FOVs the smallest (first) tube is covering fewer pixels as the FOV increases, meaning that there are fewer data points available for parameter estimation. However, for the fourth tube, the results are skewed because of the presence of some potato flakes inside the tube at the time of the scanning, which does not agree with the model and the assumption that the lumen is uniformly filled with air. The tilt angle estimation method seems to overestimate the degree of tilt for thin-walled tubes and underestimates the tilt for the thick-walled tube (*tube 5*).

### FWHM Method

Figure 6, *A* and *C*, Fig. 8, *A* and *C*, and Fig. 9, *A* and *C* show the mean error ± SE in estimating the inner and outer radii of the phantom tubes listed in Table 1 using the modified FWHM method. There is a clear underestimation bias in measuring the inner radius of the tubes by this method. Figure 6*A* shows the mean error of the estimation of the inner radius by using images acquired at different tilt angles and reconstructed with a normal kernel at a 15-cm FOV. As the tilt angle increases, the error in estimating the inner radius increases, and in the worst case it becomes approximately one pixel (pixel size = 0.29 mm). Similarly, the same trend (underestimation of inner radius) is shown when considering the different FOV data and to a lesser degree in the case of different reconstruction kernels (see Figs. 8*A* and 9*A*). It is important to note that the estimation error is approximately the same for all the FOVs as measured in millimeters. However, as the FOV increases, the pixel size increases and the error value becomes smaller when taken relative to the pixel size. When considering the images with different reconstruction kernels, the mean error decreases as the image sharpness increase except for the thick-walled fifth tube.

Considering the results for the estimation of the outer radius of the tubes (Figs. 6*C*, 8*C*, and 9*C*), the tilt-adjusted FWHM method overestimates the outer radius of thin-walled tubes and appears to underestimate the outer radius in the case of the thick-walled tube for all the tilt angles, FOVs, and normal and sharp reconstruction kernels considered. An exception to this bias pattern is images reconstructed with the very sharp reconstruction kernel. In this particular case, one cannot conclude that the tilt-adjusted FWHM method generates any specific bias pattern. As mentioned earlier, although the numerical value of the error is the same for all FOVs, the relative-to-pixel value of error is smaller as the pixel size increases.

Quantitatively, the tilt-adjusted FWHM method is able to estimate the outer radius to within 0.13 mm (∼1/2 of the smallest pixel) for all the cases. The only exception was in the case of the smallest tube analyzed with its long axis tilted at an angle of 45° to the scan plane. However, if errors were reported relative to tube size instead of absolute error values, then a trend of smaller errors occurring with larger tubes is seen (Fig. 7 in the case of different tilt angles). It is important to note that the radii reported here are calculated by using the tilt-adjusted FWHM method and not raw measurements. These results are estimations of the inner and outer radii along the horizontal and vertical axes in the image adjusted with the estimated tilt angle resulting from the ellipse-fitting method. Although the raw FWHM measurements were comparable to the tilt-adjusted FWHM measurements when the tilt was 0°, the radius estimation error becomes as large as 2 mm when measured with other tilt angles (22.5 and 45°).

### Model-Based Method

The mean estimation error is shown in Figs. 6*B*, 8*B*, and 9*B* for the inner radius and Figs. 6*D*, 8*D*, and 9*D* for the outer radius of the phantom tubes listed in Table 1 by using the model-based method and different image parameters. This method receives its input from the previous stage (FWHM and ellipse fitting), and, by minimizing the difference between the generated model and the actual image, final estimates for the inner and outer radii of the tubes are generated.

Figure 6*B* shows the mean estimation error of the inner radius of the tubes obtained by using the model-based method imaged at three different tilt angles and reconstructed with a normal reconstruction filter and a 15-cm FOV. The results show a weak pattern of overestimation of the inner radius with the exception of the large-tilt case. In addition, the effect of the tilt angle shows no specific pattern on the results. However, the error estimates are in general smaller than one-third of a 0.29-mm pixel (<0.10 mm), with the exception of the thick-walled tube at all tilts and the smallest tube at 45° tilt, where the errors increase to about one-half of a pixel.

In examining the outer radius estimation results, error reported does not follow a specific pattern and is distributed on both sides of the horizontal axis. There was no specific influence for the tilt angle on the magnitude of the error except for the smallest tube, where error increases with tilt angle. Quantitatively, the error is within one-third of a 0.29-mm pixel (<0.1 mm), except for the smallest tube at the greatest tilt considered.

Applying the model-based method to images reconstructed with different FOVs, Fig. 8*B* shows the estimation error of the inner radius. These images were reconstructed with the normal reconstruction filter and no tilt. The error values show that no pattern of bias exists; the thick-walled tube is overestimated as previously found and only the smallest tube is underestimated. Overall, the error is estimated to be within one-half a 0.29-mm pixel for all the tubes at all the FOVs. The error pattern does not appear to be affected by the change in the FOV. However, because the increase in FOV causes the pixel dimension to increase, therefore there would be a pattern such that as FOV increases the relative error decreases.

Similarly, Fig. 8*D* shows the results of the model-based approach applied to images of the phantom reconstructed with different FOVs, normal reconstruction kernel, and no tilt to the direction of table motion. The results show a somewhat weak pattern of overestimation of the outer radius of the tubes, especially the thin-walled ones. In addition, the model-based method seems to perform equally well with all FOVs and generates errors that are within a half a pixel relative to the pixel size.

Figure 9*B* shows the inner radius estimation error when analyzing images reconstructed with three different reconstruction filters, a 15-cm FOV, and no tilt. The error values show a clear overestimation of the inner radius for all the tubes. In addition, the bias pattern exhibits an increase in error with the increase of the sharpness of the reconstruction filter to the extent of tripling the radius estimation error. Numerically, the error reported is within approximately one 0.29-mm pixel.

Lastly, Fig. 9*D* shows the outer radius estimation error. In the case of the model-based method with different reconstruction filters, a weak underestimation bias is noticeable, with the smallest tube as an exception. The results do not suggest any dependence of the error values on the kernel used. In all the cases, the estimation error is within one-half a 0.29-mm pixel. In general, the error estimates of the outer radii are higher when FOVs other than 15 cm or kernels other than the normal are used.

The model-based results above have been reported in terms of the absolute error. However, when these mean error values are reported relative to the tube size (see Fig. 7 for the tilt cases and in similar figures for the different FOVs and reconstruction kernels), a pattern of error is seen that demonstrates that, as the tube radius increases, the percent error value decreases. The error measured for all tilt angles slightly deviates from this pattern. In this case, with the exception of the smallest tube at the greatest tilt, the error is within 10% of the tube inner radius and within 6% of the outer radius, irrespective of the tube size (see Fig. 7).

### Method Comparison

For each combination of imaging parameters (Figs. 6, 8, and 9), the tilt-adjusted FWHM and model-based methods were found to be significantly different (*P* < 0.05) when using both the standard paired *t*-test and the more robust nonparametric Wilcoxon signed-rank test when used in estimating the inner radii. This showed that the absolute measurement errors when using the model-based approach were statistically smaller than those when using the modified FWHM method. However, the absolute errors from estimating the outer radii with the two methods were not found to be statistically different.

## DISCUSSION

Using plastic tubes to simulate cylindrical structures, we were able to accurately estimate arbitrary tube geometry relative to the scan plane from 2D HRCT slices. Luminal dimension, wall thickness, and tilt angle are estimated by using a proposed model-based method that incorporates information about the PSF of the scanner, and the results are compared with those of the tilt-compensated FWHM method. The tilt angle estimation was found to have an accuracy of within 4° from the true tilt angle regardless of the image parameters or the actual degree of tilt (0-45°). The only exception was the smallest tube, especially at the largest tilt. This is to be expected because, owing to partial volume effects, the tilt causes the lumen of the tube to disappear and therefore makes it impractical to fit the actual data to the model. When applying this method to images obtained by using scanners with improved spatial resolution, one would expect to improve the measures obtained from smaller sized structures. It should be noted that the spatial resolution of the EBCT scanner is on the order of 9-10 line pairs per centimeter, whereas the line pair resolution of some presently marketed MDCT scanners is specified to be up to 24 line pairs per centimeter in their highest resolution mode. The loss of luminal resolution is further complicated by the same effect occurring between the outer airway wall and the parenchyma. This effect is more pronounced with larger voxel size, where a voxel has a higher chance of containing more nonhomogeneous tissue that is summarized into one gray level and with thinner walled tubes, relative to the PSF, which causes the parenchyma to be partial volumed with the inner wall and the lumen. Because parenchymal density is regionally variable and can vary with body posture and pathology, it becomes more difficult to model the PSF at the outer edge of the airway wall.

It was noted that the tilt angle mean error variation between different reconstruction kernels for the same tube was smaller than the variation in the case of the other image parameters (FOVs or tilt angles). Considering the effect of these parameters on the appearance of the tube cross sections, the reconstruction kernel causes the least change in the overall shape compared with the effects of different FOVs or different tilt angles. As a result of using the ellipse method to estimate the tilt angle, the variation in tilt error estimation increases as the effect of tilt on shape increases.

The tilt-adjusted FWHM method was found to underestimate the inner radius independent of the tilt angle, reconstruction kernel, or size of the pixel. This has been shown earlier to be true by Reinhardt et al. (24) by theoretical simulation and experimental results. They found that, particularly for thin-walled tubes (wall thickness relative to the width of the scanner PSF) and using normal reconstruction kernel, the inner radius suffers more estimation bias than the outer radius. In an exception to this generalization, the thick-walled tube reconstructed with high-spatial-frequency kernels (sharp and very sharp) is overestimated. The combination of the thick wall and the high-spatial-frequency kernels reverses the effects of the FWHM method. By reducing the blurring caused by the finite response of the scanner and the partial volume effect, the profile across the airway wall becomes steeper and more sensitive to small estimation errors.

Considering the model-based method, a less defined overestimating pattern of bias is present for the inner radius. The model-based method shows smaller error values compared with the tilt-adjusted FWHM method. The model-based results for the inner radius are individually significantly different from the results generated by using the modified FWHM method at the significance level *P* < 0.05 and are more accurate when compared using Wilcoxon signed-rank test. This shows that the model devised represents the blurring of the interface of the lumen and the inner wall with good accuracy. On the other hand, the performance of the model-based method declines when analyzing images reconstructed with high-spatial-frequency kernels. This is to be expected because the sharpening kernels change the scanner system model from that of the proposed Gaussian model. Notice that the error in estimating the inner radius more than triples (only doubles for *tube 5*) when analyzing images reconstructed with very sharp compared with normal reconstruction kernel. This suggests that sharpening filters are more useful for visualization and image display rather than model-based quantitative assessment.

Considering the bias patterns for the estimation of the outer radius of the phantom tubes, the tilt-adjusted FWHM method shows an overestimation bias for all the thin-walled tubes and an underestimation bias for the thick-walled tube. In addition, the outer radius estimation bias was less pronounced than the inner radius estimation bias. This is, again, in agreement with what Reinhardt et al. (24) have shown. This bias pattern for the outer radius estimation becomes less prominent when considering images with different reconstruction kernels. The thick-walled tube is unaffected and the method underestimates the outer radius. For the thin-walled tubes, the bias cannot be presumed to follow a certain pattern.

The outer radius estimation using the model-based approach appears to produce good results for the 15-cm FOV and normal reconstruction kernel. It does not show any specific bias patterns compared with the modified FWHM method. The model-based method, however, did not perform as well when applied to images that were reconstructed with other FOVs. In addition to the fewer sampling points represented by larger pixel size, this could be attributed to the process by which the background gray scale (α_{i}) is selected in *phase 1* in combination with the compensation for the in-plane rotation. The background gray scale is taken as the mean of the pixel gray levels comprising the outermost 8% of the profile through the center of the tube. Depending on the size of the tube, the 8% is at least one pixel wide, which with the heterogeneity of the background (which is expected to be more pronounced in actual lung images) might not be an appropriate representative value. The effect of the rotation is more prominent around the edges of the image than closer to the center. Therefore, when an image is rotated, the 8% of the profile would have been interpolated, leading to somewhat altered profile. This could explain why the Wilcoxon signed-rank test shows that the tilt-compensated FWHM method is more accurate than the model-based method. Despite the fact that the error is larger for the multiple-FOV data set, it is still within one-half a pixel relative to the pixel size of the image analyzed. It whether unclear if there is a pattern of bias for the outer radius estimation using the model-based method.

To estimate the scanner characteristics, a scanner-specific calibration step was performed using sapphire 0.4-mm radius beads to estimate σ_{x} and σ_{y}, and the larger metal and zirconium beads where used to estimate σ_{z}. Simulations have shown that the error in estimating the sigma values decreases as the size of the beads becomes smaller. This is in agreement with systems theory, which requires an impulse to be infinitesimally narrow to enable system identification. The scanner resolution and partial volume effects limit the smallest practical bead size. Because the scanner does not provide a slice thickness comparable to the in-plane pixel size, the use of the sapphire beads is limited to the estimation of σ_{x} and σ_{y}. There were only a few sample points along the *z*-axis for the sapphire beads, making an estimation of a σ_{z} unreliable. Therefore, the metal and zirconium beads were used to estimate σ_{z} because they provided greater numbers of samples along the *z* direction (20). Simulations show that beads of the size used in estimating the PSF parameters overestimate in-plane σ values by ∼8%. The experimental results show that the PSF is slightly anisotropic (σ_{x}≠σ_{y}). The geometry and the design of the scanner used (Imatron EBCT C-150) cause the PSF to be anisotropic. Chang (8) has shown that the slice characteristics of the EBCT are nonuniform and vary depending on the position within the image field and on the mode under which the scanner is operated.

The model we are using assumes σ to be constant throughout the image field and is represented by its value near the isocenter, which simplifies the model and introduces error as the region of interest moves away from the center. This intrinsic spatial variability in the scanner PSF and its interaction with the effect of the finite bead size has not been studied. Therefore, the obtained values for the PSF parameters were not corrected for the finite bead size. Because of the unique construction of the Imatron scanner used, the previous assumptions are valid approximations to the actual scanner characteristics provided that target “C” is being used (8). In our case all the images obtained were collected using the C target. Another source of error in estimating the standard deviation in the *z* direction is the problem arising from not having enough samples along the *z*-axis. This is caused by the relative thickness of the slice compared with the finer size of the in-plane pixel. State-of-the-art MDCT scanners can produce isotropic image voxels, allowing a uniform sampling of the PSF (33). These nonuniformity problems might not be as apparent if a different scanner is used, which would require the calibration step to be repeated. Calibration would also be needed even if the same scanner is used but the image is reconstructed with a different reconstruction kernel.

Other factors may be affecting the accuracy of our methods. Error could be attributed, in part, to the slight eccentricity of the phantom tubes. This source of error causes the tubes to seem more elliptical than they are and will result in a higher error especially in the 0°-tilt case. The apparent ellipticity of the tubes is also affected by the error in estimating the tilt angle. To the extent that actual airways are noncircular in cross section, the same error will be present when this method is applied to the in vivo human lung. Another source of error is the approximations implied in the airway model and the scanner PSF. The simple linear model used may not accommodate for nonlinearities in the system function owing to asymmetry in the 3D PSF of the scanner or any proprietary algorithms used in reconstructing the images. The inaccuracies in the estimation of the σ parameters due to noise or the finite size of the beads used are also another possible source of error.

The model-based method is not limited to the analysis of airways that happen to be perpendicular to the scan plane. It allows us to use those airways that are scanned tilted to the scan plane. Thus one is able to evaluate a much larger sampling of airway segments in a 2D CT slice. The 2D evaluation of the airway tree remains of importance even in the presence of newer 3D methods, because of the continued and growing worry of X-ray dose to subjects entered into study protocols requiring assessments to be done repeatedly over prolonged time periods. This is of particular concern in young children. This method can be automated by adding methods for airway tree detection and centerline determination (22, 25, 30). This method not only helps with the 2D assessment of the airway tree, but it can work well in combination with volumetric image sets of the airway (30).

In conclusion, in this paper, we have presented a new method for the measurement of airway geometry from 2D HRCT images with simultaneous estimation of airway tilt angle relative to the scan plane. The new model-based method models the scanning process as a linear convolution between an airway model and a scanner model. The airways are assumed to be circular in cross section and have a uniform lumen and wall composition, and the background parenchyma is assumed to be of varying density. The scanner is modeled as a spatially invariant 3D Gaussian function and is described by its PSF. A calibration step is required to estimate the scanner characteristics. The new method was compared against a modified FWHM method. Correcting for the tilt angle and averaging the measurements along the horizontal and vertical axes of the elliptical pattern representing the airway modified the FWHM technique.

The model-based method is able to estimate airway geometry to within one-half a 0.29-mm pixel or less in size and to within 4° in tilt for all FOVs and tilt angles reconstructed with the normal kernel. One disadvantage is the longer time it takes compared with the near instantaneous results obtained with the tilt-adjusted FWHM method. However, this method eliminates the need for reslicing the data set and increases the number of airways available for evaluation by allowing us to be able to measure airways not previously possible to analyze directly. As computers become faster and with code optimization, the computational cost could be reduced greatly. A hybrid system that combines the tilt-adjusted FWHM, for its better outer radius accuracy, and our model-based method, for its greater inner radius accuracy, will enhance the overall system accuracy as well as decrease the execution time.

## DISCLOSURES

This work was supported in part by a National Institutes of Health (NIH) Bioengineering Research Partnership Grant RO1-HL-064368, Grant 0092758 from the National Science Foundation, and by the National Institute for Environmental Health Sciences (NIEHS) through the University of Iowa Environmental Health Sciences Research Center (NIEHS/NIH P30 ES-05605).

## Acknowledgments

We thank Kenneth Beck for review of the manuscript and invaluable input and Blake Robinswood for help with the statistical analysis.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2003 the American Physiological Society