Semi-automated carotid lumen segmentation in computed tomography angiography images.

Carotid artery stenosis causes narrowing of carotid lumens and may lead to brain infarction. The purpose of this study was to develop a semi-automated method of segmenting vessel walls, surrounding tissues, and more importantly, the carotid artery lumen by contrast computed tomography angiography (CTA) images and to define the severity of stenosis and present a three-dimensional model of the carotid for visual inspection. In vivo contrast CTA images of 14 patients (7 normal subjects and 7 patients undergoing endarterectomy) were analyzed using a multi-step segmentation algorithm. This method uses graph cut followed by watershed and Hessian based shortest path method in order to extract lumen boundary correctly without being corrupted in the presence of surrounding tissues. Quantitative measurements of the proposed method were compared with those of manual delineation by independent board-certified radiologists. The results were quantitatively evaluated using spatial overlap surface distance indices. A slightly strong match was shown in terms of dice similarity coefficient (DSC) = 0.87±0.08; mean surface distance (Dmsd) = 0.32±0.32; root mean squared surface distance (Drmssd) = 0.49±0.54 and maximum surface distance (Dmax) = 2.14±2.08 between manual and automated segmentation of common, internal and external carotid arteries, carotid bifurcation and stenotic artery, respectively. Quantitative measurements showed that the proposed method has high potential to segment the carotid lumen and is robust to the changes of the lumen diameter and the shape of the stenosis area at the bifurcation site. The proposed method for CTA images provides a fast and reliable tool to quantify the severity of carotid artery stenosis.

Introduction

The carotid artery has a paired structure, namely, the left and right carotid artery and each carotid artery divides into the common (CCA), the external (ECA) and the internal carotid artery (ICA). The CCA begins in the aorta and at the neck area (bifurcation site) divides into two smaller arteries, ECA and ICA. At the circle of Willis, the ICA branches into smaller arteries that supply oxygenated blood to most of the cerebrum. Hence, the ICA plays an important role in cerebral blood circulation[1-3]. Stroke may occur when the carotid arteries become narrowed or blocked by atherosclerotic plaques or emboli[4]. Contrast enhanced computed tomography angiography (CTA) is deemed superior to magnetic resonance imaging (MRI) and ultrasound (US) for visualization and evaluation of carotid artery atherosclerosis. To diagnose the severity of stenosis, the carotid artery needs to be segmented either manually or automatically. Manual delineation of the carotid lumen is a tedious task and prone to subjectivity due to complex vessel structures, narrowing segments, significant intensity signal losses and irregularities. Unlike manual segmentation, automatic measurements are objective, more reproducible and faster. However, accurate carotid artery segmentation is a challenging task in CTA images because of the large amount of slices, variability of the size and shape, and gray levels of the artery along its path, especially in the stenosis area (or bifurcation site).

Several studies reported various lumen segmentation techniques such as using circular Hough transform[5], level-set technique[6], geometric deformable model with associated energy function[7-8], path tracking method using speed function[9] and graph-based method using the edge weighting function[10]. In this paper, a novel multi-step method with less complexity is proposed for carotid lumen segmentation. First, the mean shift technique was applied to increase lumen uniformity with minimal damage to the edges. The initial candidate points for segmentation were calculated using semi-automatic centerline extraction. Then, the watershed method was applied to reduce the segmentation cost of calculations and finally graph cut algorithm was applied.

Materials and methods

Methodology

Pre-processing

A multi-stage preprocessing technique was used to prepare CTA images for segmentation purposes (). Usually, each CTA data covers head, neck and part of the thorax. Therefore, selecting regions of interest (ROIs) reduces the computational cost and increase the accuracy of the proposed method. Since most carotid plaques occur at the bifurcation point or at ICA arteries, these locations were appointed as reference points for independent board-certified radiologists to select appropriate slices for segmentation. In this study, slices with the distance of less than 40 mm from the bifurcation point in the ICA, and 20 mm from the bifurcation point in the ECA were considered as ROI volume[11].

Due to heterogeneity of the lumen, detection of vessel boundaries is an issue using this modality. In this study, mean-shift filtering was proposed to reduce the level of heterogeneity. The mean-shift is a non-parametric iterative algorithm or a non-parametric density gradient estimation using a generalized kernel approach[12]. Mean shift defines a window level around some pixels of an image and computes the mean of data points. It then shifts the center of window to the mean and the algorithm continues until it reaches the convergence. The weighted mean of the density in the window level forn data points x, i = 1,..,n is calculated as an iterative formula:

where g is the Gaussian kernel (this function determines the weight of neighbor points for re-estimation of the mean) and h is the scaling constant[12]. The following kernel was used in this study[12]:

where C is the normalization constant, , k(x) is the simple function related to ||x||, x and x are spatial and the range factors of the feature vector, respectively, and h and h are the bandwidth of kernel: h is spatial resolution parameter effects on the smoothing and connectivity of contextual regions andh is the range resolution parameter that determines the number of segments.

Centerline extraction

Carotid arteries are surrounded by many tissue structures and veins. Therefore, for accurate and fast segmentation, selecting initial points that are located in the lumen region is necessary. For this purpose, centerline extraction was applied to detect the center positions of the artery as initial points. In this paper, centerline of the artery was extracted with three pre-defined seed points in the CCA, ECA and ICA. These seed points can be easily calculated from anatomical images without any user interaction[13] and can be applied to the entire head and neck volumetric CTA images in a few seconds.

Diagram of pre-processing techniques applied on computed tomography angiography (CTA) images.

To implement centerline extraction, first, the nearest neighbor method was used to separate carotid artery structures from the adjacent structures (). This step separates CCA-ICA artery pathway from CCA and ICA seed points. Similarly, the CCA and ECA seed points were used to separate CCA-ECA artery pathway. It should be noted that this step can only provide an estimation of artery segments and cannot separate the lumen in the narrow or nearly blocked parts of lumens. This limitation can be handled by estimation of the vessel pathway in relation to the total shape of vessel. In this study, an energy based shortest path method was used to estimate the centerline based on the total vessel shape.

The rough vessel estimation in the proposed method using three seed points.

Each path can be defined by two terminal points as CCA-ECA seed points or CCA-ICA seed points which were located at the beginning and the end of the vessel. The length of distance (C) between these points can be calculated:

where s is the length of the artery, P is the potential function and w is the constant intended fixed regulation. The optimal artery curve C(s) (with unknown length L in the 2D or 3D space) occurs when the energy of the curve reaches the minimum.

Minimizing E(C (s)) over all possible paths of C is time consuming but can be resolved using fast marching techniques[13]. In this study, the potential function (P) was obtained using a multi-scale approach analysis based on the Hessian matrix eigenvalues. The Hessian matrix for a given voxel of the image (I(x,y,z)) can be expressed as follows[14]:

The eigenvalues of the Hessian matrix for each voxel can be used as the vessel enhancement. shows the behavior of Hessian eigenvalues for different structures. Unfortunately, the value of the eigenvalues depends on the size of the objects. The variability of the size and shape of the carotid arteries along the vessel are the main challenges in the carotid lumen segmentation. A multi-scale analysis based on the scale-space theory was proposed to address these issues.

The scale-space theory introduced by Lindeberg[15] showed that a Gaussian kernel can be used to remove small objects in an image:

where L(t) is the scale-space representation at scale s, G is the Gaussian kernel and * indicates convolution operator. It is shown that the 2nd order normalized derivative of L (t) can be calculated[14]:

The following ratios were used to quantify the eigenvalues of Hessian matrix.

Although RA separates tubular shaped structures from sheet shaped structures, it cannot distinguish between tubular and a blob-shaped pattern. The RB ratio is the maximum value for blob-shaped structures and the S ratio is the lowest value for the background. The following equation was used as vesselness features for each scale:

wherea, b and c are the threshold values which control the sensitivity of the vesselness features. The potential function can be obtained by integrating the vesselness measures at different scales:

where s and s are the maximum and minimum of scale range. The centerline extraction step can be applied on both 2D Maximum Intensity Projection (MIP) maps and 3D modes. The computational cost in the 2D mode is less than that in the 3D mode, but the error (missing the correct paths) in the 2D mode is greater than in the 3D mode (). To access an accurate segmentation, centerline extraction approach was applied on the 3D mode ().

The centerline extraction in 2D mode.

The centerline extraction in 3D mode.

Segmentation

Several methods were proposed for vessel segmentation in the literatures[6-8,10,16-18]. The simple methods such as thresholding and region growing models do not need prior knowledge of the image characteristics and artery location, but these methods are not accurate, especially in the presence of the noise. In this study, a graph based segmentation method using the centerline points was implemented to segment carotid lumen from the background tissues precisely.

For graph based segmentation, each image must be mapped onto a weighted graph. Usually, each pixel of an image is mapped as a node of the graph and the weight of the edge between two nodes represents the similarity between two pixels. Labeling of the nodes can be considered for carotid image segmentation on CTA image using the graph theory[19]. The aim of this method is to assign a label set of L{0,1}to each node of the graph where 0 corresponds to the background and 1. corresponds to the carotid artery pixels. Hence, some special nodes (terminals) are used to correspond to the set of labels that can be assigned to pixels. These terminals are usually called the source, s, and the sink, t. An s-t cut, c(s; t), in a graph can be defined as a set of edges where there is no path from the source to the sink when this set is removed from the graph. Usually, the summation of the cut edge weights is called the cost of each cut.

To solve the labeling problem, one cut of s-t graph with a minimum cost must be found[19]. Unfortunately, in the graph based segmentation of the carotid arteries, this step can be time consuming because of the huge number of nodes in the volumetric data of CTA. In this study a pre-segmentation step was performed to reduce the number of the graph nodes based on the iterated region merging method with localized graph cuts using the watershed technique[20]. It converts each image into many small homogenous regions where each homogenous region was considered as a node instead of a pixel.

Following data acquisition, the extracted centerline points from the previous step were considered as source terminals for each slice. Then, at each slice, four points with a distance of as much as twice the normal diameter of the carotid arteries from the centerline point were considered as the sink terminals. Finally, an iterated watershed based graph segmentation was applied to segment the carotid lumen.

Make gold standard segmentations

Subjects in this study included non-calcified and calcified plaques and/or both types. The contour of each carotid artery was drawn manually by two independent board certified radiologists, slice by slice, using an in-house toolkit software written in MATLAB (Mathwork, version 2015) to facilitate editing of the 2D lumen segments. The ROIs excluded the border or edge of the cord (approximately 1-2 voxel from outer margin of the artery) to avoid the effects of partial volume artifacts. This toolkit software allows radiologists to review data in axial, coronal and sagittal planes and draw contours in each plane. Finally, the ground truth for each patient data was then formed by setting each pixel as the carotid segments where two radiologists reached a consensus.

Statistical analysis

shows some of volume rendering by the proposed method. The left carotid artery showed no stenosis while the others had some stenosis. It should be noted that the middle case has a lower severity of stenosis than the right case. The two most common techniques currently used to assess the performance of image segmentation techniques are spatial overlap index and surface distance metrics. The performance of the proposed method was evaluated quantitatively using the following measures compared to the ground-truth images.

Volume renderings of three segmented carotid lumens without (A) and with low (B) and high (C) value of stenosis.

Spatial overlap index

Dice similarity coefficient (DSC) represents the pixel ratio of the overlapping regions, where at any given threshold DSC values would range from 0, indicating no spatial overlap between two sets of binary segmentation results, to 1 indicating complete overlap[21].

where WG is the ground truth image and Ws is the segmented image.

Surface distance metrics

DSC is only related to the size of the contours and it does not represent the stability according to the general shape of the vessels. Therefore, surface distance metrics including mean, root mean square and maximum surfaces were presented for this purpose.

(1) Mean surface distance (Dmsd):

This Dmsd metric calculates the average distance between the obtained 3D voxel surfaces of manual segmentation and the proposed method.

Here, sdmr and sdm are the signed distance maps of the manual segmentation and the proposed method, S and S are the carotid lumen boundary surfaces obtained by manual and proposed segmentation methods, respectively, and |S| is the surface area of the i surface S.

(2) Root mean squared surface distance (Drmssd):

This metric shows the surface distance of two 3D objects by means of standard deviation of their surface voxel difference.

(3) Maximum surface distance:

This metric is similar to previous metrics except that it shows the maximum difference between the corresponding voxels of two surfaces.

Subjects

Fourteen volunteers, seven healthy subjects and seven with carotid artery disease, with the age of 60±10 (mean±standard deviation) and range of 45 to 80 years, were recruited from Tehran Heart Center Hospital (). Subjects provided written informed consent and the institutional review board approved the protocol. All patients received Visipaque (Iodixanol, 320 mg/mL), a contrast agent material. It contains 0.044 mg calcium chloride dihydrate per mL and 1.11 mg sodium chloride per mL, with a sodium/calcium ratio equivalent to blood.

Imaging

The CTA data were acquired on a 128 detector row Siemens CT scanner (SOMATOM Definition Flash, Siemens Medical Systems, Germany) with standard parameters of Carotid Angio (Adult) vascular protocol. The axial images of CTA were acquired with the following parameters: tube current= 120 mAs, tube voltage= 120 kVp, matrix size= 512 × 512, filter type= flat, focal spot= 1.2 mm, spiral pitch factor= 1.2, slice thickness= 0.6 mm, scan time= 10-14 s (depending on the individual patient's size and anatomy) and the patient position was head first-supine (HFS).

Results

The difference on the number of voxels between ground truth and segmentation volumes is depicted in . shows the results for DSC and surface distance metrics on each patient.

The difference on the number of voxels between ground truth and segmentation data

illustrates the mean and standard division for DSC, Dmsd, Drmssd and Dmax metrics. As mentioned before, the number of voxels and DSC (or any volume based metrics) depends on the size of the counters and can be changed by applying morphological operations (dilation and erosion) while the surface distance metrics are resistant to the variation of the contour size. In the subject encoded 13, the values of the volume based metrics showed some reduction while surface distance metrics showed a significant difference between the segmentation result and the ground truth data (). In this subject, a portion of the neighboring vein was detected as lumen, which was marked by a yellow arrow in . The green contours showed the segmented objects. This was viewed as an error, which can be corrected by an extra preprocessing method.

Segmentation result for case No. 13 of data (A) ground truth MIP map and (B) MIP map of segmentation by the proposed method (C) and one slice of original image (D) segmentation area in the same slice.

Discussion

In this study, a multi-step method was proposed for carotid artery segmentation in the CTA data. The proposed method used three seed points that can be obtained by anatomical markers[13]. Mean shift smoothing was used as a preprocessing step to increase uniformity of the lumen regions with minimal damage to its edge. Then, centerlines of arteries were calculated by a multi-scale 3D Hessian based shortest path method. The centerline points were used as terminals for graph cut based segmentation method. To reduce the complexity of graph cut segmentation, a watershed pre-segmentation was performed before graph cut segmentation. The 3D CTA data of 14 patients were used for evaluation of the proposed method. The results showed that the method has good performance for carotid artery segmentation, which assures that it can be used in clinical case segmentation.

Fortunately, due to an energy based centerline extraction method, the algorithm does not have a dependence on seed point selection and the user can select seed points at any part of the arteries. However, selection seed points at the beginning and final slices of ROI was preferred because plaques or narrowing usually do not occur in these slices.

This method is able to detect and segment even small branches of the carotid artery. Though in this study a refinement stage branch in the centerline extraction step has been considered to neglect these small vessels, it may be still effective for segmentation and rendering of vessels in cerebrovascular studies ().

The refinement stage in centerline extraction (A) estimated centerline with extra branch (B), and centerline after, refinement step (C) segmentation result

In the carotid atherosclerotic plaque studies, the bifurcation site and ICA artery were given special attention[22-26]. Most of the stenosis plaques occur at the bifurcation site due to the pressure changes of the blood flow and in the ICA because it is thinner than ECA. The proposed method can correctly segment the carotid arteries at the bifurcation site even if the stenosis plaque occurs ().

Stenosis plaques near bifurcation site of carotid artery.

The method can estimate the centerline in the cases with small total occlusion along the vessel path precisely (). However, if there is a large total occlusion along the vessel, it seems that the method would not be able to estimate the centerline of the entire artery, possibly because the method estimates central points incorrectly. If a calculated centerline point belongs to the other structures or background positions, the segmentation step may separate non-lumen parts of images as lumens. To solve this problem, it is suggested to add a refinement step for centerline points. In this step, centerline points that have no intensities in their positions in the ranges of the normal artery intensities should be removed from centerline set points.

shows a comparison between the results of the proposed method and those of two other methods[27-28]. Although the dice similarity index in the proposed method is less than that of other methods, it shows better values for surface distance metrics. The overall shape of the vessels will be very close to manual segmentation.

To have an analysis on the spent time, the method can be divided into the three main stages that include pre-processing, centerline extraction and graph cut stages. The centerline extraction stage requires the lowest time among these stages, because this stage only uses the lower objects (from labeling and nearest neighbor) instead of all of image pixels. The preprocessing step has the highest running time due to the use of mean shift smoothing. One approach to shortening the average preprocessing time is to reduce the mean shift range resolution in every single slice, but it may cause heterogeneity which decreases the performance of the method.

In conclusion, in this paper, a novel interactive tool was proposed for the accurate segmentation of the carotid artery in the volumetric CTA images. The experimental results showed that the proposed method has a good ability to segment the narrowed carotid arteries. It is simple to use and its user does not need to adjust any parameter and the selection of three seed points is enough. These features make it an appropriate method for clinical use as an alternative to the manual contouring. In addition, it seems that it can be used for detection and quantification of carotid stenosis in the CAD systems. Our further work is to evaluate the proposed method in a larger data set and to optimize it into a tool for detection, quantification of the stenosis and plaque segmentation.

Tab.1

The behavior of Hessian eigenvalues for some 3D models[16].

structure orientation	e1	e2	e3
noise (no preferred structure)	L	L	L
bright sheet-like structure	L	L	H–
dark sheet-like structure	L	L	H+
bright tubular structure	L	H–	H–
dark tubular structure	L	H+	H+
bright blob-like structure	H–	H–	H–
dark blob-like structure	H+	H+	H+

L=Low, H– = High with negative value, H+ = High with positive value.

Tab.2

Specification of patients that were used to evaluate the proposed method

Patient ID	Age	Carotid stenosis
1	50	No
2	45	Low stenosis in ICA
3	52	No
4	63	No
5	65	High stenosis in ICA, near the bifurcation site
6	51	No
7	67	High stenosis in ICA, near the bifurcation site
8	70	High stenosis in ICA
9	63	High stenosis in ICA
10	76	High stenosis in ICA
11	59	No
12	78	High stenosis in ICA, near the bifurcation site
13	49	No
14	61	No

ICA=internal carotid artery.

Tab.3

Validation results for the proposed method.

Patient ID	DSC	Dmsd (mm)	Drmssd (mm)	Dmax (mm)
1	0.89	0.27	0.39	2.22
2	0.91	0.22	0.29	1.7
3	0.91	0.22	0.3	1.3
4	0.89	0.25	0.43	1.92
5	0.90	0.32	0.52	3.03
6	0.85	0.21	0.29	0.87
7	0.91	0.14	0.2	1.25
8	0.93	0.21	0.28	1.12
9	0.9	0.17	0.24	1.78
10	0.83	0.26	0.34	0.87
11	0.82	0.31	0.58	3.07
12	0.88	0.32	0.40	1.27
13	0.61	1.43	2.32	8.9
14	0.90	0.2	0.26	0.71

DSC= dice similarity coefficient, Dmsd=mean surface distance, Drmssd= root mean squared surface distance, Dmax= maximum surface distance.

Tab.4

Mean and standard division of validation metrics.

DSC	Dmsd	Drmssd	Dmax
average	average	average	average	std	std	std	std
0.87	0.32	0.49	2.14	2.08	0.54	0.32	0.08

DSC= dice similarity coefficient, Dmsd= mean surface distance, Drmssd= root mean squared surface distance, Dmax= maximum surface distance.

Tab.5

Comparison of the average value of metrics between the proposed method and two other methods

Methods	DSC	D msd	D rmssd	D max
Method 1 [29]	0.80	0.86	1.57	6.1
Method 2 [30]	0.81	0.83	1.54	5.98
Proposed method	0.87	0.35	0.51	2.17

DSC= dice similarity coefficient, Dmsd= mean surface distance, Drmssd= root mean squared surface distance, Dmax= maximum surface distance.