1 Introduction and Purpose
Bronchiectasis is defined as the permanent dilatation of the airways. Patients with bronchiectasis can suffer severe exacerbations requiring hospital admission and have a poorer quality of life[Chalmers2018]. Clinicians diagnose bronchiectasis on computed tomography (CT) imaging by visually estimating the diameter of the airway/bronchus and its adjacent pulmonary artery and calculating the bronchoarterial (BA) ratio. A BA ratio greater than 1 indicates the presence of bronchiectasis[Pasteur2010].
Various groups have proposed methods to automatically and semi automatically compute the BA ratio for bronchiectatic airways[Mumcuoglu2013, PerezRovira2016, Fetita2015]. However, use of the BA ratio to diagnose bronchiectasis has two major flaws. First of all, the maximum healthy range of the BA ratio can be 1.5 times size of the artery[Hansell2010]. Second, blood vessels can change size as a result of factors including altitude[Kim1997], patient age[Matsuoka2003] and smoking status[Diaz2017]. This conflicts with the assumption that the pulmonary artery is always at a constant size.
An alternative approach to diagnose and monitor bronchiectatic airways is to analyse the taper of the airways i.e. the rate of change in the crosssectional area along the airway[Pasteur2010]. In patients with bronchiectasis, the airway is dilated and so the tapering rate must be reduced. Airway tapering is difficult to assess visually and to measure manually. As described by Hansell[Hansell2010], the observer would have to make multiple crosssectional area measurements along the airway. As mentioned in Cheplygina et al[Cheplygina2016], measuring multiple lumen is a manually exhaustive task and prone to mistakes.
1.1 Related Work
There have been various strategies to quantify tapering in the airways. The initial proposed tapering measurements by Odry et al.[Odry2006] were restricted to short lengths of the airways. A segmented airway would be spilt 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.[Venkatraman2006], 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 Oguma et al.[Oguma2015] 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.[Weinheimer2017] 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[Chalmers2015a]. 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 crosssectional area is a necessary input for the algorithms. There have been various analyses attempting to validate the reproducibility and precision of measurements against dose[Hammond2017, LeutzSchmidt2017, Jia2018], voxel size[Achenbach2009] , reconstruction kernel[Wong2008, Leader2007, Zheng2007] and normal biological variation[Brown2017]. In most of the validation experiments, area measurements were taken from phantom[Wong2008, Leader2007], porcine[Hammond2017, LeutzSchmidt2017] or cadaver[Zhang2019] models. In Fetita et al.[Fetita2014], they 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.[Oguma2015], 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. In Palagyi et al.[Palagyi2006] they 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. Secondly, 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 analyse the effect of bifurcations on tapering measurements, and consider measurement repeatability using longitudinal scans.
2 Method
We first, describe in detail the steps to acquire the airway tapering measurement. The method was initially proposed by Quan et al.[Quan2018] and summarised in Figure 1. The pipeline required two inputs. First, the most distal point of each airway of interest was manually identified by an experienced radiologist (JJ). A single voxel was marked at the end of the airway centreline. The entire analysis was completed using ITKsnap^{1}^{1}1http://www.itksnap.org, last accessed September 19, 2019.
Secondly, a complete segmentation of the airway was produced. We obtained an airway segmentation by implementing a method developed by Rikxoort et al.[Rikxoort2009] The algorithm was based on a region growing paradigm. In summary, a wave front was initialised from the trachea. Voxels on each new iteration were classed as airways based on a voxel criterion. The wave front continued until a wave front criteria was met. In certain cases, the airway segmentation was unable to reach the distal points and in these cases we extended the airway segmentation manually to the distal points. Our method is designed so that it can be incorporated to any system that provides both the segmentation and distal point to the airway of interest. Once the inputs were available, we acquired the measurement using an automatic process.
2.1 Centreline
The centreline was used to identify and order the airway segments for the tapering measurement. We implemented a curve thinning algorithm developed by Palagyi et al.[Palagyi2006] At initialisation, the algorithm used the airway segmentation and distal points acquired in Section 2. The final input was the start of the centreline 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 centreline of the trachea lay on the local maximum value of the distance transform of the segmented trachea[Aylward2002]. Algorithm 1 was used to find the centreline start point.
2.2 Recentring & Spline fitting
The next task was to separate the centreline of each individual airway from the centreline tree. To this end, we modelled the centreline tree as a graphical model similar to Mori et al.[Mori2000] The nodes corresponded to the centreline voxels and the edges linked neighbouring voxels. We performed a breadth first search algorithm[Cormen2009] on the centreline 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, centreline points corresponding to the trachea were removed from further analysis.
For each path we corrected for the discretization error  a process known as recentring[Grelard2017]. We implemented a similar method to that described by Irving et al.[Irving2014] A five point smoothing was performed along each path. We modelled the centreline as a continuous model by fitting a cubic spline denoted as
(1) 
where and . The knots where taken on every smoothed point on the centreline. The spline fitting was performed using the cscvn^{2}^{2}2https://uk.mathworks.com/help/curvefit/cscvn.html, last assessed September 19, 2019 function in Matlab. The continuous model should enable computations of the arc length and tangent at subvoxel intervals along the airway.
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 on the cubic spline . At each interval , we computed the arc length from the carina to as[Kreyszig1964]
(2) 
where and is the dot product. For our work, we considered parametric intervals of 0.25 along the spline.
2.4 Plane Cross Section
We measured the crosssectional area accurately by constructing a crosssectional plane perpendicular to the airway. Using the interval
from the arc length computation, we computed tangent vector
by(3) 
From linear algebra, points on the plane can be generated by their corresponding basis vector[Leon2009]
. To this end, we generated a set of orthonormal vectors
, using the method stated in Shirley and Marschner[Shirley2009]. The method is summarised in Algorithm 2.Assuming was the origin, each point on the plane can be written as:
(4) 
We selected the scalars such that the point spacing are 0.3mm isotropically.
2.5 Lumen Cross Sectional Area
We calculated the cross sectional area using the EdgeCued SegmentationLimited Full Width Half Maximum (FWHM_{ESL}), developed by Kiraly et al.[Kiraly2005]
The method is as follows: the crosssectional 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 centre 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 1D images with the first from the binary plane
, and second from the CT plane . We then applied Algorithm 3 to find boundary point .The final output of the FWHM_{ESL} was an array of 2D 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. [Fitzgibbon1996] 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 crosssectional area was modelled by an exponential decay along its centreline. It has been shown in human cadaver studies that the average cross section area in a branch reduces at an exponential rate at each generation[Nikiforov1985]. The same observation has been noted in porcine models[Azad2016]. Using the decay assumption, we modelled the relationship between the arc length and the crosssectional area as
(5) 
where is the arc length of the spline, is the proposed tapering measurement, is the crosssectional area and is an arbitrary constant.
In terms of implementation, for each airway track, we considered the array arc length and crosssectional area computed for each individual airway. A logarithmic transform was applied only on the crosssectional 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 Evaluation
An experienced radiologist (JJ) selected a total of 74 airways from 10 scans. The CT images were analysed from 9 patients with bronchiectasis after obtaining written informed consent at the Royal Free Hospital, London. The voxel size ranged from 0.630.80mm in plane and 0.801.5mm slice thickness. The airways were classified as healthy (n = 35) or bronchiectatic (n = 39) by the same radiologist. Details including the make and model of the scanner are provided in Table
LABEL:table_of_airways. [Milliron2015] From our dataset, many of the airways affected by bronchiectasis came from two patients. We used the same airways for the simulated low dose and voxel size experiments. A subset of the same airways were used for CT reconstruction kernel and bifurcation experiments.3.1 Simulated Images
In this experiment, we simulated images with differing radiation dose and voxel size. The purpose was to analyse 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 analysed the sensitivity of the tapering measurement against the given parameters. Finally, we investigated the reproducibility of crosssectional 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 Frush et al.[Frush2002]
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
; with a range of . 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 on Algorithm 4. For our experiment we varied from 0.5 to 5 in increments of 0.5. An example of the output image are displayed in Figure 2.To relate to the physical dose from a CT scanner, we adopted the method described in Reeves et al.[Reeves2017] 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[Gonzalez2004] 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 . Table LABEL:table_of_noise, shows values of on a selection of images against a range of . Using results from Reeves et al.[Reeves2017] and Sui et al.[Sui2016], a low dose scan with a tube currenttime product 25mAs has maximum of 55HU. Thus, we assume approximately corresponds to a low dose scan. We considered higher values of to verify any correlations in the results.
We computed the taper measurement on the noisy images using the same segmented airways and labelled distal point that were identified on the respective original image. The literature has shown in low dose scans, airway segmentation software[Rikxoort2009, Wiemker2006] 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[Rikxoort2009] and ultralow[Wiemker2006] dose scans. 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[Nagatani2015, Larbi2018, Sui2016].
3.1.2 Voxel Size
We analysed the effect of voxel sizes on the tapering measurement. For each CT image, the voxel spacing , was subsampled to new spacing of , where is a scalar constant. The intensities at each new voxel position was computed using sinc interpolation with a small amount of smoothing. We chosen Sinc interpolation to preserve as much information as possible from the original image. To compute the tapering value, we resampled the segmented airway and distal point to the same coordinate system using nearest neighbour interpolation. Morphological filtering via a closing operation[Gonzalez2004] was used on segmented airways to remove artefacts caused by the resampling. For our experiment we used the parameters; with increments of .
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 on Figure 3. We acquired the airway segmentation and distal point from a single reconstruction kernel as described in Table 1. 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 LABEL:table_of_airways.
Patients  Reconstruction Kernel used for prepocessing 
bx503  Lung 
bx507  Body 
bx513  Lung 
bx515  Body 
3.3 Biological Factors
3.3.1 Effect of Bifurcations
We analysed 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 on Figure 16. Using ITKsnap, the author (KQ) started at the cross sectional plane corresponding to the carina and scrolled towards the distal point. Using visual inspection, the following protocol was developed to identify bifurcations on cross sectional planes:

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

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.

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 is the lumen wall of both airways are completely visible and separate. The entire protocol is summarised in Figure 4.
For our experiment, we selected 19 airways from Table LABEL:table_of_airways. 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. Secondly, pairs of airways that were healthy on baseline scans and became bronchiectatic on follow up scans.
For pairs of healthy airways, a trained radiologist (JJ), manually identified 14 pairs of airways across 3 patients. . The criteria were the airway track 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 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 follow up scan, an example are displayed on Figure 5. Details of the CT images are summarised on 2 and 3.
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 on Table LABEL:table_of_airways. LABEL:table_of_airways. The tapering measurements were taken from the method discussed in Section 2.
Patients  Time between scans  Airways Labelled 

bx500  9M 6D  6 
bx504  35M 6D  7 
bx510  5M 22D  1 
P1  10M 7D  5 
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 Results
Figure 6 compares the tapering measurement between healthy and diseased airways. On a Wilcoxon Rank Sum Test between the populations, .
4.1 Dose
We analysed the difference in crosssectional area measurements and the final tapering measurements at different CT radiation doses.
For the crosssectional areas, Figure 7, compares the crosssectional areas between the original image and one of the noisy images. Each graph contains approximately 30000 unique lumen measurements. The correlation coefficients between the populations was on all graphs. The 95% confidence intervals increase with the amount of noise. For the tapering measurement, Figure 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 on all values of .
We analysed 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.
4.2 Voxel Size
We analysed 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 the correlation coefficients between measurements was . Furthermore, we analysed the error difference in arclength. On Figure 11, the mean difference shows a weak correlation coefficient with with scale . 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 .
In terms of the tapering measurement, Figure 12 compares the tapering values from the scaled images with the respective originals. The correlation coefficients between the scaled and original tapering values was on all scales . 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 . Furthermore, Figure 13, shows a positive correlation with uncertainty and scale with .
4.3 CT Reconstruction
We analysed 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 . The largest 95% confidence was in patient bx515 with mm from the mean. Figure 15, compares the differences in tapering measurement. We collected tapering measurement from 4 patients. The correlation coefficient was between the reconstruction kernel.
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 Figure 16. We compared the measurements on Figure 17, the correlation coefficient was .
The uncertainty of each tapering measurement was computed using the standard error of estimate
, defined as[Spiegal1998]:(6) 
where , is the arclength , is the estimate from the linear regression from each computed area and 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, .
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[Streiner2015] . The standard deviation of the tapering difference was mm.
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 on Figure 18. The results shows bronchiectatic airways have a greater tapering change in magnitude compared to airways that remained healthy, .
5 Discussions
In this paper, we propose a tapering measurement for airways imaged using CT and validate the reproducibility of the measurement. The tapering measurement is the exponential decay constant between crosssectional area and arclength from the carina to the distal point of the airway. Unlike other proposed tapering measurements, we assess reproducibility of the tapering measurement against simulated CT dose, voxel size and CT reconstruction kernel. Finally, we assess the effect of tapering across airway bifurcations, and examine repeatability over time using longitudinal scans.
Part of the evaluations consist of analysing the difference in tapering across longitudinal scans. The timescales between scans ranges from 9 to 35 months. The motivation for a long timescale is a proof of principle demonstration that the tapering measurement is reproducible for clinical studies. Examples include, drug trails[DeSoyza2018] and investigations in exacerbations[Chalmers2018a], where the timescales in monitoring patients were 12 months and 60 months respectively.
The pipeline consists of various established image processing algorithms. We chose the centreline algorithm developed by Palagyi et al.[Palagyi2006]. Unlike other proposed methods[Aylward2002, Jin2016, Cardenes2010] the algorithm explicitly links the distal points to the carina. Furthermore, it has been shown that the algorithm of Palagyi et al.[Palagyi2006] can be used on images with nonisotropic voxel sizes[Irving2014]. By modelling the centreline as a graphical model similar to Mori et al.[Mori2000], we performed a breadth first search[Cormen2009] to avoid analyses of false airway branches. The removal of false branches is not a trivial task[Mori2000, Irving2014, Kiraly2004].
We corrected the centreline discretisation error or recentring by smoothing points on the centreline. Smoothing has been an established method in the literature[Irving2014, Xu2015]. A recentring method was proposed by Kiraly et al.[Kiraly2004] which shifts the centreline voxels in relation to a distance transform. The process is iterative compared to a single computation of smoothing.
For our pipeline, we generated the orthonormal plane based on the method of Shirley and Marschner[Shirley2009]. We set the pixel size isotopically at 0.3mm to insure that plane image to be within the resolution of the CT image and to allow the ray casting algorithm to find the lumen at subvoxel precision. Other methods have been proposed. In Kreyszig[Kreyszig1964], they generated a binormal and principle normal. However, the method is not robust as the binormal vector can become a zero vector. Grelard et al.[Grelard2015] used Voronoi cells, a method that requires two parameters whereas Shirley and Marschner[Shirley2009] is parameter free. For our work, intensities on the crosssectional plane were computed via cubic interpolation. Various papers have used linear interpolation[Odry2008b, Kiraly2005, Tschirren2005]. However, it has been shown by Moses et al.[Moses2018] that the method can create high frequency artefacts in the image[Thevenaz2000].
Various methods have been proposed to measure the area of the airway lumen[Gu2013a, Kiraly2007, Saba2003]. We used the FWHM_{ESL} because of two distinct advantages. First, the method is parameter free. Second, the method is robust against slight variations in intensities. The method can therefore be applied to images from different scanners and images acquired using different image reconstruction kernels.
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 analyse 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, timeconsuming 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 labelling, 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.[Charbonnier2017]
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
[Moses2018, Bian2018, Yang2018].In this paper, we analyse the reproducibility of all computerized components of the tapering algorithm. The paper does not address reproducibility of manual labelling of the airways. It is noted in the literature that semimanual labelling of small airways can take hours[Tschirren2009]. Future work is required to analyse the reproducibility of manual segmentation of the airways. We hypothesise, that the segmented healthy peripheral airways consist of a small number of voxels, therefore any errors in voxel labelling will be considerable 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[Zabic2013]. 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[WonKim2014, Takenaga2016]. Thus, various groups have proposed low dose simulations using reconstructed CT images. The methods involve adding Gaussian[Takenaga2016, WonKim2014] or a combination of Gaussian and Poisson noise[Naziroglu2017] on the sinogram of the forward projection of the CT image. Whilst there has been limited validation of the appearance of lung nodules against simulated low dose simulation[Li2009], 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 is similar to results taken from low dose scans from Reeves et al.[Reeves2017] and Sui et al.[Sui2016]
Similarly, with voxel size simulation, ideally one would reconstruct the images from the original sinogram, for example in Achenbach et al[Achenbach2009]. However, as the sinograms were unavailable, we simulated the voxel size through interpolation of the original CT images similar to Robins et al[Robins2018]. We believe the simulation is sufficient as it shows the robustness and precision of the centreline, recentring and crosssectional 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 centreline computation, are repeatable across voxel sizes, we avoid resampling the image to isotropic lengths. Thus, potentially avoiding a computationally expensive[Wiemker2006] preprocessing step.
We showed the tapering measurement is reproducible by measuring the same airway across longitudinal scans with a minimum 5 month interval. The time between scans were on a similar scale from a reproducibility study on airway lumen by Brown et al.[Brown2017]. 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.[Hammond2017]. 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[Dendy1999]. For our experiment, each airway was chosen by a subspecialist thoracic radiologist. The airway was inspected to ensure it was in a healthy state, for example, with no mucus present. Thus, we assume that each pair of airways is disease free and healthy.
6 Conclusions
In this paper, we show a statistical difference in tapering between healthy airways and those affected by bronchiectasis as judged by an experienced radiologist. From Figure 6, the difference between the mean and median of the two populations was 0.011mm and 0.006mm respectively. In simulated low dose scans, the tapering measurement retained a 95% confidence interval of mm up to , equivalent to a 25mAs low dose scan. In simulations assessing different voxel sizes, the tapering measurement retained a 95% confidence between mm up to . The tapering measurement retains the same 95% confidence, mm interval against variations in CT reconstruction kernels, bifurcations and, importantly, over time in evaluating sequential scans in normal airways. Importantly, we showed as a proof of principle that the magnitude change in tapering for healthy airways is smaller than those from airways that became bronchiectatic. From our previous work[Quan2018], we showed the measurements are accurate to a sub voxel level. Our findings suggest that our airway tapering measure can be used to assist in the diagnosis of bronchiectasis, to assess the progression of bronchiectasis with time and, potentially, to assess responses to therapy.
We analysed the reproducibility of the components that constitute the tapering measurements. The reproducibility of area measurements was analysed in relation to simulated radiation dose and CT reconstruction kernels. For simulated dose, we found the 95% confidence interval retains mm in noisy images under , equivalent to a dose just higher than a 25mAs low dose scan. We note on Figure 7, there is a bias towards overestimating larger lumen sizes at lower doses. As the centreline length remains constant and bias on the smaller lumen remain stable, the overestimation results in an increase in taper magnitude. For reconstruction kernels variation, we found the largest 95% confidence interval was mm. The reproducibly of arclengths was tested against voxel sizes variability and showed that arclengths have a 95% confidence interval of up to mm for scales under . The increase in the standard deviation of arclength and area against voxel size and dose respectively correlate with uncertainty in tapering.
This paper provides useful information for clinical practice and clinical trials. An accurate prediction of the noise amplitude in a particular CT scan and its distribution is a function of the limited radiation dose of the scan, scanner geometry, reconstructed voxel size, other sources of noise, the reconstruction algorithm and any pre and postprocessing used. Many of these factors are proprietary information of the CT manufacturer and hence not available to users.[Stoel2008, Parr2004] We have undertaken an experiment to assess the dependence of our measurements on a simulated noise field added to the CT scan data and have presented the results. This gives an indication of the dependence on radiation dose assuming all other factors remain the same. We recommend that the accuracy experiment presented in this paper be repeated for the particular reconstruction, scan protocol and scanner type used to make the measurements.
Bronchiectasis is often described as an orphan disease and has suffered a lack of interest and funding[Hurst2014, Chalmers2014]. We have shown that the reproducibility of automated airway tapering measurements can assist in the diagnosis and management of bronchiectasis. In addition, we show it is feasible to use our tapering measurement in large scale clinical studies of the disease provided careful phantom calibration is taken.
Disclosures
Part of the work have been presented at the 2018 SPIE Medical Imaging conference[Quan2018].
For potential conflicts of interest: Ryutaro Tanno has been employed by Microsoft, ThinkSono and Butterfly Network (the employment is unrelated to the submitted work), Joseph Jacob has received fees from Boehringer Ingelheim and Roche (unrelated to the submitted work) and David Hawkes is a Founder Shareholder in Ixico plc (unrelated to the submitted work). The other authors have no conflicts of interest to declare.