Automatic detection of the carotid artery boundary on cross-sectional MR image sequences using a circle model guided dynamic programming
© Cheng et al; licensee BioMed Central Ltd. 2011
Received: 8 December 2010
Accepted: 11 April 2011
Published: 11 April 2011
Systematic aerobe training has positive effects on the compliance of dedicated arterial walls. The adaptations of the arterial structure and function are associated with the blood flow-induced changes of the wall shear stress which induced vascular remodelling via nitric oxide delivered from the endothelial cell. In order to assess functional changes of the common carotid artery over time in these processes, a precise measurement technique is necessary. Before this study, a reliable, precise, and quick method to perform this work is not present.
We propose a fully automated algorithm to analyze the cross-sectional area of the carotid artery in MR image sequences. It contains two phases: (1) position detection of the carotid artery, (2) accurate boundary identification of the carotid artery. In the first phase, we use intensity, area size and shape as features to discriminate the carotid artery from other tissues and vessels. In the second phase, the directional gradient, Hough transform, and circle model guided dynamic programming are used to identify the boundary accurately.
We test the system stability using contrast degraded images (contrast resolutions range from 50% to 90%). The unsigned error ranges from 2.86% ± 2.24% to 3.03% ± 2.40%. The test of noise degraded images (SNRs range from 16 to 20 dB) shows the unsigned error ranging from 2.63% ± 2.06% to 3.12% ± 2.11%. The test of raw images has an unsigned error 2.56% ± 2.10% compared to the manual tracings.
We have proposed an automated system which is able to detect carotid artery cross sectional boundary in MRI sequences during heart cycles. The accuracy reaches 2.56% ± 2.10% compared to the manual tracings. The system is stable, reliable and results are reproducible.
Adaptations of the arterial structure and function were associated with blood flow-induced changes of wall shear stress which induced vascular remodelling via nitric oxide delivered from the endothelial cell . In order to assess functional changes of the common carotid artery (CCA) over time, a precise measurement technique was necessary. In , only two static MR images representing the end-diastole and systole were taken for the measurement of the lowest and the highest arterial diameter during the heart cycle. However, this has been shown to have a higher variability than the measurement along the complete heart cycle. Furthermore the measurement on continuous images by hand-hold tracing was extremely time-consuming. Only a limited number of publications focused on the carotid arterial structure and function using MRI with advanced imaging technologies in healthy subjects were found [3–5]. The purpose of this study was to establish a novel automatic common carotid arterial wall detection algorithm in MRI sequences over several heart cycles in order to precisely determine carotid diastolic and systolic diameter changes along time and the CCA local compliance. For this purpose we have collected some MRI data from participants of the multistage ultra-aahn"Trans Europe Foot Race" in 2009 (TEFR09).
Regarding related researches on engineering aspect, the similar work to ours was found in . The images they analyzed had plaques in the artery lumen. This was one additional problem they encountered more than ours. The other problems such as: contrast variations among blood, vessel wall and surrounding tissues, image artifacts caused by blood flow and random patient motion were similar to ours. Their system needed three user interactions: giving the system the artery's center point, a seed point inside the lipid core, and a circle surrounds the vessel. With the help of the prior knowledge combined with an elliptic model and fuzzy clustering, their system was able to identify the arterial boundary and plaque boundary.
Another previous study  was a measurement on arterial wall using discrete dynamic contour (DDC) . Their method was semi-automatic because the system needed an initial contour of the inner wall contour. Moreover, their images were black blood vessel so that they were able to detect both the inner and outer wall boundaries of the carotid artery.
Another related article but focused on the coronary artery boundary detection was found in . They proposed a surface-based method to detect the coronary artery boundary in the poor quality X-ray angiography based on a 3D generalized cylinder model. Since their application was on the X-ray angiogram, therefore, the proposed method was not suitable for the application on MRI sequences.
Our contributions are to develop an automatic method to measure the arterial boundary in MR images. It is able to detect the carotid artery center position in the first stage. In the second stage, the cross sectional arterial wall boundary can be detected via Hough transform and a circle model dynamic programming. The circle model dynamic programming lets the system control the output boundary to be somewhat round but having the ability to detect the fine structure.
The paper is organized as follows. The image sources and MRI protocol are introduced in Section 2. Section 2.1-2.2 describes how to detect the artery lumen center position. In Section 2.3, the circle model guided dynamic programming is issued in details to solve our problems. Afterwards, results are given in Section 3. We then discuss properties of the proposed scheme in Section 4. Finally the conclusion is given in Section 5.
The MR-measurement of the maximal systolic and distal vessel-diameter of the CCA with additional blood pressure information leaded to local compliance of the arterial wall. After the approval of the local ethics committee in accordance to the Declaration of Helsinki, 12 participants in the TEFR09-project of the German Research Foundation (DFG-Project GZ: SCHU 2514/1-1, AOBJ: 565344) have been collected for vascular studies based on MR image sequences. One MRI sequence of one subject was randomly chosen from these 12 subjects for the validation of the novel detection algorithm of the CCA lumen presented in this study.
To acquire the vascular MRI sequences, a mobile 1.5-T MR imager (Siemens - Avanto™, Model Mob. MRI 02.05, Siemens Ltd., Erlangen, Germany) and a custom-designed four-channel phased dual mode neck matrix coil with 4 integrated preamplifiers (Siemens Ltd.) were used. The movement artifact was minimized via using a dedicated head restraint system (head coil, Siemens Ltd.) to fix the head and neck in a comfortable position (patient position: supine, head to feet).
After an initial coronal localizer, three fast localizers (triplanar TRUFI: "true fast imaging with steady state precision"; Siemens Ltd.) were used to identify the axial perpendicular acquisition location at the right CCA 10 mm inferior the carotid bifurcation. Contrast media could not be used in this study because the athletes did not accept it.
To increase the spatial resolution of the measurement (cross section view of CCA) a T2*-weighted gradient-spoiled gradient-echo cine-sequence (FLASH: "fast low angle shot", Siemens Ltd.) with prospective two dimensional ECG gating (cardiac triggering) was used. Parameters were set to be: flip angle 15°, echo time 5.41 ms, repetition time 34.74 ms, slice thickness 6 mm, field of view 289 cm2, matrix size 320×320, pixel size 0.53125 mm ISO, pixel bandwidth 250, number of images per sequence: 50 images for one RR-cycle. The imaging acquisition time was approximately 5 minutes for each sequence.
2.1 Carotid artery position detection
The carotid artery position detection is the first procedure because it can reduce the following computation time. This procedure detects only the rough artery's center position but not the artery boundary. Here we propose an easily implemented but fast algorithm to perform this work.
where T2 is the threshold obtained by Otsu's method which is a value between the gray value of the vessel lumen and other tissues. Notably the computation of T2 is based on the precondition of ignoring the background marked as -1 in the first stage. After this process vessel lumens are segmented out.
where 's' is the structure element with a disk shape (radius is 1). More clearly, it has a central pixel (the reference point) and the four-neighborhood pixels.
Afterwards, the rest features (area size and shape) are utilized to identify the carotid artery from other vessels. This process is divided into left and right parts. Assume we are processing one part of them, two largest areas are chosen and they are the internal jugular vein and the carotid artery. Their boundary points are obtained by using a simple Sobel operation . Afterwards, the PCA (principal component analysis ) is applied to get the major axis and minor axis. The length ratio of these two axes is a feature to indicate if the shape is round. Via this scheme, we are able to identify the carotid artery from other vessels. Once its center position is estimated, a region of interest (ROI) denoted as R s is extracted from the source image R to the following procedures while artery center centered at R s. The size of region (R s) depends on the pixel size of MR images and the anatomic knowledge how large the carotid artery can be. This can be calculated in prior.
2.2 Carotid artery boundary detection
2.2.1 Directional gradient
The results obtained by the method addressed in Section 2.1 do not offer accurate artery boundaries. This is because the Otsu's thresholding technique never offers good results in case that the intensity is not consistent for each object (here the carotid artery) to be measured. Especially it is possible that the morphological opening operation shrinks the artery's actual size. Therefore, we propose a method to detect the accurate artery boundary. Since the artery boundary has intensity different to its surrounding area, the gray level gradient is useful information. However, the internal jugular vein is very close to the artery in images so that it makes the boundary detection difficult if we consider only the intensity gradient only. This case is worse if the intensity gradient on the vein boundary is stronger than that on the artery boundary. We therefore take the direction into consideration and name this gradient to be directional gradient. In literatures we do not find any similar publication except the directional gradient vector flow in . Our consideration is that: since the artery lumen is brighter than its surrounding areas, gradients resulted from bright pixels to dark pixels are of interest. If the artery center is known, then the directional gradient can be found that is parallel to the radiation lines centered from the artery. Actually the directional gradient is a special case of multi-directional gradient.
Where t = 4(θ + π)/π. The calculations when is located in other regions are similar. Via this scheme an edge map representing the gradient intensity can be obtained and denoted as R e. Notably, we are looking for gradients resulted from bright pixels to dark pixels. Therefore, the edge map we look for is a minimal value. The negative gradient denotes the boundary gray level changing from bright pixel to dark pixel along vector which is we want. Via this way, the positive gradient resulted from the jugular vein very close to carotid artery will not affect the searching of the artery boundary. Afterwards, positive values in R e are set to zeros and the Otsu's thresholding technique is applied again to find a threshold value. The binarisation technique is applied to R e: values under the threshold are set to ones. The resultant binary image is denoted as R b.
2.3 Circle model and dynamic programming
The round shape information of artery is important. It is used to avoid possible errors caused by local noises. These errors include the heterogeneous gradient obtained in the artery lumen and at the boundary. To alleviate this problem we apply Hough transform  to detect round objects in R b. The Hough transform is a feature extraction technique used in image analysis, computer vision, and digital image processing. There are three variables to be determined, i.e., the circle center position (x and y coordinates) and the radius. Although the radius of artery is unknown in advance, its range of variance can be estimated. Therefore, we have to calculate its Hough transform by varying the radius from r 1 to r n ∈ N. After transformations each radius obtains an accumulator matrix. We search the maximum value in each accumulator matrix and find out the one which has the maximum value among all accumulator matrices. The corresponding radius and position is the artery's radius and center position, respectively.
Although the artery lumen is in general round, however, it is not the case from the pixel's view point. In addition, not all artery lumen is pure round during the heart beat cycle. Some of them are elliptic. A fine tuning is then necessary to obtain a sub-pixel accuracy. For this reason, we propose a method to identify the artery boundary based on a circle model.
Dynamic programming is a method of solving complex problems by breaking them down into simpler steps used in mathematics and computer science . It is applicable in image processing to solve optimal problems such as finding a minimum (or maximum) with some given constraints[16, 17]. However, the limitation of using this technique in images is that it cannot solve the closed form contour. One solution is to transform the image from Cartesian coordinate to polar coordinate  and then apply dynamic programming on the polar coordinate. Note that this procedure applies only on a region of interest (ROI). However, two preconditions have to be satisfied: 1) the rough center position is known; 2) the contour has a star-like shape. Our problem meets these two preconditions.
Curve continuity: A variable for continuity is considered. Let d r denote the maximal range that nodes in column x - 1 are allowed to jump onto the next column x in either up or down directions. Therefore, each node has maximum (2d r + 1) possible link paths to its previous column. If d r is set larger, the curve's roughness and the computation time are both increased. The smoothness of the curve is quantified by j = y x-1- y x which is embedded into the cost function.
Circle model: The circle model having a known radius is embedded into the structure to guide the dynamic programming. This is based on the fact that the artery boundary is near round and the radius is estimated by Hough transform in prior. A Gaussian model is used to generate the strength how strong the dynamic programming is guided by the circle model. Let r denote the known circle radius, the strength is formulated as , where σ is a variable controlling the strength of guide. If σ is getting smaller, the Gaussian has a thin but sharp shape and the circle model has a larger effect on the result, i.e. it is a more circle-like boundary. If σ is getting larger, the Gaussian term vanishes and it works like a normal dynamic programming without the circle model. Since y and r are both integers, a look-up table can be set to reduce the computation time.
Directional gradient: The directional gradients are very important information to detect the artery boundary accurately. Gradients having negative values denote the carotid artery boundary, while positive gradients denote other boundary which we treat them as noises. Thus, the boundary detection problem is then transformed to an optimization problem which searches an optimal contour:
subject to 2 ≤ x ≤ N, 1 ≤ y ≤ M;
Therefore, the index can be stored in the coordinate matrix X(x, y) = y + j*. In this construction, small cost values indicate higher likely boundary information. The position with the minimum cost value in the cost map C(x, y) is searched. With a backward search from N to 1 in X, the complete coordinates (p 1 p 2 p 3...p N ) of the artery boundary can be determined, which is the optimal solution to this problem. Notably, the processes addressed in Section 2.3 are applied only on the extracted ROI (R s, R e, and R p) to reduce the computation time tremendously. Results obtained by dynamic programming are integers. A moving average technique  is applied to make the boundary smoother. In order to reduce the computation time, we simply average the neighbouring 4 points and the center point. The resultant polar coordinates are then transformed back to Cartesian coordinates.
2.4 System reliability
The proposed system has been studied for the reliability against different noise levels and different contrasts.
where AAutomated (i) and AManual (i) are areas calculated by the automated and the manual drawing on image number i, respectively.
where M × N denotes the image dimension, g (x, y) and n(x, y) are the intensities of image and noise at image coordinate (x, y), respectively. The MR image at (x, y) has an intensity g(x, y)≥0 and the noise intensity might be negative. To calculate the SNR we define if g (x, y) = 0 or n (x, y) = 0.
The detected artery lumen center position is used to predict the center in the next image. Similarly, the detected artery lumen radius is used to set the size of ROI in the next image. Here we expend 1.5 times radius from the center to each side (left, right, up, and down) to define the size of ROI. For the reason of explanation, the ROIs shown in Figure 4 are larger than 1.5 times.
where εi has been defined in equation (7), N = 50.
The computer system has Intel® Core™ 2CPU T5600, 1.83 GHz, with 2 GB RAM. The programs are based on the Matlab platform . The computation time for 50 images processing is around 30 seconds. More results are downloadable under our website.
In this study we use a circle model in Hough transform and the dynamic programming instead of the ellipse Hough transform because of the consideration on the computation time. Full ellipse detection is rarely performed because it is very slow. It requires a 5 dimensional parameter space - as opposed to 2 for straight line detection and 3 for circle detection. Although the gradient direction can be taken into consideration to reduce the computation time , it still needs much more time than that in our design. Moreover, the artery's shape can be changed if a plaque exists. Our design has the advantage to detect boundary which is not a circle or an oval.
Since our algorithm uses area, shape, and intensity as features to identify the carotid artery position, the conditions (prerequisites) that the carotid artery is able to be identified are: 1) there are less or no plaques in the artery; 2) the blood flow maintains in a level so that the intensity in artery lumen in MRI is distinguishable from neighbouring tissues. Our subjects are healthy runners, although some are old people, there are less plaques in the artery lumen. Therefore all carotid artery lumens can be modelled by a circle model. The blood flow in the carotid artery is different from that in the femoral artery. It does maintain at a level so that the intensity in the lumen is distinguishable from the artery wall and other tissue nearby. Therefore, there is no problem to identify the carotid artery centers using our proposed algorithm.
Regarding the chosen of parameters, there are three parameters in our scheme: α, d r , and σ . Normally, the discontinuity is prevented by setting d r = 1. We suggest the range of α to be (0,1]. If α becomes smaller, the output curve is more rough. On the contrary, the curve becomes smooth if α is set near to 1. The standard deviation σ is a control by the circle model. The range we suggest is [0.4, 4]. If σ ≤ 0.5, then the output curve is nearly round. If σ ≥ 1 then it is able to detect the fine structure such as plaques in the artery lumen.
Based on the study on different contrast resolutions (from 50% to 100%) and different noise levels (SNR ranges from 16 dB to 20 dB), the proposed method has shown its robustness and reliability against contrast resolution and noise.
In conclusion, we have proposed a fast and robust scheme to detect the carotid artery boundary in MR image sequences fully automatically. This scheme is divided into two stages: (1) detect the center of the carotid artery (2) detect the boundary of the carotid artery. We combine the circle model with the dynamic programming so that the resultant boundary is circle-like shape. The accuracy control shows that the averaged relative error of the automated results compared to the manual results is 2.56% and the standard deviation is 2.10%. Via this system we are able to analyze tremendous amount of images and all results are repeatable.
This work was supported in part by the German Research Foundation (DFG: "Deutsche Forschungsgemeinschaft"), under Grants SCHU 2514/1-1 and SCHU 2514/1-2 and by the National Science Council (NSC), Taiwan, under Grant NSC 99-2221-E-039-012.
- Kingwell BA, Sherrard B, Jennings GL, Dart AM: Four weeks of cycle training increases basal production of nitric oxide from the forearm. Am J Physiol Heart Circ Physiol 1997, 272: H1070-H1077.Google Scholar
- Mohiaddin RH, Underwood SR, Bogren HG, Firmin DN, Klipstein RH, Rees RS, Longmore DB: Regional aortic compliance studied by magnetic resonance imaging: the effects of age, training, and coronary artery disease. Br Heart J 1989, 62(2):90–96. 10.1136/hrt.62.2.90View ArticleGoogle Scholar
- Chow TY, Cheung JS, Wu Y, Guo H, Chan KC, Hui ES, Wu EX: Measurement of common carotid artery lumen dynamics during the cardiac cycle using magnetic resonance TrueFISP cine imaging. J Magn Reson Imaging 2008, 28(6):1527–1532. 10.1002/jmri.21527View ArticleGoogle Scholar
- Lin AP, Bennett E, Wisk LE, Gharib M, Fraser SE, Wen H: Circumferential strain in the wall of the common carotid artery: comparing displacement-encoded and cine MRI in volunteers. Magn Reson Med 2008, 60(1):8–13. 10.1002/mrm.21621View ArticleGoogle Scholar
- Keenan NG, Locca D, Varghese A, Roughton M, Gatehouse PD, Hooper J, Firmin DN, Pennell DJ: Magnetic resonance of carotid artery ageing in healthy subjects. Atherosclerosis 2009, 205(1):168–173. 10.1016/j.atherosclerosis.2008.11.018View ArticleGoogle Scholar
- Adame IM, van der Geest RJ, Wasserman BA, Mohamed MA, Reiber JHC, Lelieveldt BPF: Automatic segmentation and plaque characterization in atherosclerotic carotid artery MR images. MAGMA 2004, 16: 227–234. 10.1007/s10334-003-0030-8View ArticleGoogle Scholar
- Ladak HM, Thomas JB, Mitchell JR, Rutt BK, Steinman DA: A semi-automatic technique for measurement of arterial wall from black blood MRI. Medical Physics 2001, 28(6):1098–1107. 10.1118/1.1368125View ArticleGoogle Scholar
- Lobrget S, Viergever MA: A discrete dynamic contour model. IEEE Trans. on Medical Imaging 1995, 14: 12–24. 10.1109/42.370398View ArticleGoogle Scholar
- Kayikcioglu T, Gangal A, Turhal M, Kose C: A surface-based method for detection of coronary vessel boundaries in poor quality X-ray angiogram images. Pattern Recognition Letters 2002, 23(7):783–802. 10.1016/S0167-8655(01)00156-8View ArticleGoogle Scholar
- Otsu N: A threshold selection method from gray-level histograms. IEEE Trans Sys Man Cyber 1979, 9: 62–66. 10.1109/TSMC.1979.4310076View ArticleGoogle Scholar
- Gonzalez RC, Woods RE: Digital Image Processing. Prentice Hall; 2002.Google Scholar
- Jolliffe IT: Principal component analysis. Springer-Verlag; 1986.View ArticleGoogle Scholar
- Cheng J, Foo SW, Krishnan SM: Directional gradient vector flow for snakes. Proceedings of the Fourth IEEE International Symposium on Signal Processing and Information Technology 2004.Google Scholar
- Shapiro L, Stockman G: Computer vision. Prentice-Hall; 2001.Google Scholar
- Dasgupta S, Papadimitriou CH, Vazirani UV: Algorithms. McGraw-Hill; 2008.Google Scholar
- Cheng DC, Schmidt-Trucksaess A, Liu CH, Liu SH: Automated Detection of the Arterial Inner Walls of the Common Carotid Artery Based on Dynamic B-Mode Signals. Sensors 2010, 10(12):10601–10619. 10.3390/s101210601View ArticleGoogle Scholar
- Cheng DC, Jiang X: Detections of Arterial Wall in Sonographic Artery Images Using Dual Dynamic Programming. IEEE Trans. On Information Technology in Biomedicine 2008, 12(6):792–799. 10.1109/TITB.2008.926413View ArticleGoogle Scholar
- Domínguez AR, Nandi AK: Improved dynamic programming based algorithms for segmentation of masses in mammograms. Medical Physics 2007, 34(11):4256–4269.View ArticleGoogle Scholar
- Arce GR: Nonlinear Signal Processing: A Statistical Approach. New Jersey: Wiley; 2005.Google Scholar
- Matlab The MathWorks 2008. [http://www.mathworks.com/]
- Ballard DH: Generalizing the Hough transform to detect arbitrary shapes. Pattern Recognition 1980, 13(2):111–122. 10.1016/0031-3203(81)90009-1View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.