Lumen diameter (LD) and Intima-Media thickness measured using automatic and manual methods (n=30). The difference (), the coefficient of variability () and the correlation () between both measurements are also presented.
Ultrasound imaging is widely used to assess carotid, brachial, femoral, as well as other arteries. There are major advantages of using ultrasound in comparison to other imaging techniques, such as its non invasiveness and its capability to produce real-time visualization of the arterial lumen and vessel wall that is not possible with any other imaging modality . Recent clinical studies have benefited from continuous improvements in ultrasound image quality, new imaging techniques and signal processing algorithms with the aim of identifying the vulnerable carotid plaque based on the mechanical wall motion behavior [2,3].
The vulnerable arterial plaque may cause atherothrombotic events, myocardial infarction and stroke, which are responsible for approximately 35% of the total mortality in the western world, and are the leading causes of morbidity world-wide . The first indication of cardiovascular disease is a thickening of the intimal and medial layers of the arterial wall. It involves lipid accumulation and the migration and proliferation of many cells in the sub-intimal and medial layers, which results in the formation of plaques. It is the rupture of such plaques that causes myocardial infarcts, cerebrovascular events, peripheral vascular disease and kidney infarcts. The impact of the intima-media thickness (IMT) on the incidence of cardiovascular events in the Rotterdam study by B-mode ultrasound indicates that the risk of myocardial infarction increases 43% per standard deviation increase (0.163 mm) in common carotid IMT . The main conclusions resulting from this study were supported by other independent investigations which reveal that an IMT higher than 0.9-1.0 mm indicates a potential atherosclerotic disease, which translates into an increased risk of a cardiovascular event. Hence, the robust segmentation and measurement of IMT by B-mode ultrasound has a considerable impact in the early diagnosis of atherosclerosis, prognosis prediction, and in the monitoring of responses to lifestyle and prescribed pharmacological treatments.
Ultrasound signals have also been extensively used in clinical sites, by exploiting Doppler effect to measure vascular blood velocity and flow, among other applications [5,6]. Typically, a spectrum of frequencies related to the different velocities of the blood cells is presented as a curve of velocity versus time. The analysis of this curve can reveal important relationships between the frequency spectral pattern along the cardiac cycle and the presence of cardiovascular diseases [7,8], among other examples [9,10].
2. Arterial vessels image analysis
2.1. Problem statement
The current clinical practice in the assessment of the early cardiovascular diseases involves acquisition of B-mode ultrasound data from large superficial arterial vessels such as the common carotid artery (CCA). The image acquisition process generates sequences of two dimensional ultrasound images along the time (2D+time) that are currently interpreted using either manual annotation procedures or commercially available semi-automatic image processing environments. While the manual annotation generally results in accurate IMT measurements, it is subject to intra and inter-observer variability. Moreover, this procedure requires annotating multi-frame ultrasound data, which is not only labor intensive but also highly dependent on the experience of the observer. All these factors stimulated the investigation of automatic segmentation techniques, which can greatly support the clinical practitioners in their evaluation and may have substantial benefits in the quality of the medical act. As a consequence of this clinical interest, a large number of studies were focused on the development of automatic IMT segmentation algorithms in order to provide an accurate analysis of IMT measurements [11-31].
The IMT complex is best visualized in longitudinal sections of the CCA. Fig. 1 shows a representative B-mode ultrasound image of the CCA and a schematic illustration of the relevant leading edges of echo responses. Previous studies [10-12] have shown that the leading edges can be mapped to the following interfaces: near-wall media-adventitia, far-wall lumen-intima and far-wall media-adventitia. The lumen diameter (LD) is defined as the distance between the intima-lumen interface of the near-wall and the lumen-intima interface of the far-wall. The far-wall IMT is defined as the distance between the far-wall lumen-intima (FWLI) and the far-wall media-adventitia interfaces (FWMA).
The determination of ultrasonic measurement of the artery becomes equivalent to accurately detecting the echo boundaries presented in Fig. 1. However, the existence of ultrasonic imaging artifacts such as speckle, reverberations and dropouts make the accurate definition of a boundary very difficult.
2.2. IMT segmentation algorithms
Since the IMT complex is defined by two distinguishable interfaces in the ultrasound image data, the majority of studies that addressed the IMT segmentation were built on several assumptions regarding the intensity profiles associated with each interface. The automatic detection of such interfaces requires a priori knowledge that implies a certain amount of user intervention. In recent studies, substantial efforts have been dedicated to reduce the level of user intervention during the IMT segmentation process. A review of the methods in this area is presented by Molinari et al. , which address the main directions of research in the field of IMT segmentation. The published methods in this area can be classified in three classes: edge-based, dynamic programming and probabilistic IMT segmentation methods. The edge based segmentation schemes aim to reconstruct the IMT complex from gradient data and in this process prior knowledge relating to the intensity profiles is enforced in the process of selecting the FWLI and FWMA interfaces. The early segmentation schemes that addressed the identification and measurement of the IMT were based on semi-automatic methods where the user intervention was critical to obtain accurate results. De Groot et al.  analyzed the contribution of the patient and observer variability in the process of the IMT measurement in serial carotid and femoral B-mode ultrasound scans. Through the use of variance components analysis they found that 75% of the variance in the measurement of the far-wall thickness could be attributed to the differences among patients and a 7% variability recorded in the analysis of the mean IMT thickness was caused by the ultrasound equipment. In a study presented by Selzer et al. , the initial IMT boundary position was manually selected and this information was used to guide an IMT edge-based segmentation process. Dwyer et al.  developed a semi-automatic IMT segmentation algorithm that was applied in B-mode ultrasound images. In the algorithm the average distance between the FWLI and FWMA interfaces was used to approximate the IMT. A selection of the frames used in this study was carried out by the user to ensure that there would be a clear identification of the interfaces.
Liguori et al.  proposed a multi-step semi-automatic IMT segmentation algorithm that has been developed for the analysis of single frames B-mode ultrasound data. The first step of their method involves the manual selection of the region of interest (ROI) followed by a threshold procedure that is applied to set all pixels with intensity values lower than a pre-defined threshold. The IMT detection entails an analysis of the intensity profiles associated with the gradient data under the assumption that all pixels corresponding to the lumen are anechoic and the image areas that define the tunica intima and tunica adventitia are the most reflective arterial layers.
The two IMT interfaces are selected by analyzing the strength of the gradient in the direction of the ultrasound beam, the FWLI interface corresponds to the first relative maximum, while the FWMA is given by the second one. A similar approach was employed in [18,19], where the authors also analyzed the pixel intensity profiles to detect the salient intensity transitions that are characteristic for the two IMT interfaces. Ilea et al.  adopted a multiscale approach that embeds a statistical shape model with the aim of identifying the two interfaces that form the IMT without any user intervention. The developed algorithm was validated on 49 single frame B-mode ultrasound images and results were compared against manually annotated data. Delsanto et al.  implemented a hybrid algorithm where active contours were applied to refine the initial FWLI and FWMA estimates. The reported results indicate the efficiency of this approach in reducing the level of outliers, but several problems that are caused by the gaps in the IMT structure started to surface. The assumption that in longitudinal images of the carotid artery the IMT is defined by a pair of active contours was used in subsequent contributions [22-26]. However, all contributions require a user interaction in the process of initializing the active contours.
Gutierrez et al.  proposed a different semi-automatic active contour-based IMT segmentation algorithm where the edge information is extracted using a multiresolution approach. The major objective of this paper was to measure the lumen diameter and the IMT and in their experiments the authors assessed the performance of their semi-automatic algorithm against the manually segmented ultrasound data using metrics such as the coefficient of variability and correlation.
Dynamic programming algorithms were proposed [28,29] as a computationally efficient alternative to the standard heuristic search methods when applied in the context of boundary tracing, and due to their intrinsic properties they positioned as an attractive approach for IMT segmentation. An example is represented by the work of Liang et al.  where the authors addressed the IMT segmentation problem by adopting a multiscale approach. Rocha et al.  applied a related segmentation scheme to more challenging ultrasound images that exhibit arterial plaques. Their approach starts with the detection of the media adventitia layer by searching for the best fit of a cubic spline to the edge data by taking into account the anatomical characteristics of the adventitia. This is followed by the segmentation of the lumen boundary by applying a hybrid dynamic programming-based active contour technique.
A distinct category of algorithms relied on probabilistic schemes to identify the IMT interfaces such as the work of Destrempes et al. . Their algorithm was developed based on the assumption that the echogenicity of the region of interest where the IMT is located can be modeled using a mixture of three Nakagami distributions and the parameters of the distributions are estimated using an expectation maximization (EM) algorithm. The proposed method proved accurate but it requires user interaction for the ROI initialization.
The discussion will be continued with a technical presentation of an artery boundary segmentation method to measure lumen diameter and IMT.
2.3. Artery boundary enhancement
One of the first steps in to segment artery leading edges is to enhance such interfaces from B-mode ultrasound images. To enhance border detection accuracy, a multiscale border identification can be implemented using filters in the form of scaled convolution operators [32,33]. The scale space of an image is constructed through convolution of the image with a two-dimensional (2D) Gaussian density kernel with zero mean and standard deviation:
where denotes the dimension of the input domain. A blurred replica of the original image is obtained by convolution with for a specific . The stack of images as a function of increasing scale parameter is coined a linear scale space. Hence, as increases the detailed object structures vanish while gross structures persist. Fig.2 shows the 2D gradient magnitude calculated for a carotid vessel image (B-mode ultrasound) in three different scales.
Based on these features a scaled artery image is used to identify the approximated position of the near and far walls. Two complementary images are obtained based on the gradient value in y-direction: one that enhances pixel values transitions from high to low echoes, such as edges encountered in near wall tissue interfaces, and other that enhances pixel values transitions from low to high echoes (such as edges encountered in far wall tissue interfaces). Fig. 3 shows the boundary enhancement of the near and far wall.
2.4. Contour modeling
The contour of each wall can be modeled following the Geometrically Deformed Model proposed by Lobregt and Viergever . In this model, a set of vertices connected by straight line segments or edges forms the basic contour structure (Fig. 4).
In Fig. 4, the position of a vertex is represented by a vector , and the edge between and by a vector . The contour deformation is caused by a combination of forces which act on the vertices. The resulting acceleration in vertex is denoted by a vector .
The contour local curvature at a vertex is defined as the difference between the directions of the two edge segments that join at that location:
The local tangential unit vector is defined as the normalized sum of the unit vectors of two joining edge segments:
The local radial direction at a vertex is obtained from by a rotation over radians:
2.5. Dynamic force formulation
In the model definition, the dynamic in each vertex must satisfy the Newton´s second law,
Where is a coefficient that has a mass unit, and and are the damping (or viscous), the internal and the external forces, respectively.
The internal force can be estimated from the local contour curvature along the local r-axis :
The external force acting in each vertex can be approximated by some image feature. In this paper we used the information obtained from the local image gradient as the external force.
The damping force is proportional to the velocity of the vertex and points in opposite direction:
The total force acting on a vertex is a weighted combination of damping, internal and external forces :
where and are the weighting factors.
The deformation process over the contour is implemented as a numerical time integration process in which the complete state of the contour is calculated at a sequence of discrete positions in time . A set of state equations controls the deformation process in terms of position, velocity and acceleration of each vertex on the contour:
Where , and define the position, velocity and acceleration, respectively, of the vertex in a incremental time . The vertex mass, , is setting constant for all vertices and the resulting force, , is calculated using equation (8).
Fig. 5 shows the result of the deformation of the active contours during the process of detection the FWLI and FWMA in two patients in systole and diastole. The vessel’s lumen diameter and IMT can be obtained by automatic measurement of the distance between segments estimated from the linear regression of the contours in each interface.
2.6. Vascular blood velocity and flow
As presented in the previous sections, B-mode ultrasound is capable of reliably and accurately imaging peripheral arteries and can be used for vessel diameter measurement. Ultrasound signals have also been extensively used in clinical sites by exploiting Doppler effect to measure vascular blood velocity and flow.
The Doppler effect refers to an increase or decrease in the frequency of a wave that is traveling toward or away from the observer, respectively. This principle is applied in Doppler ultrasound to measure the direction and velocity of blood flow in the vessels. The relationship between the Doppler shift and blood flow velocity is:
where is the speed of sound in blood (1540 m/s), is the angle between the ultrasound beam and the direction of blood flow, and and are the frequencies of the transmitted and returned signals, respectively. Due to the angle dependency, alignment of the ultrasound beam parallel with the direction of blood flow is essential to obtain accurate flow estimates.
In vascular studies using commercial ultrasound equipment, a spectrum of frequencies related to the different velocities of the blood cells is presented as a curve of velocity versus time (Fig. 6). This information can reveal important relationships between the frequency spectral pattern along the cardiac cycle and the presence of cardiovascular diseases [7-10].
However, commercial ultrasound systems are primarily dedicated to get instantaneous data for individual patients, and they have, in general, low flexibility to perform large-scale researches. Thus, to make this kind of study easier in clinical protocols involving hundreds of patients, computational tools have to be developed to extract quantitative data from spectral display of Doppler ultrasound images.
After two steps defined by the user: calibration and selection of the region of interest (ROI), a Gaussian filter (=1 pixel, precision 90%) can be applied to the grayscale input image to attenuate high frequency noise.
The detection process of the time axis (‘X’) considers the smallest variation of the pixels intensities occurs in the horizontal direction. Thus, equations (11), (12) and (13) calculates the ordinate ‘y’ expected for the axis ‘X’, which will be the reference (0 m/s) to the velocity calculation.
where I(i, j) is the image intensity (grayscale level) at i and j coordinates; m is the image width (in pixels); ymin and ymax are, respectively, the ordinates of the superior and inferior lines that delimit the rectangular ROI defined by the user; z is an empirical pre-defined threshold to reject the graphic background;
The image can be transformed to a binary image depending on a threshold level that can be defined by the user. After image binarization, a median filter (size: 3 x 3 pixels) can be applied to edge smoothing and spurious suppression.
An envelope detection step can be initialized with horizontal lines at the top and at the bottom of the ROI. Each point of these lines is moved down or up to the border of the binary curve. Then, the algorithm holds either the superior or the inferior contour (Fig. 7), assuming that, in general, the desired one has higher amplitude variation.
Finally, after automatic detection of the peaks (Fig. 7), the algorithm can compute the mean peak velocity (), the mean envelope velocity () and the velocity time integral (), according to the equations (14), (15) and (16), respectively.
where is the total number of peaks detected in the curve; is the total number of pixels in the curve; and are the calibrations in pixels for the time and the velocity axes, respectively; is the amplitude in pixels for each point in the curve.
In addition, if B-Mode images are available, an arterial wall interface detection module can determine the vessel diameter and the mean peak blood flow () and the mean blood flow () can be estimated using (16) and (17).
3. Clinical applications
3.1. Automatic IMT measurement
The method presented in the previous section can be used in clinical applications to assess large sequences of 2D+time B-mode ultrasound images of the CCA. We performed an experiment that comprised the analysis of 180 images from 30 patients (3 images in diastole and 3 images in systole for each patient), all of which included manually defined interfaces for reference. The minimum and maximum artery diameters were measured for each patient using the manual and the automatic procedure.
In order to study the variability between the automatic and manual definition of artery boundaries, the pooled mean, , and the standard deviation, , for the difference between automated and manual measurements of lumen diameter were computed. The coefficient of variation () was calculated according equation (19).
The strength of the relationship between automated and manual methods is indicated by the correlation, , between the two measurements:
where is the covariance between the automated and manual measurements. and are the standard deviation of automated and manual measurements, respectively.
Means () and standard deviations () for the differences between the automatic and manual methods were calculated for the population (n=30). The coefficients of variability () and the correlation () between both methods were also obtained.
The results obtained for the parameters , , and are summarized in Table 1.
|Intima Media Thicknes|
A comparative analysis between commercial ultrasound systems operated by specialists and the method presented in the previous section to measure blood velocity and flow automatically can be performed. A new experiment was performed using systolic mean peak velocities (102 samples) and velocity time integrals (75 samples) of common carotid and brachial arteries under basal condition, brachial arteries in the reactive hyperemic response and echocardiographic exams. Table 2 shows the number of images and samples used in this analysis.
|Artery||Number of images||Number of samples|
|Peak velocity||Velocity time integral|
|Common carotid arteries||30||39||31|
|Brachial arteries under basal condition||23||35||24|
|Brachial arteries in the reactive hyperemic response||10||15||15|
According to the procedure described in the previous section, the peak velocities measured from the carotid arteries were negative (average: -0.59 m/s), while from the brachial arteries under basal condition, as well as in the reactive hyperemic response, the peak velocities were positive (averages: 0.63 m/s and 1.18 m/s, respectively). In the echocardiographic images the measurements were either positive or negative and the average of the absolute values was 1.48 m/s.
Similarly, positive or negative velocity time integrals were dependent on the exam type. However, the number of cardiac cycles used to get these measurements was not standardized, leading to a range of the samples, from -150 to 91 cm.
Fig. 8 shows Bland-Altman’s  and Linear Regression graphics for the systolic peak velocity analysis, where ‘A’ refers to the measurements done with a commercial ultrasound system and ‘B’ refers to the methodology described in section 2.6. Bias, standard-deviation, correlation coefficient, and linear equation results are presented in Table 3. Like peak velocity, Figure 9 and Table 4 were obtained for velocity time integral analysis.
|Systolic Peak Velocity (N=102)|
|Bias (m/s)||sd (m/s)||rAB|
|Linear Regression Equation|
|0.02||0.05||0.9985||A = 0.9938*B – 0.0190|
|Bias (cm)||sd (cm)||rAB|
|Linear Regression Equation|
|1.25||3.86||0.9984||A = 1.030*B – 0.9287|
Measurements of lumen diameter (LD) and intima-media thickness (IMT) of carotid, brachial and femoral arteries from B-mode ultrasound are defined as the average distance of interfaces between vessel tissue layers. In order to determine the interfaces location, manual tracing is commonly used. However, this approach is a time consuming procedure and based on subjective operator assessment. Besides, it inevitably results in inter and intra-observer variability due to the complex nature of the echogenic zones, especially at the lumen-intima interface, which frequently present weak echoes, echo dropouts and irregularities caused by scattering.
In this chapter, we have reviewed some methods for automatic or semi-automatic interface detection and presented in detail an approach that uses the active contour technique. In this technique, external forces are proportional to the local image gradient obtained from a multiscale analysis. The automated measurements, when compared to those obtained by manual tracing, are equally accurate and the coefficients of variability between both methods are below 1,0% for Lumen Diameter and 6,5% for IMT thickness measurements.
Vascular blood velocity can be measured by using Doppler effect, and if lumen diameter is available the blood flow can also be estimated. However, commercial ultrasound systems are primarily dedicated to get instantaneous data for individual patients, and they have, in general, low flexibility to perform large-scale researches. Thus, to make this kind of study easier in clinical protocols involving hundreds of patients, computational tools have to be developed to extract quantitative data from spectral display of Doppler ultrasound images.
We briefly presented a methodology to extract automatically blood velocity and flow from Doppler ultrasound images that permits extensive clinical studies. The small bias and high correlation for both, peak velocity and VTI, indicate the reliability of this methodology and these findings are better than those presented by Tschirren et al.  (bias: 0.40 m/s for peak velocity and 7 cm for VTI), though their results refers to a dilatation study of the brachial artery, where data values were about ten times higher than the present study.
It is important to note that for VTI statistics shown in Table 4, the threshold used to get the binary image was 60, instead of the default 40 used to extract the systolic peak velocity. This change was motivated by the higher bias (1.70 cm) and standard deviation (6.78 cm) obtained with the default value for VTI. Despite these numerical results, it is not possible to conclude that the threshold of 60 is more appropriate than 40, since the human operation to get measurements using the ultrasound equipment may also be subject to systematic errors and deviations. For instance, visual results showing the envelopes drawn on the blood velocity graphics point that, by using the proposed methodology (Fig. 7), the envelope line is more refined than that obtained by manual operation of an ultrasound system (Fig. 6). In the latter case the user does not notice or simply disregards small image brightness variations, which means that this procedure is highly dependent on the user’s subjective evaluation and it is hardly reproducible.
By processing a diversity of common carotid, brachial and echocardiographic Doppler image samples, collected from four different commercial ultrasound systems, the proposed methodology to measure peak velocity and VTI was validated by the Bland-Altman’s analysis and by the correlation coefficient. Visual analysis also confirmed that the spectrum envelope detection is very satisfactory.
The methodology was implemented in a user friendly graphical interface that has a semi-automatic characteristic. The delivery of this tool was intended to help clinicians in their studies based on Doppler ultrasound images with the following advantages: to save operational time, to lower subjective results, and to support measurement reproducibility.
However, ultrasound is still an observer-dependent modality in which the image quality depends on an experienced observer, appropriated technique and equipment. Automated systems and algorithms which can improve measurement accuracy and reproducibility, without the observers input, still remains an open area of study and research.
This work was supported in part by the National Institute of Science and Technology—Medicine Assisted by Scientific Computing INCT MACC, and the Zerbini Foundation.