Automated Detection of the Arterial Inner Walls of the Common Carotid Artery Based on Dynamic B-Mode Signals

In this paper we propose a novel scheme able to automatically detect the intima and adventitia of both near and far walls of the common carotid artery in dynamic B-mode RF (radiofrequency) image sequences, with and without plaques. Via this automated system the lumen diameter changes along the heart cycle can be detected. Three image sequences have been tested and all results are compared to manual tracings made by two professional experts. The average errors for near and far wall detection are 0.058 mm and 0.067 mm, respectively. This system is able to analyze arterial plaques dynamically which is impossible to do manually due to the tremendous human workload involved.


Introduction
Carotid atherosclerosis may appear as intima-media thickening or plaque formation. Common carotid artery intima-media thickness (CCA-IMT) measurements have been confirmed to be an early OPEN ACCESS marker of atherosclerosis [1] and have been associated with a higher risk of stroke [2] and myocardial infarction [3]. The non-invasive sonographic examination has demonstrated its potential in early prediction of cardiovascular diseases [4]. Carotid plaque formation is a different form of atherosclerotic damage of the carotid artery. An increasing grade of stenosis caused by plaques is associated with higher incidence of stroke and major cardiovascular events. Further, total plaque area is a major determinant of cardiovascular events and mortality.
B-mode sonography is widely used to scan the carotid artery in the diagnosis of atherosclerosis. It is noninvasive and low cost. The carotid IMT can be measured either manually [5], semi-automatically or automatically [6][7][8][9][10][11][12][13][14]. Carotid plaques are usually measured manually by measurement of the largest diameter in length and height, respectively. This can only be seen as a rough estimation of the plaque area. Measurement of stenosis degree is mostly performed via few points placed manually at the nearest point of the stenosis.
The automatic methods of outlining the outer boundaries of IMT and plaques have the potential to make results reproducible and eliminate the strong variations observed in manual tracings by different observers. Particularly, it is almost impractical to trace the boundaries of the inner wall on image sequences manually since this requires a tremendous amount of time. By using an automatic system the processing time can be considerably reduced. The motivation of this study was to develop a fiable system which is able to automatically detect the intima and adventitia of both near and far artery walls along a section, with or without plaques, using dynamic B-mode sonographic image sequences even under strong speckle noise conditions. This system can identify not only the IMT but also the lumen diameter (LD) during systolic and diastolic cycles, from which the artery's volume change with respect to time can be outlined. A precise measurement of the IMT over the heart cycle in areas without plaque formation is the basis for the measurement of the local arterial stiffness. The same measurement in the area of a plaque offers the possibility of measuring the deformation of a plaque over the heart cycle, giving the perspective to analyze the deformability of a plaque. This may be a step forward of the detection of vulnerable plaques. Via this system, the dynamic processes of the carotid artery (CA) can be represented by some parameters such as IMT variation, LD variation, and the frequency analysis on the LD changes.
The carotid artery stiffness (or elasticity, CS) is one of the important factors for predicting cardio-vascular events [15,16]. In order to measure CS one has to calculate the cross-sectional area changes of the CCA during the systolic and diastolic cycle, i.e., CDist = ∆A/(A•∆P), where A is diastolic lumen area and ∆A is the maximum area change during a systolic-diastolic cycle, and ∆P is the local blood pressure change measured by a sphygmometer on branchial artery. The CDist can be converted to CS by giving CS = (CDist) −½ [15]. Our method is an automated scheme which is able to calculate the diameter changes of CCA during the heart beat cycle. It provides direct information for calculating ∆A and A.
In previous studies CS was calculated by measuring the LDs in systolic and diastolic states, i.e., only two points. The contribution of this study is to propose a fiable technique for detecting the sectional IMT and LD changes along a section of CCA with respect to time, with and without plaques and even having speckle noises. This method is in general different from previously published works and it provides good results. Technically, we use a novel dual dynamic programming combined with some anatomic knowledge incorporated into the scheme which makes the detection more robust towards the presence of speckle noises.

Materials
After at least 15 minutes of rest in the supine position, ultrasonic examinations of the right and left CCA were performed. An Esaote Mylab 25 (Florence, Italy) ultrasound scanner with a transducer (LA523) having variable-band linear array (13)(14). The neck of the study subjects was turned slightly to the left or right side for scanning the right or left carotid artery, respectively. The transducer was positioned at the anterior-lateral side of the neck. The lumen was maximized in the longitudinal plane with an optimal image of the near and the far vessel wall of the CCA or carotid bifurcation. Thus, typical double lines could be seen as the intima-media layer of the artery. Plaques were scanned in a longitudinal and cross-sectional plane showing the highest diameter. At least three heart cycles of every subject were acquired for measurement of the IMT or plaque, respectively. All sequences were stored digitally in the ultrasound device and transferred to a personal computer for further image analysis. This study complies with the Declaration of Helsinki and was approved by the ethics committee of the Technical University Munich, Germany. All patients signed an informed consent before being examined.

Problems Statement
The problems encountered in this study can be itemized as follows:  Speckle noises;  Calcification deposition on intima;  Moving target;  Parameter determination in different sonographic modalities.
The first problem is very common in the near arterial wall. Due to the jump of echo intensity between intimal surface and lumen, echoes (speckle noises) are transferred to the lumen which cause blurring of the intima-lumen interface. The second problem is that calcification deposition generally causes a shadow area under it. This results in a problem in detecting the far wall adventitia or the near wall intima depending on the location of the plaque. The third problem states that the artery changes in lumen diameter, IMT and position with systolic and diastolic phases. Since sequential images are processed, the tracking problem appears. The last problem is very common in many algorithms. Due to different properties of different ultrasonic machines, one has to determine the optimal parameter set for each specific machine.

Image Feature
A typical B-mode image is shown in Figure 1. There are plaques on both near and far walls. On the contrary, a typical CCA image without plaque having thin and thick IMT is shown in Figures 2 and 3, respectively. Outlining the boundaries of intima and adventitia is not easy, especially on the near wall.
Some features have been studied in previous studies [7,11], in which the most frequently used feature is the gradient of gray-values. From our previous studies [7,14,17], it is found that MacLeod operator combined with some specific enhancement filters are very suitable for feature extraction in intima and adventitia detections.   The dynamic B-mode sonographic sequences are recorded and converted to a series of static images and saved in BMP format. A graphic user interface (GUI) is designed for the user to select a region of interest (ROI) for the following automatic analysis. The selected ROI is a rectangular area which is large enough to cover the possible area during a complete systolic and diastolic cycle. Similarly, when the IMT is analyzed, the rectangular area covers only the intima, adventitia, and part of artery lumen area. However, since the IMT does not change too much during the heart cycle, the rectangular area for IMT measurement is much smaller than the one for lumen diameter measurement. After the ROI is selected, the following analysis and measurement processes are fully automatic.
The feature image (F) is obtained from the ROI sub-image by means of the following filtration: where I denotes the ROI sub-image; M f is the MacLeod operator; and x f is the enhancement filter. The design of x f for the IMT and LD detection is different. The MacLeod filter is defined as follows [18]: x where: The MacLeod filter is used to detect the gradient in different orientations. In total we used four angles, i.e., 4 3 and , 2 , . However, due to the noise interference in the sub-intima area, a CCA lumen single MacLeod filter cannot finish this work. Therefore, we add another enhancement filter to reduce the effect of noises. This is because the gradient direction is considered in the enhancement filter to degenerate the noises effect in the sub-intima computed by MacLeod filter. There are two kinds of f x , one for IMT detection and another one for LD detection. Each one has its own duality, one for near wall and another one for far wall. The enhancement filter for far wall LD detection has been published in [17]. For simplicity we do not repeat it again because it is a large matrix. The enhancement filter for far and near wall IMT detection is as follows: Figure 4 demonstrates the noise reduction effect of the enhancement filter. In cases without plaques, the thickness of sub-intima does not change on a large scale. If it has plaques, the noise layer existing in the sub-intima is not so strong and therefore, it can be ignored. The size of the filter is based on the consideration of the sensor's frequency and the anatomic structure. Since the normal IMT is around 0.5 mm for young persons and the sensor is 10 MHz, therefore, the pixel size is about 0.1 mm and the IMT is around 5 pixels. That is the reason why the filter size is 7 × 7. Via this design, the IMT can be detected with less interference from the noise layer.

(c)
On the task of lumen diameter detection, detecting the intima of both near and far walls precisely is very important. To achieve this goal, we separate the detection into two phases, i.e., the near wall IA (intima-adventitia) and far wall IA detection. From the two detections the intima layers of both walls can be obtained. The reason for using DDP (Dual Dynamic Programming) two-times instead of using the 4-dynamic programming is because it requires less computation time.
Most previous studies use the adventitia of the near wall instead of the initima to calculate LD. This is an approximation when there is no plaque. Since some patients have plaques, this approximation might cause significant faults in LD calculation.
According to our experience the enhancement filter in (2) combined with MacLeod filter produces good results for the far wall IMT detection. However, it is not the case for the near wall. Figure 5 demonstrates one typical problem. This image is a sub-image extracted from a raw image in a sequence. Some noises existed under the intima when the artery is in the systolic cycle because its larger moving speed. This will cause ambiguity in IA detection because most of the rest images have no such kind of noises. We therefore propose a method to alleviate this problem. The ambiguity results from the fact that the directional gradient of gray values cannot differentiate the interface between the lumen-intima interface and lumen-noise interface. If the gradient of the interface between the noise and artery lumen is larger than that between media and adventitia ( Figure 5 shows the case), then the following DDP will have a faulty result. The problem then becomes how to enhance the interface between media and adventitia. We observed that there are large gray values under both lumen-intima and media-adventitia interfaces in the near wall. We can make good use of this information. Nevertheless, we do not hope to enhance the outer wall simultaneously, since it will result CCA lumen CCA lumen in ambiguity in searching the dual lines, i.e., the intima and the adventitia. To solve this problem we apply dynamic programming to find out the outer wall first. Then an additional feature image (H) is calculated as follows: where I is the sub-image of the size M × N and r is determined by the knowledge and the image resolution; y outer is the coordinate of the outer wall. Notably, I is normalized so that 0 ≤ I(y,x) ≤ 1.

Dual Dynamic Programming (DDP)
In our previous study we have developed dual dynamic programming (DDP) for IMT detection.
Here we briefly describe this scheme as follows. Let H(x, y) denote the feature image grid of size M × N, where M and N are the number of rows and columns, respectively. Assume the DDP running from left to right in order to find dual curves such as intima and adventitia (or intimal layers of both near and far walls in case of detecting the lumen diameter). The DDP intends to find the global maximum, which is the summation of the feature values on the grid where the dual curves go through. The dual curves can be denoted as a point set , and its corresponding feature values are Finally, the cost function can be defined to be: matrix. Take k = 1 as an example as follows: where the row and column index pair is ) , ( 2 1 j j . Since the minimization problem depends on ) , ( 2 1 j j only, the process can be accelerated via finding the minimum in the following matrix: . The indices pointing backwards are stored in coordinate matrices C1 and C2 as follows: After the DDP is processed, the cost values are computed and saved in the cost map is searched, which indicates the last points of the dual lines. Having a backwards tracking from N to 1 in the coordinate matrices C1 and C2, the complete coordinates of the dual lines can be obtained.

Artery Movement Tracking (AMT)
The artery has extension and possible movement during the heart cycle. The movement is mostly in the y-direction, i.e., orthogonal to the artery flow. This is because the pressure from the upper skin is less than the pressure from the tissues below. Therefore, the ROI is changing in each image of a sequence. A simple automatic tracking scheme is embedded in order to track a correct ROI for the following DDP process. Let be the rectangular ROI selected by the user in the first image of a given sequence, where is the size of the ROI, and ) , ( T and ) (s T i , respectively; ∆s is a small integer covering a possible range of artery movement in the y-direction. The optimal value is the maximal NCC. Thus,

Procedure
We summarized the whole process in Section 3 to be a procedure as follows: Step 1: The user selects a directory in the hard disk where images of an image sequence are unpacked and stored.
Step 2: The image number (L) in the directory is counted. The first image is automatically loaded and shown on the monitor.
Step 3: The user selects a rectangular area where the artery should be analyzed including the near and the far wall of the CA.
Step 4: The selected rectangular area is equally divided into two parts: the upper and lower parts.
Record the coordinates. To analyze near wall, go to Step 5. To analyze far wall, go to Step 6. . Normalize I so that , apply DDP to find out the outer wall coordinate outer y .

Accuracy
An accuracy study is an assessment of the reliability of a system. In order to compare results made by the proposed automated system, we developed software for manual drawings. The manual drawings were performed by well-trained radiology technicians and then confirmed and modified by a medical physician. Four layers were drawn, i.e., intima and adventitia of both near and far walls. The linear interpolation is applied to fill out points between two points in manual drawings.
Assume a layer goes from left to right ( 1 x to N x ) and the i-th corresponding y-coordinate of the manual drawing and automated result are i y and i ŷ , respectively. The averaged error of one layer is

Results
The B-mode sonographic image videos were captured by an experienced physician of the Institute of Exercise and Health Sciences Vascular Lab at Basel University. Each video has more than 78 static images after decoding. The image resolution is 652 × 800 and the pixel size is about 0.106 mm. We took three test videos and show some of them in the paper.
The IMT detection has duality, i.e., near and far wall IMT detections. Traditionally, the far wall is taken for IMT measurement because the near wall intima is insufficiently visible due to the noise impact in many cases. Our method has better robustness against noise impact on the near wall and it is able to detect IMT on both near and far walls with and without plaques. Moreover, the measurement of LD is the distance between both intima of near and far walls.
We have tested three image sequences: S1, S2, and S3, they include plaques, thin, and thick IMT, respectively. All sequences have at least three heart cycles. Figures 6 to 8 show the test results of IMT detection on the far wall, whereas Figures 9 to 11 show the test results of IMT detection on the near wall of the three sequences.    In Figure 6 there is a plaque on the right-hand side which is successfully detected. A total of 78 images are in sequence S1. Another example shows the ability to detect a thin IMT as shown in Figure 7. A total of 111 images are in sequence S2. Figure 8 shows the thick IMT test result. A total of 86 images are in sequence S3. Figure 9 shows that our method has the ability to detect the plaque on the near wall. Figure 10(b) shows the problem using the method in [7] to detect the dual lines when the near wall is not clear. However, the proposed method can have a good result, as shown in Figure 10(c). As well, it produces good results in thicker IMT detection (Figure 11).    It can be noted that there are two -hills‖ during this period. It reveals that intima has a lower intensity than the adventitia, which results in a longer blurring (image 40-80) than in the adventitia (image 40-60). However, this error can be treated as low because its averaged error during these two -hills‖ is around 0.8 pixels. Compare this to the total average error (0.6 pixel for near wall and 0.55 pixel for far wall in Table 1) the difference (between 0.8 pixel and 0.55 pixel) is only 0.25 pixel (0.0265 mm). Absolutely it is in the range of 0.048 mm, thus within 5% of the IMT. We summarize all errors of the automated results comparing to the expert's manual tracing in Table 1. The average error on near wall detection (0.63 pixel or 0.067 mm) is larger than the far wall detection (0.55 pixel or 0.058 mm). This is because the near wall is in general vaguer than the far wall, especially in the systolic cycle. The intima detection is in general easier than the adventitia detection if there is a plaque. In thin or thick IMT cases the difference is not significant.

S2
In order to test our system, we analyzed more patients' image sequences (S4-S7) and calculated the unsigned (absolute error) and signed error. In this experiment we detect intima on both near and far walls. This is because we want to analyze the artery's lumen diameter change with respect to time during the heart cycle. The unsigned (Equation 8a) and signed (Equation 8b) error in j-th frame, similar to the equation in section 3.6, are calculated using the equations: where y Automatic (i) and y Manual (i) denote the i-th y-coordinate of automatic result and manual drawing. However, N is limited to 10 points in each image because the manual drawing is a tremendous human work. The averaged error in a sequence is calculated by: The results are shown in Table 2.  Table 2 it is noted that the signed error at the near wall is larger than that for the far wall. The far wall has almost no bias between the automated results and the manual drawings. This confirms that the detection on near walls is in general more difficult than the detection on far walls. Figure 14 depicts our manual drawing GUI. The solid lines denote the boundary of the automated analysis. The dashed lines denote where the expert should give their judgment on the lumen-intima interface. Manual drawing involves a tremendous amount of human work. Although we are not able to provide many manual drawings for comparison, however, the automated results were controlled visually by experts and all results are similar to those shown in Tables 1 and 2. The system is built on the Matlab platform. Some kernel functions are written in C language in order to reduce the computation time.

Discussion
According to our experience, the most difficult part is to detect the intima as well as the adventitia of the near wall. This is due to the fact that the near wall has a greater moving speed than that on the far wall during the systolic and diastolic cycle.
The dynamic analysis of the IMT and lumen of the artery faces problems which do not occur with static images. The static image is a single shot providing a clear view of the intima-lumen and adventitia-media interfaces. During the heart cycle the near wall usually moves more than 1 mm and might thus be associated with the blurring of the IM layers. This usually affects the intima-lumen interface more than the adventitia-medial interface due to the lower echogeneity of the intima. In CCA near walls with a thinner IMT this problem might also happen in the adventitia, as shown in Figure 10. However, the tracing error is relatively low compared with the overall movement of the artery. Notably, the accuracy of the detection on intimal layer does not show significant difference between the cases which have or have not plaque. Therefore, the existence of plaque does not affect the accuracy in calculation of the lumen diameter.
The artery adventitia is difficult to recognize if a plaque exists. For instance, the formation of a vulnerable plaque might be due to a severe microbial invasion. The microbial products and autoantibodies stick to the lipoproteins, which aggregate and obstruct the capillary lumen, leading to local ischemia and inflammation [19]. This may result in a blurred boundary between media and adventitia. This might be one of the reasons that the accuracy of automated identification on adventitia in the sequence S1 is worse than the other two sequences (Table 1). Precise plaque detection is a prerequisite for the calculation of total plaque area. For automatic detection, noises within a plaque do seriously disturb the tracing process.
The main contributions of this study are: (1) we propose a scheme which is able to detect the near and far wall intima and adventitia with and without plaque; (2) we solve the ambiguity problem caused by artery movement. These are the two major problems in the artery wall dynamic analysis in the B-mode RF sensors. Dynamic analysis of the artery lumen area might be a direct way to measure the local compliance and distensibility of the common carotid artery. Due to the heterogeneity of the arterial wall, the measurement at a single line-as it is in traditional M-mode imaging-may not be representative for the elastic properties of the arterial across a segment. The measurement of IMT as well as the LD in B-Mode image series might 1. Be more representative and 2. Offer the opportunity to show differences within a certain segment which stand for different atherosclerotic processes.
This may not be true for the plaque. However, to show changes of plaque thickness over time might be an indicator for vulnerability of the plaque, which is of significant clinical importance with respect to future cardio-vascular events.

Conclusions
A new method of analyzing IMT on both near and far wall of CCA in dynamic B-mode sonographic image sequences is proposed. This method is able to detect thin or thick IMT, with or without plaques. We have tested three image sequences and all results are compared to the manual tracing made by two experts. We have compared this method with the previous study and the results show that this method has better performance. The average error is of a sub-pixel level, i.e., 0.55 pixel (0.0583 mm) for the far wall and 0.63 pixel (0.0668 mm) for the near wall detection.