1 Introduction
In fibrosing lung disease, contraction of the lung interstitium pulls on airway walls, and causes a dilatation. Airway dilatation is calculated using crude visual scores have been shown to be powerful predictors of outcome in idiopathic pulmonary fibrosis (IPF) [8]. Computed Tomography (CT) coupled with image analysis algorithms to measure the airway lumen can provide an excellent tool to quantitatively assess these dilatation. CT gives a very high image contrast between air within the airway lumen and the surrounding airway wall. However, airway lumen measuring algorithms are subjected to sources of noise. Variations in voxel sizes [1] and CT reconstruction algorithms [6] have been shown to affect the precision of the crosssectional area measurements of the airways. In addition, there can be variations in airway crosssectional area due to normal biological processes such as points of airway bifurcation where the airway transiently dilates [12, 3]. These sources of noise makes postprocessing of measurements difficult to locate dilatation on an airway track. Identifying the point of dilatation and subsequently calculating the volume of abnormal airways in the lung potentially provides a sensitive means for quantifying subtle worsening of disease in IPF across longitudinal scans.
In this work, we introduce a variant of Bayesian Changepoint Detection (BCPD) models to automatically identify the location of airway dilatation in the presence of strong measurement noise, given two measurements of the same airway acquired at different time points. BCPD model, typically used in DNA sequencing [13] and stock market data analysis [10], aims to capture abrupt variations in the underlying distributions of a given signal or a time series. The method processes a series of airway crosssectional area changes between baseline (first) and followup (second) CT scans, and generates the posterior distribution over multiple possible points of abnormal variations, whose the mode is taken as the final prediction. We test the efficacy of the method on (1) CT images of healthy airways with simulated dilatation, and (2) pairs of real images of IPFaffected airways acquired approximately 1 year intervals. For the simulated dataset, we measured the detection accuracy with respect to the commonly used naive thresholding and maximum likelihood based methods [9]. For the longitudinal IPF dataset, we compare the predictions of our model to measurements from radiologists based on two different protocols.
2 Method
In this section, we introduce our method for quantifying progression of IPF on a patients serial CT scans. The method proceeds as follows. Firstly, we fit a tubular shape model to the airways in both baseline and followup CT scans, and acquire estimates of the crosssectional areas along them (Sec.
2.1). We then treat the difference in the crosssectional areas as a 1D signal, and employ the proposed Bayesian changepoint model to estimate the posterior distribution over locations of abrupt airway dilatation (Sec. 2.2 and 2.3). Lastly, we postprocess this posterior distribution to determine the region of dilatation (Sec. 2.4).2.1 Airway Preprocessing
In this work, for each airway track, we acquire a series of crosssectional area measurements using the method proposed by Quan et al [12] (Fig. 1). Following a semimanual segmentation of the airway, the method computes the airway centreline and the corresponding normal planes at contiguous intervals, each of which is then used to estimate the crosssectional area. The final output is a 1D function of crosssectional area along the airways for baseline and the followup scans.
The next task was to align the signals on both baseline and followup
scans. To this end, we resample the signal to 1mm using cubic interpolation. We considered the first 50 points on both signals
from the start of the carina. To register the two airways together, we apply the transformation where(1) 
The longest of the two signal were truncated from the right hand side such that both signals were of the same length. For the rest of the methodology, we will only consider the series difference defined as .
2.2 Bayesian Changepoint Model
The progression of fibrosing lung diseases can manifest itself as dilatation of the internal airways. Thus, we hypothesise that at the start of dilatation, the series undergoes an abrupt variation, which we refer to as a changepoint. More formally, given signal of length , we define a changepoint as the location where there exists a change in parameters in the underlying distribution . In other words, at changepoint , the observations can be separated at such that:
(2) 
where . This definition can be naturally extended to the scenario with
changepoints; we denote the changepoint location vector by
, with parameters for each respective segment. For ease of notation, we also denote and . Assuming statistical independence between segments, the likelihood factorises as:(3) 
where We also specify prior distributions on the the number of changepoints , the locations of the changepoints , and the parameters of the corresponding segments where , and represent the hyperparameters.
Given the likelihood and the prior distributions above, we would like to estimate the posterior distributions over the number of changepoints , the locations of changepoints and the parameters of the respective segments
. Commonly in Bayesian changepoint analysis, likelihoods with a conjugate prior are chosen to allow calculation of analytical posteriors for changepoints and parameters, enabling faster convergence and better accuracy
[4]. However, in our work, we wish to relax this conjugacy assumption in order to allow flexibility in the choices of prior and likelihoods. In the next section, we describe the inference scheme in such highly general setting, viable for arbitrary forms of the likelihood distribution and the prior. In this paper, was defined as the Student t distribution with parameterswith degrees of freedom
, meanand variance
. We chose tdistribution for its robustness to outliers
[11] caused by noise from the area measurements.2.3 Reversible Jump MCMC for Posterior Inference
Posterior inference with our model possesses two challenges. Firstly, without the conjugacy assumption, computing the posterior distributions is intractable. Secondly, the dimensionality of the posterior distribution over the changepoints is given by and varies during inference. To combat the first problem, we use the MetropolisHasting (MH) algorithm [2], a variant of Markov Chain Monte Carlo (MCMC) methods that can sample from the posterior, with or without conjugacy. Given that the number of changepoints is known, MH can be used to sample from the posterior distributions over the changepoints and segment parameters . To address the second problem of varying posterior dimensionality , we extend the above sampling scheme to the Reversible Jump MCMC framework [5]. Taken all together, the method is capable of traversing the full posterior distributions for and we refer to this as Reversible Jump Metropolis Hasting (RJMH) algorithm.
Overview of RJMH:
The RJMH proceeds by randomly executing one of four possible moves at each iteration (see Algorithm 1 for the pseudocode):

The parameters of every segment is updated.

One of the changepoints is selected randomly and updated, and the parameters of its neighbouring two segments are updated.

Randomly choose a location in the time series and add it as a changepoint. The parameters of the resultant two new segments are also sampled.

One of the changepoints is randomly removed and the parameters of the new segment are sampled
We denote these moves as: respectively. We also define the maximum number of changepoints and at the boundary cases for , we impose restrictions such that is not possible when and is not possible when . Each move updates the appropriate subset of parameters by sampling from the corresponding proposal distributions and , and is only executed if it passes the associated acceptance criteria .
We can subdivide these four moves into two sets: fixed (MH Steps) and changing (RJ Steps) i.e. resampling parameters and moving changepoints does not affect the changepoint model whereas adding or deleting a changepoint changes the value of .
MH Steps:
For , we set N. This step resamples parameters of each segment by proposing Gaussian perturbations around the current values of parameters for all segments.
For , we set Poi where Binary. This step selects a changepoint at random and shifts it with a Poisson perturbation. The segments neighbouring this new changepoint location have parameters resampled as in move using the current segment parameters.
Note that each proposal in the MH steps will be symmetrical and therefore, cancel out with the reverse proposal term (See Algorithm 1).
RJ Steps:
For the changing moves, the algorithm moves in dimensionality and we need to introduce alternative proposal terms for the new segments produced as there are no current segment parameters to base our proposals on.
For i.e. , we proposed random new changepoints over our data, U. The proposed split an existing segment into a new left segment and new right segment . Our proposal for are defined by a Gaussian perturbation on empirical values of the respective segments (Fig. 2). The proposal for , is Gaussian perturbation of the previous update . Due to dependence of the proposal, a Jacobian term is introduced .
For i.e. , we remove a changepoint . As before, proposals for , are defined using empirical values of the segments and Gaussian perturbation (Fig. 2). The proposal for is the mean of the previous . The move introduces Jacobian term .
Initialization:
The overall algorithm structure can be seen in the pseudocode Algorithm 1. With the model defined, we mentioned the priors used and address some common implementation techniques. The priors used are given as follows:
(4)  
(5)  
(6)  
(7) 
The hyperparameters for were chosen to be noninformative and within plausible ranges. In order to reduce the number of changepoints we could detect (in the case of acquisition noise from the CT), we set the expectation for to be sufficiently low. In terms of implementation, we follow a standard procedure of setting a burnin for the number of iterations to ignore any potential issues with intialisation at approximately 25% of the total iteration count. To remove any chance of autocorrelation, we thin the number of samples by only storing the iteration, after the burnin period. Additionally, we choose the number of iterations, , to be sufficiently large in order to achieve convergence.
2.4 Locating airway dilatation
For airways affected by IPF, dilatation starts at the distal point and progresses in the proximal direction [7]
. Therefore, we topologically can assume that each affected airway undergoes a single changepoint from which dilatation starts. To locate such unique changepoint, we consider the posterior probability of the changepoint
, and perform the following postprocessing steps.On each airway track, the proximal region is surrounded by cartilage [14]. As the airway track loses cartilage support, the geometry of the airway changes and results in a changepoint. Since such biological changepoint is independent of the disease state and occurs prior to dilatation, we eliminate it by discounting the most proximal peak in the posterior distribution . We then selected highest peak on the modified posterior as the final estimation for the starting point of dilatation.
3 Evaluation & Results
We evaluated our proposed method with two experiments. Firstly, we used a set of healthy airways, in which we had simulated dilatation to quantify the accuracy of our method and compared the results with conventional changepoint methodologies. Secondly, we assessed the clinical utility of our method on a dataset of airways affected by IPF. We compared our labels with those from two experienced thoracic radiologists. The properties of the images used in both experiments are displayed on Tab. 1.
Experiment  Patients  BL Voxel Size  FU Voxel Size  Airways  Time between scans 
1  1  0.67, 0.67, 1.00  0.56, 0.56, 1.00  6  9M 6D 
1  2  0.63, 0.63, 1.00  0.78, 0.78, 1.00  7  35M 6D 
1  3  0.72, 0.72, 1.00  0.72, 0.72, 1.00  1  5M 22D 
2  1  0.72, 0.72, 1.00  0.64, 0.64, 1.00  3  12M 5D 
2  2  0.67, 0.67, 1.00  0.87, 0.87, 1.00  1  10M 24D 
3.1 Disease simulation
To quantitatively assess accuracy, a ground truth is required. To this end, we applied our point detection algorithm on augmented healthy airway series to simulate the airway dilatation caused by IPF. After obtaining written informed consent from 3 patients, a trained radiologist (R1) selected 14 pairs of healthy airways in both baseline and followup scans. The image properties are displayed on Tab. 1. They were acquired from different scanners and used different reconstruction kernels. The airways were preprocessed as described in Sec. 2.1 to produce a function of area change along the length of the airway. We interpret this output signal as the overall noise caused from normal biological changes and acquisition.
We modelled the change in dilatation with a logistic function;
(8) 
We interpreted as the magnitude of dilatation and as the point of dilatation. The parameters are set such that the dilatation starts 1040mm from the distal point in increments of 5mm. In addition, we set to range from 0.33 in increments of 0.25. Finally we set mm, in order to create an abnormal increase in cross sectional area. To simulate the dilatation on the airway; the logistic function was added to the area change of the healthy airway, as shown on Fig. 3. We applied every permutation of and on each of the 14 healthy airways. After the signal augmentation, we applied our proposed Bayesian changepoint detection algorithm.
The proposed method was compared against two conventional methods. First, we use a basic thresholding method. For a given signal, we applied a 5mm moving mean average and thresholded the point at which the signal reached above the upper quartile from the right hand side. Secondly, we implemented the method based on Lavielle
[9], in summary we consider changepoint and these changepoints , minimize the function:(9) 
where is modified such that the function finds less than changepoints. For our paper, different were evaluated and we found
(10) 
gives the most accurate results. To replicate the post processing of our proposed method, we consider possible changepoints. This takes into account the changepoint caused by the support cartilage. Finally, we set a minimum distance of 20mm. Once we acquired the changepoints, the most peripheral point was chosen as the point of dilatation. The implementation was performed through Matlab inbuilt function; findchangepts^{1}^{1}1https://www.mathworks.com/help/signal/ref/findchangepts.html last accessed on August 24, 2019..
Fig. 4 shows the accuracy for each individual method as a heatmap. The metric used to quantify accuracy was the displacement between the ground truth and the changepoint given by the proposed method. A positive displacement (mm) corresponds to an overestimation of the ground truth towards the distal point. Each entry on the heatmap corresponds to the median average of all displacements over 14 airway pairs.
When the magnitude of dilation is larger than , our proposed method achieves consistently higher accuracy than Lavielle [9]. We note that the accuracy gain in the peripheral regions of the airways at 1030mm from the distal point are the most clinically relevant in IPF as lung parenchymal damage begins in the lung periphery and progresses proximally [7]. Furthermore on the same peripheral regions 1030mm, the baseline method showed systematic bias in accuracy towards the central airways. This was due to the baseline method being influenced by outliers from the longer expanses of normal airway regions. The proposed method uses the tdistribution as the likelihood thus making it robust to possible outliers within the data [11]. On the other hand, the proposed method suffers from poor accuracy below magnitudes of dilatation . However, in physical terms a dilatation of corresponds to a percentage increase in cross sectional area of . This is within the range of normal biological change of the airways [3].
3.2 Application to Airways affected by IPF
We applied our method to airways affected by IPF. The purpose was to compare our measurement to the labels provided by a radiologist and compute the volume change of the identified diseased airway regions. For our dataset, we acquired 4 airway pairs from 2 patients after obtained a waiver for consent from the local Research Ethics Committee. All airways were judged by the radiologist R1 to be dilated as a consequence of IPF on baseline and to have visually worsened on followup imaging. Image properties are displayed in Tab. 1.
We compared the performance of our method against two trained thoracic subspecialist radiologists. Two radiologists R1, R2 identified the point of at which a given airway was seen to demonstrate increased dilatation on the followup CT scan. To assessed the reproducibility of manual labeling, each radiologist labelled the same airway twice through two different protocols. In the first method, the radiologists interrogated axial CT images. Using 2 separate workstations and the airway centreline, the radiologists identified the point on the centreline (on the followup scan) where the airway demonstrated definitive worsened dilatation. For the second method, the radiologist compared the aligned reconstructed crosssectional planes on baseline and followup scans. The radiologist then selected the slice where the airway had worsened when evaluated against the baseline scan. An example of both protocols are displayed on Fig. 5.
Fig. 6 compares the predictions of our method with the labels from radiologists obtained in two different protocols. The results indicate that the predictions for Airway 1, 3 and 4 are within the range of the radiologists’ labels. In the case of Airway 2, although our method based on the maximum peak overestimates with respect to the radiologists’ predictions, the posterior distribution contains another equally probable peak that underestimates the radiologists’ labels (see the second highest peak at 70mm), potentially indicating a more proximal point of dilatation. To test this, we delineated the boundary of the lumen on the reconstructed cross sectional slices at baseline in the neighbouhood of this peak, and Fig. 7 shows the initial few slices (6264mm). We observed that, when the delineated boundary from baseline was superimposed on the followup scan, the boundary is contained inside of the followup lumen with several lumen pixels are consistently outside from the boundary. Thus, this result indicates that the starting point of dilatation is more proximal than the labels from the radiologists.
To demonstrate the clinical utility of locating the starting point of airway dilatation, we compared longitudinal airway volume changes in diseased and healthy regions of each of the airway tracks. To find the volume of the airway, we considered the aligned signals, ,. These signals are measurements of area against the airway arc length. Thus volume can be computed via the area under the curve. For our work, we used the trapezium rule to find the volume. Three volumetric regions were considered: (i) the entire airway i.e. from the carina, to the distal point, , (ii) carina to the dilatation point , denoted as . (iii) the dilatation point to the distal point, denoted by . Note that and does not overlap with . Tab. 2, shows the results of the percentage volume change. The volume change in had greater sensitivity for selecting progressive airway dilatation in IPF than the volume change in the entire airway .
Airway  PVC of  PVC of  PVC of 

1  2.6 %  32.9%  1.1% 
2  3.2 %  129.7%  2.2% 
3  2.6 %  47.4 %  0.3% 
4  7.4 %  48.4 %  7.1% 
4 Discussion & Conclusion
In this paper, we modelled changes in area along airway tracks as a time series with the purpose of detecting abnormal dilatation caused by IPF using Bayesian changepoint detection. Experiments on simulated data show that our model is able to detect the starting location of airway dilatation with superior accuracy than the relevant baseline methods. The results on the IPF longitudinal dataset display reasonable agreement with radiologists, while in one case indicating a more plausible location of dilatation, potentially missed by the experts.
Identifying changepoints and thereby calculating a change in airway dilatation over time could become a sensitive measure of IPF aggravation. This would form an important secondary endpoint in drug trials, where our measurements could indicate whether a new drug ameliorates disease progression better than existing medications. In the future, we hope to evaluate such utility of our proposed method on a larger cohort of IPF subjects. Furthermore, our method is applicable to other lung airway diseases, characterised by dilatation (such as cystic fibrosis) or other geometrical deformations, and such extensions remain valuable future work.
5 Acknowledgements
Kin Quan is supported by the EPSRCfunded UCL Centre for Doctoral Training in Medical Imaging (EP/L016478/1) and the Department of Health NIHRfunded Biomedical Research Centre at University College London Hospitals. Ryutaro Tanno is supported by Microsoft Research Scholarship. Micheal Duong is funded by the Centre for Doctoral Training in Financial Computing and Analytics. Arjun Nair is funded by National Institute for Health Research Biomedical Research Centre. Mark Jones and Christopher Brereton are supported by National Institute for Health Research. Joseph Jacob is a recipient of Wellcome Trust Clinical Research Career Development Fellowship 209553/Z/17/Z.
References
 [1] T Achenbach et al. Influence of pixel size on quantification of airway wall thickness in computed tomography. Journal of Computer Assisted Tomography, 33(5):725–730, 2009.
 [2] S Chib et al. Understanding the MetropolisHastings Algorithm. The American Statistician, 49(4):327–335, 1995.
 [3] L Gazourian et al. Quantitative computed tomography assessment of bronchiolitis obliterans syndrome after lung transplantation. Clinical Transplantation, 31(5):e12943, 2017.
 [4] A Gelman et al. Bayesian Data Analysis, Third Edition, volume 1. CRC Press, 2014.
 [5] P J Green et al. Reversible jump MCMC. Genetics, 155(3):1391–1403, 2009.
 [6] S Gu et al. Computerized identification of airway wall in CT examinations using a 3D active surface evolution approach. Medical Image Analysis, 17(3):283–296, 2013.
 [7] J Jacob et al. HRCT of fibrosing lung disease. Respirology, 20(6):859–872, 2015.
 [8] J Jacob et al. Mortality prediction in idiopathic pulmonary fibrosis: evaluation of computerbased CT analysis with conventional severity measures. European Respiratory Journal, 49(1):1601011, 2017.
 [9] M Lavielle. Using penalized contrasts for the changepoint problem. Signal Processing, 85(8):1501–1510, 2005.
 [10] T Mikosch et al. Changes of structure in financial time series and the GARCH model. REVSTAT Statistical Journal, 2(1):41–73, 2004.
 [11] S J D Prince. Computer Vision: Models, Learning, and Inference. Cambridge University Press, 2012.
 [12] K Quan et al. Tapering analysis of airways with bronchiectasis. Proceedings of SPIE, 10574:105742G, 2018.
 [13] D Siegmund et al. Detecting simultaneous varient intervals in aligned sequences. The Annals of Applied Statistics, 5(2A):645–668, 2011.
 [14] E R Weibel. Morphometry of the Human Lung. Springer, 1963.