Automated Detection Framework of the Calcified Plaque with Acoustic Shadowing in IVUS Images

Intravascular Ultrasound (IVUS) is one ultrasonic imaging technology to acquire vascular cross-sectional images for the visualization of the inner vessel structure. This technique has been widely used for the diagnosis and treatment of coronary artery diseases. The detection of the calcified plaque with acoustic shadowing in IVUS images plays a vital role in the quantitative analysis of atheromatous plaques. The conventional method of the calcium detection is manual drawing by the doctors. However, it is very time-consuming, and with high inter-observer and intra-observer variability between different doctors. Therefore, the computer-aided detection of the calcified plaque is highly desired. In this paper, an automated method is proposed to detect the calcified plaque with acoustic shadowing in IVUS images by the Rayleigh mixture model, the Markov random field, the graph searching method and the prior knowledge about the calcified plaque. The performance of our method was evaluated over 996 in-vivo IVUS images acquired from eight patients, and the detected calcified plaques are compared with manually detected calcified plaques by one cardiology doctor. The experimental results are quantitatively analyzed separately by three evaluation methods, the test of the sensitivity and specificity, the linear regression and the Bland-Altman analysis. The first method is used to evaluate the ability to distinguish between IVUS images with and without the calcified plaque, and the latter two methods can respectively measure the correlation and the agreement between our results and manual drawing results for locating the calcified plaque in the IVUS image. High sensitivity (94.68%) and specificity (95.82%), good correlation and agreement (>96.82% results fall within the 95% confidence interval in the Student t-test) demonstrate the effectiveness of the proposed method in the detection of the calcified plaque with acoustic shadowing in IVUS images.

lumen, the vessel wall, and the plaque components manually. However, the manual analysis is a very time-consuming process and can easily lead to the inter-and intra-observer variability.
Therefore, it is highly desired to develop the computer-aided algorithm for analyzing IVUS images in order to increase diagnostic efficiency since well-designed computer-aided analysis could generate more objective conclusions. One important application of the computer-aided analysis in IVUS images is the determination of the stenosis degree and the area of blood flows inside the lumen by detecting the lumen border and the media-adventitia border. A large number of computer-aided algorithms have been developed in this area. Most of these algorithms were implemented based on the image descriptor and the contour evolution: (1) among the image descriptor based algorithms, the gradient is one widely used image descriptor [5]. However, it is susceptible to speckle noise and varieties of artifacts in IVUS, such as stents and the guide wire. Therefore, some other image descriptors, such as edge pattern [6], textural information [7,8] and statistical properties of the pixel intensity [9,10], have been proposed to overcome the disadvantages of the gradient descriptor; (2) in contour evolution based algorithms, the vessel borders are usually detected by evolving the per-defined initial contour until convergence. These algorithms could be divided into snake-based [4,11] and level-set-based algorithms [12,13]. The snake-based algorithms always treat the evolving contour as a parametric curve, and compute the desired border by minimizing the curve energy. Unlike the snake-based algorithms, the level-setbased algorithms often apply implicit functions to format the contour and evolve the contour of the zero level by solving the partial differential equations [14]. Another important application of the computer-aided analysis of IVUS images is the detection of the calcified plaque with acoustic shadowing. In the clinical expert consensus documents from the American College of Cardiology and the American Heart Association [15,16], the calcified plaque can be used to detect obstructive coronary artery disease, and moreover it can also be considered as an accurate indicator of atherosclerotic disease and the risk of acute myocardial infarction [17], especially in diabetics [18], smokers [19] and elderly [20]. The methods for detecting the calcified plaque with acoustic shadowing could be divided into two categories: (1) the first kind of approaches will look for the calcified plaque inside the region between the lumen and the vessel wall, after the lumen border and the media-adventitia border are both detected. Because the identification of these two borders could narrow down the region of interest (ROI) for facilitating the detection of the calcified plaque, some classification algorithms have been used to extract the calcified plaque from the ROI, such as Bayesian classifier [13] or fuzzy k-means [21]. However, the performance of these classification algorithms relies on the accurate location of the lumen and media-adventitia borders, and the calcified plaque with acoustic shadowing can likely damage the accuracy because it can be considered as an artifact to obstruct the detection of the two vessel borders. For the media-adventitia border, it may be shaded by the acoustic shadowing because the most of the ultrasonic wave can be reflected by the calcified plaque, and for the lumen border, it may be represented as a weaker border than the border between the calcified plaque and the soft plaque; (2) the other kind of approaches is to detect the calcified plaque directly, which always adopted a coarse-to-fine strategy to locate the calcification. In the first step, the classification algorithms, such as Otsu method [22], thresholding method based on 2-D Renyi's entropy [23] or kmeans [24], were applied to the IVUS image s in order to coarsely specify the locations of calcified plaques. And in the second step, the acoustic shadowing would be the recognized by comparing the gray intensity of the calcified plaque with that of the acoustic shadowing [25,26]. However, in the coarse detection, the ROI detection of the calcified plaque did not depend on its acoustic shadowing, and thus some bright non-calcified regions might be considered as the ROI of the calcification, such as collagen and guide-wire artifacts. Furthermore, the refinement relied on the image contrast between the calcified plaque and the acoustic shadowing, which is controlled by the adjustments of gamma curve in image acquisition [27], and thus the contrast may vary among different IVUS images. So the refinement of the calcified plaque should not rely on the gray intensity since it would reduce the performance of the algorithm.
In this paper, we develop a coarse-to-fine method to automatically detect the calcified plaque with acoustic shadowing in IVUS images in three steps: (1) the coarse detection of calcified plaque is implemented by Rayleigh mixture model (RMM) and Markov random field (MRF); (2) the refinement of calcified plaque is carried out by five predefined constraints; (3) the border of calcified plaque is detected by the graph searching algorithm. The validation of our method is performed over in-vivo IVUS images (DICOM format) acquired from eight patients. The results of our method are compared with the manual drawing results from one experienced cardiologist. In the experiments, the ability of our method to distinguish between the IVUS images with and without the calcified plaque is firstly evaluated by the sensitivity and specificity. Next, in the IVUS images with calcified plaque, two different evaluation methods, linear regression and Bland-Altman analysis, are performed in order to validate the correlation and the agreement between our method and the manual detection. In addition, we compute how much overlapping between two calcified plaques separately produced by our method and the manual detection, in order to evaluate the location accuracy of the calcified plaque detected by our method.
The rest of this paper is organized as follows: The second section discusses the proposed method, including the pixel classification by the RMM, the detection of the angular location of the calcified plaque, and the refinement and the border detection of the calcified plaque. The third section shows the experiments and results. The fourth section demonstrates the discussion, and the last section draws the conclusions of this study.

Methodology
In this section, we develop a method to detect the calcified plaque with acoustic shadowing in IVUS images. The flowchart of the proposed method is shown in Figure 1. First, a RMM is applied to cluster pixels in IVUS images in order to distinguish between the hyperechoic and hypoechoic regions, and solved by the EM algorithm. In the solution of the RMM, a new prior probability is produced based on the spatial relationship among neighboring pixels. Second, a MRF is employed to coarsely detect the location of the calcified plaque, and solved by the BP algorithm. In the solution of the MRF, the results of the RMM and a curve called maximum intensity curve (MIC) defined in the IVUS image are combined to compute the relationship between the observed variables and hidden variables, which can be considered as the prior information of the MRF. Third, the pseudo calcified plaques are removed by five predefined constraints. At last, the border of the calcified plaque is detected by the graph searching method algorithm [28] with cost function formatted by MIC and the image gradient.

Pixel Classification in IVUS Images
Different vascular tissues appear in different gray levels and shapes in IVUS images, so we can use the RMM to describe the distribution of the image intensity and cluster pixels in IVUS images [29]. The RMM formulas are shown in Appendix A. There are three important issues in the RMM: (1) the selection of the number of the mixture components; (2) the construction of the prior probability; (3) the estimation of the parameters.
First, the number of the mixture components K is empirically predefined as 5 in our study, different from K~3 in previous approaches [30] [31], because the purpose of the pixel classification in the proposed method is to label pixels belonging to the acoustic shadowing into the same group, which can help the following steps for detecting the calcified plaque, rather than distinguish among three layers of the blood vessel [30], or among different types of atherosclerotic plaques [31]. Here we give a comparison of the RMM classification between K~3 and K~5, shown in Figure 2(b) and 2(c). For K~5, 97.74% of pixels in the acoustic shadowing are labeled into the same group, while for K~3, the percentage is only 82.35%. Therefore, the RMM with K~5 can more appropriately extract the acoustic shadowing in the IVUS image.
Second, different from many previous methods [30] [31], we add the neighboring relationship in the prior probability p ij , rather than assume the independence between pixels. The computation of p ij is shown in Appendix A.
Third, the parameters of the RMM are estimated by the maximum likelihood method (ML) [24], and the most direct way to solve the ML method is to compute the derivative of the likelihood function with respect to every unknown parameter. However, the computation is intractable due to the complexity of the class-conditional probability density P(Y i Dv j ) in Equation (8) and the prior p ij in Equation (13). Therefore, we apply the expectation-maximization algorithm (EM) [32] to solve the ML method. The EM algorithm is an iterative method for estimating the parameters of the maximum likelihood method. In each iteration, the EM algorithm can be divided into two step: E-step and M-step. The E-step is used to determine the objective function in this iteration, and the M-step is used to maximize the objective function in order to update the parameters. The details of the Estep and the M-step are shown in Appendix B.
From above, we present the detailed steps of the method for estimating RMM parameters in Appendix C. After the parameters of the RMM are determined, the pixels in the IVUS image can be clustered by the maximum a posterior estimation method (MAP) [24]: the ith pixel can be assigned to the jth class if where The results of the pixel classification by the RMM are shown in Figure 2(c).

Detection of Angular Location of the Calcified Plaque
Because the growth of the calcified plaque is spatially continuous along the vessel wall, the pixels belonging to the same calcified plaque have the relationship along the angular direction in polar coordinates of the IVUS image. In order to model the relationship, we firstly transform the IVUS image and its corresponding image classified by the RMM both from the Cartesian domain to the polar domain. The IVUS polar image is denoted by I and its corresponding image after pixel classification is denoted by I R , shown in Figure 2. Then every column of I can be described by a random variable, and these random variables generate a MRF, which contains this kind of relationship [33]. Let the numbers of rows, columns and gray-scale value at the coordinate (x,y) in I be H, W and I(x,y), respectively, where x is the column coordinate and y is the row coordinate. Then the hidden variables in the MRF are denoted by Z 1 ,Z 2 ,:::,Z W , where Z i (i~1,2,:::,W ) corresponds to the ith column of I. The domain of Z 1 ,Z 2 ,:::,Z W is -1 and +1, where Z i~z 1 means the ith column of I contains the calcified plaque, and Z i~{ 1 means it does not contain the calcified plaque. The neighborhood of Z i includes Z i{1 and Z iz1 . In addition, Y 1 ,Y 2 ,:::,Y W denote the observed variables corresponding to Z 1 ,Z 2 ,:::,Z W , respectively. Through the MRF, the possibility of every column I containing the calcified plaque with acoustic shadowing can be calculated. Therefore, the columns in I can be divided into two parts: one part contains the calcified plaque with acoustic shadowing and the other does not, where the columns including the calcified plaque can be considered as the angular location of the calcified plaque. There are two important issues in the MRF: the determination of the prior and the solution of the MRF.
1. The prior represents the relationship between the observed variables Y 1 ,Y 2 ,:::,Y W and the hidden variables Z 1 ,Z 2 ,:::,Z W . In our method, the prior F can be computed by the gray-scale information, the pixel classification by the RMM and the maximum intensity curve in I, shown in Appendix D. The maximum intensity curve is defined in Definition 1.
Definition 1 (Maximum intensity curve). For V i[f1,2,:::, W g, the pixel with the maximum gray-scale value in the ith column of I can be found, and its row coordinate value is denoted by p i . Let L m~f p 1 ,p 2 ,:::,p W g, and then L m is called the maximum intensity curve, abbreviated as MIC.
2. The exact solution of the MRF is usually performed by the maximum likelihood (ML) or MAP parameter estimation [34]. In these methods, the parameter estimation can make the free energy of the MRF reach the extremum, and moreover it can be considered as a training process, because it requires the known observation (training data) to optimize the parameters. However, in many cases, there is no closed form solution of the MRF for the ML or MAP parameter estimation [34]. Thus, it is needed to use the approximate method to solve the MRF, and the training process aforementioned is implicit in the approximation process. The Bethe free energy is a common approximation to the free energy of the MRF [35], and has been proved to be equivalent to the belief propagation (BP) [36]. Therefore, we can use the BP algorithm to approximately solve the proposed MRF. The BP algorithm is an iterative ''sum-product message passing'' algorithm and applies a kind of variables called ''belief '' to approximate the probability of the random variables [37] [38]. The prior of the BP algorithm equals F in Equation (28). The process of the proposed BP algorithm is shown in Appendix E.
The final beliefs can be denoted by b Ã (Z 1 )~b (tz1) (Z 1 ), b Ã (Z 2 )~b (tz1) (Z 2 ),:::,b Ã (Z W )~b (tz1) (Z W ) if the BP algorithm is converged in the tz1 iteration. Then the belief of every column with the calcified plaque b Ã (Z i~z 1) is binarized in order to find which columns of I contain the calcified plaque, formulated as At last, we can find the angular location of the calcified plaque in I as follow: for two column d l and d r ( are all equal to 1, and b Ã (Z d l {1~z 1), b Ã (Z d r z1~z 1) are both equal to 0, then the angular range between the d l th column and the d r th column of I can be considered as the angular location of one calcified plaque, where d l and d r are called the left border and the right border of the calcified plaque, respectively.

Refinement of the Calcified Plaque
The detection of the angular location of calcified plaques is affected by the sensitivity of MIC to the noise and the classification error from the RMM, which leads to the pseudo calcified plaques. Therefore, the refinement of the calcified plaque is needed. Here we define five constraints in order to refine the calcified plaque, that is, if a calcified plaque does not satisfy these constraints, then it can be considered as a pseudo calcified plaque. The details of the five constraints are shown in Appendix F. These constraints describe different aspects of the calcified plaque. Firstly, the Constraint 1 considers the brightness of the calcified plaque in IVUS images because the minimum threshold T 1 of pixel values limits the overall brightness of the calcified plaque. Then the Constraint 2 and the Constraint 3 aim to reduce the influence from other tissues on the recognition of the calcified plaque, because some tissues on adventitia, such as collagen, have highlevel ultrasonic echo and thus can be represented as bright regions with high gray intensity, which may be mistaken as the calcified plaque. The Constraint 2, concerning the growth direction of the calcified plaque, limits the slope of two endpoints of the calcified plaques in order to recognize the bright region with too large slope in polar coordinate as the non-calcified region, because the calcified plaque grows along the vessel wall in histology. However, only relying on the Constraint 2, some non-calcified tissues on adventitia along the vessel wall still cannot be recognized. Therefore, the length of the calcified plaque must be considered. The limitation of the angular length of the calcified plaque along the polar axis in the Constraint 3 can exclude the bright regions with too long angular length. At last, the Constraint 4 and the Constraint 5 focus on the acoustic shadowing behind the calcified plaque. Because the calcified plaque can reflect most ultrasonic signals, the acoustic shadowing behind the calcified plaque is darker than its nearby tissues. Therefore, in the Constraint 4, the limitation of the minimum difference between the brightness of the acoustic shadowing and the nearby tissues can facilitate to recognize the shadowed regions on adventitia, ensuring that every detected calcified plaque has acoustic shadowing behind it. The Constraint 5 aims to recognize the pseudo acoustic shadowing through the shape of the acoustic shadowing, where Equation (47) represents the closeness between the calcified plaque and its acoustic shadowing because the amplitude of the ultrasonic signal behind the calcified plaque decreases sharply. Moreover, the ultrasound wave is omnidirectionally emitted from the catheter, and thus the acoustic shadowing is rectangle-like in polar coordinates. Equation (48) is used to describe the similarity between the acoustic shadowing and the rectangle-like shape in polar coordinates.

Border Detection of the Calcified Plaque
After the refinement, we detect the border of the calcified plaque, which can be divided into four part: the left border, the right border, the upper border and the lower border. The left border d l and the right border d r represent the angular location of the calcified plaque, and have been computed in the above. Therefore, only the upper border and the lower border of the calcified plaque should be detected in this part. In order to trace the two borders, we employ the graph searching algorithm [28], which is a searching algorithm to find the optimal path in the weighted graph (cost function) connecting two endpoints. Because MIC crosses through the calcified plaque, we consider the upper border is inside a region, denoted by R 1 , above MIC and between the d l th column and the d r th column, and consider the lower border is inside a region, denoted by R 2 , below MIC and between the d l th column and the d r th column. R 1 and R 2 are shown in Figure 3(a). Then the graph searching algorithm is separately applied to R 1 and R 2 to extract the upper border and the lower border of the calcified plaque. There are two important issues in the graph searching algorithm: the selection of the two endpoints and the formulation of the cost function.
(1) In order to generate a traced border, two endpoints of the border, called the beginning point and the ending point, should be acquired in advance for the graph searching algorithm. With respect to R 1 or R 2 , the column coordinates of the two endpoints are d l and d r . In addition, the upper border and the lower border of the calcified plaque are the junctions between the calcified plaque and non-calcified tissues, which may be influenced by the noise and tissues with high level echo. Therefore, in order to reduce above influence, the row coordinates of the endpoints are not specified beforehand, that implies the optimal border traced by the graph searching algorithm has the global minimum cost.
(2) Because the graph searching algorithm is separately applied to R 1 and R 2 , two cost functions should be computed applied to R 1 and R 2 , respectively. With respect to R 1 , the upper border can be considered as the transition between the bright region and the dark region in the IVUS image. Therefore, the cost function in R 1 can be constructed based on the gradient information, formulated as With respect to R 2 , the acoustic shadowing occupies most of its area, and moreover, the region of acoustic shadowing far from the catheter cannot provide useful information for the detection of the lower border; on the contrary, the noise in the region may affect the tracing of the lower border. Therefore, we extract a subregion R' 2 from R 2 for tracing the lower border of the calcified plaque in order to reduce the influence from the noise. R' 2 is defined as the region between the d l th column and the d r th column, and between MIC and the p'th row computed in Equation (46). R 1 , R 2 and R' 2 acquired from I are shown in Figure 3(a) and Figure 3(b). Because the lower border can also be considered as a dark-bright transition, the cost function in R' 2 can be defined similarly to Equation (4) as The cost function M in Figure 3(b) is shown in Figure 3(c).
The results of the graph searching algorithm is shown in Figure 3(d), where the calcified plaque is the region inside the yellow contour. In order to reconstruct the border of calcified plaque in original IVUS images, it is transformed from polar domain to Cartesian domain, shown as the yellow contour in Figure 3(e).

Results
The data acquisition was performed with the approval of the Health Science Research Ethics Committee of the Department of Cardiology at Guangdong General Hospital, and the participants provided written informed consent before beginning the acquisition. All the IVUS data were collected over a Volcano machine (Volcano Corp., Rancho Cordova USA) with a pullback speed of 0.5mm/s using a 20MHz solid-state IVUS catheter. The IVUS images were reconstructed and processed at a commercially available IVUS console (In-Vision Gold, Volcano Corp., Rancho Cordova USA) and saved as DICOM format into CDs for off-line analysis. All the codes in this study were implemented by Matlab R2012a on a desktop computer with Intel(R) Xeon(R) CPU E5-2650(2.00 GHz) and 32GB memory. In our experiments, calcified plaques detected by our method were compared with those manually detected by one cardiologist. Totally 996 in-vivo IVUS images are used in our study, where 498 images contain the calcified plaque with acoustic shadowing (considered as the experimental group), and other 498 images do not contain any calcified plaque with acoustic shadowing (considered as the control group), and there are three kinds of imaging diameter (8mm, 10mm and 12mm) in these IVUS images. The selection protocol of the IVUS images contains: (1) every calcified plaque in IVUS images should have the acoustic shadowing behind it; (2) all IVUS images are collected at different times from eight patients.
The results of detecting the calcified plaque in some IVUS images are shown in Figure 4. The performance of our method was evaluated in two aspects: one was to evaluate the ability of our method to correctly distinguish between IVUS images with and without the calcified plaque; the other was to evaluate the location accuracy of the calcified plaque. Firstly, the ability to correctly identify IVUS images with the calcified plaque or without was evaluated by the sensitivity and specificity. The sensitivity represents the percentage of correctly identified IVUS images with the calcified plaque in the experimental group, and the specificity represents the percentage of correctly identified IVUS images without the calcified plaque in the control group, which can be computed by true positive (TP), false positive (FP), true negative (TN) and false negative (FN):  where TP is the number of IVUS images with the calcified plaque correctly identified, TN is the number of IVUS images without the calcified plaque incorrectly identified, FP is the number of IVUS images with the calcified plaque incorrectly identified, and FN is the number of IVUS images without the calcified plaque correctly identified. Table 1 shows the sensitivity and the specificity tested over all patients are 94.68% and 95.82%, respectively. Secondly, the location accuracy of correctly identified calcified plaques was quantified by five measurements: plaque area (PA), angular length (AL), angular center (AC), plaque thickness (PT) and the distance to catheter (DC), shown in Figure 5. PA can measure the size of the calcified plaque. AL and AC can measure the angular location (the location along the angular direction) of the calcified plaque between two angles h 1 and h 2 , formulated as PT and DC can measure the radial location (the location along the radial direction) of the calcified plaque, where DC equals the distance between the leading edge of the calcified plaque and the catheter's center. These measurements calculated from our method and the manual drawing method were compared separately by two different evaluation methods: linear regression and Bland-Altman analysis [39], in order to analyze their correlation and agreement, respectively. In every evaluation method, we firstly analyze the comparative results on all IVUS images in order to evaluate the overall performance of our method, and then on IVUS images from every patient in order to evaluate the between-patient difference of our method. In linear regression, the correlation coefficient r and the root-mean-square error (RMSE) are used to investigate the overall correlation between our method and the manual drawing method. Figure 6 shows the overall correlation. For AL, AC, PA, PT and DC, the values of r are 0.9557, 0.9747, 0.8311, 0.4461 and 0.9828, respectively, and the values of RMSE are 9.5055, 20.5260, 0.1964, 0.0674 and 0.0758, respectively. Figure 7 shows the betweenpatient difference for the five measurements. The standard deviation of r with respect to the five measurements is smaller than 0.15. Table 2 shows the corresponding numerical results. In the Bland-Altman analysis, four indices, r 1 , r 2 , r 3 and CI, are used to investigate the agreement between our method and the manual drawing method. Let X be the difference between same two measurements computed separately by our method and by manual, and let Y be the mean value of the two measurements. CI represents the 95% confidence interval of Y , formulated as (d{2s,dz2s), where d and s are the mean value and standard deviation of Y , respectively. r 1 represents the ratio of d to the mean value of X . r 2 represents the ratio of the maximum absolute value of Y to the mean value of X . r 3 represents the frequency of the points without CI in the Bland-Altman plot. Figure 8 shows the overall agreement between our method and the manual drawing method. With respect to the five measurements, v7.1% of results fall without CI in the Student t-test. Figure 7 shows the between-patient difference of our method. The standard deviation of the number of results within the 95% confidence interval with respective to the five measurements is smaller than 2.61%. Table 3 presents the corresponding numerical results. Additionally, the overlap between two calcified plaques detected by our method and drawn by manual separately can also represent the location accuracy of the calcified plaque, which is measured by the sensitivity and specificity. The sensitivity and specificity can also be calculated by Equation (6), where TP is the number of pixels inside the calcified plaque correctly identified, FP is the number of pixels inside the calcified plaque incorrectly identified, FN is the number of pixels outside the calcified plaque incorrectly identified, and TP is the number of pixels outside the calcified plaque correctly identified. Table 4 shows that the mean values of the sensitivity and specificity in the measurement of the overlap are 83.79% and 99.74%, respectively.
In addition, Table 5 compares the computational costs of our method separately applied to the experimental group and the control group. In the experimental group, the mean value and standard deviation of the computational cost are 6.81s and 0.21s, respectively. And in the control group, the mean value and standard deviation are 5.26s and 0.08s, respectively.

Discussion
In this section, the accuracy and robustness of our method, the computational cost, the limitation and the future work will be discussed.

Accuracy
The accuracy of our method is evaluated through comparing the results from the manual drawing by one cardiologist and from the computer-aided automatic detection by our method. The evaluation of accuracy can be divided into two parts: one is to test the accuracy of distinguishing between IVUS images with and without the calcified plaque; the other is to evaluate the location    accuracy of the calcified plaque drawn by manual and detected by our method. Firstly, the identification of the IVUS image with or without the calcified plaque can be considered as a binary classification, and thus the sensitivity can be used to represent the correct classification rate of the images with the calcified plaque and the specificity can be used to represent the correct classification rates of the images without the calcified plaque. In Table 1, the sensitivity and specificity computed from all IVUS images are 94.68% and 95.82%, respectively, representing that our method has good ability to distinguish between the normal vascular cross section and the pathological vascular cross section with respect to the calcification. Secondly, the location accuracy of the calcified plaque is evaluated separately by linear regression and Bland-Altman analysis. The linear regression focuses on the correlation between our results and manual drawing results by the correlation coefficient r and RMSE. In Table 2, we can see that our results and the manual drawing results have high statistical relationship in detecting the angular location of the calcified plaque because r equals 0.9557 for AL and equals 0.9747 for AC. And our method also has high correlation with the manual drawing method in detecting radial location of the calcified plaque because r equals 0.9828 for the distance DC from its leading edge to the catheter (shown in Figure 5). However, the correlation for PT has low level (r = 0.4461), implying that the detection of the trailing edge of the calcified plaque do not have enough accuracy (shown in Figure 5). That is because the ultrasonic signals transmitted from the catheter are almost completely blocked by the calcified plaque, forming a shadow behind the calcified plaque that makes it difficult to estimate the real thickness of the calcium. Moreover, the high correlation for PA (r = 0.8311) represents the sizes of the same calcified plaque detected by manual and our method are close to each other. In addition, the values of b 0 are all larger than zero for the five measurements, which means that the five measurements computed by our method is the underestimation of the manual drawing results. In the Bland-Altman analysis, we evaluate the agreement between our method and the manual drawing method. The values of r 2 for AL, PA, PT are 0.3030, 0.7590, 0.7184, respectively, representing that the maximum difference between our method and manual drawing method is not small, and however the values of r 1 for the three measurements (0.0630, 0.2113, 0.2295) show that the variation in measuring AL, PA and PT are small, implying that r 2 for AL, PA, PT are large only in a small amount of IVUS images. Moreover, r 3 ƒ7.03% for AL, AC, PA, PT and DC in the Bland-Altman analysis shows most of our results fall into 95% confidence interval with high reliability in the Student t-test. In addition, the overlapping rate between the calcium regions detected separately by manual and by our method can be measured by the sensitivity and specificity, which can represent how many image pixels belong to the calcified plaque. In Table 4, the mean values of the sensitivity and the specificity are respectively 83.79% and 99.74%, representing that the difference between the calcification regions obtained by manual and by our method is at a low level, and furthermore the number of correctly identified pixels inside the calcium region is smaller than the number of correctly identified pixels outside the calcium region, implying that the calcified plaque detected by our method is a little underestimation of the plaque manually drawn. In conclusion, we can consider that the accuracy of the detection of the calcified plaque by our method is at a high level.

Robustness
With respect to the robustness of our method, it is investigated from two aspects: the process of our method and the betweenpatient performance of our method. Firstly, the kernel of our method includes the pixel classification by the RMM, the angular location detection by the MRF, and the refinement of the calcified plaque by the predefined constraints. Due to the superposition of ultrasonic signals in the vessel and the influence from the speckle noise, the pixel intensity of different tissues in the IVUS image can be described by the probabilistic model. Comparing with the widely used Gaussian mixture model [40] [41], the RMM is more appropriate to describe the distribution of pixel intensity in the IVUS image because the statistics of the pixel intensity relies on the speckle induced by a lot of randomly located scatters, and can be well modeled by the Rayleigh distribution. The pixel classification by the RMM facilitates to the subsequent step of the proposed method for the location of the calcified plaque, because it can make the class number of pixels decrease from the level of the gray intensity to the component number of RMM, which can also be considered to improve the between-class difference and reduce the within-class difference of the image pxiels. Moreover, unlike the RMM applied in IVUS previously [9] [10] [42], we add the spatial relationship into the RMM classification. Because the vascular tissues have the spatial continuity in anatomy, the pixel values of the neighboring pixels in IVUS images are statistically dependent. The dependence can help to improve the robustness of the RMM classification to the noise because the prior probability generated by the dependent relationship in Equation (10) -Equation (13) has an effect on the local smooth. Similarly, the calcified plaque can also be considered to have the angular spatial relationship because it grows along the vessel wall. Therefore, we consider that the calcified plaque has the Markovianity along the angular axis of the IVUS image in polar coordinate, and the MRF can be used to determine the angular location of the calcified plaque based on the results of the RMM classification. However, the angular location of the calcified plaque is still influenced by the noise and error in the RMM classification. Therefore, the refinement of the calcified plaque is needed. In the refinement, five predefined constraints are proposed by considering different aspects of the calcified plaque: the gray intensity, the growth direction, the angular length and the corresponding acoustic shadowing. And any calcified plaque not satisfying these constraints can be considered as the pseudo calcification and then removed. In addition, different some previous works [13], our detection method does not rely on the beforehand detection of the lumen border and/or the media-adventitia border, which can save the computation cost and moreover avoid the error in the detection of the calcified plaque resulting from the error of the vessel border detection.
Secondly, the between-patient difference of our method is also investigated to further examine the robustness of our method. In the linear regression and the Bland-Altman analysis, the performances of our method tested separately from Patient 1 to Patient 8 are compared. Table 2 presents the standard deviation of the correlation coefficient r in linear regression for every measure (AL, AC, PA, PT and DC) in different patients are smaller than 0.15. Table 3 shows the standard deviation of r 1 , r 2 and r 3 in the Bland-Altman analysis for every measure in different patients are smaller than 0.11, 0.23, and 2.61%, respectively, which means that there is no significant between-patient difference of our method with respect to the correlation and the agreement. In addition, the overlap between the calcified plaques detected by our method and by manual drawing in different patients is shown in Table 4, where the standard deviation of the mean values of the sensitivity Table 5. The computing time of the proposed method is formulated as mean + standard deviation measured by the second.  Table 4. The results of the overlapping rate between the calcified plaques detected by our method and drawn by manual. The sensitivity represents the percentage of the pixels correctly identified inside the calcified plaque, and the specificity represents the percentage of the pixels correctly identified outside the calcified plaque. doi:10.1371/journal.pone.0109997.t004 and specificity (3.85% and 0.09%) represents that overlapping rate has small difference in different patients. Therefore, we can conclude that our method has similar performance in different patients.

Computational Cost
The computational cost of our method is presented in Table 5. Firstly, the computation cost over the experimental group is compared with the control group. The result shows that our method runs faster about 1.6s in the control group than in the experimental group, because the refinement of the calcified plaque will spend less time and the border tracing by the graph searching will not be executed when the IVUS image does not contain the calcified plaque with acoustic shadowing. Secondly, by comparing different parts in our method, we can see that the pixel classification by the RMM and angular location detection of the calcified plaque by MRF take the two most computational costs in the process of our method (35% and 49% in the experimental group and 39% and 53% in the control group), because the two parts of our method contain plenty of iteration operations. Moreover, the standard deviations of the cost in each parts of our method are small (v0.14s), and that implies the computational cost is hardly affected by the variation of the plaque size. In addition, the cost of our method (''others'' in Table 5) excluding the RMM classification, the angular location detection by MRF, the refinement and the border detection spends 13% and 9% of the total cost in the experimental group and the control group, respectively, and it might be reduced by optimizing the programming code in order to improve the computational efficiency.

Future Work
In the future, we will improve the efficiency of our method in two aspects: one is to improve the model for IVUS; the other is to improve the prior knowledge about the plaque. Firstly, we will try to apply different probabilistic mixture models for more appropriately describing the pixel intensity in the IVUS images, and use more complex graphical models to instead of the MRF in order to improve the detection accuracy of the calcified plaque. Secondly, we will add some high-level prior knowledge into the RMM in order to consider each kind of tissues as a whole for facilitating the pixel classification. In addition, in order to combine the above two aspects, the supervised scheme can be used by training the plaque models based on the more complex predefined features of the calcified plaque for better performance of our method.

Conclusion
The computer-aided detection of the calcified plaque with acoustic shadowing is vital in the quantitative analysis of IVUS images. In this paper, a method based on the Rayleigh mixture model, the Markov random field and the graph searching algorithm is proposed for the automatic detection of the calcified plaque. Additionally, the proposed method includes a refinement strategy to remove pseudo calcified plaques based on the predefined constraints about the calcification. The performance evaluation between the proposed automatic detection method and the manual drawing method (by one cardiologist) was performed over 996 in-vivo coronary IVUS images acquired from 8 patients with two aspects. The first aspect is the evaluation of the accuracy to correctly identify the IVUS images with or without the calcified plaque, and it was measured by the test of the sensitivity and specificity. The results shown that there is high correct rate of identifying the IVUS image including or excluding the calcified plaque with sensitivity 94.68% and specificity 95.82%. The other aspect is the evaluation of the accuracy to locate the calcified plaque by the linear regression, the Bland-Altman analysis and the overlapping rate between the calcified plaques separately detected by our method and by manual drawing. The linear regression shown that the purposed automatic method is correlated very well with the manual drawing method over four measurements (angular length, angular center, plaque area and the distance to catheter) with high correlation coefficients, and the Bland-Altman analysis demonstrated a high level of agreement between our method and the manual drawing method with w93% results falling within the 95% confidence interval. In addition, the sensitivity (83.79%) and specificity (99.74%) used to measure the overlapping rate shown that the calcified plaques detected by our method and drawn by manual measured are matched well. These results demonstrated the effectiveness of our method in the detection of the calcified plaque with acoustic shadowing in IVUS images.

A Constructing the RMM model
The RMM can be formulated as, where Y i is the ith pixel in the IVUS image, P(Y i ) is the probability density function of Y i in the RMM, p ij is the prior probability, K is number of mixture components in the RMM, v j is the parameter vector of the jth component, and P(Y i Dv j ) is the probability density function of Y i in the jth component. P(Y i Dv j ) can be formulated by the Rayleigh distribution with parameter v j [43], where s j is the mode, a j is the translation component and v j~f a j ,s j g.
The construction of the prior p ij in the RMM considers the neighboring relationship among pixels in IVUS images, inspired by [44] and computed as follows: 1. Define the weight of the ith pixel in the jth component of the RMM as j j (Y i ), where c j and b j (j~1,2,::,K) are parameters, and Y i is the mean value of pixels in the neighborhood N i of the ith pixel, formulated as where M is the number of pixels in the neighborhood N i . 2. Define the weight of the ith pixel and its neighborhood N i in the jth component of the RMM as q j (Y i ), where a is used to control the function shape of q j (Y i ).
3. Define the prior probability p ij as In Equation (13), the normalization of q j (Y i ) aims to translate the range of p ij into ½0,1.

B Formulating the E-step and the M-step in the EM algorithm
In the E-step, the objective function of the EM algorithm for the RMM can be formulated as where H (t) is the known parameters in the tth iteration, H is the unknown or updating parameters in the tth iteration, P(v j DY i ,H (t) ) is the posterior distribution and P(Y i Dv j ) is class-conditional probability density. Because the first term in the righthand of Equation (14) does not contain the parameters fv j g, and the second term of that does not contain the prior probabilities fp ij g, the maximization of the objective function Q(H,H (t) ) can be simplified to separately maximize the first term in the righthand of Equation (14) with respect to the parameter set fv j g, and to maximize the second term of that with respect to the prior fp ij g. In the M-step, we employ the steepest-descent method [45] to maximize the objective function Q(H,H (t) ) in order to update the parameter set H, where g~0:01 is the learning rate. Let Q be the abbreviation of  (16) can be calculated as follows: 1. The partial derivative of Q with respect to a j is 2. The partial derivative of Q with respect to s j is 3. The partial derivative of Q with respect to c j is 4. The partial derivative of Q with respect to b j is 5. The partial derivative of Q with respect to a is LQ La~X

C Estimating the RMM parameters by the EM algorithm
The estimation of the RMM parameters by the EM algorithm can be divided into the following steps: Step 1: initialize the parameter set H.
(1) The class number K is set at 5 and a is initialized as a (0)~2 . Then the k-means algorithm [24] is used to cluster pixels in the IVUS image based on the gray intensity in order to find the mean gray-scale values of K classes m (3) The relationship among the mean value m, the mode s and the translation component a are deduced from the Rayleigh distribution [46] is So the mode s (0)~½ s (0) 1 ,s (0) 2 ,:::,s (0) K can be computed from Equation (22) as (4) Let c (0)~a(0) and b (0)~s(0) .
Step 4: reach the convergence.
After the tth iteration, if DH (tz1) {H (t{1) Dv0:0001, the iteration process can reach the convergence. If not, return to the Step 2 and continue the iteration.

D Computing the prior of the MRF
The prior of the MRF can be computed as follows. Firstly, the dark region in I below MIC is located by the results of RMM classification, which contains the acoustic shadowing behind the calcified plaque. In I R , we can compute the mean gray-scale values of pixels from the Class 1 to the Class K, denoted by m Ã 1 ,m Ã 2 ,:::,m Ã K , respectively, where the two smallest values m Ã i ,m Ã j can be found, and the set of the pixels in the two classes corresponding to m Ã i and m Ã j can be considered as the dark region in I, denoted by R D . Then we define a binary image I D , as the same size as I, to indicate the location of R D , formulated as where I D is shown in Figure 2(f). The upper border of R D can be defined as the boundary between R D and the other region in I D , denoted by L u~f q 1 ,q 2 ,:::,q W g, where q i is the row coordinate value of the pixel on the upper border and in the ith column of I D . L u can be computed as where L u is shown in Figure 2(f). Secondly, the relationship between Z 1 ,Z 2 ,:::,Z W and Y 1 ,Y 2 ,:::,Y W can be influenced by three factors: the distance from the upper border L u to the edge of the imaging area in IVUS, the distance from L u to MIC, and the mean gray-scale value in the region between L u and MIC, which are denoted by

E Solving the MRF by the BP algorithm
The BP algorithm can be divided into three parts: the initialization of the MRF and belief propagation messages, the iterative process and the convergence criterion.
(1) Initialization For i~1,2,:::,W , the marginal probabilities of Z i~z 1 and Z i~{ 1 are firstly initialized to be equal, that is, the marginal probabilities are formulated as P(Z i~{ 1)~P(Z i~z 1)~0:5 ð29Þ and the belief at Z i is initialized as the marginal probabilities, Then the local evidence w(Z i ,Y i ) can be represented as the prior of the BP, that is, quantified by F~ff 1 ,f 2 ,:::,f W g in Equation (28) as where fw i (Z i ,Y i )g are invariant in the iteration. Next, if Z j is in the neighborhood of Z i , the relationship between Z i and Z j can be represented by the compatible function y(Z i ,Z j ), formulated as where y i,j (Z i ,Z j ) are not updated in the iteration. In addition, the message from Z i to Z j , denoted by m i,j (Z j ), is initialized as (2) Iteration If Z j is in the neighborhood of Z i , we can compute the message m (t) i,j (Z j ) from Z i to Z j at the tth iteration as Then the belief b (t) i (Z i ) of Z i can be computed as where N i is the neighborhood of Z i . The superscript ''(t)'' in Equation (34) (35) represent the corresponding variables are in the tth iteration. (

3) Convergence
The convergence criterion is formulated as follow, where e is set to be 0.001 in our method. If Equation (36) is satisfied, the iteration can be terminated. where T 3~1 in our method. If Equation (40) is not satisfied, the calcified plaque can be considered pseudo and removed out of R c .
(3) Constraint 3: the length of the calcified plaque along the angular axis in I can be represented as W c~d r {d l , and then the length should be smaller than the threshold T 4 , formulated as where T 4~3 00 in our method. If Equation (41) is not satisfied, the calcified plaque can be considered pseudo and removed out of R c . (4) Constraint 4: two columns s l and s r in I are firstly computed as Then the left region, right region and shadowing region of the calcified plaque can be defined based on s l and s r , shown in Figure 3 (2) similarly, the right region is surrounded by four borders: the left border, right border and lower border are three lines between (d r ,p d r ) and (d r ,H), between (s r ,p s r ) and (s r ,H), between (d r ,H) and (s r ,H), respectively, and the upper border is a part of MIC between (d r ,p d r ) and (s r ,p s r ); (3) the shadowing region is also surrounded by four borders: the left Figure 9. The left and right yellow points represent the points with coordinate (e,p e ) and (g,p e ), respectively. The red curve and the green curve are L m and L u , respectively. The region below L u is R D . N a is the number of pixels inside the intersection between the blue dashed rectangle and R D . N u is the number of pixels inside the intersection between the red dashed rectangle and R D . doi:10.1371/journal.pone.0109997.g009 border, right border and lower border are three lines between (d l ,p d l ) and (d l ,H), between (d r ,p d r ) and (d r ,H), between (d l ,H) and (d r ,H), respectively, and the upper border is a part of MIC between (d l ,p d l ) and (d r ,p d r ). Next, the mean gray values of pixels within the shadowing region, the left region and the right region of the calcified plaque are computed, denoted by m m , m l and m r , respectively, and the three mean values should satisfy m l wm m , m r wm m 1 2 (m l zm r ){m m wT 5 ð43Þ where T 5~2 0 in our method. If Equation (43) is not satisfied, the calcified plaque can be considered pseudo and removed out of R c .
Similarly, with respect the right border d r , another pixel (g,p g ) on MIC satisfying (gwd r ) Vi,j[½d r ,g, if jwi, then p i §p j ð45Þ Then the maximum value between p e and p g can be found, formulated as p'~max (p e ,p g ) ð46Þ Next, let N u be the number of pixels in the intersection between R D and a rectangle with vertices (e,p'), (g,p'), (e,H), (g,H), and N a be the number of pixels in the intersection between R D and a rectangle with vertices (e,p'), (g,p'), (e,1), (g,1). The regions correspond to N a and N u are shown in Figure 9, and the following conditions should be satisfied: N a =(g{e)wT 6 ð47Þ where T 6~0 :3 and T 7~5 in our method. If Equation (47) and Equation (48) are not both satisfied, then the calcified plaque can be considered pseudo and removed out of R c .