Segmentation of tree structures comprising of vessels, neurons, airways etc. are useful in extraction of clinically relevant biomarkers [1, 2]. The task of extracting trees, mainly in relation to vessel segmentation, has been studied widely using different methods. A successful class of these methods are based on techniques from target tracking. Perhaps the most used tracking strategy is to proceed from an initial seed point, make local-model fits to track individual branches in a sequential manner and perform regular branching checks [3, 4]. Such methods are prone to local anomalies and can prematurely terminate if occlusions are encountered. The method in  can overcome such problems to a certain extent using a deterministic multiple hypothesis testing approach; however, it is a semi-automatic method requiring extensive manual intervention and can be computationally expensive. In 
, vessel tracking on 2D retinal scans is performed using a Kalman filter. They propose an automatic seed point detection strategy using a matched filter. From each of these seed points vessel branches are progressively tracked using measurements that are derived from the image data. A gradient based measurement function is employed which fails in low-contrast regions of the image, which are predominantly regions with thin vessels. Another major class of tracking algorithms are based on a stochastic formulation of tracking[5, 6] using some variation of particle filtering. Particle filter-based methods are known to scale poorly with dimensions of the state space .
In spirit, we propose an exploratory method like particle filter-based methods, with a salient distinction that the proposed method can track branches from several seed points across the volume. We use linear Bayesian smoothing to estimate branch states, described using Gaussian densities. Thus, the method inherently provides an uncertainty measure, which we use to discriminate true and false positive branches. Further, unlike particle filter-based methods, the proposed method is fast, as Bayesian smoothing is implemented using the RTS (Rauch-Tung-Striebel) smoother  involving only a set of linear equations.
We formulate tracking of branches in tree structures using probabilistic state-space models, commonly used in target tracking and control theory . The proposed method takes image data as input and outputs a collection of disconnected branches that taken together forms the tree structure of interest. We first process the image data to obtain a sequence of measurements and track all possible branches individually using Bayesian smoothing. We then use covariance estimates of individual branches to output a subset of the most likely branches yielding the tree structure of interest. Details of this process are described below.
2.1 Tracking individual branches
We assume the tree structure of interest, , to be a collection of
independent random variables, where individual branches are denoted . Each branch of length is treated as a sequence of states, . These states are assumed to obey a first-order Markov assumption, i.e.,
The state vector has seven random variables,
describing a tubular segment centered at Euclidean coordinates , along an axis given by the direction vector with radius .
The observed data, image , is processed to be available as a sequence of vectors. We model the measurements as four dimensional state vectors consisting only of position and radius. This is accomplished using a multi-scale blob detector . The input image with voxels is transformed into a sequence of measurements, with position and radius information, denoted , where each . This procedure applied to the application of tracking airway trees is described in Section 2.5.1.
2.2 Process and Measurement Models
Transition from one tracking step to another within a branch is modelled using the process model. We use a process model that captures our understanding of how individual branches evolve between tracking steps and has similarities with the model used in . We assume first-order Markov independence in state transitions from (1), captured in the process model below:
where is the process model function and is the process noise. is assumed to be a zero mean Gaussian density, i.e, , with process covariance, , acting only on direction and radius components of the state vector,
where only the non-zero part of the matrix is shown and
is the process variance. The parametercan be seen as step size between tracking steps. As (3) is a recursion, the initial point (seed point), , comprising of position, scale and orientation information is provided to the model. Seed points are assumed to be described by Gaussian densities, , with mean and covariance . We present an automatic strategy to detect such initial seed points in 2.5.2.
The measurement model describes the relation between each of the 4-D measurements, in the sequence, , and the state vector, , as shown in Figure 1. A simple linear measurement model captures this relation,
where are observations generated by true states of the underlying branch at step , is the measurement function. is the measurement noise with covariance that is a diagonal matrix with entries, , which correspond to variance in the observed position and radius, respectively. All possible measurement vectors obtained from the image are aggregated into the measurement variable .
2.3 Bayesian Smoothing
The state-space models presented above enable us to estimate branches using the posterior distributions, , using standard Bayesian methods. We employ Bayesian smoothing as all the measurements are available at once, when compared to sequential observations that are more common in object tracking applications. Due to a linear, Gaussian process and measurement models, Bayesian smoothing can be optimally performed using the RTS smoother . RTS smoother uses two Bayesian filters to perform forward filtering and backward smoothing. Forward filtering is identical to performing Kalman filtering and consists of sequential prediction and update with observed information of the state variable. Once a branch is estimated using forward filtering, the saved states are used to perform backward smoothing using a Kalman-like filter which improves state estimates by incorporating additional information from future steps. Standard equations for an RTS smoother are presented below .
2.3.1 Forward Filtering
Equations in the first column of Table 1 are used to perform prediction and update steps of the forward filtering. In the prediction step, process model is used to predict states at the next step. Mean and covariance estimates of the predicted Gaussian density, i.e, of state conditioned on the previous state, denoted with subscript , are computed in (6),(7). In the update step, described in (8) – (12), predicted density is associated with a measurement vector to obtain posterior density. First, the new information from measurement is computed using (8) and is aptly called the “innovation”, denoted as . Uncertainty in the new information, innovation covariance , is computed in (9). Then, predicted mean is adjusted with weighted innovation and predicted covariance is adjusted with weighted innovation covariance to obtain the posterior mean and covariances, in (11) and (12), respectively. The weighting computed in (10), denoted as , is the Kalman gain which controls the extent of information fusion from process and measurement models.
We continue estimation of the posterior density (described by posterior mean and covariance) in a sequential manner for the branch until no new measurements exist for updating. After the final update step, a sequence of posterior mean estimates and posterior covariance estimates , obtained from the forward filter are saved, for further use by the backward smoother.
2.3.2 Backward smoothing
The smoothed estimates are obtained by running a backward filter starting from the final tracked state of the forward filter. The intuition behind backward smoothing is that the uncertainty in making predictions in the forward filtering can be alleviated using information from future steps. It is implemented using the equations in the second column of Table 1.
When performing the RTS smoother recursions, the forward filter expects a single measurement vector for the update step. We employ rectangular and ellipsoidal gating to reduce the number of measurements handled during the update step .
First, we perform simple rectangular gating which is based on excluding measurements that are outside a rectangular region around the predicted measurement in equation (8) using the following condition:
where is the covariance of the predicted measurement in equation (9). The rectangular gating coefficient, , is usually set to a value . Rectangular gating localises the number of candidate measurements relevant to the current tracking step. To further narrow down on the best candidate measurement for update, we follow rectangular gating with ellipsoidal gating . With ellipsoidal gating we accept the measurements within the ellipsoidal region of the predicated covariance, using the following rule:
where is the rectangular gating threshold, obtained from the gating probability , which is the probability of observing the measurement within the ellipsoidal gate,
2.4 Tree as a Collection of Branches
Once a branch is smoothed and saved using Bayesian smoothing described previously, we process new seed points and start tracking branches until no further seed points remain to track from. This procedure yields a collection of disconnected branches. The next task is to obtain a subset of likely branches that represent the tree structure of interest by discarding false positive branches.
2.4.1 Validation of Tracked Branches
An advantage of using Bayesian smoothing to track individual branches is that apart from estimating the branch states from the image data (using the smoothed posterior mean estimates), we can also quantify the uncertainty of the estimation at each tracking step (using the smoothed posterior covariance estimates). Thus, we have the possibility of aggregating this uncertainty over the entire branch to validate them. We explore this notion to create a criterion for accepting or rejecting branches.
By aggregating variance for all tracking steps in each branch, we obtain a measure of the quality of branches. A straightforward approach is to use total variance, obtained using the trace of each of the smoothed posterior covariance matrices. We average the sum total variance over the length of each branch, , to obtain a score, , which is then thresholded by a cut-off to qualify the branches,
2.5 Application to Airways
The proposed method for tracking tree structures can be applied to track airways, vessels or other tree structures encountered in image processing applications. We focus on tracking airways from lung CT data and present the specific strategies used to implement the proposed method.
2.5.1 Multi-scale representation
The measurement model discussed in Section 2.2
assumes a 4-D state vector as measurements to the RTS smoother. This is achieved by first computing an airway probability image using a k-Nearest Neighbour voxel classifier trained to discriminate between airway and background, described in. Blob detection with automatic scale selection  for different scales, , is performed on the probability image to obtain the 4D state measurements as blob position and radius. Indistinct blobs are removed if the absolute value of the normalized response at the selected scale, , is less than a threshold . This makes the representation sparse, , and the tracking more efficient than if performed at voxel-level. An example of the sparse representation can be found in Figure 2.
2.5.2 Initialisation of Branches
The multi-scale representation of the image data discussed above also provides a response corresponding to the best scale. As this response is normalised for scales, we incorporate this information in selecting the initial seed point for every branch. We start tracking from the seed point with the largest scale and the largest response. The initial direction information is obtained from eigen value analysis of the Hessian matrix computed at the corresponding scale provided in the measurement vector. Once a branch is tracked along the initial direction, we track from the same seed point but in the opposite direction. Thus, if a seed point is obtained from the middle of a branch we can track it bidirectionally. After tracking in both directions, all the involved measurements including the seed point are removed from the measurement vector, and the next best candidate seed point is chosen and tracking commences from there. The tracking procedure on the entire image is complete when no more seed points are available.
3 Experiments and Results
The evaluation was carried out on 32 low-dose CT chest scans from a lung cancer screening trial . Training and test sets comprising of 16 images each were randomly obtained from the data set. All scans have a resolution of approximately 1mm 0.78mm 0.78mm. The reference segmentations consist of expert verified union over the results of two previous methods [11, 12]. The proposed method is compared with region growing on the probability images.
3.2 Error Measure, Initial Parameters and Tuning
We use an error measure defined as the average of two distances, . The first distance, , captures the false positive error and is the average minimum Euclidean distance from segmentation centerline points to reference centerline points. similarly defines the false negative error, as the average minimum Euclidean distance from reference centerlines points to segmentation centerline points.
There are several parameters related to the RTS smoother that need to be initialised. These parameters were tuned using the training set and fixed for the evaluation on the test set to: standard deviations of the process noise,, measurement noise on radius mm and measurement noise on position mm. The initial covariance, across branches was set to . The most crucial parameter in the proposed method is the threshold parameter presented in Section 2.4.1. The threshold to validate branches is tuned to be . The gating probability was set to a high value, .
Figure 3 illustrates features of the proposed method by visualising centerlines overlaid on the reference segmentation. Influence of the threshold parameter is illustrated with the segmentation results for a single volume without any threshold (seen in Figure 2(a)) and after applying the tuned threshold (seen in Figure 2(b)). Evidently, thresholding the average total variance of a branch eliminates false positive branches.
The final output obtained from the method is a collection of disconnected branches. While such collection of branches are still useful in extracting biomarkers, for evaluation purposes we merge the results obtained with the segmentations from region growing on probability images and extract centerlines from the merged segmentation using 3D thinning, as seen in Figure 2(a) and 2(b). This also allows us to demonstrate the improvement our method provides by extracting peripheral airway branches, which are typically the challenging ones. One such combined result is shown in Figure 2(c), where the yellow centerlines correspond to region growing and blue one is the combined result.
Performance on the test set for two different scenarios of the proposed method is reported in Table 2 along with the numbers for region growing on probability images. The result for the best performing region growing on probability images is denoted with RG and those obtained by combining the proposed method with region growing are denoted as RTS+RG. We first combine the proposed method with the best performing region growing case (with minimum ) results and it is denoted as (RG+RTS). We observe an improvement of about 36% on . It is to be noted, there is substantial reduction in , indicating that many branches missed by region growing are now segmented. There is a very small increase in false positives which could also be due to the missing branches in the reference segmentation; however, the net result is a large improvement. To test whether the proposed method can simultaneously reduce the number of false positives and false negatives compared to region growing, we merge the proposed method with the region growing result that yields non-optimal , and do observe a reduction in both and when compared to the best performing RG as seen in the entries for (RG+RTS).
The computational expense for running the proposed method is small. The largest chunk of it is used in generating the multi-scale representation of the images, which is in the range of 10-15s per volume. Tracking using the RTS smoother and obtaining the segmentation takes about 4s on a laptop with 8 cores and 32 GB memory running Debian operating system.
4 Discussion and Conclusions
We presented an automatic method for tracking tree structures, in particular airways, using probabilistic state-space models and Bayesian smoothing. We demonstrated that branches can be tracked individually from across the volume, starting from several seed points. This approach of tracking branches from across the volume has the advantage that even in the presence of occlusions, such as mucous plugging or image acquisition noise, the chances of detecting branches beyond the occlusions are higher. An inherent measure of uncertainty in the branch estimates has been presented due to the Bayesian nature of the method. We demonstrated the use of thresholding this uncertainty measure to discriminate detected branches. The use of sparse representation of voxels in the image using blob detection makes the method computationally efficient.
A possible limitation with the proposed method is that it yields a disconnected tree structure. For applications where this is an issue, one can enforce a global connectivity constraint on the disconnected set of branches to obtain fully connected tree as done in  or similar. It is also possible to derive biomarkers directly from the disconnected branches, as shown in .
We performed an evaluation of the results obtained from the proposed method by combining it with the results from region growing on probability images. We showed that there is substantial improvement in the segmentation results, indicating that the exploratory approach taken up in our method has potential in improving tree segmentations.
This work was funded by the Independent Research Fund Denmark (DFF) and Netherlands Organisation for Scientific Research (NWO).
-  Lesage D, et.al. A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes. Medical image analysis. 2009 Dec 31;13(6):819-45.
-  Lo P, et.al. Extraction of airways from CT (EXACT’09). IEEE Transactions on Medical Imaging. 2012 Nov;31(11):2093-107
-  Friman O, et.al. Multiple hypothesis template tracking of small 3D vessel structures. Medical image analysis. 2010 Apr 30;14(2):160-71.
-  Yedidya T, et.al. Tracking of blood vessels in retinal images using Kalman filter. InComputing: Techniques and Applications. Digital Image 2008 (pp. 52-58).
-  Florin, C., et.al. Particle filters, a quasi-monte carlo solution for segmentation of coronaries. MICCAI. 2005 (pp. 246-253). Springer Berlin Heidelberg.
Lesage, D., et.al. Adaptive particle filtering for coronary artery segmentation from 3D CT angiograms. Computer Vision and Image Understanding. 2016.151, pp.29-46
-  Särkkä, Simo. Bayesian filtering and smoothing. Cambridge University Press, 2013.
-  Lindeberg, Tony. Feature detection with automatic scale selection. International journal of computer vision 30.2. 1998.
-  Bar-Shalom Y, Willett PK, Tian X: Tracking and data fusion. YBS publishing; 2011
-  Pedersen, Jesper H et. al. The Danish randomized lung cancer CT screening trial-overall design and results of the prevalence round, Journal of Thoracic Oncology, 2009.
-  Lo, Pechin, et.al. Vessel-guided airway segmentation based on voxel classification. First International Workshop on Pulmonary Image Analysis. MICCAI. 2008
-  Lo, Pechin, et.al. Airway tree extraction with locally optimal paths. MICCAI. 2009
-  Graham, Michael W., et al. Robust 3-D airway tree segmentation for image-guided peripheral bronchoscopy. IEEE transactions on Medical Imaging (2010)
-  Sørensen, Lauge, et al. Dissimilarity-based classification of anatomical tree structures. Information Processing in Medical Imaging. Springer Berlin/Heidelberg, 2011.