Multiple Hidden Markov Model for Pathological Vessel Segmentation

One of the obstacles that prevent the accurate delineation of vessel boundaries is the presence of pathologies, which results in obscure boundaries and vessel-like structures. Targeting this limitation, we present a novel segmentation method based on multiple Hidden Markov Models. This method works with a vessel axis + cross-section model, which constrains the classifier around the vessel. The vessel axis constraint gives our method the potential to be both physiologically accurate and computationally effective. Focusing on pathological vessels, we reap the benefits of the redundant information embedded in multiple vessel-specific features and the good statistical properties coming with Hidden Markov Model, to cover the widest possible spectrum of complex situations. The performance of our method is evaluated on synthetic complex-structured datasets, where we achieve a 91% high overlap ratio. We also validate the proposed method on a real challenging case, segmentation of pathological abdominal arteries. The performance of our method is promising, since our method yields better results than two state-of-the-art methods on both synthetic datasets and real clinical datasets.


Introduction
Automatic vessel segmentation in three-dimensional (3-D) medical computed tomography (CT) images plays a fundamental role in many clinical fields: study of anatomical structure [1], quantification of vascular diseases (stenosis, occlusion, and calcification) for clinical diagnosis [2], surgery planning [3], and patient-specific flow simulations [1]. Vessel segmentation can help clinical workers to have an intuitive impression of vessels and blood supply. Based on the vessel segmentation, clinical workers can also establish the patients' response to treatment and determine the stage of diseases, to further plan a minimally invasive surgery. All these applications ask for a competent segmentation technique, which has the capability of segmenting vessels accurately, not only for normal vessels but also for vessels with the presence of pathologies.
In recent years, many methods [1][2][3] have been proposed to segment vessels with the presence of pathologies. These methods can be roughly classified into three categories: (a) feature based segmentation approaches [4,5], which have been proven to be efficient in detecting vessels at different scales; (b) tracking based segmentation approaches [6][7][8], which have the capability to be robust against noise; and (c) model based approaches [9][10][11][12], which resolve subsequent two-dimensional (2-D) slices of vessels using tubular shape priors.
By targeting the curvilinear structures, the feature based segmentation approaches are efficient in identifying vessels from other tissues [5,[13][14][15][16][17][18]. They have several common detection procedures: firstly, assuming that vessels have identifiable curvilinear structures and then calculating the eigenvalues of the Hessian matrix to traverse the whole 3-D image at different scales, by convolving with 3-D Gaussian filters. The standard deviation of the Gaussian filter defines the detecting scale. Third, a response function based on these eigenvalues is constructed, which can detect the curvilinear structures locally at a certain scale. Since the curvilinear structures signed to vessels are significantly different from other tissues, the feature based segmentation approaches can recognize vessels from other nonvessel structures, such as planar structure, blob, noise, or no structure. The response 2 BioMed Research International consisting of eigenvalues can represent the local structure, when it comes to its maximum over different scales [5,16].
However, the feature based segmentation approaches do not produce direct vessel segmentation, they output responses instead of segmentation. As a result, their results can be used in structure analysis, where many researchers consider it as the advantage of the feature based segmentation approaches [19,20], especially in combining the tracking based segmentation approaches. Based on the responses of the feature based segmentation approaches, the tracking based segmentation approaches resolve subsequent twodimensional (2-D) slices of vessels, by using the tubular shape priors [19,20]. Based on the tubular shape structure, Tyrrell et al. [21] propose another 3-D structure, 3-D cylindroidal superellipsoids, combining with the local regional statistics method. This tracking methodology extracts topological information from microvasculature networks.
The tracking based segmentation approaches are shown to be robust against noise [19][20][21]. However, the shape priors they use are too exclusive [22], which lead to false detection with the complex vessel boundaries. Targeting at this challenge, the statistical mixture models have been proposed [22][23][24]. They combine the statistical mixture models with expectation-maximization algorithm, to detect complex vessel boundaries. However, these approaches ask for an accurate parametrical estimation or nonparametric modeling. As a result, they all rely on the second-order derivative information, which is associated with the principal curvatures of image intensities. They suffer from sensitivity caused by local deformations, which are considered as pathologies in this work.
In the past decade, after many sophisticated vascular segmentation algorithms have been proposed, there are still several challenges remaining unsolved, such as detecting vessel radius with the presence of pathologies (stenosis and occlusion), delineating vessel accurately with low contrast vessel boundary, and distinguishing vessels from vessellike tissues. All these challenges are caused by pathologies, especially by the vessel calcification. As far as we know, the accurate delineation of pathological vessel is valuable in clinical fields. However, the proposed challenges may lead to false segmentation, which makes the vessel segmentation useless. For further information, several good reviews can be found in [1][2][3].
Addressing these unsolved challenges, we propose a Multiple Hidden Markov Model (MHMM) for pathological vessels, by taking the advantages of the redundant information embedded in multiple vessel-specific features and the good statistical properties coming with the Hidden Markov Model (HMM). As we target at pathological vessels, the shape complexity of obscure vessel boundary may not be captured by one single existing vessel-specific feature. Thus, we dig the redundant information embedded in multiple vessel features automatically by using multiple HMMs and try to combine multiple vessel features optimally by employing the good statistical properties of HMM. Different from most of the other researchers, we build the MHMM to describe the transition course from the inner vessel to the outside of vessel.
By doing so, the MHMM has the capability to delineate the obscure vessel boundaries accurately. This paper is organized as follows. Section 2 introduces the proposed method MHMM. Section 2.1 details the training samples and the testing samples we use in this work. The multiscale vessel-specific feature set is followed in Section 2.2. Then we present the details of training a single HMM in Section 2.3. The MHMM combining multiple vessel features with multiple HMMs is given in Section 2.4. The evaluation datasets, the evaluation methods, and the experimental results are shown in Section 3. In Section 4, we discuss the advantages and disadvantages of our method, possible improvements, and future work. Finally, we conclude the paper in Section 5.

Method
This paper describes a Multiple Hidden Markov Model (MHMM) for abdominal artery with the presence of pathologies, which takes two advantages: (1) the redundant information included in multiple vessel features and (2) the good statistical properties of Hidden Markov Model (HMM) [25]. First, we extract vessel axis from a gray-scale 3-D CT scan, by using our previous method [9]. Then, along this extracted vessel axis, voxel series are generated on the cross-sections for training and testing purposes. Subsequently, multiple vessel features are calculated for each voxel in the voxel series. Finally, a MHMM is trained based on these vessel features. The trained MHMM can detect vessel boundary with the presence of pathologies, by expressing the similarity marching degree (posterior probability).

Training Samples and Testing Samples.
We randomly selected 50 CT datasets of abdominal artery from 120 CT datasets, which were collected by our collaborators (the clinical workers from the 2 affiliated Hospital of Harbin Medical University). These 120 CT datasets were collected from 120 patients, 73 males and 47 females with mean age 65 years (min 58, max 91), who entered our study nonconsecutively from October 2014 to October 2017. These 50 selected CT datasets are known with aortoiliac stenosis, calcification, and occlusion.
Two expert raters were asked to segment the datasets. We understand that the ground truth segmentations of those clinical datasets do not exist, even segmented by multiple experts. Thus, we developed a graphical user interface (GUI) in our lab [9] to help the experts mark different tissues in 3-D CT images, including pathologies. With this GUI, the user can do the marking job in two stages: propagating stage and checking stage. In the propagating stage, the GUI can propagate to mark the target tissue in 3-D CT image automatically, with seed points given by experts. In the checking stage, the GUI helps the user to cine-page through the CT slices, scroll in and out of individual tissue, adjust window setting, and zoom to improve visualization. As a result, the user can check the propagating result and set new seed points. By using this GUI, two expert raters were asked to mark vessels and pathologies (calcification and iliac stent).

Vessel axis
Cross-section We collect voxel series as training samples and testing samples, based on the segmentations provided by the expert raters. Firstly, vessel axis is extracted by using our previous method [9], see Figure 1(a) for a snapshot. Then, rays are spaced every 5 ∘ from the axis point on the cross-section, shown in Figure 1(b). Finally, voxel series are collected along the rays (Figure 1(c)). Along the vessel axes, the voxel series are gathered cross-section by cross-section, to form the training samples and the testing samples. Let us denote all the voxel series as = { 1 , . . . , , . . . , }, where = {V ,1 , . . . , V , , . . . , V , } is a voxel series for one ray, V , is a voxel in the voxel series , and is the number of voxel series. Please notice that one voxel series is considered as a whole by the proposed MHMM, to describe the tissue transition from the inner vessel to the outside of vessel.
We target at accurate vessel segmentation with the presence of pathologies. However, the pathology obscures the vessel boundary and makes it difficult to identify vessel from pathology. Thus, we divide the voxels into 4 classes: lumen (vessel), intima (vessel boundary), pathology (pathological tissue), and adventitia (outer tissue). By doing so, we can train the MHMM to focus on obscure vessel boundary. To form the vessel boundary class, we strip the outermost layer away from vessel along the ray. Since the vessel wall thickness changes with vessel scale, we change the size of the outermost layer from 2 voxels to 4 voxels, to cover the vessel scale changes, according to the CT scan resolution in this work.

Feature Set for Multiple Hidden Markov
Model. Due to the presence of pathologies, it is challenging to segment vessels accurately. And the variance of vessel intensity makes this even worse. In our previous work [26], we found that no vessel feature is optimal in this challenging case. Thus, we employ multiple vessel-specific features here and try to explore the redundant information included in these features, which can cover the widest possible spectrum of situations. Since the diameter (scale) changes in vessels, all the vesselspecific features we chose are multiscale vessel features. These features are extracted by using Hessian matrix, which uses Gaussian convolution to capture the gradient based shape information. And the scale of vessel, , is represented by the standard deviation of Gaussian convolution. In the following, we introduce the multiscale vessel-specific features employed in this paper.
(1) Sato feature ( 1 ) [5]: this feature is designed for healthy vessel, with bright tubular structure in a dark environment. For one scale, , the Sato vessel feature 1 can be described as follows: where 1 , 2 , and 3 are the eigenvalues of Hessian matrix, subject to 1 > 2 > 3 [5]. These eigenvalues are calculated by using Gaussian convolution, and the scale is the standard deviation for Gaussian convolution; 1 and 2 are two preset parameters, 1 < 2 . In this work, 1 and 2 were fixed to 0.5 and 2, respectively.
(2) Frangi feature 2 [14]: this feature promotes the enhancement of tubular structure and smoothes other structures out. In the original work, this feature gives promising performance to the vessels with low contrast. For a single scale, , the Frangi vessel feature is given as follows: where 1 , 2 , and 3 are defined in (1). The scale is also the standard deviation for Gaussian convolution; = | 2 |/| 3 | (subject to ) identifies plate-like structures from tubular structures, while = | 1 |/√| 2 3 | (subject to ) identifies blob-like structures from tubular structures; = √ 2 1 + 2 2 + 2 3 (subject to ) smoothes background noise out. By following the original research, we set the parameters as follows: = 0.5 and = 0.5. Here is equal to half of the maximum Frobenius norm of the Hessian over all Frobenius norms computed on the whole image.
(3) Shikata feature 3 [27]: based on Frangi feature, Shikata feature assumes the cross-section of tubular structure has 2-D Gaussian intensity distribution. This assumption gives Shikata feature the capability of detecting small vessels. The feature 3 is given as follows: where 2 is the second eigenvalues of Hessian matrix, defined in (1) and ( → ) is the intensity at voxel → .
(4) Li feature 4 [28]: this feature enhances the tubular structures by using curvature analysis, which gives Li feature the potentiality to distinguish tubular structures from other structures. This Li feature 4 can be described as follows: (5) Manniesing feature 5 [29]: focusing on background noise and small vessels, this feature improves the Frangi feature, by using a nonlinear anisotropic diffusion approach. This feature, 5 , can be given as follows: where = | 2 |/| 3 | is designed to identify plate-like structures from line-like structures, = | 1 |/√| 2 3 | has the capability to discriminate blob-like structures from line-like structures, and = √ 2 1 + 2 2 + 2 3 eliminates background noise. By following the original research, we preset = 0.5, = 0.5 in this study. And is defined in (2).
By following our previous work [26], we set the vessel scale as = 0.6 × 2 ( −1)/2 , = 1, 2, . . . , 6, for each of the five features, to cover both large and small vessels (the diameter of vessel changes from 0.7 mm to 6.0 mm). The response of one vessel feature reaches its maximum, when the scale of the feature matches the size of local tubular structures. The maximum response can be collected as the local vessel feature.

HMM-Based Recognition Model.
The proposed MHMM is the combination of multidimensional HMMs. Before presenting the MHMM, we would like to introduce one individual HMM first. This HMM is trained by using one of the five features.
HMM is considered to be the best for modeling variations in observations of similar data series and for finding dissimilarities across different series. One HMM ( ) can be expressed as a five item array as = ( , , , , ), where is the number of invisible tissue states. In this paper, we set the number of invisible tissue states as 4 to represent four tissue states: lumen (vessel), intima (vessel boundary), pathology (pathological tissue), and adventitia (outer tissue); is the number of observation values, which are decided by the vessel features. We will describe this observation value later; is the initial state distribution = { }, = 1, . . . , . Since we start detecting lumen and intima from the axis point on the cross-section, the initial state distribution is 1 = 1, = 0, = 2, . . . , ; is the state transition matrix with size × . And is the emission matrix with size × . These two matrixes are the key parameters of HMM, which we are going to estimate in the training process.
The HMM, , describes the state transition from the inner vessel to the outside of vessel, as shown in Figure 2. This state transition is a visible random process, which gives the similarities between states and voxels. Let us denote the four tissue states (lumen, intima, pathology, and adventitia) as = { 1 , . . . , 4 }. As shown in Figure 2, we model the pathologies only existing outside the intima. The pathologies we target at are calcification, stenosis, occlusion, and low contrast vessel. Within these pathologies, calcification is a main pathophysiological process in vasculature, and the other pathological tissues come with it. More specifically, calcification determines the location of the other pathological tissues [30]. The location of calcification can be divided into superficial or deep. In the calcification process, deep calcification is followed by plaque shrinkage and arterial remodeling, in which a kind of internal elastic lamina covers the calcification plaque [30]. Based on this finding, many clinical trials consider deep or multiple calcifications as superficial calcification, for blood supply prediction [31]. Thus, we model the pathologies to be outside of the intima.
Different from other algorithms (such as SVM and AdaBoost classifier), the vessel feature values cannot be the input of HMM directly, due to the computational time issue and the convergence issue in the training process [25]. For one vessel feature, the feature values vary in a large area due to different tissues. Associating each feature value with a state will cause a heavy burden for the HMM, which may lead to slow convergence or nonconvergence. Thus, the HMM asks for quantization of feature space into a smaller number of discrete space, called observation space. This quantization is more like a binning process. Firstly, we divide the feature space into several bins. Let us denote the bin set as = Dividing method One bin Average Average One bin Number of bins 1 The method of dividing bins is the key for HMM. In our experiment, we found that the capability of HMM correlates with the number of bins. More specifically, the number of bins determines the capability of HMM in identifying different tissues, in one given area of the feature space. The larger the number of bins is, the more powerful the HMM is. Here, we take the Sato feature, for example, to demonstrate the method we use in dividing bins. Firstly, we look into the distribution of feature values along the ray (Figure 3(a)). Let us denote the distributions of lumen, intima, pathology, and adventitia as , , ℎ , and V , respectively. And the distributions of different tissues along the ray are shown in Figure 3(b). It can be easily observed that the distribution of pathology ( ℎ ) overlaps the distributions of lumen ( ) and intima ( ). Since we target delineating the vessel boundary (intima) accurately with the presence of pathologies, we want the HMM to be powerful in these two areas ( and ). Assume that ℎ overlaps and in intervals [ 1 , 2 ) and [ 2 , 3 ), respectively. We give more bins to the overlapping intervals, to make the HMM more powerful in these areas. The dividing method we use here is shown in Table 1. In this work, we set 1 and 2 as 1 = 5, 2 = 3 for all the features. The five features we use in this work have different ranges of overlapping areas, which are shown in Table 2. Please notice that the number of bins should not be large, due to the convergence issue. The computation time will be discussed in Section 4. After turning the feature series into observation series, the HMM, = ( , , , , ), can be trained by using Baum-Welch algorithm [32]. The initialization condition used for HMM training can be summarized as follows: (1) the number of invisible states = 4; (2) the number of observation values = 10; (3) the initial state distribution 1 = 1, and = 0, = 2, . . . , 4; BioMed Research International 7  (4) the state transition matrix : uniform distribution by using general principle [25]; (5) the emission matrix : uniform distribution by using general principle [25].

Multiple Hidden Markov
Model. Inspired by redundant information (multiple vessel features) [26] and good statistical properties (HMM random process), we formulate the problem of vessel segmentation as the joint segmentation of multidimensional HMMs. Since we target accurate vessel segmentation with the presence of pathologies, covering the widest possible spectrum of challenging situations is the key. We try to combine multidimensional HMMs by using multiple vessel features, as shown in Figure 4. The combination of multidimensional HMMs is the Multiple Hidden Markov Model (MHMM), Φ, which can be given as follows: where is one HMM trained by using the training subset , = 1, . . . , 5; , = 1, . . . , 5, are nonnegative terms, which give the relative strengths of the respective HMMs. These nonnegative terms are calculated in the training process, which we will describe later.
In the detecting (testing) process, after the testing voxel series = { 1 , . . . , 72 } (the rays space every 5 ∘ ) are collected on the cross-section, feature series = { 1 , . . . , 5 } can be calculated for the voxel series, by using the five proposed vessel features. With the feature series , the trained MHMM Φ can detect vessel boundary (intima) on this cross-section by calculating the posterior probability ( | Φ, ), where = { 1 , . . . , 4 } are the four tissue states. ( | Φ, ) is a similarity map between the four tissue states and the voxels in the testing voxel series, which is defined as follows: where ( | , ) is the similarity map estimated by the th HMM. These similarity maps are calculated by the forwardbackward procedures [32], in the same way as for training a HMM. A similarity map for one voxel series, calculated by ( | 1 , 1 ), is shown in Table 3 for example. As the maximums of the similarities are shown in bold, the state transition from the axis point to the outside of vessel can be easily obtained.

BioMed Research International
For each testing voxel series, we can obtain a state series by using maximum likelihood estimation as follows: Based on the state series (label series), the lumen (vessel) and intima (vessel boundary) can be segmented on the cross-section along the rays. Then along the vessel axes, we can obtain the whole vessel segmentation. Since we collect voxel series along the rays every 5 ∘ for training and testing, a given voxel can appear in multiple series and then gets different labels from different series potentially. To handle this multiple labels issue, mean similarity values of the four states are calculated for one given voxel, and the state (label) with maximum value will be signed to this voxel. In our work, we observe a few voxels have equal values in the lumen state and the intima state. Since we consider voxels in these two states as vessel finally, it does not affect the final segmentation result.
As we described in Section 2.3, multidimensional HMMs are trained separately. For feature (training subset ), the HMM ( ) can be trained by adjusting the parameters of to get the local maximum of the posterior probability ( | , ), by using the Baum-Welch algorithm [32]. Then, the trained HMM ( ) can describe the state transition from the inner vessel to the outside of vessel with the similarity map. . Finally, their combination (Eq. (6)), the proposed MHMM, can be obtained.

Results
In order to evaluate the performance of the proposed MHMM, we test it in two challenging cases: (1) synthetic vascular experiment with different levels of Gaussian noise; (2) abdominal artery segmentation with the presence of pathologies. Moreover, we compare our method with two methods to illustrate the capability of our method in delineating accurate vessel boundary with the presence of pathologies.

Two Methods for Comparison.
The two chosen methods are state-of-the-art methods: a vessel tractography method proposed by Cetin [33] and a learning-based regression method introduced by Schaap [12]. The vessel tractography method (Cetin's method) is a vessel axis-based method, which combines a vessel tensor with a vessel axis tracing method. This vessel tensor is the key for Cetin's method, which models the vessel as a cylinder. And the learning-based regression method (Schaap's method) can learn the geometry and the appearance of the vessels from annotated data. Notably, these two methods target at geometrical features of vessel. Since the pathologies do not have an explicitly geometrical feature, we believe that these two methods have the capability to identify vessels from pathologies. That is the reason we choose these two state-of-the-art methods for comparison. We implemented these two algorithms ourselves, with changes needed to utilize our clinical datasets, by following the original researches.
Cetin's method and Schaap's method ask a single seed point and an approximate centerline for initializations, respectively. With the vessel axis extracted by using our previous method, we set the start point of the vessel axis as the seed point for Cetin's method and set the vessel axis as the approximate centerline for Schaap's method.
Although Cetin's method is an unsupervised clustering method, some key parameters are estimated by using the training samples as follows: the calcium threshold = 1500; the height of the cylinder ℎ = 4 ; the number of directions for tensor fitting = 24; the number of directions for branch detection = 256; the initial radius range [0. 7 6]; the termination ratios 0 = 0.2 and 1 = 0.9; and the angle parameters 1 = 5 /3 and 2 = /9.
For the implementation of Schaap's method, several parameters are optimized by following the original optimization sequence [12]: the calcium removal parameters = 20, = 300, and = 250; the radius of the region of interest 25 voxels; the expected error of the -nearest neighbors (KNN) estimate = 0.25; the number of shape parameters = 12; the number of principal components = 32; and the number of neighbors = 16. The remaining parameters, the regression matrix̂and the shape parameters * , can be trained by using our training samples.

Evaluation Metrics and Statistical Tests.
In order to present the comparison results between the proposed MHMM and the two state-of-the-art methods, quantitative analysis is employed to quantify the performance of these three methods, by calculating the value of overlap and accuracy. Here, the manual reference segmentations (Section 2.1) are used as ground truth. Two volume and surface based metrics are borrowed from literatures: Dice Overlap Coefficient (DOC) [34] and Average symmetric Surface Distance (ASD) [35]. The values of DOC and ASD are calculated by comparing the segmented vessels with the manual reference segmentations (ground truth). For DOC (ASD), the larger (smaller) the value is, the better the segmentation result is. And the DOC (ASD) is given in percent (millimeters).
In order to acquire obvious results in comparison, we employ the paired -test, which can assess the differences in segmentation accuracy between the MHMM and the other two methods. The significant differences in comparison are marked with symbol * . And the differences with a level set at < 0.05 will be considered to be statistically significant.

Synthetic Vascular Datasets.
The synthetic vascular datasets we use here are released by Hamarneh and Jassi [36] in March 2013, which contain 120 datasets. Since the vascular trees in the synthetic vascular datasets are simulated, they can be used as ground truth segmentations. As we discussed previously, the ground truth segmentations of real clinical datasets do not exist, which makes the comparison on synthetic vascular datasets valuable. Furthermore, we simulate the challenging case, the vessel boundary obscured by pathologies, by adding different levels of Gaussian noise. Figure 5 shows the synthetic vessel tree example before and after adding Gaussian noise. All the 120 datasets available online are used in our synthetic validation experiments. Gaussian noise with different levels, 2 = 20, 2 = 40, 2 = 60, and 2 = 80, is added to the datasets to form the training datasets (240 datasets) and the testing datasets (240 datasets). 2 = 20 and 2 = 40 are used to simulate the normal noise levels in CT scans [37]. And 2 = 60 is for the low dose scans [37]. To our knowledge, noise level changes with body size, and higher noise level can be found in larger body size [38]. Based on this finding, 2 = 80 is added to the synthetic datasets, to simulate the real challenging case (low dose scans with large body size). The three methods are tested on the testing datasets. The comparison results are summarized in Figure 6.
Let us focus on the DOC results first. As we can easily observe, the three methods yield similar results, with the presence of low levels of noise ( 2 = 20 and 2 = 40). They all produce a high average DOC, as high as 92.91 ± 7.98%. When the level of noise comes to 2 = 60, the average DOC of Cetin's method and Schaap's method goes down dramatically to 85.28 ± 6.49% and 84.01 ± 7.44%, respectively. With the highest level of noise, the situations for these two methods become even worse. However, the performance of our method maintains relatively stable, with the four levels of noise. Our method yields a high average DOC of 92.37 ± 6.94%, 92.22 ± 4.35%, 90.06 ± 6.07%, and 85.08 ± 7.66%.
The results summarized in Figure 6 demonstrate that the MHMM can delineate vessel boundary accurately, when the boundary is blurred by high level of noise. Moreover, the MHMM is more resistant to Gaussian noise, compared  Figure 6: Summary of comparison evaluation on the testing datasets (240 synthetic vascular volumes). Dice overlap coefficients (DOC) using the three methods are plotted. The averages of DOC and average ASD, as well as the statistical significance, are shown in the following table, where the best results are presented in bold, and * indicates the statistically significant differences between our method (MHMM) and the other two methods at a significance level of 0.05. to Cetin's method and Schaap's method. The differences of DOC and ASD between our method and the other two are statistically significant ( < 0.05) with high levels of noise, which are denoted by the symbol * . Figure 7 shows the 3-D views of segmentation results on two randomly selected synthetic vascular trees, with the Gaussian noise of 2 = 20 and 2 = 80. It can be easily observed that the segmentations obtained by our method are more accurate than those of Cetin's and Schaap's methods.

Abdominal Artery Datasets.
Recently, there has been a trend towards minimally invasive surgery in the treatment of abdominal artery disease, which means less pain, a shorter hospital stay, and fewer complications. This minimally invasive surgery asks for accurate segmentation of abdominal artery, which can help the surgeon operate with precision, flexibility, and control. However, the extensive vessel calcification in abdominal artery disease obscures the vessel boundary (see Figure 8 for a snapshot) and makes it really challenging to segment vessel accurately. Thus, we collected our abdominal artery datasets from patients with known abdominal artery diseases. We obtained 50 abdominal CT images (120 kV, 40 mAs) by using a Philips Brilliance 64  CT scanner. For all the images, the in-plane resolution varies from 0.65 to 1.00 mm, and the slice thickness varies from 0.5 to 1.0 mm. For training and testing purpose, we randomly divide the 50 abdominal CT images into two datasets. Each of them contains 25 abdominal CT images. Please note that all the 50 abdominal CT images were segmented by two expert raters.
The comparison results of the three methods on the abdominal artery datasets are summarized in Figure 9. With the symbol * , we can easily observe that the differences between our method and the other two are statistically significant, in both DOC and ASD. The DOCs are 84.82 ± 3.91%, 85.10 ± 5.14%, and 90.15 ± 3.19% for the Cetin's, Schaap's, and our methods, respectively, while ASDs are 2.33± 1.51mm, 1.97±1.87mm, and 0.88±0.71mm. These significant differences indicate our method achieves better segmentation accuracy than the other two state-of-the-art methods.
As far as we know, the abdominal artery with the presence of pathologies should be segmented accurately for the surgery plan, since the blood supply prediction is critical for the minimally invasive surgery. Here, in order to give an intuitive impression of our method, we present the comparison results in 3-D, Figure 10. Two abdominal arteries shown in Figure 10 are with the challenges in Figure 8. We use the white rectangles to show the comparison areas with significant differences between different segmentations. In comparison area (a) (Figure 10(a)), the vessel kisses a vessel-like calcification belt. Cetin's method and Schaap's method suffer obvious oversegmentation. In the stenosis area ( Figure 10(b)), Cetin's method and Schaap's method underestimate or overestimate the radius of vessel. In the low contrast areas (Figures 10(c) and 10(d)), these two methods fail to identify vessel boundary from pathologies. However, our method segments vessel accurately in these challenging cases. These comparison results show that our method has the capability in dealing with the challenging cases caused by pathologies, while the other two methods do not. Moreover, we show several examples of 2-D crosssections with pathologies in Figure 11. In the calcification areas, the first and second rows in Figure 11, Cetin's and Schaap's methods have more oversegmentations than our method. And in the stenosis area, the third row in Figure 11,   Figure 9: Summary of accuracy evaluation on 30 carotid artery datasets for the three methods. DOCs using Cetin's, Schaap's, and our methods are plotted. Average DOC and ASD values and statistical significance are also shown just below the box plot in a table. In the table, the best results are presented in bold, and symbol * indicates the statistically significant differences between our method and the other two.
Cetin's and Schaap's methods suffer from undersegmentation and oversegmentation, respectively. It can be observed that our method produces the most similar results with the manual reference segmentations, which are considered as ground truth in this paper.

Discussion
At the very beginning, we tried to model a single HMM with multiple vessel features, by concatenating features to form a large feature vector. However, one single HMM with multiple features suffers from poor performance due to large number of parameters to learn, even nonconvergence [39]. As a result, we failed to train a single HMM with five features by using our datasets. Then, inspired by joint segmentation, we turned to MHMM, which combines multidimensional HMMs. The proposed MHMM was implemented in Matlab R2016b, with Core i7-8700 (6 cores at 3.2 GHz) and 16 GB RAM, as well as the other two methods. The computation time of the three proposed methods is shown in Table 4 for comparison.
It can be observed that the computation time of our method is much larger than that of the other two methods. Although the proposed MHMM produces accurate segmentations with the presence of pathologies, it suffers from slow convergence. That is the reason we did not apply it to large-scale scans with complex branches, such as lung datasets and brain datasets. In the future, we plan to design a cascade-AdaBoost-MHMM classifier, which organizes the AdaBoost classifiers and MHMM classifiers in a cascade way. This cascade-AdaBoost-MHMM classifier can wipe healthy vessels away, just leaving the problem areas to MHMM.
Although vessel segmentation plays a fundamental role in many clinical fields, the low segmentation accuracy issue caused by pathologies prevents it from widespread use. Several methods [10,12] have been proposed to deal with some particular challenges, such as low contrast vessel, by using minimal radius based method or area based method. However, these methods fail to identify vessels from pathologies, when facing stenosis and complicated vessel cross-sections. Thus, we propose the MHMM, combining the redundant information (multiple vessel features) and the statistical properties (HMM), to cover the widest possible spectrum of challenging situations.
The absence of a systematic evaluation workflow [40] makes the quantitative validation challenging in this work. It is challenging to compare our method with others, because of different databases, different target populations, and different quality metrics. In order to deal with this challenge and present our method in a more intuitive way, we compare our method with two state-of-the-art algorithms on the same database with the same metrics. We show the comparison

Conclusion
Targeting at vessels with the presence of pathologies, we propose a posterior probability classifier by combining multiple HMMs for accurate vessel segmentation. This classifier is an axis-based method, which works with a vessel axis + cross-section model. This model gives our method two advantages: (1) constraining the classifier around the vessel axes and (2) making the classifier focus on identifying vessels from vessel-like structures. The redundant information that comes with the multiple vessel features helps our method to delineate the obscure vessel boundary accurately. In the synthetic experiment, we use 240 complex-structured datasets (after adding noise) to evaluate the performance of our method. Our method achieves both good DOC and ASD scores, even with the highest level of noise. And the comparison experiment on real clinical datasets makes our method more promising, since the scores achieved by our method are much higher than two state-of-the-art methods.

Data Availability
(1) The synthetic vascular data used to support the findings of this study are openly available from the Vascusynth Team at http://vascusynth.cs.sfu.ca/Data.html. (2) The crotch artery data used to support the findings of this study have not been made available because the data are owned by a third party (the 2nd Affiliated Hospital of Harbin Medical University).

Conflicts of Interest
The authors declare that they have no conflicts of interest.