Carotid Artery Wall Motion Estimated from Ultrasound Imaging Sequences Using a Nonlinear State Space Approach

Abstract

It is very challenge to investigate the motion of the carotid artery wall in ultrasound images, because of the high nonlinear dynamics of this motion. In our study, the nonlinear dynamics of carotid artery wall motion is first approximated by our nonlinear state-space approach driven by a mathematical model of the mechanical deformation of carotid artery wall. Then, the two-dimensional motion of carotid artery wall is computed by solving the nonlinear state-space approach using the unscented Kalman filter. We have then evaluated the performance of our approach by comparing it with the manual tracing method (the correlation coefficient equals 0.9897 for the radial motion and 0.9703 for the longitudinal motion) and three other state-of-the-art methods for 73 subjects. The results indicate the reliable applicability of our approach in tracking the motion of the carotid artery wall and its potential usefulness in routine clinical diagnosis.

Z. Gao and Y. Sun—Contributed equally to this work.

1 Introduction

Early diagnosis of atherosclerotic disease is very important since it can direct therapies to prevent its complications. Traditional measurements (such as intima-media thickness), however, cannot effectively detect the preclinical and subclinical lesions due to the insignificant morphological changes in the vessel during this period. On the other hand, the motion of carotid arteries (Fig. 1) can be used to identify the atherosclerotic disease in the early stage by characterizing the artery stiffness [1]. In particular, it has been recently proved that the longitudinal motion of carotid artery wall is associated with atherosclerotic disease, and could be sensitive to show its early lesions [2]. For instance, the magnitude of longitudinal motion has been shown to be reduced in subjects with carotid plaques, type 2 diabetes and periodontal disease [2]. More importantly, in recent literatures, the longitudinal motion of the carotid artery wall was considered as a potential predictor of cardiovascular events in a large population screening, as it was associated with traditional risk factors while being independent of traditional risk markers [3].

However, manual tracing of the motion trajectories of the carotid artery wall is labor-intensive and time-consuming. Furthermore, it is very challenging to extract the longitudinal motion of the carotid artery wall from ultrasound image because of the poor longitudinal spatial resolution of ultrasound probe. Different block matching (BM) methods have been frequently applied to track the motion of the carotid artery wall. In the previous BM approaches, many similarity measures, including normalized cross correlation [4] and sum of square difference [3], have been proposed to locate the target block in the searching region. However, if those similarity measures are used directly without considering the noises in the reference and targeting images, the conventional BM method might not be able to locate the best-matched block in the targeted image [4]. Therefore, a recent study adopted the Kalman filter to estimate the motion of carotid artery wall by treating the noises as the uncertainties in the state space [3]. However, the movements of the block between two neighboring imaging frames are modeled as a constant process, and this may provide a biased motion of the carotid artery wall because of its nonlinear dynamics.

In order to properly handle the nonlinear arterial wall dynamics, we have developed a nonlinear periodic function with a deterministic input from a mathematic model of the arterial wall deformation, in order to appropriately match the dynamics of the carotid artery wall [5, 6]. Then, we have built a nonlinear state-space approach based on this nonlinear periodic function, and solved this problem by combining the BM method and unscented Kalman filter (UKF). We have then validated the performance of our approach by comparing the results of our method with the manual tracing method and three state-of-the-art methods [3, 4, 6] on monitored ultrasound imaging sequences from 73 subjects.

2 Methodology

In this study, a nonlinear state space approach have been developed to track the two-dimensional (2D) motion of the carotid artery wall. In this state space approach, the nonlinear dynamics of carotid artery wall is approximated by using a nonlinear periodic function [7], and the biased motion trajectories are corrected by means of a mathematical model of carotid artery wall [5]. Then, we have solved the state space problem by combined the UKF and the BM method together in this work. In order to reduce the influences from the out-of-plane motion, weak echoes, and motion artifacts during the ultrasound imaging, three user-defined reference blocks were selected in the first frame of each ultrasound sequence before the motion tracking, shown in Fig. 1. The size of these reference blocks were all 0.5 mm $\times$ 1.7 mm. Then, three search regions were selected by making their centers the same as the corresponding reference blocks. All search regions were set at 1.3 mm $\times$ 2.5 mm.

Fig. 1.

(a) Sample carotid ultrasound image. The blue curve and the yellow curve are respectively the lumen border and the media-adventitia border of the carotid artery. The green rectangles are user-defined reference blocks, and the red rectangles are the tracked blocks corresponding to the green rectangles. (b), (c), (d) and (e) the examples of the blocks’ locations tracked by BM, KBM, OP and our approach in the same ultrasound frame. (Color figure online)

Prediction model. In tracking the 2D motion of carotid artery wall, the dynamics of the reference block in the ultrasound sequence was described by the discrete-time nonlinear state space equations, as

$$\begin{aligned} \begin{aligned} \mathbf {x}_{n+1} = f(\mathbf {x}_n,\mathbf {u}_n) + \mathbf {w}_n,\quad \mathbf {y}_n = \mathbf {x}_n + \mathbf {v}_n \end{aligned} \end{aligned}$$

(1)

where n is the time index or frame index. $\mathbf {w}_n$ and $\mathbf {v}_n$ are the noise obeyed by gaussian distribution with covariance matrices $\mathbf {Q}$ and $\mathbf {R}$, that is $\mathbf {w}_n \sim \mathcal {N}(0,\mathbf {Q})$ and $\mathbf {v}_n \sim \mathcal {N}(0,\mathbf {R})$. $\mathbf {Q} = 0.01\mathbf {I}$ and $\mathbf {R} = 100\mathbf {I}$, where $\mathbf {I}$ is the identity matrix. $\mathbf {x}_n = [x^1_n,x^2_n]^\mathrm {T}$ is the center point’s location of reference block $\mathcal {B}^{ref}_n$ in the nth frame, and $x^1_n, x^2_n$ are the y-coordinate and the x-coordinate respectively. Similarly, $\mathbf {x}_{n+1} = [x^1_{n+1},x^2_{n+1}]^\mathrm {T}$ is the center point’s location of reference block $\mathcal {B}^{ref}_{n+1}$ in the $n+1$th frame, and $\mathbf {y}_n = [y^1_n,y^2_n]^\mathrm {T}$ is the center point’s location of reference block $\mathcal {B}^{best}_n$. $\mathbf {u}_n = [u_1,u_2]^\mathrm {T}$ is the input signal inspired by a mathematical model of carotid arterial wall’s mechanical deformation [5, 6]. Therein, T denotes the number of ultrasound frames sampled in a cardiac cycle. $\mathbf {u}_n = \mathbf {g}_1(n)$ when $n\leqslant 0.4T$, and $\mathbf {u}_n = \mathbf {g}_2(n)$ when $0.4T<n\leqslant T$, where

$$\begin{aligned} \mathbf {g}_1(n) = \left[ \begin{array}{l} 0.25\tanh (1.22(n-0.4T))(1-\tanh (1.22n))\sin ^2\frac{\pi n}{3.5T} \\ 0.25\tanh (2.07(n-0.4T))(1-\tanh (2.07n))\sin ^2\frac{\pi n}{0.65T} \end{array}\right] \end{aligned}$$

(2)

and

$$\begin{aligned} \mathbf {g}_2(n) = \left[ \begin{array}{l} 0.25(1+\tanh (1.22(n-T)))(1+\tanh (1.22(0.39T-n)))(1.15-0.4n) \\ 0.25(1+\tanh (2.07(n-T)))(1+\tanh (1.22(0.39T-n)))(7-0.2n) \end{array}\right] \end{aligned}$$

(3)

Fig. 2.

Bland-Altman plot for the radial motion (left) and the longitudinal motion (right). The x-axis and y-axis represent the average and the difference between the motion positions separately tracked by our approach and manual tracing method. The blue solid line shows the bias between our approach and manual tracing method. The intervals between green dotted lines are the 95 % confidence interval. The black dotted lines are linearly fitting lines of all the scatter points.

The nonlinear function f in Eq. (1) is defined, based on [7], as below,

$$\begin{aligned} \begin{aligned}&x^1_{n+1} = x^1_n - \frac{T}{5} x^2_n + Tu_1 \\&x^2_{n+1} = (1 - \frac{T}{5}) x^2_n + Tx^1_n - \frac{T}{373}x^1_n x^2_n + Tu_2 \end{aligned} \end{aligned}$$

(4)

Then the state space equations in Eq. (1) are solved to obtain the best estimate of $\mathbf {x}_{n+1}$ (denoted by $\hat{\mathbf {x}}_{n+1}$) in the $n+1$th frame by using the unscented Kalman filter (UKF) [8]. The UKF samples a set of points (called sigma points) around the mean and then approximates the nonlinear transformation by the evolution of the sigma points with time. Thus, it can reduce the error from the linearized transformation of the nonlinear function. The process of UKF can be divided into two steps: state updating and observation updating. In this process, the state of the reference block $\mathbf {x}_n$ in the nth frame and $\mathbf {x}_{n+1}$ in the $n+1$th frame are respectively considered to be the posterior state $\hat{\mathbf {x}}^+_{n}$ and $\hat{\mathbf {x}}^+_{n+1}$ in UKF.

State updating. State updating estimates the priori state $\hat{\mathbf {x}}^-_{n+1}$ and priori error covariance $\hat{\mathbf {P}}^-_{n+1}$ of the reference block $\mathcal {B}^{ref}_{n+1}$ in the $n+1$th frame from the information in the last frame, formulated as

$$\begin{aligned} \begin{aligned} \hat{\mathbf {x}}^-_{n+1} = \frac{1}{2M} \sum _{s=1}^{2M} \hat{\mathbf {x}}^s_{n+1},\quad \hat{\mathbf {P}}^-_{n+1} = \frac{1}{2M} \sum _{s=1}^{2M} \left( \hat{\mathbf {x}}^{s}_{n+1} - \hat{\mathbf {x}}^-_{n+1}\right) \left( \hat{\mathbf {x}}^{s}_{n+1} - \hat{\mathbf {x}}^-_{n+1}\right) ^\mathrm {T} + \mathbf {Q} \end{aligned} \end{aligned}$$

(5)

where M is the dimension of the state $\mathbf {x}$. $\hat{\mathbf {x}}^s_{n+1}$ are the sigma points in the $n+1$th frame, and can be obtained from the sigma points $\hat{\mathbf {x}}^s_{n}$ in the nth frame by the nonlinear function f, i.e. $\hat{\mathbf {x}}^s_{n+1} = f(\hat{\mathbf {x}}^s_{n},\mathbf {u}_n)$. Therein, $\hat{\mathbf {x}}^s_{n}$ can be computed by the posterior state $\hat{\mathbf {x}}^+_{n}$ and posterior error covariance $\hat{\mathbf {P}}^+_{n}$ in the nth frame

$$\begin{aligned} \begin{aligned}&\hat{\mathbf {x}}^{s}_n = \hat{\mathbf {x}}^{+}_{n} + \check{\mathbf {x}}^{s}_n,\quad s=1,2,...,2M \\&\check{\mathbf {x}}^s_n = \left( \sqrt{M\mathbf {P}^+_{n}}\right) _s,\quad \check{\mathbf {x}}^{M+s}_n = -\left( \sqrt{M\mathbf {P}^+_{n}}\right) _s,\quad s=1,2,...,M \end{aligned} \end{aligned}$$

(6)

where $\left( \sqrt{M\mathbf {P}^+_{n}}\right) _s$ is the ith row of the matrix $\sqrt{M\mathbf {P}^+_{n}}$.

Fig. 3.

Errors of our approach for the radial motion (left), the longitudinal motion (middle) and the 2D motion (right) in different noise levels. The x-axis represents the signal-noise ratio of the Rayleigh noise added into the ultrasound dataset. The red lines show the errors without adding the noise. (Color figure online)

Observation updating. Observation updating estimates the posterior state $\hat{\mathbf {x}}^+_{n+1}$ and priori error covariance $\hat{\mathbf {P}}^+_{n+1}$ of the reference block $\mathcal {B}^{ref}_{n+1}$ in the $n+1$th frame, formulated as

$$\begin{aligned} \begin{aligned}&\mathbf {K}_{n+1} = \mathbf {P}_{xy}\mathbf {P}_y^{-1},\quad \hat{\mathbf {x}}^+_{n+1} = \hat{\mathbf {x}}^-_{n+1} + \mathbf {K}_{n+1}(\mathbf {y}_{n+1} - \hat{\mathbf {y}}_{n+1}) \\&\hat{\mathbf {P}}^+_{n+1} = \hat{\mathbf {P}}^-_{n+1} - \mathbf {K}_{n+1} \mathbf {P}_y \mathbf {K}_{n+1}^\mathrm {T} \end{aligned} \end{aligned}$$

(7)

where $\mathbf {P}_y$ is the covariance of the measurement noise, $\mathbf {P}_{xy}$ is the cross covariance between the state variable and the measurement variable and $\hat{\mathbf {y}}_{n+1}$ is the estimate of the measurement variable. They can be obtained as:

$$\begin{aligned} \begin{aligned}&\hat{\mathbf {y}}_{n+1} = \frac{1}{2M} \sum _{s=1}^{2M} \hat{\mathbf {y}}^s_{n+1},\quad \hat{\mathbf {y}}^s_{n+1} = \hat{\mathbf {x}}^s_{n+1} \\&\mathbf {P}_y = \frac{1}{2M} \sum _{s=1}^{2M} \left( \hat{\mathbf {y}}^s_{n+1} - \hat{\mathbf {y}}_{n+1} \right) \left( \hat{\mathbf {y}}^{s}_{n+1} - \hat{\mathbf {y}}_{n+1}\right) ^\mathrm {T} + \mathbf {R} \\&\mathbf {P}_{xy} = \frac{1}{2M} \sum _{s=1}^{2M} \left( \hat{\mathbf {x}}^{s}_{n+1} - \hat{\mathbf {x}}^-_{n+1}\right) \left( \hat{\mathbf {y}}^{s}_{n+1} - \hat{\mathbf {y}}_{n+1}\right) ^\mathrm {T} \end{aligned} \end{aligned}$$

(8)

Measurement acquisition. The measurement $\mathbf {y}_{n+1}$ in Eq. (7) is the center point’s location of the best-matched block $\mathcal {B}^{best}_{n+1}$ in the $n+1$th frame. Then, $\mathcal {B}^{best}_{n+1}$ can be computed from the reference block $\mathcal {B}^{ref}_{n+1}$ in the $n+1$th frame of the ultrasound sequence, by the block matching (BM) method with normalized cross-correlation [4].

3 Experiments and Results

A total of 14 healthy subjects (7 males, 62.0 ± 8.7 age years; 7 females, 59.7 ± 7.8 age years) and 59 unhealthy subjects (37 males, 59.0 ± 13.6 age years; 22 females, 60.6 ± 12.2 age years) were involved in this study. The inclusion criterion for the unhealthy subjects was the presence of one of the three diseases diagnosed: type 1 or 2 diabetes, hypertension and heart disease. The data collection was implemented by an ultrasound physician with more than 10-year experience by a high-resolution ultrasound system (iU22, Philips Ultrasound, Bothell, WA, USA) and a 7.5 MHz liner array transducer.

Fig. 4.

Example of the motion trajectories of three reference blocks on the carotid artery wall from a healthy subject (left three columns) and an unhealthy subject (right three columns) for manual tracing method (green), as well as for our approach (red), CBM (blue), KBM (black) and OP (magenta). The x-axis of each subfigure represents the time of the ultrasound sequence (second). The y-axis of each subfigure represents the displacement of the block (mm) with respect to its position in the first frame. (Color figure online)

Accuracy. The accuracy of our approach was evaluated by comparing the results computed by our approach with those manually performed by an experienced ultrasound physician (considered as the ground truth). The physician identified the center points for each of the three blocks in the first frame and then visually detected the three center points in every single subsequent frame. After the tracking process, two measurements were computed to show the results of our approach: the positions of the radial motion (RP) and the longitudinal motion (LP). The radial motion denotes the dynamics of the best-matched block along the radial direction across the ultrasound sequence, and the longitudinal motion denotes that in the longitudinal direction. The positions of radial/longitudinal motion show the location of the best-matched block in every frame. Three kinds of errors (defined in [9]) were used to measure the accuracy of our approach: errors of radial motion, longitudinal motion and 2D motion. The results of linear regression showed that our approach was highly correlated with the manual tracing method for different subjects (RP: $r=0.9897\pm 0.0140$, p-value$\,<\,$0.001; LP: $r=0.9703\pm 0.0539$, p-value$\,<\,$0.001). The Bland-Altman plots shown in Fig. 2 illustrate that the bias of our approach with respect to the manual tracing method were 0.049 mm for RP and 0.0069 mm for LP. Then, we linearly fitted the scatter points in the Bland-Altman plots, and the slopes of the fitting straight lines (black dotted lines) were 0.0054 for RP and 0.0062 for LP. In addition, the coefficient of variation (CV) of the difference between our approach and the manual tracing method was 0.8681 for RP and 1.1046 for LP. These results indicate that difference between our approach and the manual tracing method is small, and moreover not much influenced by their averages. Finally, we added different levels of Rayleigh noise (1 dB $\sim$ 30 dB signal-noise ratio) in the ultrasound dataset and retested the proposed approach. Figure 3 shows that our approach has good performance to resist the noise disturbance.

Fig. 5.

Comparative results of the tracking error between our approach and the other methods (CBM, KBM and OP) for the radial motion (left), longitudinal motion (middle) and the 2D motion (right). The digits above the boxes in each figure show the p-values of the difference between the computer-aided method and the manual tracing method.

Comparison with other methods. Our approach was compared with three state-of-the-art methods on all subjects: conventional block matching method (CBM) [4], Kalman-based block matching method (KBM) [3] and optical flow method (OP) [6], with all parameters shown in their own publications except the sizes of the reference block and search region. The values of the two sizes used in the three methods are same as our approach. Figure 4 illustrates the motion trajectories of the carotid artery wall computed by our approach and other methods in a healthy subject and an unhealthy subject; Fig. 1(b)–(e) show examples of the blocks’ location obtained by our approach and other methods in the same ultrasound frame. Figure 5 displays the comparative results between our approach and other methods. The results shows that the accuracy of our approach (RP: 0.0895 ± 0.0777 mm, p-value$\,<\,$0.001; LP: 0.1619 ± 0.1789 mm, p-value = 0.939) was higher than those of other methods. This indicates that the use of the nonlinear periodic function and the mathematical model of the mechanical deformation of carotid artery wall in the prediction model was effective. Moreover, we compared the computational cost of our approach with those of other methods. The results showed that our approach (0.0257 ± 0.0046 second per frame) was moderate compared to the other methods (CBM: 0.0249 ± 0.0049 second per frame, KBM: 0.0772 ± 0.0048 second per frame, OP: 0.0101 ± 0.0011 second per frame).

4 Conclusion

In this study, we have developed an approach to automatically track the dynamics of the carotid artery wall. Our methodology involves developing a nonlinear periodic function with input from the mathematic model of mechanical deformation of the carotid artery wall, for the purpose of properly handling the nonlinear components of the dynamics. The performance of our approach was tested on a dataset with 73 carotid ultrasound sequences, and the results were compared with those manually performed by one experienced ultrasound physician and three state-of-the-art methods. The results demonstrated the effectiveness of our approach in the motion tracking of carotid artery wall in ultrasound image sequences.