1.1.

Related Work

There have been various strategies to quantify tapering in the airways. The initial proposed tapering measurements by Odry et al.¹¹ were restricted to short lengths of the airways. A segmented airway would be split into four equal parts. Each segment had an array of computed lumen diameters. The tapering was measured as the linear regression of the lumen diameters along the branch. The method shared similarity to Venkatraman et al.,¹² but the diameter measurements were taken across the central half of each branch. Various analyses attempted to measure the taper of airways containing multiple branches. In the work of Oguma et al.,¹³ the region of interest from the carina to the fifth generation airway was measured; however, this was only performed in patients with COPD. Finally, Weinheimer et al.¹⁴ used a graphical model of the airways for their proposed tapering measurement. The graphical model was based on a graphical tree originating at the trachea and extending into distal branches, depending on airway bifurcations. A tapering measure was assigned to the edge of the graph depending on the lumen area and generation. They also proposed a regional tapering measurement based on the segments of the lobes of the lungs.

The described tapering measurements have two key limitations. First of all, there is no detailed quantification of reproducibility when considering differences in specifications of the CT scanner, or reconstruction kernel, making it difficult to compare tapering statistics from different machines or from the same CT scanner employing different scanning parameters. Second, the region of interest for the tapering measurement was restricted to airways that were segmented using the respective airway segmentation software. Bronchiectasis is a heterogeneous disease—it can affect any area in the lung including the peripheral regions.¹⁵ Thus to encapsulate the disease in the tapering measurement, one would need to consider the region of interest as the entire airway, from the trachea to the most distal point.

In all the proposed tapering measurements, obtaining the cross-sectional area is a necessary input for the algorithms. There have been various analyses attempting to validate the reproducibility and precision of measurements against dose,^16–18 voxel size,¹⁹ reconstruction kernel,^20–22 and normal biological variation.²³ In most of the validation experiments, area measurements were taken from phantom,^20,21 porcine,^16,18 or cadaver²⁴ models. Fetita et al.²⁵ used synthetic models of the lung. None of these experiments were explicitly performed on scans with bronchiectasis. Furthermore, the area measurements were not taken at contiguous intervals along the lumen thus missing possible dilatations from a bronchiectatic airway.

For our work and the method from Oguma et al.,¹³ the tapering measurement involves the computation of the arc length of the airway at contiguous intervals. In the literature, investigations of the reproducibility of arc length computation in airways are limited. The work of Palágyi et al.²⁶ used simulated rotation of in vivo scans. The assessment of the reproducibility was based on the lengths of a single branch rather than multiple generations of branches, thereby precluding estimations of reproducibility of airway quantitation from the carina to an airway’s most distal point.

1.2.

Contributions of the Paper

To our knowledge, there is no detailed analysis on the reproducibility of a global tapering measurement of airways using CT. Thus the purpose of this paper is as follows. First of all, we will discuss in detail a tapering measurement of the airways on CT imaging. Second, we quantify the reproducibility of the measurement against variations in simulated dose and voxel sizes. In addition, we compare the variability of the tapering measurements across different CT reconstructions kernel. Finally, we analyze the effect of bifurcations on tapering measurements and consider measurement repeatability using longitudinal scans.

2.1.

Centerline

The centerline was used to identify and order the airway segments for the tapering measurement. We implemented a curve thinning algorithm developed by Palágyi et al.²⁶ At initialization, the algorithm used the airway segmentation and distal points acquired in Sec. 2. The final input was the start of the centerline at the trachea. The shape of the trachea was assumed to be tubular, with an approximate constant diameter and orientated near perpendicular to the axial slice. Thus the centerline of the trachea lay on the local maximum value of the distance transform of the segmented trachea.³⁰ Algorithm 1 was used to find the centerline start point.

Algorithm 1

Locating the start of centerline on the trachea.

2.2.

Recentering and Spline Fitting

The next task was to separate the centerline of each individual airway from the centerline tree. To this end, we modeled the centerline tree as a graphical model similar to Mori et al.³¹ The nodes corresponded to the centerline voxels and the edges linked neighboring voxels. We performed a breadth first search algorithm³² on the centerline image. Starting from the carina, we iteratively found the next set of sibling branches. When a distal point was found at the end of a parent branch, the path leading to the distal point was saved. The output was an array of ordered paths describing the unique route from the trachea to the distal point. The proposed tapering measurement started at the carina. Thus centerline points corresponding to the trachea were removed from further analysis.

For each path, we corrected for the discretization error—a process known as recentering.³³ We implemented a similar method to that described by Irving et al.³⁴ A five point smoothing was performed along each path. We modeled the centerline as a continuous model by fitting a cubic spline $F:[0,k_n]\to {\mathbb{R}}^3$ denoted as

2.3.

Arc Length

The tapering measurement required an array of arc lengths at contiguous intervals from the carina to the distal point. For our pipeline, we considered small parametric intervals $t_i$ on the cubic spline $\mathit{F}(t)$. At each interval $t_i$, we computed the arc length from the carina to $t_i$ as³⁶

2.4.

Plane Cross Section

We measured the cross-sectional area accurately by constructing a cross-sectional plane perpendicular to the airway. Using the interval $t_i$ from the arc length computation, we computed tangent vector $\mathit{q}\in {\mathbb{R}}^3$ by

From linear algebra, points on the plane can be generated by their corresponding basis vector.³⁷ To this end, we generated a set of orthonormal vectors ${\mathit{v}}_1$, ${\mathit{v}}_2\in {\mathbb{R}}^2$ using the method stated by Shirley and Marschner.³⁸ The method is summarized in Algorithm 2.

Algorithm 2

Constructing the basis for the plane reconstruction, adapted from Ref. 38.

Assuming $\mathit{F}(t_i)$ was the origin, each point $\mathit{u}\in {\mathbb{R}}^3$ on the plane can be written as

We selected the scalars ${\alpha }_1,{\alpha }_2\in \mathbb{R}$ such that the point spacing is 0.3 mm isotropically.

2.5.

Lumen Cross-sectional Area

We calculated the cross-sectional area using the edge-cued segmentation-limited full-width half-maximum (${\mathrm{FWHM}}_{\mathrm{ESL}}$), developed by Kiraly et al.³⁹ The method is as follows. The cross-sectional planes were aligned on both the CT image and airway segmentation. The intensities of the plane were computed for both images using cubic interpolation. Fifty rays were cast out in a radial direction, from the center of the plane. Each ray sampled the intensity of the two planes at a fifth of a pixel via linear interpolation. Thus each ray produced two 1-D images with the first from the binary plane $r_b$ and second from the CT plane $r_c$. We then applied Algorithm 3 to find boundary point $l$.

Algorithm 3

Summary of the FWHMESL, adapted from Ref. 39. The purpose of the algorithm was to find the point of the ray that crossed the lumen.

The final output of the ${\mathrm{FWHM}}_{\mathrm{ESL}}$ was an array of 2-D points corresponding to the edge of the lumen. Finally, we fitted an ellipse based on the least square principle. The method was developed and implemented in MATLAB by Fitzgibbon et al.⁴⁰ We considered the cross-sectional area as the area of the fitted ellipse.

2.6.

Tapering Measurement

We assumed for a healthy airway that the cross-sectional area was modeled by an exponential decay along its centerline. It has been shown in human cadaver studies that the average cross-sectional area in a branch reduces at an exponential rate at each generation.⁴¹ The same observation has been noted in porcine models.⁴² Using the decay assumption, we modeled the relationship between the arc length and the cross-sectional area as

In terms of implementation, for each airway track, we considered the array arc length and cross-sectional area computed for each individual airway. A logarithmic transform $\log (x)$ was applied only on the cross-sectional area array. We fitted a linear regression on the signal; the tapering measurement is defined as the gradient from the line of best fit.

3.1.

Simulated Images

In this experiment, we simulated images with differing radiation dose and voxel size. The purpose was to analyze the reproducibility of the tapering measurement against various properties of the CT image. Furthermore, we varied the noise and voxel sizes at regular intervals. Thus we also analyzed the sensitivity of the tapering measurement against the given parameters. Finally, we investigated the reproducibility of cross-sectional area and airway length measurements with changes in dose and voxel sizes, respectively.

3.1.1.

Dose

To simulate the images acquired with different radiation doses, we used the method adapted from Ref. 44. We performed a Radon transform on each axial slice of the original CT image. The output is a sinogram of the respective axial slice. To simulate different radiation doses, Gaussian noise was added on each sinogram with standard deviation $\sigma ={10}^{\lambda }$, with a range of $\lambda$. The noisy sinograms were then transformed back into physical space using the filtered back projection. The final output is a noisy CT image in Hounsfield units in integer precision. A MATLAB implementation is displayed in Algorithm 4. For our experiment, we varied $\lambda$ from 0.5 to 5 in increments of 0.5. An example of the output image is displayed in Fig. 2.

Algorithm 4

Adapted MATLAB code to simulate noise from differing doses.

function noisySlice = AddingDoseNoise (axialSlice, lambda)

%Creating the sinogram

sinogram = radon(axialSlice,0:0.1:179);

%Adding the Gaussian Noise

noisySinogram = sinogram + randn(size(sinogram))*10^(lambda);

%Converting the noisy sinogram into physical space using Filter Backprojection.

noisySlice = iradon(noisySinogram,0:0.1:179,length(axialSlice));

%Converting into integer precision intensities

noisySlice = int16(noisySlice);

end

Fig. 2

(a) CT images with simulated noise against varying $\lambda$ and (b) an image subtraction of the simulated noisy image with the original. The units are in HU.

To relate $\lambda$ to the physical dose from a CT scanner, we adopted the method described by Reeves et al.⁴⁵ This paper quantified the dose of an image with a homogeneous region in the chest CT scan. To this end, we used the homogeneous region inside the trachea. Using the airway segmentation, we considered the first 60 axial slices of the segmented trachea. To avoid the influence of the boundary, the tracheas were morphologically eroded⁴⁶ with a structuring element of a sphere of radius 5. All segmentations were visually inspected before further processing. Finally, we computed the standard deviation of the intensities inside the mask, denoted as $T_n$. Table 2 shows values of $T_n$ on a selection of images against a range of $\lambda$. Using results from Reeves et al.⁴⁵ and Sui et al.,⁴⁷ a low-dose scan with a tube current-time product 25 mAs has maximum $T_n$ of 55 HU. Thus we assume $\lambda =3.5$ approximately corresponds to a low-dose scan. We considered higher values of $\lambda$ to verify any correlations in the results.

Table 2

Table of standard deviation of intensity Tn (HU) in the inner lumen mask for a selected image against differing λ.

λ	bx500	bx503	bx504	bx505	bx507	bx508	bx511	bx512	bx513	bx515
Ground truth	15	28	15	21	16	28	20	33	21	17
0.5	15	28	15	21	16	28	20	33	21	17
1.0	15	28	15	21	16	28	20	33	21	17
1.5	15	28	15	21	16	28	20	33	21	17
2.0	15	28	15	21	16	28	20	33	21	17
2.5	16	29	16	22	17	28	20	33	22	17
3.0	21	32	21	26	22	31	25	36	26	22
3.5	49	54	49	51	49	54	51	57	51	49
4.0	146	150	148	148	148	148	149	149	148	145
4.5	461	467	467	463	466	462	467	459	463	457
5.0	1456	1473	1477	1464	1474	1458	1475	1449	1463	1445

We computed the taper measurement on the noisy images using the same segmented airways and labeled distal point that were identified on the respective original image. The literature has shown in low-dose scans, airway segmentation software^29,48 cannot segment airways to the lung periphery as well as standard dose scans of the same patient. But these methods can still segment a large number of branches in low- and ultralow-dose scans.^29,48 Furthermore, research has shown that there are minor differences in the performance of radiologists when attempting to detect features from standard and low-dose CT scans.^47,49,50

3.2.

CT Reconstruction

On a subset of images, four patients were scanned using the Toshiba Aquilion ONE scanner. On the same scan, two different images were computed. The images were reconstructed using the lung and body kernels, respectively. An example of the reconstruction kernels is displayed in Fig. 3. We acquired the airway segmentation and distal point from a single-reconstruction kernel as described in Table 3. The tapering measurement was computed on both reconstruction kernels using the same airway segmentation and distal points. We used the same airways as described in Table 1.

Fig. 3

Images from the same CT scan with the (a) body kernel and (b) lung kernel. Both images are displayed in the same intensity window.

Table 3

The images used for the reconstruction kernel experiment. This table lists which reconstruction kernel was used to generate the airways segmentation and distal point labeling. The make, model, and voxels size of the images are displayed in Table 1.

3.3.

Biological Factors

3.3.1.

Effect of bifurcations

We analyzed the effect of airway bifurcations on the tapering measurement. To this end, we manually identified regions of bifurcating airways. On a selected subset of airways, we considered the reconstructed airway image described in Fig. 16. Using ITK-SNAP, the author (K. Q.) started at the cross-sectional plane corresponding to the carina and scrolled toward the distal point. Using visual inspection, the following protocol was developed to identify bifurcations on cross-sectional planes:

• 1. The scrolling stops when the airway is almost or at the point of separation.

• 2. The author scrolls back until the airway stops decreasing in diameter. An alternative interpretation is when the airways are about to enlarge due to the bifurcation.

• 3. Starting at the point of enlargement and scrolling forward, each slice is delineated as a bifurcating region until complete separation of bifurcating airways is reached. The criteria for a complete separation are the lumen wall of both airways that are completely visible and separate. The entire protocol is summarized in Fig. 4.

Fig. 4

(a) A region of bifurcation along the reconstructed slices. The green, blue, and red regions are the slices corresponding to the enlargement, break, and separation slices, respectively. The labeled region consists of slices from green to red. (b) A cross-sectional plane where the airway is at the point of bifurcation, indicated by the blue arrow. (c) First slice of the bifurcation region. (d) The final slice of the bifurcation region. The slides are chronologically ordered with the protocol described in Sec. 3.3.1.

For our experiment, we selected 19 airways from Table 1. The data consisted of 11 healthy and 8 bronchiectatic airways. The entire analysis was performed on ITK-SNAP..

3.3.2.

Progression

We examined possible changes in tapering of airways in patients over time. In this experiment, we consider two sets of longitudinal scans. First, pairs of airways that were healthy on both baseline and follow up scans. Second, pairs of airways that were healthy on baseline scans and became bronchiectatic on follow up scans.

For pairs of healthy airways, a trained radiologist (J. J.), manually identified 14 pairs of airways across 3 patients. The criteria were the airway track that must have a healthy appearance on both baseline and follow up scans. For the second population, the same radiologist manually identified 5 pairs of airways from a single patient P1. The scans were obtained from the University College London Hospital and acquired with written consent. The criteria for selection were airways that appear healthy on baseline scans and became bronchiectatic at the follow up scan, an example is displayed in Fig. 5. Details of the CT images are summarized in Tables 4 and 5.

Fig. 5

The same pair of airways from longitudinal scans: (a) initial healthy airway and (b) the same airway at the same location becoming bronchiectatic.

Using two separate work stations, the airways were visually registered between the longitudinal scans. Airways were taken from various regions of the lungs and were different to the airways displayed in Table 1. The tapering measurements were taken from the method discussed in Sec. 2.

Table 4

List of the images for progression experiment. This table includes time between scans in months (M) and days (D). The airways on this table are different to Table 1.

Patients	Time between scans	Airways labeled
bx500	9M 6D	6
bx504	35M 6D	7
bx510	5M 22D	1
P1	10M 7D	5

Table 5

List of make, models, and voxel sizes of CT images for progression experiment. The voxel sizes are displayed as x, y, z and in mm units. Abbreviations: GEMS, GE medical systems; SS, Siemens SOMATOM.

Patients	Date 1 CT scanner	Date 1 voxel size	Date 2 CT scanner	Date 2 voxel size
bx500	Toshiba Aquilion ONE	0.67, 0.67, 1.00	Toshiba Aquilion ONE	0.56, 0.56, 1.00
bx504	Toshiba Aquilion ONE	0.63, 0.63, 1.00	Toshiba Aquilion ONE	0.78, 0.78, 1.00
bx510	GEMS Lightspeed Plus	0.70, 0.70, 1.00	Philips Brilliance 64	0.72, 0.72, 1.00
P1	GEMS Discovery STE	0.85, 0.85, 2.50	SS Definition AS	0.66, 0.66, 1.50

4.1.

Dose

We analyzed the difference in cross-sectional area measurements and the final tapering measurements at different CT radiation doses.

For the cross-sectional areas, Fig. 7 compares the cross-sectional areas between the original image and one of the noisy images. Each graph contains $\sim 30,000$ unique lumen measurements. The correlation coefficients between the populations were $r>0.99$ on all graphs. The 95% confidence intervals increase with the amount of noise. For the tapering measurement, Fig. 8 displays the measurements from all the noisy images compared to their respective original images. The correlation coefficient between noisy and original tapering measurements was $r>0.98$ on all values of $\lambda$.

Fig. 7

A series of Bland–Altman⁵⁴ graphs comparing area measurements from a simulated low-dose scan and the original image. On all graphs, the correlation coefficient was $r>0.99$.

Fig. 8

A series of Bland–Altman⁵⁴ graphs comparing tapering measurement between simulated dose and the original image. On all graphs, the correlation coefficient was $r>0.98$.

We analyzed the tapering difference between the original images and simulated images. We interpret the mean and standard deviation of the error difference as the bias and uncertainty, respectively. Figure 9 shows an overestimation bias with an increase in noise and a positive correlation between uncertainty and dose.

Fig. 9

(a) Mean and (b) standard deviation of the difference in tapering between original images minus the simulated lower dose.

4.2.

Voxel Size

We analyzed the computed spline and tapering for all the scaled images. We used the arclength of the spline as the metric for comparison for the computed spline. Figure 10 compares the arclengths computed from the scaled splines with the respective originals. On all scales $\sigma$, the correlation coefficients between measurements were $r>0.98$. Furthermore, we analyzed the error difference in arclength. In Fig. 11, the mean difference shows a weak correlation coefficient with $r=0.55$ with scale $\sigma$. The mean difference shows both an overestimation and underestimation bias with the arclength measurement. Figure 11 shows a weak correlation between standard deviation and scale with $r=0.51$.

Fig. 10

A series of Bland–Altman⁵⁴ graphs comparing arc lengths between scaled images and the original images. On all graphs, the correlation coefficient was $r>0.98$.

Fig. 11

(a) Mean and (b) standard deviation of the difference in arclength between original images minus the scaled images. The correlation coefficient of the mean and standard deviation against scale are $r=0.55$ and $r=0.51$, respectively.

In terms of the tapering measurement, Fig. 12 compares the tapering values from the scaled images with the respective originals. The correlation coefficients between the scaled and original tapering values was $r>0.97$ on all scales $\sigma$. In addition, we examined the error difference of the original minus the scaled tapering. Figure 13 shows a negative correlation with both overestimation and scale with $r=-0.98$. Furthermore, Fig. 13 shows a positive correlation with uncertainty and scale with $r=0.94$.

Fig. 12

A series of Bland–Altman⁵⁴ graphs comparing tapering between original images and scaled images. On all graphs, the correlation coefficient was $r>0.97$.

Fig. 13

(a) Mean and (b) standard deviation of the difference in tapering between original images minus the scaled images. The correlation coefficients of the mean and standard deviation against scale are $r=-0.98$ and $r=0.94$, respectively.

4.3.

CT Reconstruction

We analyzed the difference in cross-sectional area and tapering measurement between reconstruction kernel. Figure 14, compares the difference in area measurements. On all patients, in cross-sectional area measurements, the correlation coefficient between the two measurements was $r>0.99$. The largest 95% confidence was in patient bx515 with $\pm 1.98{\mathrm{mm}}^2$ from the mean. Figure 15 compares the differences in tapering measurement. We collected $n=44$ tapering measurement from 4 patients. The correlation coefficient was $r=0.99$ between the reconstruction kernels.

Fig. 14

Bland–Altman⁵⁴ graphs comparing the cross-sectional area between the lung and body kernels. On all four images, the correlation coefficient was $r>0.99$.

Fig. 15

Bland–Altman⁵⁴ graph comparing tapering measurements $(n=44)$ between lung and body kernels $r=0.99$.

4.4.

Clinical Results

4.4.1.

Bifurcations

We compared tapering measurements with and without points corresponding to bifurcations. On the first dataset, the tapering measurements were computed using all area measurements. The second dataset had tapering measurements computed without area measurements from the bifurcating regions as described in Fig. 16. As we compared the measurements in Fig. 17, the correlation coefficient was $r=0.99$.

Fig. 16

(a) A signal of area measurement with bifurcation regions (red) and tubular regions (blue). (b) The same signal with tubular regions (blue) only. On both graphs, the black line is the linear regression of the respective data. The gradient of the line is the proposed tapering measurement. (c) The reconstructed bronchiectatic airway of the same profile. The blue-shaded and red-shaded regions correspond to the tubular and bifurcating airways, respectively. A reconstructed healthy airway has been discussed by Quan et al.²⁷ Similar reconstructed cross-sectional images of vessels have been discussed by Oguma et al.,¹³ Kumar et al.,⁵¹ Alverez et al.,⁵² and Kirby et al.⁵³

Fig. 17

(a) Bland–Altman⁵⁴ graph showing the relationship of the taper rates $(n=19)$ with and without bifurcations $r=0.99$ and (b) comparison of the standard error from linear regression between airways with and without bifurcations. On a Wilcoxon rank-sum test between the two populations, $p=7.1\times {10}^{-7}$.

The uncertainty of each tapering measurement was computed using the standard error of estimate $s$ defined as⁵⁵

Eq. (6)

$$s=\sqrt{\frac{\sum _{i=1}^N{(Y_i-y_i)}^2}{N}},$$

where $y_i$ is the arclength, $Y_i$ is the estimate from the linear regression from each computed area $x_i$, and $N$ is the number of points in the profile. Figure 16 compares the uncertainty between the two populations. There was a statistical difference between the populations, on a Wilcoxon rank-sum test, $p=7.1\times {10}^{-7}$.

4.4.2.

Progression

For healthy airways, we grouped tapering values between the baseline and follow up point. Figure 18 compares the measurements between the two time points. The results demonstrated good agreement with an intraclass correlation coefficient⁵⁶ $\mathrm{ICC}>0.99$. The standard deviation of the tapering difference was $1.45\times {10}^{-3}{\mathrm{mm}}^{-1}$.

Fig. 18

(a) Bland–Altman⁵⁴ graph comparing tapering measurement on the same airways from the first and subsequent scan $\mathrm{ICC}>0.99$ and (b) The change in tapering between healthy airways (H to H) and airways that became bronchiectatic on the follow up scans (H to D), $p=0.0072$.

For airways that became bronchiectatic, we considered the change in tapering, i.e., tapering value at follow up minus tapering value at baseline, the results are displayed in Fig. 18. The results show that bronchiectatic airways have a greater tapering change in magnitude compared to airways that remained healthy $p=0.0072$.

5.1.

Limitations

In this study, we compared the tapering measurement for healthy and diseased airways using a Wilcoxon rank-sum test. The test assumes the data points are independent. However, we used a variety of airways from the same lung. Thus the tapering profiles of the same patients will have a degree of overlap. Future work is needed to analyze data points that are not dependent on each other.

A key limitation of the tapering measurement is the requirement of having a robust airway segmentation. In this paper, the airway segmentation software was often unable to reach the visible distal point of an airway. Thus time-consuming manual delineation was needed to extend the missing airways. The distal point is usually located at the periphery of the lungs. Thus to avoid manual labeling, a segmentation algorithm would need to automatically segment the airways past the sixth airway generation. From the literature, the state-of-the art software developed by Charbonnier et al.⁷¹ using deep learning could still only consistently segment airways to the fourth generation. The segmentation of small and peripheral airways is not a trivial task.^66,72,73

In this paper, we analyze the reproducibility of all computerized components of the tapering algorithm. This paper does not address reproducibility of manual labeling of the airways. It is noted in the literature that semimanual labeling of small airways can take hours.⁷⁴ Future work is required to analyze the reproducibility of manual segmentation of the airways. We hypothesize that the segmented healthy peripheral airways consist of a small number of voxels; therefore, any errors in voxel labeling will be considerably smaller then a dilated peripheral airway affected by bronchiectasis.

In this work, we simulated low-dose scans through performing Radon transforms on existing CT images, adding Gaussian noise on the sinogram and using backprojection to reconstruct noisy CT images. There are proposed methods to simulate a low-dose scans by adding a combination of tailored Gaussian and Poisson noise on the sinogram.⁷⁵ These methods assume the original high-dose sinogram are available for simulation; however, it has been acknowledged that sinograms are generally not available in the medical imaging community.^76,77 Thus various groups have proposed low-dose simulations using reconstructed CT images. The methods involve adding Gaussian^76,77 or a combination of Gaussian and Poisson noise⁷⁸ on the sinogram of the forward projection of the CT image. Although there has been limited validation of the appearance of lung nodules against simulated low dose simulation,⁷⁹ there has been no validation on the efficacy of these methods on the appearance of airways. We believe that our low-dose simulation is sufficient because the measured standard deviation of the trachea mask $T_n$ is similar to results taken from low-dose scans from Reeves et al.⁴⁵ and Sui et al.⁴⁷

Similarly, with voxel size simulation, ideally one would reconstruct the images from the original sinogram.¹⁹ However, as the sinograms were unavailable, we simulated the voxel size through interpolation of the original CT images similar to Robins et al.⁸⁰ We believe that the simulation is sufficient as it shows the robustness and precision of the centerline, recentering, and cross-sectional plane algorithms in the pipeline. Changes in voxel sizes will change the combinatorics or arrangement of the binary image. By showing steps in the pipeline like centerline computation are repeatable across voxel sizes, we avoid resampling the image to isotropic lengths. Thus potentially avoiding a computationally expensive⁴⁸ preprocessing step.

We showed that the tapering measurement is reproducible by measuring the same airway across longitudinal scans with a minimum 5-month interval. The time between scans was on a similar scale from a reproducibility study on airway lumen by Brown et al.²³ An ideal experiment to assess reproducibility of the same airway from different scans would be to acquire follow up scans immediately after baseline scans similar to Hammond et al.¹⁶ However, that work was performed on porcine models. Due to considerations of radiation dose, it is difficult to justify the acquisition of additional scans of no clinical benefit.⁸¹ For our experiment, each airway was chosen by a subspecialist thoracic radiologist. The airway was inspected to ensure it was in a healthy state, e.g., with no mucus present. Thus we assume that each pair of airways is disease free and healthy.