Reliable and Accurate Calcium Volume Measurement in Coronary Artery Using Intravascular Ultrasound Videos

Abstract

Quantitative assessment of calcified atherosclerotic volume within the coronary artery wall is vital for cardiac interventional procedures. The goal of this study is to automatically measure the calcium volume, given the borders of coronary vessel wall for all the frames of the intravascular ultrasound (IVUS) video. Three soft computing fuzzy classification techniques were adapted namely Fuzzy c-Means (FCM), K-means, and Hidden Markov Random Field (HMRF) for automated segmentation of calcium regions and volume computation. These methods were benchmarked against previously developed threshold-based method. IVUS image data sets (around 30,600 IVUS frames) from 15 patients were collected using 40 MHz IVUS catheter (Atlantis® SR Pro, Boston Scientific®, pullback speed of 0.5 mm/s). Calcium mean volume for FCM, K-means, HMRF and threshold-based method were 37.84 ± 17.38 mm³, 27.79 ± 10.94 mm³, 46.44 ± 19.13 mm³ and 35.92 ± 16.44 mm³ respectively. Cross-correlation, Jaccard Index and Dice Similarity were highest between FCM and threshold-based method: 0.99, 0.92 ± 0.02 and 0.95 + 0.02 respectively. Student’s t-test, z-test and Wilcoxon-test are also performed to demonstrate consistency, reliability and accuracy of the results. Given the vessel wall region, the system reliably and automatically measures the calcium volume in IVUS videos. Further, we validated our system against a trained expert using scoring: K-means showed the best performance with an accuracy of 92.80 %. Out procedure and protocol is along the line with method previously published clinically.

Appendix

Automatic threshold-based algorithm

Otsu et al. [Otsu-IEEE-1979] developed a nonparametric and unsupervised method based on the maximization of the separability of the resultant classes. The procedure use the 0^th and the first-order cumulative moments of the gray level histogram. The procedure is very simple and effective, hence used as a threshold estimator in our work.

Let the pixels of any given image be represented in M gray levels [1, 2, ⋯, M]. The number of pixels at gray level i is denoted by n _i and the total number of pixels by N = n ₁ + n ₂ + ⋯ + n _M . After normalization the gray-level histogram may be regarded as a probability distribution:

$$ \begin{array}{cc}\hfill {p}_i=\frac{n_i}{N},\hfill & \hfill\ {p}_i\ge 0,{\displaystyle \sum_{i=1}^M}{p}_i=1\hfill \end{array} $$

(1)

Pixels are separated into two classes, C ₀ and C ₁ (object and background or vice versa), by a threshold at gray level k; C ₀ denotes pixels with grey levels [1, ⋯, N]. Then, the probabilities of class occurrence and the class mean gray levels, respectively, are given by:

$$ {w}_0={P}_r\left({C}_0\right) = {\displaystyle \sum_{i=1}^k}{p}_i=w(k) $$

(2)

$$ {w}_1={P}_r\left({C}_1\right)={\displaystyle \sum_{i=k+1}^M}{p}_i=1-w(k) $$

(3)

and

$$ {\mu}_0={\displaystyle \sum_{i=1}^k}i{P}_r\left(i\Big|{C}_0\right)={\displaystyle \sum_{i=1}^k}i{p}_i/{w}_0=\mu (k)/w(k) $$

(4)

$$ {\mu}_1={\displaystyle \sum_{i=k+1}^M}i{P}_r\left(i\Big|{C}_1\right)={\displaystyle \sum_{i=k+1}^M}i{p}_i/{w}_1 = \frac{\mu_T-\mu (k)}{1-w(k)} $$

(5)

where

$$ w(k)={\displaystyle \sum_{i=1}^k}{p}_i $$

and

$$ \mu (k)={\displaystyle \sum_{i=1}^k}i{p}_i $$

are the 0^th and first order cumulative moments of the histogram up to the k ^th level, respectively, and

$$ {\mu}_T=\mu (M)={\displaystyle \sum_{i=1}^M}i{p}_i $$

The optimum threshold is defined as the value that maximizes the between class variance

$$ {\sigma}_B^2(k) = \frac{{\left[{\mu}_Tw(k)-\mu (k)\right]}^2}{w(k)\left[1-w(k)\right]} $$

(6)

Thus the optimal threshold k * is given by:

$$ {\sigma}_B^2\left(k*\right)=\underset{1\le k<M}{ \max }{\sigma}_B^2 $$

(7)

K-means algorithm

Hartigan and Wong [Hartigan-JRSS-1979] proposed a K-means algorithm to divide n points into m clusters so that the within-cluster sum of squares is minimized. K-means algorithm is an unsupervised clustering algorithm that classifies the input data points into multiple classes based on their inherent distance from each other and hence used in our work. The algorithm assumes that the data features form a vector space and tries to find natural clustering in them. The points are clustered around centroids μ _i  ∀ i = 1 ⋯ m which are obtained by minimizing the objective

$$ V={\displaystyle \sum_{i=1}^m}{\displaystyle \sum_{x_j\in {S}_i}}{\left({x}_j - {\mu}_i\right)}^2 $$

(8)

where there are m clusters S _i , i = 1, 2, . . ., m and μ _i is the centroid or mean point of all the points x _j  ∈ S _i . Various steps of the algorithm are as follows:

• 1. First we compute the intensity distribution of the intensities.

• 2. Then we initialize the centroids with m random intensities.

• 3. We have to repeat the following steps until the cluster labels of the image do not change anymore.

• 4. Now, we cluster the point based on distance of their intensities from the centroid intensities.

$$ {c}^{(i)}= \arg \min {\left\Vert {x}^{(i)} - {\mu}_j\right\Vert}^2 $$

(9)

• 5. Finally, we compute the new centroid for each of the clusters.

$$ {\mu}_i = \frac{{\displaystyle {\sum}_{i=1}^m}I\left\{{c}^{(i)}=j\right\}{x}^{(i)}}{{\displaystyle {\sum}_{i=1}^m}I\left\{{c}^{(i)}=j\right\}} $$

(10)

Where m is a parameter of the algorithm (the number of clusters to be found), I iterates over all the intensities, j iterates over all the centroids and μ _i are the centroid intensities.

Fuzzy c-means algorithm

Zadeh [Zadeh-I&C-1965] proposed a fuzzy set theory and gives an idea of uncertainty of belongings, which is described by a membership function. One of the most widely used approaches to fuzzy prototype-based clustering is Fuzzy c-Means (FCM) clustering. The FCM algorithm is an extension of the K-means algorithm and provide a generalization of the FCM algorithm and hence used in our work. Given the data set X, we first choose the number of clusters c, the weighting exponent m, and the termination tolerance ε > 0. Now we compute the cluster centers:

$$ {v}_i^{(l)} = \frac{{\displaystyle {\sum}_{k=1}^N}{\left({\mu}_{ij}^{\left(l-1\right)}\right)}^m{X}_k}{{\displaystyle {\sum}_{k=1}^N}{\left({\mu}_{ij}^{\left(l-1\right)}\right)}^m},\ 1\le i\le c $$

(11)

Then we compute the distances:

$$ {D}_{ik}^2 = {\left\Vert {x}_k-{v}_i\right\Vert}^2 = \sqrt{{\left({x}_k-{v}_i\right)}^T\ \left({x}_k-{v}_i\right)},\ 1\le i\ \le c,\ 1\le k\ \le N $$

(12)

Now we update the partition matrix:if

$$ {D}_{ik}>0\ \mathrm{f}\mathrm{o}\mathrm{r}\kern1.75em 1\le i\le c,\ 1\le k\ \le N $$

(13)

$$ {\mu}_{ik}^{(l)} = \frac{1}{{\displaystyle {\sum}_{j=1}^c}{\left(\frac{D_{ik}}{D_{jk}}\right)}^{\frac{2}{m-1}}}\kern2em \mathrm{otherwise}\kern2em {\mu}_{ik}^{(l)}=0 $$

(14)

until

$$ \left\Vert {U}^{(l)} - {U}^{\left(l-1\right)}\right\Vert < \varepsilon $$

Where, U is the cluster center.

Hidden markov random field algorithm

Wang et al. [Wang-arXiv-2012] implement a MATLAB toolbox named HMRF-EM-image for 2D image segmentation. This tool can implement an edge-prior-preserving image segmentation and hence used in our work. HMRF-EM algorithm is given below:

• 1. First we start with initial parameters set Ø⁽⁰⁾.

• 2. Now we calculate likelihood distribution P ^(t)(y _i |x _i , θ _xi ).

• 3. Current parameter set Ø^(t) is used to estimate labels by MAP estimations.

• 4. Now posterior distribution is calculated for all l ∈ L and all pixel y _i

$$ {P}^{(t)}\left(l\Big|{y}_i\right) = \frac{G\left({y}_i;{\theta}_l\right)P\left(l\Big|{x}_{N_i}^{(t)}\right)}{P^{(t)}\left({y}_i\right)} $$

(15)

Where ${x}_{N_i}^{(t)}$ is the neighbourhood configuration of x ^(t) _i and ${P}^{(t)}\left({y}_i\right)={\displaystyle \sum_{l\in L}}G\left({y}_i;{\theta}_l\right)P\left(l\left|{x}_{N_i}^{(t)}\right.\right)$

Here we have

$$ P\left(l\left|{x}_{N_i}^{(t)}\right.\right)=\frac{1}{Z} \exp \left( - {\displaystyle \sum_{j\in {N}_i}}{V}_c\left(l,{x}_j^{(t)}\right)\right) $$

P ^(t)(l|y _i ) is used to update the parameters

$$ {\mu}_l^{\left(t+1\right)} = \frac{{\displaystyle {\sum}_i}{P}^{(t)}\left(l\Big|{y}_i\right){y}_i}{{\displaystyle {\sum}_i}{P}^{(t)}\left(l\Big|{y}_i\right)} $$

$$ {\left({\sigma}_l^{\left(t+1\right)}\right)}^2 = \frac{{\displaystyle {\sum}_i}{P}^{(t)}\left(l\Big|{y}_i\right){\left({y}_i - {\mu}_l^{\left(t+1\right)}\right)}^2}{{\displaystyle {\sum}_i}{P}^{(t)}\left(l\Big|{y}_i\right)} $$

For minimizing the total posterior energy we need to solve x*

$$ \begin{array}{l}{x}^{*}=\mathrm{argmin}\ \left\{U\left(y\Big|x,\O \right)+U(x)\right\}\\ {}\kern1.68em \mathrm{x}\in x\end{array} $$

(16)

where, the likelihood energy is

$$ U\left(\left(y\Big|x\right),\O \right)={\displaystyle \sum_i}U\left(\left({y}_i\Big|{x}_i\right),\ \O \right) = {\displaystyle \sum_i}\left[\frac{{\left({y}_i - {\mu}_{xi}\right)}^2}{2{\sigma}_{xi}^2}+ \ln {\sigma}_{xi}\right] $$

In the image domain, we assume that one pixel has at most 4 neighbours, than the clinic potential defined on the neighboring pixel is

$$ {V}_c\left({x}_i,\ {x}_j\right) = \frac{1}{2}\left(1-{I}_{xi,\kern0.5em xj}\right) $$

(17)

where

$$ {I}_{xi,\kern0.5em xj} = \left\{\begin{array}{c}\hfill 0\kern1.5em \mathrm{if}\ {\mathrm{x}}_{\mathrm{i}}\ne {x}_j\hfill \\ {}\hfill 1\kern1.5em \mathrm{if}\ {\mathrm{x}}_{\mathrm{i}}={x}_j\hfill \end{array}\right. $$

An iterative algorithm is designed to solve Eq. number 16:

• 1. From the previous loop of EM estimator we have initial estimate x ⁽⁰⁾.

• 2. Provided x ^(k), for all 1 ≤ i ≤ N, we find

$$ {x}_i^{\left(k+1\right)}= argmin\ \left\{U\left({y}_i\Big|l\right)+{\displaystyle \sum_{j\in {N}_i}}{V}_c\left(l,{x}_j^{(k)}\right)\right\} $$

(18)

• 3. Repeat the above step until U((y|x), Ø) + U(x) converges or a maximum k is achieved.