Reproducibility of an airway tapering measurement in CT with application to bronchiectasis

09/16/2019 ∙ by Kin Quan, et al. ∙ UCL 7

Purpose: This paper proposes a pipeline to acquire a scalar tapering measurement from the carina to the most distal point of an individual airway visible on CT. We show the applicability of using tapering measurements on clinically acquired data by quantifying the reproducibility of the tapering measure. Methods: We generate a spline from the centreline of an airway to measure the area and arclength at contiguous intervals. The tapering measurement is the gradient of the linear regression between area in log space and arclength. The reproducibility of the measure was assessed by analysing different radiation doses, voxel sizes and reconstruction kernel on single timepoint and longitudinal CT scans and by evaluating the effct of airway bifurcations. Results: Using 74 airways from 10 CT scans, we show a statistical difference, p = 3.4 × 10^-4 in tapering between healthy airways (n = 35) and those affected by bronchiectasis (n = 39). The difference between the mean of the two populations was 0.011mm^-1 and the difference between the medians of the two populations was 0.006mm^-1. The tapering measurement retained a 95% confidence interval of ±0.005mm^-1 in a simulated 25 mAs scan and retained a 95 up to 1.5 times the original voxel size. Conclusion: We have established an estimate of the precision of the tapering measurement and estimated the effect on precision of simulated voxel size and CT scan dose. We recommend that the scanner calibration be undertaken with the phantoms as described, on the specific CT scanner, radiation dose and reconstruction algorithm that is to be used in any quantitative studies. Our code is available at https://github.com/quan14/AirwayTaperingInCT

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 13

page 16

page 17

page 18

page 30

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 broncho-arterial (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, Perez-Rovira2016, 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 cross-sectional 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 cross-sectional 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 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[Hammond2017, Leutz-Schmidt2017, 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, Leutz-Schmidt2017] 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 ITK-snap111http://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.

Figure 1: Summary of steps in our pipeline

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.

Data: Distance image on the th axial slice
Result: Start point of trachea
First slice at the top of the trachea.
while  do
      
end while
Algorithm 1 Locating the start of centreline on the trachea

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 re-centring[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 cscvn222https://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 sub-voxel 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 cross-sectional area accurately by constructing a cross-sectional 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.

Data: Unit tangent vector of the spline
Result: Basis of the orthogonal plane
Arbitrary vector such that and are not collinear
Algorithm 2 Constructing the basis for the plane reconstruction, adapted from Shirley and Marschner [Shirley2009].

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 Edge-Cued Segmentation-Limited Full Width Half Maximum (FWHMESL), developed by Kiraly et al.[Kiraly2005]

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

Data: The rays: , where is the length from the centre to the border of the plane.
Result: The position of the lumen edge, .
The first index of the ray such that
Local maximum intensity in nearest to s
The index such that
Minimum intensity in from to
The index such that
The index such that and
Algorithm 3 Summary of the FWHMESL, adapted from Kiraly et al.[Kiraly2005] The purpose of the algorithm was to find the point of the ray which crossed the lumen.

The final output of the FWHMESL 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 cross-sectional 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 cross-sectional area as

(5)

where is the arc length of the spline, is the proposed tapering measurement, is the cross-sectional area and is an arbitrary constant.

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 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 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.63-0.80mm in plane and 0.80-1.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 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 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.

function noisySlice = AddingDoseNoise(axialSlice,lambda)
%Creating the sinogram
sinogram = radon(axialSlice,0:0.1:179); %Adding the Gassian 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
Algorithm 4 Adapted Matlab code to simulate noise from differing doses
Figure 2: TOP: CT images with simulated noise against varying . BOTTOM: An image subtraction of the simulated noisy image with the original. The units are in HU.

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 current-time 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 ultra-low[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
Table 1: The images used for the reconstruction kernel experiment. The table lists which reconstruction kernel was used to generate the airways segmentation and distal point labelling. The make, model and voxels size of the images are displayed in Table LABEL:table_of_airways.
Figure 3: Images from the same CT scan with the body kernel (LEFT) and lung kernel (RIGHT). Both images are displayed in the same intensity window.

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 ITK-snap, 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:

Figure 4: FAR LEFT: 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 labelled region consists of slices from green to red. CENTRE LEFT: A cross sectional plane where the airway is at the point of bifurcation, indicated by the blue arrow. CENTRE RIGHT: First slice of the bifurcation region. FAR RIGHT: The final slice of the bifurcation region. The slides are chronologically ordered with the protocol described in Section 3.3.1
  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 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.

Figure 5: The same pair of airways from longitudinal scans. LEFT: Initial healthy airway. RIGHT: 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 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
Table 2: List of the images for progression experiment. The table includes time between scans in months (M) and days (D). The airways on this table are different to Table LABEL:table_of_airways
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
Table 3: 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. Abbreviation: GEMS - GE Medical Systems, SS - Siemens SOMATOM.

4 Results

Figure 6 compares the tapering measurement between healthy and diseased airways. On a Wilcoxon Rank Sum Test between the populations, .

Figure 6: Comparing the proposed tapering measurement with labelled healthy and bronchiectatic airways. On a Wilcoxon Rank Sum Test between the populations, .

4.1 Dose

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

For the cross-sectional areas, Figure 7, compares the cross-sectional 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 .

Figure 7: A series of Bland-Altman[Bland1986] graphs comparing area measurements from a simulated low dose scan and the original image. On all graphs, the correlation coefficient was
Figure 8: A series of Bland-Altman[Bland1986] graphs comparing tapering measurement between simulated dose and the original image. On all graphs the correlation coefficient was .

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.

Figure 9: Mean (LEFT) and standard deviation (RIGHT) of the difference in tapering between original images minus the simulated lower 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 .

Figure 10: A series of Bland-Altman[Bland1986] graphs comparing arc lengths between scaled images and the original images. On all graphs, the correlation coefficient was
Figure 11: Mean (LEFT) and standard deviation (RIGHT) of the difference in arclength between original images minus the scaled images. The correlation coefficient of the mean and standard deviation against scale are and respectively.

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 .

Figure 12: A series of Bland-Altman[Bland1986] graphs comparing tapering between original images and scaled images. On all graphs the correlation coefficient was
Figure 13: Mean (LEFT) and standard deviation (RIGHT) of the difference in tapering between original images minus the scaled images. The correlation coefficient of the mean and standard deviation against scale are and respectively.

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.

Figure 14: Bland-Altman[Bland1986] graphs comparing the cross-sectional area between the Lung and Body kernels. On all four images the correlation coefficient was
Figure 15: Bland-Altman[Bland1986] graph comparing tapering measurements between Lung and Body kernels, .

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 .

Figure 16: TOP LEFT: A signal of area measurement with bifurcation regions (red) and tubular regions (blue). TOP RIGHT: 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. BOTTOM: The reconstructed bronchiectatic airway of the same profile. The blue-shaded and red-shaded regions corresponds to the tubular and bifurcating airways respectively. A reconstructed healthy airway have been discussed in Quan et al.[Quan2018] Similar reconstructed cross sectional images of vessels have been discussed in Oguma et al.[Oguma2015], Kumar et al.[Kumar2015], Alverez et al.[Alvarez2017]. and Kirby et al.[Kirby2018]

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

Figure 17: LEFT: Bland-Altman[Bland1986] graph showing the relationship of the taper rates with and without bifurcations, . RIGHT: Comparison of the standard error from linear regression between airways with and without bifurcations. On a Wilcoxon Rank Sum Test between the two populations, .

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

Figure 18: RIGHT: Bland-Altman[Bland1986] graph comparing tapering measurement on the same airways from the first and subsequent scan, . LEFT: The change in tapering between healthy airways (H to H) and airways that became bronchiectatic on the follow up scans (H to D), .

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 cross-sectional 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 non-isotropic 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 sub-voxel 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 cross-sectional 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 FWHMESL 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, 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 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 semi-manual 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 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 centreline computation, are repeatable across voxel sizes, we avoid resampling the image to isotropic lengths. Thus, potentially avoiding a computationally expensive[Wiemker2006] pre-processing 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 post-processing 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.

Acknowledgements.
Kin Quan would like to thank Prof Simon Arridge and Dr Andreas Hauptmann for their helpful conversations on the low dose simulations. This work is supported by the EPSRC-funded UCL Centre for Doctoral Training in Medical Imaging (EP/L016478/1) and the Department of Health NIHR-funded Biomedical Research Centre at University College London Hospitals. Ryutaro Tanno is supported by Microsoft Research Scholarship. Joseph Jacob is a recipient of Wellcome Trust Clinical Research Career Development Fellowship 209553/Z/17/Z

References