Automated Segmentation of Common Carotid Artery in Ultrasound Images

We propose a basis splines-based based active contour method for the segmentation of lumen boundary and media adventitia boundary from transverse and longitudinal B-mode ultrasound images. Basis-spline has M knots in the shape template and five free parameters i.e., a pair of center coordinates, scaling in the horizontal and vertical directions, and the rotation angle. The segmentation of the region of interest in ultrasound images is done by minimizing the local energy function using gradient descent technique. Further optimization is carried out using Green’s theorem. Automatic localization poses an equal importance as segmentation and is achieved using sum of absolute difference method. The result of experimental validation has been done on the Signal Processing Lab, Brno University’s database of transverse and longitudinal B-mode ultrasound images consisting of 971 and 84 images, respectively. We attained accurate segmentation of the lumen boundary from transverse and longitudinal B-mode ultrasound images with an accuracy of 99.68% and 96.98% and Dice index of 93.33% and 91.70%, respectively. In addition, we attained an accuracy of 99.29% with Dice index of 91.78% for the segmentation of media adventitia boundary from transverse B-mode ultrasound images, which is the highest among the previously proposed methods.


I. INTRODUCTION
A S of 2020, the World Health Organization (WHO) reported cardiovascular diseases to be the leading cause of death worldwide, accounting for 31% of lives lost [1]. Among which 5.5 million deaths are due to carotid atherosclerosis, commonly known as stroke. Every year 13.7 million new cases of strokes are reported with one in every four people above the age of 25 expected to suffer from a stroke in their lifespan [2]. Screening of atherosclerosis can be performed with the help of ultrasound (US) imaging of the common carotid artery, a popular, affordable and nonintrusive method [3]. In general population, patient's risk of cardiovascular events can be predicted using parameters such as artery stiffness, lumen diameter, and intima-media thickness of the arteries. Localization and segmentation of media adventitia boundary (MB) and lumen boundary (LB), as shown in Fig. 1 is important for the measurement of intima-media thickness on B-mode US images. The intimamedia thickness is regarded as an important indicator in the management and treatment of atherosclerotic risk. Carotid intima-media thickness progression is generally associated with cardiovascular risk [4], but analyzing the US images for screening of cardiovascular risks can be an arduous task for clinicians.
There are several challenges associated with the careful and comprehensive manual/automatic segmentation of carotid MB and LB in US images. The manual segmentation is laborious, time consuming, and requires a high degree of expertise in identifying and delineating carotid MB and LB. It is also challenging to identify complex atherosclerotic plaque and lesions in complex carotid arterial geometries. There exists lack of expert radiologists with extensive experience in the evaluation of carotid vessels. In addition, the B-mode US images suffer from speckle noise and low signalto-noise ratio. The quality of the image is operator dependent. Varied echo characteristics on MB and LB for different US imaging instruments makes the segmentation task difficult. There is a need for an automated technique to segment the carotid MB and LB from B-mode US images with high accuracy and robustness. The outcome of the segmentation can be used to measure the carotid wall thickness and other clinically relevant measurements such as lumen area, plaque area, and percent stenosis. This paper focuses on the important first step of carotid MB and LB segmentation.

A. OUR CONTRIBUTION
We present a fully automated Basis spline-based active contour model, which was inspired by the work of Thevenaz et al., Unser et al., Delgado-Gonzalo et al., and Pediredla et al. [5]- [10]. Automatic detection of LB, in transverse and longitudinal B-mode US images, and MB, in transverse B-mode US images, was achieved by using contrast limited adaptive histogram equalization (CLAHE)-based sum of absolute differences technique [11]. The shape specificity, ellipse in transverse and rectangle in longitudinal B-mode US images, is preserved for coarse segmentation to save computations. For fine segmentation, the converged shape will deform using B-spline-based active contour by minimizing a normalized local contrast function. Further optimization is achieved by the use of Green's theorem and gradient descent technique [12]. The use of a single method for the segmentation of LB from both transverse and longitudinal B-mode US images sets it apart from the rest of the reported techniques.

II. PRIOR WORK
Several researchers have addressed the segmentation of various layers and parameters associated with the common carotid artery such as lumen cross-section, lumen boundary, and media adventitia layer from longitudinal and transverse B-mode US images. The previous contributions can mainly be grouped into methods based on contour fit models, mor-phological operations, machine learning techniques, and clinical studies.
In the contour fit-based carotid artery segmentation category we discuss the following contributions. Mao et al. described a deformable contour model with manual initialization for the segmentation of lumen cross-section in transverse B-mode US images. Their algorithm used entropy mapping for initiation of the deformable contour and morphological operations based on the local statistical values, such as mean gray level, standard deviation of gray level and so on [13]. Loizou et al. utilized snakes technique for the segmentation of intima-media layer of the far wall of the common carotid artery in longitudinal B-mode US images [14]. Golemati et al. introduced Hough transform-based boundary extraction technique for the segmentation of arterial lumen from longitudinal and transverse B-mode US images [15]. Stoitsis et al. suggested an improvement over their previously introduced method [15] by the use of active contour-based methodology initialized by Hough transforms for arterial lumen boundary segmentation in B-mode US images [16]. Ukwatta et al. devised a two-level-set-based semi-automated method for the segmentation of LB and MB from a set of 3D US images. The two-level functions were coupled using a boundary separation-based energy [17]- [19]. Dhupia et al. developed an active oblong-based method coupled with post processing methods of median filtering and Canny edge detection for the segmentation of the LB from longitudinal mode US images [20]. Luo et al. introduced a modified level set method with probability density function-based double adaptive thresholding for semi-automatic segmentation of lumen from time-of-flight magnetic resonance angiogram images. It outperformed the traditional level set algorithm on computing time, stability, and accuracy of segmentation. The algorithm was validated on 283 cases achieving a Dice index of 0.88 ± 0.07 [21]. In later years, Yang et al. also presented an active shape model-based technique for the segmentation of media adventitia boundary and lumen boundary on transverse view slices extracted from three-dimensional US images. They validated the proposed technique on 340 US images [22]. Santos et al. presented an image contrast characteristics-based identification of lumen followed by morphological operators and anisotropic diffusion filter for noise reduction. They performed the segmentation of lumen boundary using a modified Chan-Vese level set model, where a local-based framework was used, allowing them to move the segmentation contour according to the internal energy defined in the Chan-Vese model. They performed the validation on a set of 15 B-mode US images [23]. Kumar et al. introduced a circular active disc-based formalism for the segmentation of the lumen and media adventitia boundary from transverse mode US images [24]. In later years they proposed an improved elliptical active disc-based model for lumen boundary segmentation from transverse B-mode US images, which outperformed the previously introduced model [25]. They validated their results on Brno university's signal processing (SP) lab database, which consisted of 971  US images [26]. Kumar [27], [28]. Smitha et al. developed a greedy snake segmentation method with stopping criterion and active resampling for the segmentation of intima media. They were also able to delineate the plaque well with the developed method. They also measure the intima media thickness using the modified snake algorithm and validated with the measurements obtained form an expert [29]. Nagaraj et al. proposed a support vector machine based technique for the segmentation of intima media complex. They conducted a quantitative evaluation on 49 carotid US images [30].
We discuss the following morphological operations-based methods for the segmentation of LB and MB from US images. Hamou et al. presented a boundary extraction-based model for the segmentation of the common carotid artery. Their model was aided by the use of Canny edge detection, histogram equalization, and morphological operations for the suppression of speckle noise, generally found in US images [31]. Dayem et al. described a thresholded fuzzy region growing approach for LB segmentation of the common carotid artery in B-mode US images [32]. Dayem et al. devised a watershed algorithm for the segmentation of the carotid artery, consisting of preprocessing, watershed segmentation, region merging, and boundary extraction as the major stages [33]. Yangmesh et al. described a self-adaptive histogram equalization, Canny's edge detection, non-linear filtering, and morphology-based algorithm for the segmentation of common carotid artery lumen from transverse B-mode US images [34]. Yang et al. also proposed a mathematical morphology approach and gradient vector field snake method for the segmentation of lumen and adventitia contours. They validated the method on 110 B-mode US images [35]. Zhao et al. has developed a state space framework to segment the intima-media boundary from US images throughout the cardiac cycle. They have evaluated their model by comparing it to three synthetic models and 156 subjects from four medical centers [36]. Nagaraj et al. proposed a structured random forest based technique for the segmentation of carotid artery wall in longitudinal mode US images. They conducted a quantitative evaluation on 70 carotid US images [37].
In recent year, numerous machine learning-based methods have been developed for the segmentation of carotid artery from US images. Saxena et al. discussed an active dynamic thermography-based method for the classification and screening of carotid artery stenosis. A difference in the recovery rate of the skin surface temperature is used to classify the carotid artery stenosis [38]. Meshram et al. reported standard U-Net and dilated U-Net architecture-based carotid plaque segmentation algorithm. Algorithms were validated on 101 severely stenotic carotid plaque patients. U-Net algorithm achieved a Dice index of 0.48 and 0.83 for automatic and semi-automatic segmentation, respectively. In comparison, the dilated U-Net algorithm achieved a Dice index of 0.55 and 0.84 for automatic and semi-automatic segmentation, respectively [39]. Zhou et al. proposed a deep learning-based semi-automated method for the segmentation of media-adventitia boundary and lumen-intima boundary. The media-adventitia boundary was segmented using a dynamic convolutional neural network for a pixel-by-pixel classification. The lumen-intima boundary was segmented using a U-Net model. The method was validated on 144 3-dimensional US images. They achieved a Dice index of 0.9646 and 0.9284 for media-adventitia boundary and lumenintima boundary, respectively [40]. Menchon-Lara et al. proposed a fully automatic technique based on machine learning and pattern recognition for the measurement of intimamedia thickness. The technique was validated on a set of 55 longitudinal US images of common carotid artery [41]. Benes et al. devised a non-linear model of artery-based RANSAC method for the localization of common carotid artery in longitudinal B-mode US images. They carried out the validation on 30 US images [42]. Zhou et al. has proposed a U-Net based pixel-by-pixel classification algorithm for the segmentation of carotid plaque from 3D US images. The algorithm was validated on 26 3D US images, containing 34 VOLUME 4, 2016 plaques, with a Dice similarity index of 90.72 ± 6.2% [43]. Abd-Ellah et al. have developed a computer-aided diagnosis system for common carotid artery disease diagnosis. The developed system consists of four parts namely, localization, segmentation, intima-media thickness measurement, and a two-stage classification of the common carotid artery. The localization of the common carotid artery from transverse mode US images was achieved using the faster R-CNN. They attained a Jaccard similarity index of 90.86% on a database containing 196 images [44].
Loizou et al. reported a study investigating the differences in intima-media thickness and diameter measurements of the common carotid artery. The investigation was conducted on 73 normal subjects and 83 renal failure disease patients. They carried out the statistical analysis on the obtained measurements and concluded that the mean intima-media thickness and diameter measurements were considerably higher for the renal failure disease patients compared to that of normal subjects [45]. Rizi et al. has analysed the longitudinal motion in US clinically and technically. They concluded that their results enables the use of methodological tools to quantify the longitudinal motion of the carotid artery [46].

III. PROPOSED METHODOLOGY
In this section, we develop a novel and unified Basis splinebased active contour technique for the segmentation of LB and MB from transverse and longitudinal B-mode US images. Cubic B-spline kernel is used to parameterize the concentric shape template. Inherent properties of cubic B-spline such as non-negativity, compact support, continuous secondorder derivative, and minimum curvature make it suitable for the application [47], [48]. Additionally, satisfying the Riesz basis property ensures the stability of representation between the continuous contour and the discrete coefficients used in the representation [49].
The active contour energy is a normalized contrast function, which measures the difference in average intensities between the inner contour and the annular region between the inner and outer contours. The normalization with respect to the area of the active contour ensures that tight-fitting outlines are obtained. The optimization is carried out with respect to the affine transformation parameters such as scaling, translation, and rotation.

A. MODELING OF THE SHAPE TEMPLATE AND ACTIVE CONTOUR
Let r 0 (t) = (u 0 (t), v 0 (t)) T and r 1 (t) = (u 1 (t), v 1 (t)) T be a pair of parameterized inner and outer contours, respectively. Here t is an independent variable. 0 and 1 are subscripts used to differentiate between the inner and outer contour, respectively. Kernel function ζ(t) and its integer shifted version is used to represent the function u 0 (t), u 1 (t), v 0 (t), and v 1 (t). The parameterization has been achieved in accordance with [50], and vector representation form of the results are: where and v i (t) are periodic for closed curves which leads to the following equivalent representation. where c ui,k = c ui,k+M , c vi,k = c vi,k+M , i = 0, 1, are the periodized coefficient sequences with period M (M represents the number of knots). An example of a shape template designed using cubic B-spline is shown in Fig. 2.
Outer contour (r 1 (t)) of an origin-centered convex shape can be scaled variant of the inner contour (r 0 (t)), i.e., r 1 (t) = αr 0 (t), where α is a scalar variable that is greater than 1. As a result, the outer contour coefficients get scaled up by a factor of α times the inner contour coefficients. The local contrast function's computation region and the control over the annular region width are achieved by the parameter α. The active contours are derived from the shape template in (1) according to the affine transformation below: where i = 0, 1 corresponds to the inner and outer contours, respectively, of the active contour. X 1 and X 2 are the scale parameters (horizontal and vertical direction, respectively); θ is the angle of rotation (clockwise); u c and v c are the translation parameters (center coordinates).

B. NORMALIZED LOCAL CONTRAST FUNCTION
For an active contour initialized on an image p, let ℜ 0 and ℜ 1 be the regions enclosed by the inner and outer contours, respectively. We define the active contour energy as a normalized local contrast function: where a 0 and a 1 are the areas enclosed by the inner and outer contour and E 0 and E 1 represents the energy of the inner and outer contour of the shape template, respectively. The term X 1 X 2 is the normalizing factor (net scaling in area) of the active contour. Minimizing E enables the active contour to lock on to objects that are brighter than their immediate surroundings, and maximizing it would have the opposite effect.

C. OPTIMIZATION AND PARAMETERIZATION OF B-SPLINE BASED ACTIVE CONTOUR
The optimization of the energy function E is desired for the convergence of active contour on the region of interest. We use gradient descent technique [51] to minimize the energy function. The equations associated with gradient descent are Here, P n denotes the parameter (X 1 , X 2 , θ, u c , and v c ) at iteration 'n'. The computation of partial derivatives of the energy E with respect to the free parameters X 1 , X 2 , θ, u c and v c , are essential for gradient descent-based optimization [12]. Energy E consists of double integrals, and the free parameters define the boundaries. Green's theorem [52] is used to reduce the computations required to find the partial derivatives with respect to the free parameters. It reduces the surface integral to an ordinary double integral in finding the area of a closed contour in every iteration until its convergence on the region of interest. For an image p, final expressions for partial derivatives of E with respect to the free parameters in terms of u 0 (t), v 0 (t), u 1 (t), v 1 (t) are given by

∂E ∂θ
The weighted sum of B-splines, as in (1), can be used to express the active contour template. B-spline coefficient adaptation from the partial derivatives of parameters shown in (4) is done as in (5) and (6). Convergence towards the LB and MB is accelerated due to the use of pre-computed kernels b k−l and ζ ′ (t). The partial derivatives of the parameter are bound to be finite-length sequences as kernels b k−l and ζ ′ (t) are finite.
The iterative computation of kernels Q p,1 and Q p,2 reduces the overall computations required. Recalculation of free parameters and change of the corresponding energy terms in the cost function is required after every iteration. The Locality property of B-splines means that its coefficients often remain unchanged.

D. AUTOMATIC LOCALIZATION
The localization of the B-spline contour is achieved by the sum of absolute differences method. We compute the absolute difference between the US image and the template at where b k−l (t) = ζ(t)ζ ′ (t + k − l) and Mi denotes the number of knots.
each point on the image using the following equation where A represents the US image; B represents the template; x and y, m and n represents the indices of the US image and the template, respectively. The coordinates (x c , y c ) is the point in the US image where the sum of absolute differences value is the least, that point is chosen for the localization of B-spline contour.

IV. EXPERIMENTS
The developed algorithms were validated on the Artery data set that is made publicly available by SPLab of Brno Univer-sity. The data set consists of images captured with two US devices, an Ultrasonix and a Toshiba device. They captured 50 video sequences from 15 adults using the Ultrasonix device. They set aside every 20th frame of the video sequence to form a data set of 538 US images. They also recorded 22 video sequences from the Toshiba device from two adults and one child. From these video sequences every 10th frame was included in the Toshiba data set totalling to 433 US images. Together the Artery data set consisted of 971 transverse Bmode US images. The Longitudinal mode US images were also publicly made available by SPLab of Brno University. This data set consists of 84 B-mode US images in the longitudinal section. The images were captured by a Sonix OP US scanner from 10 volunteers. About eight images were captured per volunteer   The ground truth for the segmentation of the MB from the transverse mode US images is provided by SPLab, Brno University. The ground-truth for the segmentation of LB in both transverse and longitudinal mode US images was obtained by the outlines provided by two expert radiologists with 35 and 10 years of experience. The manual annotations were done anonymously using drawing tools of ImageJ [54]. Various methods have been proposed for the segmentation of LB and MB in the past. Some methods were based on traditional image processing methods, while the recently proposed methods are mostly based on machine learning. Machine learning methods clearly have the upper hand when  the methods are well trained. There is a clear lack of data for training machine learning models for the segmentation of LB or MB from US images. Due to which the proposed method based on traditional image processing method performs very well. The proposed method has been validated on 971 transverse and 84 longitudinal B-mode US images, the highest number of images used for the validation of the method among the previously reported methods. The proposed method also segments MB and LB from transverse and longitudinal mode US images, which is a first among the previously reported methods. The ground truth provided for the segmentation was oval or circular. This can work against the proposed method as it is a amorphous segmentation method. Hence, the previously reported shape-specific segmentation methods are bound to perform better than the proposed method. Due to which the Dice index values may appear to not be the best among the reported methods. But, it can be seen that even with this handicap the proposed method out performs the previously proposed methods for the segmentation of the MB, and is comparable to the previously proposed methods for the segmentation of LB.
The term agreement-of-experts (AOE) is the similarity between the ground truth provided by the two experts. This value is the dice index obtained by comparing the ground truth of the two expert radiologists. We obtained an average AOE of 0.9628 and 0.9427 for the 84 longitudinal mode US images and 971 transverse mode US images, respectively. These values show that even two certified experts' annotation are similar to only about 94.28% to 96.27% for the obtained ground truth. The AOE value reveals that the two experts only agree with approximately 95% each other's annotations. This shows that the proposed method with its dice index scores of 91.70% and 93.33% is comparable to the experts' annotations.

V. RESULTS
The proposed method was implemented as an ImageJ [54], [55] plugin. Adaptive histogram equalization is used as a standard preprocessing step to improve the detection and segmentation. The localized and cropped carotid artery region only is then inverted and used for further processing. This reduces the computations required in further steps as the optimization deals with only cropped region of interest instead of the entire US image. The technique is validated on Brno University's SP Lab US image databases which are publicly available [26], [28], [56]- [60]. The results of LB segmentation from transverse and longitudinal B-mode US images are compared with the ground truth provided by two expert radiologists. The results of MB segmentation are compared with the ground truth provided by SPLab.
Over the years various parameters have been used to quantify the performance of a segmentation method. We have used parameters such as Sensitivity (S e ), Specificity (S p ), Accuracy (Ac), Dice (D i ), and Jaccard (J i ) index to better quantify the performance of the proposed technique Fig. 7 is used to better visualize the parameters. The red-colored circle represents ground truth and the yellow-colored circle represents the algorithm output.
The performance metrics are then computed and compared as shown in Tables  The high values in the performance metrics show that the results obtained from the proposed method have a high degree-of-agreement among the experts, even with the intraobserver variability. The distribution of Dice index throughout the database along with the agreement of experts (AOE) VOLUME   are shown in Fig.6. Fig.6 (a) and (b) shows the distribution of Dice index on comparison of the results obtained by the algorithm with the ground truth. We can conclude from the plot that the results obtained from the algorithm is accurate and reliable. There are two sets of ground truth available for the segmentation of LB, so the Dice index obtained for each set of ground truth is represented by a different colored line plot. The AOE is also provided to observe intraobserver variability. Fig.6 (c) represents the distribution of the Dice index when the results obtained from the algorithm is compared with the ground truth provided. As only one set of ground truth is available for the segmentation of MB the graph consists of just one line plot. It shows that the proposed method maintains a Dice index above 0.9 for multiple experts throughout the database.
The distribution of performance metrics for the segmentation of LB, from transverse and longitudinal, and MB from transverse B-mode US images are shown in Fig. 8. Distribution of each parameters for each individual set of ground truth vs. algorithm segmentation is provided in Fig. 8 (a)−(e). From the Box plots its evident that S e for the segmentation of LB, from longitudinal mode and transverse mode, and of MB transverse mode have values predominantly distributed above 0.80, meaning that the true positive rate of the algorithm is very high. The distribution of the S p throughout the result was very close to 1. This shows that the true negative rate of the algorithm is very high. The A c of the results are also distributed among the values above 0.95, showing that the algorithm predominantly segments the correct region. The distribution of Dice and Jaccard indices are in agreement with the statistics represented in the Tables 1, 2 and Fig.6.

VI. CONCLUSIONS
We proposed a method for the segmentation of lumen boundary and media adventitia boundary from both transverse and longitudinal B-mode ultrasound images. The proposed technique uses a B-splines based active contour for the segmentation. The computational cost involved in obtaining the best fit contour on the region of interest is optimized by the using gradient descent technique along with Green's theorem. Validation of the proposed method has been done on Brno University's SPLab transverse and longitudinal B-mode ultrasound images containing 971 and 84 images, respectively. Validation for the segmentation of lumen boundary using the proposed method, in transverse and longitudinal Bmode ultrasound images, has been carried out on the ground truth offered by two expert radiologists. Validation for the segmentation of media adventitia boundary using the proposed method was done using the publicly available ground truth provided by SPLab. The proposed method achieved an average accuracy of 99.68% and 96.98% and average Dice index of 93.33% and 91.70% for the segmentation of lumen boundary in transverse and longitudinal B-mode ultrasound images, respectively, and average accuracy of 99.29% and average Dice index of 91.78% for media adventitia boundary segmentation from transverse B-mode ultrasound images. The method portrays very competitive segmentation accuracy over the existing methods. The method also excels at being able to segment areas of interest from both transverse mode and longitudinal mode ultrasound images.