Automatic Coronary Artery Segmentation Using Active Search for Branches and Seemingly Disconnected Vessel Segments from Coronary CT Angiography

We propose a Bayesian tracking and segmentation method of coronary arteries on coronary computed tomographic angiography (CCTA). The geometry of coronary arteries including lumen boundary is estimated in Maximum A Posteriori (MAP) framework. Three consecutive sphere based filtering is combined with a stochastic process that is based on the similarity of the consecutive local neighborhood voxels and the geometric constraint of a vessel. It is also founded on the prior knowledge that an artery can be seen locally disconnected and consist of branches which may be seemingly disconnected due to plaque build up. For such problem, an active search method is proposed to find branches and seemingly disconnected but actually connected vessel segments. Several new measures have been developed for branch detection, disconnection check and planar vesselness measure. Using public domain Rotterdam CT dataset, the accuracy of extracted centerline is demonstrated and automatic reconstruction of coronary artery mesh is shown.


Introduction
Three-dimensional (3-D) reconstruction of the coronary arteries from coronary computed tomographic angiography (CCTA) would lead to higher accuracy and reproducibility in the diagnosis and to better precision in the quantification of severity of the coronary arteries diseases. It is also essential for the 3-D reconstruction and post-processing tools such as curved multi-planar reformatted (MPR) images through the lumen of each coronary artery. Furthermore, it is also one of the prerequisite steps in subsequent analysis, such as detection of lesions [1] and image fusion [2,3]. Also, full reconstruction of coronary artery tree with lumen boundary is important.
Since advances in computational fluid dynamics and image-based modeling now permit determination of rest and hyperemic coronary flow and pressure from CCTA, these techniques have been used to non-invasively compute fractional flow reserve (FFR), which is the ratio of maximal coronary blood flow through a stenotic artery to the blood flow in the hypothetical case that the artery was normal, using CTA images [4]. It is pre-assumption that the entire artery tree structure with accurate lumen boundary has to be obtained for the use of such technology.
A significant amount of research has been done on the segmentation of vascular structures and of the coronary arteries in particular [5,6]. However, many of these vessel segmentation methods require user interaction: manual definition of the start and end points and vessel directions [7], or manual insertion of the intermediate points to bridge gaps [8][9][10]. Under these interactive segmentation methods, the accuracy of segmentation result is strongly dependent on the interaction from individual users implying a low reproducibility in segmentation results. A trade-off between robustness and automation has to be carried out whenever the segmentation algorithm is to be embedded in a commercial package for clinical applications. Moreover, a certain amount of interaction usually leads to the increment of the processing time.
Several approaches were proposed for automatically tracking the coronary arteries. They appear in various forms. Tek et al. [11] detected ostium locations on the surface of aorta and then started the multi-scale medialness-based vessel tracking from the ostium locations. Medialness filter is their base of vessel finding. Kitslaar et al. [12] uses connected component and morphological filter to detected vessel regions and refines centerlines. Bauer et al. [13] presented an automatic approach that consists of generic methods for detection of tubular objects, extraction of their centerlines, and grouping of these centerlines into tree structures. Zambal et al. [14] calculated candidates for coronary artery seeds and tracks the vessel segments based on vessel surface gradients by approximating the position of the heart and then cylindrical sampling patterns were fitted. Kitamura et al. [15] proposed a method to train a classifier of a tubular 3-D object with a dimension reduction approach using Hessian analysis. Zhou et al. [16] segmented the vascular structures within the heart region using a multi-scale coronary response and then tracked the coronary arteries by a 3-D dynamic balloon tracking algorithm.
The purpose of this work is to extract the whole coronary artery tree segmentation system by actively searching for branches and disconnected vessels. It utilizes the statistical branch occurrence model and checking disconnection model by image appearance. While tracking vessels with the stochastic model for curvature penalized vessel geometry, it finds the branches, a stenosis lesion and seemingly disconnected vessels by atherosclerotic plaque. The whole tree structure is robustly constructed with a probabilistic branch detection method. Vessel Tracking as MAP Estimation Problem. Our general framework for vessel tracking and segmentation is to estimate the geometry of the vessel by formulating the problem as MAP (maximum a posteriori) estimation. We look for that vessel for which Pðhypothesized vessel model j image dataÞ is a maximum. Since this posteriori likelihood of a hypothesized vessel model given the image data can be written as

Bayesian Formulation of Vessel Tracking Problem
The vessel geometry at i is denoted X i , X i ¼ ðx i ;r i Þ 0 .x i andr i represent the coordinate of the center point in 3-D and radius at i'th position respectively. Denote by X the vector having all the x i , r i in a vessel as its components and by Y the vector having as its components the CT intensities y l,m,n at all the voxels (l, m, n) in the vessel region. Then the proposed system seeks the vessel geometry maximizing the following equation: where PðYjXÞ is the likelihood and PðXÞ is the prior distribution. Geometric Vessel Model. The vessel model is generated based on the following assumptions and statistically trained. The cross-section of a vessel is closed curve with an arbitrary form satisfying the flowing shape restrictions: 1. The vessel radius variation is small and the vessel radius change is likely to be slow.
2. The vessel direction changes are likely to be slow.
3. The local average of gray level in a vessel is likely to vary slowly. 4. The gray level variation between a vessel and background is likely to be large.
Meir [17] uses stochastic model for winding road from satellite image. Their stochastic model reflects faithfully the behavior of long and snaky shape object which narrows and widens along its centerline with the geometric constraint. This model can be easily extended to three dimensional problems. Vessels in human body have similar characteristics as the roads in satellite images. They both change the propagating directions and widths slowly along its centerline with occasional branches. For our problem, a three dimensional stochastic process model is built exhibiting the preceding vessel geometry.
Specifically, autoregressive processes are designed to model vessel center line, vessel width, gray level within the vessel, edge strength at the vessel boundary, and gray levels outside the vessel and adjacent to the boundaries. These regions are referred as background. Note that the vessel geometry processes are hidden, i.e., they are not observed directly in the data. The stochastic processes are functions of a discrete parameter i which can be thought of as time or distance, in voxels, along three dimensional axis. Vessel center-curve at i isx i , which takes values in 3-D. This variable is not quantized, it takes arbitrary real values. The fx i g process is given by Eq (4), where ε x i , is a zero mean, white, Gaussian driving noise. This process will generate a straight line if ε x i is zero.
In Eq (5), vessel radius is modeled such thatr i is perpendicular circle to center curve fx i g. The stochastic processes ε x i and ε r i are independent Gaussian white noise sequences with zero means and variances σ d i , and σ r i , respectively. The vessel obtained by the un-forced solution(ε x i = 0, ε r i = 0) will be straight cylinder with radius r i as shown in Fig 2. Vessel boundary location is uniquely determined by thex i andr i .
A second order Markov Process is employed to model the mean intensity of the CT image in sequential spherical data of the vessel to be consistent with our vessel model. For foreground and background model,f f i andb i are the mean intensities inside and outside of sphere i, and ε f i , ε b i are their Gaussian white noise sequences with zero mean and variance σ f i and σ b i , respectively.
Here is the list of state variables and equations used in vessel propagation model. General state equation of markov random process can be written for foreground and background as follows.
In Eq (9), the first row gives geometric constraint which comes from Eq (4), and the second and third rows of the equation come from Eqs (6) and (7). They all together sum up to the following equation where maximum-likelihood estimate is found. Furthermore, MAP estimatê x i;r can be also obtained if proper prior is found for vessel geometry.
Curvature and Torsion. In differential geometry, the fundamental theorem of curves states that any regular curve in three dimensions with non-zero curvature has its shape completely determined by its curvature and torsion. Curvature measures how much a curve bends on a plane and torsion measures how much a curve deviates from its osculating plane. The shape of a vessel does not change abruptly. Therefore, for curve reconstruction problems, it makes sense to regularize the solution with curvature and torsion prior [18]. If the reconstructed curve should follow closely the shape of an underlying surface which is locally planar. In such situations, it may be advantageous to penalize torsion. As shown in Eq (10), the vessel trajectory is estimated using the paths with two previous points f i−1 and f i−2 . In each step of tracking procedure, three points are used for computing likelihood. In other words, the curvature of vessel trajectory is obtained with three points and the curvature constraint is used in our geometric model depending on ε x i in Eq (4). If ε x i is small, it will prefer smoother trajectory and if big, it will allow more complex trajectory. Curvature and torsion regularization method for curves is extensively performed in [19,20]. However, torsion is not currently used in our vessel tracking process. The vessel trajectory model may be improved by employing torsion constraint.
Usage of Region Growing in our Method. Region growing is an approach to image segmentation in which neighboring pixels are examined and added to a region class. It is one of the practicable ways to achieve image segmentation and many conventional algorithms in medical image analysis employ region growing method. Region-growing algorithms turn out to be very practicable way to achieve a time-saving vessel segmentation even with the disadvantages due to their inherent sequential nature. However, in vessel segmentation problem, by just region growing implementation, highly accurate vessel segmentation cannot be achieved since many other regions next to vessel have similar HU level, e.g., left ventricle or other vascular organs. It has to be incorporated with other processes to obtain successful vessel segmentation.
To overcome such disadvantages, we merely use region growing as a rough guide to reduce the number of computation. Later our geometric vessel model will separate vessels from background correctly.
Adaptive Region growing first selects the candidate region where vessel may lie in. This practical procedure defines the region where the probability of vessel is almost zero and reduces the amount of computation. Fig 3 shows the region growing and vessel geometrical model together. The region growing selects the candidate area (in red) in advance.

Branch Detection
For the reconstruction of complete coronary artery tree, the tracking or segmentation system requires to find the branches of vessels and track the sub-branch vessels. When a new branch is encountered while tracking a vessel, the tracker must initiate a new vessel seed for tracking and branch detection mechanism is required. However, patients with stenosis occasionally show disconnected vessels in CT due to narrowing or blockage of the coronary arteries. It often occurs near branch points [21]. Hence we devise a method named "active search" which can search and find such disconnected arteries and branches.
For vessel branch detection, conventional methods can be grouped into two categories; region growing based methods and statistical analysis based methods of the vessel propagation distribution. And there are also variants of two categories. For region growing methods, Tek et. al. [11] proposed a method that after extracting centerlines of the three major coronary arteries, the algorithm starts to trace side branches. First, the bifurcation of a side branch is detected on a major centerline using region growing based lumen segmentation. Starting from a centerline point, bright voxels connected to the current point are added iteratively. The growing front is traced. If a side branch presents, the region growing procedure goes into this side branch. A side branch is detected when it finds a front with a distance to the existing major centerline larger than a threshold. At each detected bifurcation point, a data-driven centerline tracing process is initialized. Jiang [22] adopts a sphere propagation method to determine the coronary artery branches. By using different size of sphere and uniform ray casting they find the branch points. They produce incremental size spheres R+1, R+2 and R+3 and propagate uniformly distributed rays. If rays reach outermost shell satisfying pre -specified condition and concentrate at separate two clusters, they can find side branches. The progressive region growing method has been developed and skeletonization method is used for tree structure finding [23]. Such analysis methods are supported by extracting the skeletons or medial axes of the segmented voxel sets.
For statistical method for branch detection, Dip test has been used [24,25]. A Dip statistic test is to test the multi-modality of the re-sampled particles against the null hypothesis of a unimodal distribution. The Dip test, a measure of departure from uni-modality in 1D, measures the maximum distance between the empirical distribution and the best fitting unimodal distribution. Then, if a junction is detected by using dip test, k-means clustering (k = 2) is performed to find the direction of the daughter branches.
Our Branch Detection Procedure. We devise a robust branch detection system incorporating conventional statistical method and newly developed active search method. As shown in Fig 1, while tracking artery, the system continuously detects bifurcation or multifurcation by checking if there exists separate clusters for high probability candidate points. If a cluster or multiple clusters are found, those points are denoted as branch points. For active search scheme, the branches are actively searched in three categories. They are discussed in next subsection.

Active Vessel (Branch) Search
In CCTA images, vessel sometimes can be seen to be disconnected due to plaque build-up or noise and sequential tracking methods will stop since there is a gap between vessel segments. In such cases, the active search method will overcome the gap.
The active search are used: to check if a vessel is locally disconnected (type 1), to find disconnected branch (type 2) and according to branch occurrence model (type 3). Their usage is shown in the workflow chart in Fig 1. When a vessel is abruptly disconnected or there is probable branch area, active search is employed. Conventional methods for branch detection described in the above sub-section utilize the fact that sub-branch pipelines are locally detectable and well separated. They test if there is a separation of highly contrasted regions using either region growing or statistical test for probability distributions of tracking mode. However, when local discontinuity exists such methods will fail. This active search method is also applied for branch detection. The examples of both cases follow. In Fig 4 which is one of Rotterdam test data (testdata 26), a vessel is disappearing and re-appearing along several axial slices of images. In the figure, RCA is clearly shown in (a-b) then in (c-e) it is lost and re-appeared again in (f-h). In this case, coronary artery tracking method will stop due to low intensity of invisible vessel. Hence an active search will start when abrupt termination is detected. In other hands, Fig 5 shows a discontinuity at LAD branch point. Near a bifurcation or branching point, discontinuity occasionally occurs due to complex fluid turbulence. However, the active search is an exhaustive method and brute force search will require of a large number of computation. To reduce the number of computation, the branch occurrence model has been devised. The active search comes from three different manners as described above. The following is the 1. Eq (10) gives the likelihood of each direction (pointx i ). Many points will have high probabilities and cluster into vessel propagation directions. From the cluster of such points, the branch can be detected. If there exists two significant clustering, bifurcation is found. (So far, this is closely related to conventional branch detection DIP test method).
2. If no branch point is found, active search is employed at highly probable region. (i:e:; if the neighboring volume S has higher intensity above some threshold.) 3. Active search will investigate equally spaced pointss in the neighboring volume S (s 2 S), where dist(s,x i ) T th and dotðs Àx i ;x i Àx iÀ1 Þ T y where dist is distance measure of two points and dot is dot product of two vectors.
4. If new points are found, they will be queued as new seed points.
Branch Occurrence Model using Possion Distribution. This type of active search (type 3) is used for ensuring that left main (LM) artery's first bifurcation is not missing. For the branch occurrence model, the first branch location of left main (LM) artery is measured from the aorta. LM coronary artery usually bifurcates into left anterior descendent artery (LAD) and left circumflex artery (LCX) branches. The distances of this first branch of left main coronary arteries from ostium are shown in Table 1. The training set of Rotterdam data is used for the distance measure. Automatic Coronary Artery Segmentation Using Active Search from Coronary CT Angiography  The occurrence based process can be modeled using Poisson process. Poisson is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time and/or space if these events occur with a known average rate and independently of the time since the last event. Poisson distribution is written as In Eq (11), x is the distance from the ostium and the λ is the expected distance value. Using medical expert annotated distances of first branch from ostium in Rotterdam training data [26,27], the Poisson distribution model has been trained as shown in Table 1. which shows the distance of first branch from the ostium.l is the estimate of Poisson parameter and 17.3325. In the experiments, all of the first branch (LM toward LCX and LAD) were correctly found in 32 dataset by utilizing the active search method together with the Poisson branch occurrence model. As described before, highly stenotic coronary arteries appear sometimes as disconnected as shown in Figs 4 and 5 even though they are connected. In such cases, we use active search scheme which can detect disconnected vessel and connect them. While MAP estimate is found via maximizing using Eqs (2) and (10), active searching scheme runs at the same time as shown in Fig 1. With two previous center point x i−2 , x i−1 and vessel geometry constraint, the prior probability distribution of vessel geometry is assigned. Search candidates points x i 's are uniformly generated on the hemisphere with radius S (step size), and each point x i is assigned with the prior probability depending on geometric constraint. Also the likelihood from foreground and background model and others is computed. The posterior probability is computed from the prior and the likelihood. For branch detection, the candidate points with relatively high posterior probability on the hemisphere is now clustered into small group and the distance of clusters are measured. If the distance is over some threshold the branch is detected. At the same time, active search with Poisson process for the branch occurrence model and active search for high probable region from the HU level is also investigated.

Planar Vesselness Measure
It is necessary to analyze the shape of coronary artery for determining if detected or segmented object is coronary arteries and not others. While the tracking system tracks down the coronary arteries and their branches, the tracker may stray into wrong region like left ventricle. Tracker needs to check if the tracked object is a vessel. Hence, detection and calculation of topological features for such 3-D tube-like objects are required via description of shapes and other characteristics of complexity. This should allow to represent interior structure of 3-D objects.
The following planar vesselness measure is used in our system. We use a fast and efficient planar version of vesselness measure. It can tell when a vessel tracker must stop or if searched point is a vessel-like point in the active search process.
Analysis of Cross Sectional Shape. To speed up the tracking process, simple planar Hessian based measure is employed in our system. To detect foreground objects solely, Frangi Automatic Coronary Artery Segmentation Using Active Search from Coronary CT Angiography proposed a method that can find elongated, i.e., tube-like objects [28] and the method shows successful results. However it requires an extensive 3-D computation along the scale space in conventional Frangi's vesselness measure. A simple planar version is devised here. If the planar measure scores low in consecutive cross sections. One can tell it is not a vessel. The underlying assumption is that the cross section of coronary artery or its area is elliptic and has higher HU value compared to myocardium or surrounding materials. Let the linear scalespace representation of the cross sectional image I 0 (x, y) at scale σ given by: where G(x, y; σ) is the Gaussian kernel with scale factor σ. It starts with iso-intensity surface definition which is where σ is omitted for simplicity. Let δ = (Δx, Δy) and using Taylor's expansion in vector notation, it can be simplified as: where Hðx; y; zÞ ¼ is the Hessian matrix. This quadratic equation represents a general quadric surface including ellipsoids. Eq (14) can be divided into two parts: If H is positive or negative definite, it represents convex or concave region, respectively (i.e., dark area or bright area). Since we are looking for coronary arteries in CCTA, negative definite will be a good indicator. Therefore, in 2-D cross sectional views of coronary arteries, the Hessian will be H ¼ I xx I xy I yx I yy ! and the bright elliptic feature like cross section of a vessel is equvallently represented as the set of fðx; yÞjHðx; yÞ < 0g ¼ fðx; yÞjI xx þ I yy < 0 \ I xx I yy > I 2 xy g ð17Þ Hence the planar vesselness measure v is computed from Eq (16). It can be written as the function of I xx +I yy and I xx I yy À I 2 xy . For detail, please refer to [29]. The performance of our planar vesselness measure is compared to conventional 3-D vesselness filter in Fig 8. Both measure have significantly high value along the tracked vessel (especially proximal and mid part of a vessel) and seems to be correlated in the figure. When there is no vessel, both value will be almost zero. The following subsection summarizes the overall branch detection procedure.
Initiation of Seed Point for Tracking. In the experiment, the start points (seed points) are fixed and given in rotterdam coronary artery algorithm evaluation framework. The initial direction is found via ML estimators like the cylinder in Fig 9. At a given seed point SP and searching direction in 3-D, the cylindrical mask is generated using Gaussian shaped kernel along the long axis. The most probable direction is found. When active search are used, the initial direction of the vessel is found the same way.

Experiment
The accuracy of the centerline extracted from the suggested method is tested on 32 CCTA public Rotterdam challenge dataset [26,27]. These datasets were provided by the organizers of the "Coronary Artery Tracking Challenge". On the evaluation framework, it measures the accuracy of the centerlines of four significant vessels for each dataset and in total the number of tracked vessel is 128. While the evaluation method only assesses the estimated centerlines, our system finds the centerlines and automatically reconstructs the vessel lumen boundary mesh at the same time. Our centerline extraction experiment using Rotterdam dataset will be categorized as method #2 following the Rotterdam coronary artery algorithm evaluation framework rules. One single point S (starting point at ostium) is pre-defined at the challenge site [26,27] for each left main artery (LM) and right coronary artery (RCA). The detailed experiment   (Tables 3-6). Finally, all together with 32 dataset, the result is OV 83.6%, OF 65.3%, OT 88.2%.
The running time is from 4 to 6 minute depending on complexity of each dataset. The proposed method finds the entire tree structure including small vessel branches as shown in Figs 10 and 11. In the first figure, the colored spheres represent reference points for the identification of selected vessels for the Rotterdam challenge. Fig 11 shows the reconstructed lumen boundary mesh which is to be used for CFFR (computed fractional flow reserve). For the generation of meshes, the lumen boundary surface is automatically optimized via surface constraint and image gradient as a post process.
In Tables 2-6, the scores and ranks for individual dataset are shown in detail. The rank is given from the score which is computed from non-linear increasing function of overlap ratio [26]. In Table 3, our Test dataset result is shown. In Table 4, it shows that the mean percentage values of OV and OT are significantly different (p < 0.0001).
For the comparison of challenge participating methods, the overlap scores of all the challenge participants (OV, OF, OT and average of them) are shown in Table 7. For the average overlap score, the median is 53.7 and ours is 51.8. The first and third quartiles are 47.1 and 63.2, respectively. Hence, interquartile range is 16.1. These statistics are obtained from the

Discussion
A new vessel tracking method based on stochastic geometric processes with active branch search is developed for 3-D coronary artery segmentation. The suggested method reconstructs the whole coronary artery tree structures with the estimation of the vessel lumen boundary which can immediately generate mesh of the coronary structure (Fig 11). The result of developed method is compared with the results of other methods as the centerline accuracy problem using Rotterdam coronary artery algorithm evaluation framework [26]. It shows high Such methods employ complicated hybrid methods using prior information like the left ventricle segmentation [30] (ModelDrivenCenterline in Table 7) or requires multiple processes [31] (GFVCoronaryExtractor). However, since our system is based on pure MAP (maximum a posteriori) estimation via vessel geometry model, we believe it is easy to adopting other techniques and can be combined for improvement. With regards to the accuracy distance measure number in Tables 5 and 7, our method shows rather large distances compared to other methods. We believe the reason is that our system concentrates on the reconstruction of the lumen wall instead of finding accurate centerlines and it does not optimize or adjust the centerlines after the first detection. For the summary of newly developed active search scheme, it comes in three ways. The first one is for disconnected vessel (type 1) and the second is for disconnected branches (type 2). The third one is using branch occurrence model (type 3). For the disconnected vessel the search is straightforward and just tries to find the vessel segment over the gap. Then it connects them. And for the missing branches, both high probabilistic region investigation (type 2) and the branch occurrence model (Poisson) are employed (type 3). Such active search method plays an important role as follows. Especially when there is a local discontinuity in vessel image, newly suggested active search method overcomes difficulties due to the nature of sequential vessel tracking techniques. In CCTA images, vessel sometimes can be seen to be disconnected due to plaque or noise in the image. In such cases, the active search method will overcome the seemingly disconnected gap. In addition, such discontinuity frequently occurs near branch points. When there are bifurcations, there can be blood turbulence at the bifurcation area and plaque can be easily built-up [21]. Such plaque can makes vessel branch seem to be disconnected. This is why the type 2 active search is employed. This is an improvement to conventional branch detection and vessel finding algorithm since it provides an efficient local search scheme. Benefiting from the new active search scheme, the proposed method has high capability of finding seemingly disconnected vessel and branches.
For the termination of vessel tracking, a newly devised practical planar vesselness measure is used. It is computationally efficient and comparable to popular vesselness measure [28] as shown in Fig 8. More recently, a new deep neural network learning based vessel classifier is devised and tested. It is based on a convolutional neural network (CNN) and a preliminary test shows a high classification rate in distinguishing vessels from non-vessel substances [32].
As for statistical branch occurrence location model, it is used only for the first branch of the left main artery for now but the complete statistical branch location model for the entire coronary artery system is under study. More analysis for the active search is required. After the entire tree structure is found, the lumen mesh is optimized via surface constraint and gradient information and this mesh is to be used for other tasks (see Fig 11). Also, automatic aorta and ostium segmentation is under development now and the system will be full automatic without requiring starting points.