Tracking Carotid Artery Wall Motion Using an Unscented Kalman Filter and Data Fusion

Analyzing the motion of the common carotid artery (CCA) wall yields effective indicators for atherosclerosis. In this work, we propose a state-space model and a tracking method for estimating the time-varying CCA wall radius from a B-mode ultrasound sequence of arbitrary length. We employ an unscented Kalman filter that fuses two sets of measurements produced by an optical flow algorithm and a CCA wall localization algorithm. This fusion-and-tracking approach ensures that feature drift, which tends to impair optical flow based methods, is compensated in a temporally consistent manner. Simulation results show that the proposed method outperforms a recently proposed optical flow based method.


I. INTRODUCTION A. BACKGROUND, MOTIVATION, STATE OF THE ART
According to the World Health Organization (WHO), over 31% of all deaths in 2016 were caused by coronary events, most of which can be attributed to a single progressive disease called atherosclerosis [1], [2]. It is well known that the progression of atherosclerosis is associated with increasing arterial stiffness [3]. Several diagnostic methodologies exploit this association in order to estimate the risk of coronary events. These methodologies are generally based on measurements of quantities that are intrinsically associated with arterial stiffness [4], [5]. One of these quantities is the motion of the arterial wall in response to variations in blood pressure or blood flow, which occurs naturally in the human body.
This article presents a new methodology for the analysis of the motion of the wall of the common carotid artery (CCA) from a B-mode ultrasound sequence. In screening tests, methods based on the analysis of the motion of the CCA wall have been shown to produce effective indicators for atherosclerosis [6], [7]. The main challenge here is an accurate and reliable estimation of CCA wall motion from an ultrasound sequence. Typically, computer-aided methods for this task explicitly or implicitly involve a procedure called speckle tracking [8]. One of the first applications of speckle tracking to the The associate editor coordinating the review of this manuscript and approving it for publication was Halil Ersin Soken . estimation of CCA wall motion was based on block matching [9], [10]; however, the performance of block matching methods is negatively affected by a phenomenon known as speckle decorrelation [9]. Subsequent studies addressed this issue by using a state-space model for the evolution of the reference block [11], [12] or for the movement of the artery [7], [13].
As an alternative to block matching, several authors proposed speckle tracking methods based on optical flow [14]. For estimating CCA wall motion in particular, the Lucas-Kanade algorithm and its extensions were used [15], [16]. A comparative study of different speckle tracking methods in [17] demonstrated that a modified Lucas-Kanade algorithm outperforms other optical flow algorithms as well as block matching algorithms. In optical flow methods, speckle decorrelation translates into a phenomenon known as feature drift [18], which again negatively affects the performance. Indeed, the position of a selected point (feature) is tracked through cumulative summation of the estimated displacements between consecutive ultrasound frames, and hence the errors in the overall displacement estimates are accumulated as well. This may result-especially in the case of long ultrasound sequences-in a progressive divergence of the estimated point trajectory from the true trajectory, i.e., the estimated point ''drifts away'' from its true position. In [18], a Lucas-Kanade algorithm based method with explicit feature drift compensation was proposed; this method outperforms an earlier method described in [15]. All the mentioned studies analyze longitudinal scans of the CCA, with the exception of [15] and [18], where transverse scans of the CCA are analyzed using the Lucas-Kanade algorithm. In a transverse scan, the shape of the CCA wall cross section can be approximated by a circle, which simplifies the modeling.

B. CONTRIBUTION
In this work, we propose a new model-based method for estimating CCA wall motion. The CCA wall cross section is modeled as a circle, and our goal is to estimate the time-varying radius of that circle from-possibly long-B-mode ultrasound sequences showing the CCA transverse section. Based on a new state-space model for CCA wall motion, and using the result of an optical flow algorithm, we employ the unscented Kalman filter (UKF) [19] to track the time-varying center point and radius of the CCA wall circle.
As mentioned earlier, the performance of optical flow methods is impaired by the phenomenon of feature drift. In order to counteract feature drift, our implementation of the UKF fuses measurements from two different sources. One source is the pyramidal implementation of the Lucas-Kanade optical flow algorithm [20], which tracks the positions of a set of feature points (FPs) around the CCA wall. Here, an FP is a point in the ultrasound frame that is assumed to represent the local tissue movement, as manifested by a high response of the Harris detector [21]. The other source is the circle localization algorithm proposed in [22], which provides preliminary estimates of the center point and radius of the CCA circle. In this new fusion approach, feature drift in the tangential direction (with respect to the CCA circle) is reflected by a deviation from the state-space model, while feature drift in the radial direction is reflected by a deviation from the estimates provided by the circle localization algorithm. In contrast to the method of [18], the proposed method compensates for the feature drift in a temporally consistent manner that takes into account the natural smoothness of CCA wall motion. This inherent temporally consistent feature drift compensation is a major advantage of our method.
We assess the accuracy and robustness of the proposed method in several different synthetic scenarios with various levels of additive and multiplicative noise and various degrees of speckle decorrelation. We also demonstrate experimentally that our method outperforms the method proposed in [18] in all considered scenarios. We chose the latter method as a reference method since, to the best of our knowledge, it is the only existing method performing an explicit feature drift compensation, which is required for good performance in the case of long ultrasound sequences. The reference method was already benchmarked to a conventional artery wall tracking method in [18]. Finally, we present results obtained for a real ultrasound sequence and discuss their physiological validity.
The proposed method differs from the method previously presented in our conference publication [23] in that it estimates (tracks) the fundamental frequency-corresponding to the heart rate-used in the underlying state-space model.
Thus, in contrast to [23], the fundamental frequency does not have to be precisely known in advance and is allowed to vary with time. Furthermore, compared to [23], our validation of the performance of the proposed method is based on more realistic simulated sequences with a time-varying fundamental frequency, and we present a detailed description of the generation of these sequences. Finally, the technical presentation is more detailed.

C. PAPER ORGANIZATION
The remainder of this paper is organized as follows. An overview of the proposed method is presented in Section II. The calculation of the input to the UKF is described in Section III. The signal and state-space models on which the proposed method is based are developed in Section IV. The operation of the UKF is discussed in Section V. Finally, experimental results obtained for synthetic and real data are presented in Section VI.

II. METHOD OVERVIEW
For tracking a circle representing the CCA wall cross section, the proposed method fuses the results produced by the pyramidal Lucas-Kanade optical flow algorithm of [20] and by the circle localization algorithm of [22]. This fusion is performed by an UKF [19], which tracks the circle parameters and related quantities in a frame-sequential, frame-recursive manner. The input to our method is a B-mode ultrasound video sequence consisting of frames I t that show temporally successive transverse scans of the CCA. One recursion of the method, corresponding to frame (or time) index t, is visualized in Figure 1. The basic geometry of the CCA wall representation by a circle and of the surrounding FPs is depicted in Figure 2. We denote the center point of the circle by c t , the circle radius by r t , and the two-dimensional Cartesian position vector of the n-th FP by y (n) t . In Figure 3, VOLUME 8, 2020 we visualize each step of the proposed method (as presented in Figure 1).
For each frame I t with t = 1, 2, . . . , our method calculates a listŶ t of estimates of N t FP positions,ŷ (n) t for n = 1, 2, . . . , N t , as well as an estimate of the CCA circle center point,ĉ t , and of the CCA circle radius,r t . For the first frame, t = 1, the FP position listŶ 1 is empty and the CCA circle parameter estimatesĉ 1 andr 1 are calculated by applying the circle localization algorithm of [22] to I 1 .
For the further frames, t = 2, 3, . . . , the FP position list Y t and the circle parameter estimatesĉ t andr t are calculated frame-sequentially as follows (see Figures 1 and 3). In a first stage, as explained in Section III, we determine quantities that serve as the input to the UKF. First, the previous FP position listŶ t−1 is updated. To this end, we perform an FP inclusion step in which new FPs are identified in an annular search region in frame I t−1 that depends on the previous circle parameter estimatesĉ t−1 andr t−1 , and the positions of these new FPs are estimated and included in the list. The resulting augmented list will be denoted byŶ t−1 . Next, in the FP migration step, the pyramidal Lucas-Kanade optical flow algorithm is used to move the FP positions inŶ t−1 such that they conform to the current frame I t ; this results in an updated FP position listỸ t and a list of corresponding error metrics E t . Finally, in the FP exclusion step, FPs with a low confidence of successful localization by the optical flow algorithm, as determined from E t , are excluded from the FP position list. This yields a preliminary FP position listỸ t consisting of preliminary FP position estimatesỹ (n) t , n = 1, 2, . . . , N t . We note that the FP inclusion and exclusion steps are performed because of speckle decorrelation. Indeed, an FP may become progressively harder to localize (i.e., with less confidence) via the optical flow algorithm as the ultrasound image evolves. Such FPs need to be excluded from the tracking process and possibly replaced by new FPs.
In parallel to the FP inclusion, migration, and exclusion operations, preliminary circle parameter estimatesc t andr t are calculated by applying the circle localization algorithm of [22] to the current frame I t .
Then, in a second stage, the UKF fuses the information provided by the preliminary FP position listỸ t and by the preliminary circle parameter estimatesc t andr t . This information constitutes the input (i.e., the ''observed measurements'') of the UKF. The UKF then calculates final FP position estimateŝ y (n) t , n = 1, 2, . . . , N t , which make up the final FP position list Y t , as well as final circle parameter estimatesĉ t andr t . This is described in Section V. The UKF is based on a stochastic ''system model'' that describes the temporal variation of the CCA circle parameters, the evolution of the FP positions, and the extraction of the UKF input. This system model will be developed in Section IV.

III. CALCULATION OF THE UKF INPUT
The input to the UKF at frame time t consists of the preliminary FP position listỸ t comprising preliminary FP position estimatesỹ (n) t , n = 1, 2, . . . , N t , as well as of the preliminary circle parameter estimatesc t andr t (see Figure 1). These quantities are calculated as described next.

A. FP INCLUSION
First, the previous FP position listŶ t−1 is updated by performing the FP inclusion, migration, and exclusion steps. In the FP inclusion step, we use the Harris detector [21] to detect new FPs and estimate their positions within an annular search region in frame I t−1 . The annular search region is defined by the center pointĉ t−1 , the inner radiusr t−1 − r , and the outer radiusr t−1 + r , whereĉ t−1 andr t−1 are the previous circle parameter estimates and r > 0 is a fixed parameter. From the set of detected new FPs, we discard all FPs whose Harris response [21] is lower than a threshold t−1 , n = 1, 2, . . . , N t−1 , which will be denoted asŶ t−1 . By the above construction, the FP positionsŷ

B. FP MIGRATION
The augmented FP position listŶ t−1 was calculated from frame I t−1 , and thus it does not yet take into account the current frame I t . Therefore, in the next step of our methodthe FP migration step-we adaptŶ t−1 to I t . To this end,  Figure 1-for two successive ultrasound frames I t −1 (in (a)) and I t (in (b)-(e)). Bullets represent FPs. In particular, the green bullets in (a) depict FPs newly included by the FP inclusion step, while the red bullets in (c) depict FPs to be excluded by the FP exclusion step. The dotted circles in (a) represent the boundary of the annular search region used in the FP inclusion step. The dashed circle in (d) represents the preliminary CCA circle calculated by the circle localization algorithm of [22]. The solid circle in (e) represents the final CCA circle estimate produced by the UKF. The width of the annular search region in (a) and the length of the FP motion vectors in (b) are exaggerated for better visualization.
we use the pyramidal implementation of the Lucas-Kanade optical flow algorithm [20] to move each FP positionŷ t is calculated as the 1 metric between the two patches, i.e., where I t−1 (k, l) denote corresponding pixels of the two patches. The list of error metrics ε (n) t , n = 1, 2, . . . , N t−1 will be denoted as E t . The use of this metric is based on the assumption that the motion of the examined tissue in the neighborhood of an FP can be approximated sufficiently well by a simple translation. This assumption can be relaxed, as discussed in Section VII.

C. FP EXCLUSION
The FP migration step described above is followed by the FP exclusion step. An FP positionỹ (n) t is excluded from the FP position listỸ t if the corresponding error metric ε where ε t is the arithmetic mean of all the error metrics in E t and α > 1 is a fixed parameter. This rule is motivated by our empirical observation that there is usually a small proportion of FPs whose error is much higher than the average, and excluding these ''outliers'' improves the tracking accuracy. The FP positions surviving this exclusion step are denoted asỹ (n) t , n = 1, 2, . . . , N t , and the corresponding list asỸ t . This list constitutes one of the inputs to the UKF.

D. PRELIMINARY CIRCLE PARAMETER ESTIMATION
The second input to the UKF is obtained by applying the circle localization algorithm of [22] to the current frame I t . This results in ''preliminary'' circle parameter estimatesc t andr t . The circle localization algorithm of [22] differs from standard algorithms such as the Hough transform [24] in that it utilizes the pixel intensity directly rather than the image gradient. This is beneficial in our context because the blurred appearance and speckled texture of ultrasound images make accurate gradient determination challenging. The preliminary circle parameter estimatesc t andr t provide a complementary ''absolute'' position information that is used in the UKF to compensate for feature drift. Without this absolute position information-i.e., relying only on the relative optical flow information underlying the FP migration step-the FP position estimates and circle radius estimates produced by the UKF would drift away from the CCA wall (as explained earlier in Section I-A).

IV. SYSTEM MODEL
A major contribution of our work is a new state space model that constitutes a stochastic description of the dynamics of CCA wall motion (within the transverse cross section considered) and of a related measurement process. This model is inspired by the state-space model of [25], which uses a Fourier series with time-varying coefficients whose temporal evolution is modeled by a random walk. The proposed state space model consists of a state evolution model and a measurement model, and it provides the basis for the operation of the UKF. We note that although this model is designed for the CCA transverse section, the underlying approach can be used in many other speckle-tracking scenarios where a quasi-periodic motion is observed.

A. CCA CIRCLE RADIUS AND FP RADII
Our state space model is based on models for the time-dependence of the CCA circle radius and of the FP radii. VOLUME 8, 2020 The regularity of the heart beat manifests itself in an approximate periodicity of the movement of the CCA wall. This approximate periodicity suggests the use of a Fourier series to model the CCA wall movement [26,Ch. 11]. Concretely, we model the time-dependence of the CCA circle radius r t by the sum of a constant R that describes the mean CCA circle radius, a quasi-periodic function that is constructed as the superposition of M harmonic components with timevarying Fourier coefficients a (m) t and b (m) t , m = 1, 2, . . . , M and time-varying fundamental frequency, and a slowly varying component γ t that represents the effect of breathing [18], [27]. That is, where mϕ t is the instantaneous phase of the m-th component and ϕ t can be considered as the ''fundamental'' instantaneous phase (whose temporal derivative-if t were continuouswould be a time-varying fundamental frequency).
We consider N t FPs at each frame time t. The position of the n-th FP at frame time t, relative to the CCA circle center point c t , can be expressed as Here, ρ (n) t (referred to as the ''FP radius'') and θ (n) t represent the polar coordinates of FP n with respect to the center c t , as illustrated in Figure 2. The FP radius is modeled as where r t is the CCA circle radius,δ (n) denotes an initial approximation of the radial deviation of the FP from the CCA circle (which was calculated at the frame time when the FP was included in the FP position list, see Section V-A), and δ (n) t denotes a time-varying radial deviation.

B. STATE
The UKF estimates a time-varying state vector x t in a framerecursive manner, using at each frame time t the results of the previous frame time t −1 and a measurement vector z t . In our method, the state is a vector of dimension L t 2M + 2N t + 7 defined as Here, related to the CCA circle radius model in (2), a t (a (1) t ) T are the vectors of time-varying Fourier coefficients of the quasi-periodic component of the CCA circle radius, ϕ t (ϕ tφt ) T comprises the fundamental instantaneous phase and its derivative, i.e., the time-varying fundamental frequency, γ t (γ t γ t−1 ) T comprises the current and previous samples of the breathing component of the CCA circle radius, and R t is a (formally) time-varying version of the mean CCA circle radius R. Furthermore, related to the FP model in (3) and (4), t ) T is the CCA circle center point in Cartesian coordinates, and the FP position parameter vector characterizes the position of FP n in polar coordinates. Note that from the state x t , it is possible to determine the CCA circle center point c t and radius r t as well as all the FP positions (up to the initial radial deviation approximationsδ (n) ).

C. STATE EVOLUTION MODEL
The UKF relies on a state evolution model, which is a stochastic model for the one-frame evolution of the state x t , i.e., the transition from x t−1 to x t . For a basic formulation of the state evolution model, we temporarily assume that N t = N and L t = 2M + 2N + 7 = L are constant; an extension to time-varying N t and L t will be given in Section V-A. We choose a linear state evolution model given by Here, is an L × L block-diagonal matrix defined as where ϕ 1 1 and φ γ , φ δ ∈ (0, 1). Furthermore, the driving process u t is an independent and identically distributed (iid) L-dimensional zero-mean process with covariance matrix where ϕ σ 2 This model summarizes the following individual state evolution models: Here (·) i denotes the i-th element of the vector in parentheses and (·) i,j denotes the vector comprising the i-th and j-th elements of the vector in parentheses. Note that the state evolution models for a  1, 2, . . . , N ) are random walk models, ϕ t is modeled by a constant velocity model, the state evolution model for γ t is a second-order auto-regressive (AR(2)) model, R t is modeled as constant (because the variance of (u t ) 2M +5 is zero), and δ (n) t is modeled by a first-order auto-regressive (AR(1)) model. We use AR models, rather than simple random walk models, for γ t and δ (n) t because the AR model effectively limits the range of feasible values. This helps prevent the FPs from drifting away from the CCA wall, and is one reason why our method is able to compensate the feature drift resulting from the optical flow algorithm. (The other reason is the use of the results of the circle localization algorithm of [22] as an additional input to the UKF, as explained in Section III-D.) For a larger variance of (u t ) l , the respective state component (x t ) l tends to change more rapidly and less smoothly.

D. MEASUREMENT MODEL
The UKF also relies on a measurement model in addition to the state evolution model. The measurement model is a stochastic model for the dependence of the observed measurements on the state x t . We recall from Sections II and III that our measurements are given by the preliminary FP position estimatesỹ Here, ρ (n) t is given by Equation (4), i.e., ρ t , in which, in turn, r t is given by Equation (2) with R formally replaced by R t . Thus, expression (10) involves the state components (cf. t ) T for n = 1, 2, . . . , N ; note that ϕ t and θ (n) t enter in a nonlinear manner. Furthermore, the measurement noise processes v (n) t are mutually independent, iid, zero-mean, two-dimensional processes with covariance matrix C v (n) diag{σ 2 v,1 , σ 2 v,2 }. In a similar manner, the preliminary CCA circle parameter estimatesc t andr t are modeled as noisy versions of the true CCA circle parameters c t and r t , respectively, i.e., where v (c) t and v (r) t are mutually independent, iid, zero-mean processes with covariance C v (c) σ 2 v (c) I 2 and variance σ 2 v (r) , respectively. The expressions (11) involve the state components (cf. (5)) a t , b t , ϕ t , γ t , R t , and c t . Finally, the overall measurement vector z t comprises all the measurements, i.e., We can now summarize the measurement models (10) and (11) into an overall measurement model where and v t v of v t is given by

V. UKF OPERATION
The UKF [19] is a suboptimal sequential Bayes filter that calculates at each frame time t = 2, 3, . . . estimates of the posterior mean E{x t |z 1:t } and posterior covariance Cov{x t |z 1:t } of the state x t , where z 1:t (z T 1 · · · z T t ) T . These estimates will be denoted byx t and P t , respectively; note thatx t ≈ E{x t |z 1:t } provides an estimate of x t . We use a UKF, rather than simply a Kalman filter [28,Ch. 2], because of the nonlinearity of our measurement model (12). Also, in comparison to the extended Kalman filter [28,Ch. 2], the UKF usually provides a better estimation accuracy while its computational complexity is of the same order [19]. We use a slightly modified version of the UKF algorithm that includes an FP allocation/deallocation step to account for the time-varying number of state variables. Each UKF recursion thus consists of the following steps, which are described further below: FP allocation/deallocation, state prediction, calculation of sigma-points, measurement prediction, update, and estimation. These steps are performed at each frame time t = 2, 3, . . . .

A. FP ALLOCATION/DEALLOCATION
According to (5), the inclusion of a new FP n in the previous FP position listŶ t−1 -as described in Section III-A-entails a corresponding allocation (insertion) of the FP position parameter vector p (n) t−1 in the state vector x t−1 . This implies an analogous allocation of an estimate of p (n) t in the state estimate vectorx t for t = t − 1, t, . . . We initialize this estimate as (cf. (6))p t−1 is the initial FP position obtained in the FP inclusion step (see Section III-A), andĉ t−1 is the previous estimate of the circle center. Furthermore, corresponding variances are allocated as two additional elements on the diagonal of the approximate covariance matrix P t for t = t − 1, t, . . .. These variances are initialized in P t−1 as s 2 δ and s 2 θ , which are parameters expressing our a priori VOLUME 8, 2020 uncertainty about δ (n) t−1 and θ (n) t−1 , respectively. Finally, we calculate an initial approximation of the radial deviation of the new FP from the previous estimate of the CCA circle asδ (n) = ξ t−1 2 −r t−1 , where · 2 denotes the 2 (Euclidean) norm; we recall thatδ (n) is used in the FP radius model (4) for t = t −1, t, . . ..
Similar allocation operations are applied to the transition matrix in (8) and the covariance matrix C u in (9), which thereby become time-varying. More specifically, the values φ δ and 1 are inserted as additional diagonal elements of t , and the values σ 2 δ and σ 2 θ are inserted as additional diagonal elements of C u,t , in both cases for t = t −1, t, . . ..
Conversely, the exclusion of FP n at time t-as described in Section III-C-entails a deallocation (removal) of the corresponding FP position parameter vector p (n) t−1 in the state vector x t−1 . This implies analogous deallocation operations inx t , P t , t , and C u,t for t = t −1, t, . . ..

B. STATE PREDICTION
In the state prediction step of the UKF recursion at frame time t, a ''predicted'' state meanx t|t−1 and covariance P t|t−1 are calculated by propagating the previous meanx t−1 and covariance P t−1 through the linear state evolution model (7). This results in

C. CALCULATION OF SIGMA-POINTS
Next, so-called sigma-points are used to approximate the propagation of the predicted meanx t|t−1 and covariance P t|t−1 through the nonlinear measurement model (12). According to [19], for state vector dimension L t , there are 2L t + 1 sigma-pointsx (l) t and corresponding weights w (l) t , l = 0, 1, . . . , 2L t , which are calculated aŝ Here, ( 1/2 t ) l denotes the l-th column of 1/2 t , which is a square root of the matrix t L t 1−w (0) t P t|t−1 , i.e., any square matrix satisfying t can be chosen differently from 1/3, as discussed in [19].

D. MEASUREMENT PREDICTION
The sigma-points are now propagated through the (noiseless) nonlinear measurement model (12), i.e., Then, a ''predicted'' measurementẑ t|t−1 and an innovation covariance matrix S t|t−1 are calculated aŝ

E. UPDATE
Next, we perform the update step of the ordinary Kalman filter [28,Ch. 2], i.e., Here G t is a sigma-point approximation of the Kalman gain matrix, which is given by

F. ESTIMATION
Evaluating g(·) in (13) for the state estimatex t yields as the final estimates of the FP positions y (n) t , for n = 1, 2, . . . , N t ; these estimates constitute the final FP position estimate listŶ t . Similarly, we useĉ t andr t as the final estimates of the CCA circle parameters c t and r t . Note that the FP position estimate listŶ t and the center point estimateĉ t are used only internally, namely, as an input to the FP inclusion step at the next frame time t+1, as explained in Sections II and III-A. On the other hand,r t forms the output of the overall method.

G. INITIALIZATION
At frame time t = 1, no FPs are allocated yet, i.e., N 1 = 0. We initialize the meanx 1 and covariance matrix P 1 as (cf. (5)) Here,φ 1 is an initial value 1 of the time-varying fundamental frequencyφ t , andr 1 andĉ 1 are initial CCA circle parameter estimates that are calculated by the circle localization algorithm of [22]. The variances s 2 ab , s 2 ϕ , s 2 γ , s 2 R , and s 2 c are parameters expressing our a priori uncertainty about the Fourier coefficients a  1, 2, . . . , M ) as well as about ϕ 1 , γ 1 , R 1 , and c 1 , respectively. 1 The fundamental instantaneous phase ϕ t and the time-varying fundamental frequencyφ t are part of the state x t -see (5)-and thus are estimated (tracked) along with the other state components. Therefore, in contrast to [23], the fundamental frequency does not have to be precisely known in advance. However, a reasonable initialization ofφ t at t = 1 is required for the method to converge.

VI. NUMERICAL RESULTS AND DISCUSSION
In this section, we present numerical results demonstrating the performance of the proposed method in comparison to an improved version of the state-of-the-art reference method described in [18]. The latter method was chosen as a reference method because, similarly to the proposed method, it performs an explicit feature drift compensation and is thus suited to the case of long ultrasound sequences.

A. BASIC SIMULATION SETUP
We considered the following types of experimental data. Data set S was created synthetically through a time-varying transformation of a real ultrasound image, and was used to evaluate the robustness of the proposed method to additive and multiplicative noise. Data set F was created synthetically by means of the Field II simulation program [29], [30], and was used to determine the method's performance in the presence of speckle decorrelation. Finally, we validated our method on a real ultrasound sequence. More detailed descriptions of these experimental data will be provided in later subsections. The parameters of the proposed method that we used in our experiments are listed in Tables 1 and 2. In the course of our experiments, we observed that the circle localization algorithm used in the reference method [18] performs poorly on our data. Therefore, for a fair comparison of tracking accuracy, we modified the reference method such that it uses the circle localization algorithm used in our method (see Section III-D). Furthermore, in both the proposed method and the reference method, we included a correction of the radius estimation bias exhibited by the circle localization algorithm. More specifically, the preliminary circle radiusr t , used subsequently by the UKF, was obtained by subtracting a constant from the radius estimate computed by the circle localization algorithm. For data sets S and F, we estimated this constant from sequence S ∞ and F 0 , respectively, by averaging over all frames the error of the radius estimates produced by the circle localization method. Note that we were able to determine this error because for data sets S and F the ground truth is available.
The proposed method was implemented in MATLAB on an Intel(R) Core(TM) i5-7500 CPU with a base frequency of 3.40 GHz, without the use of multi-threading. The processing of one frame always took less than 40 ms. In particular, the processing of a single frame of 348 × 280 px from dataset S (discussed later) took on average 26 ms.

B. MOTION MODEL
To generate the simulated sequences, we modeled the CCA wall transverse section by a circle with a time-varying radius r t that is the sum of a pulse wave component and a breathing component. For the pulse wave component, we employed the parametric model proposed in [31], whose parameters we chose such that they conformed to radius estimates obtained from real ultrasound sequences. In order to introduce smooth variations in the pulse rate, we frequency-modulated the pulse wave by a sine function with a period of half a minute. During one period of the modulating sine function, the heart rate varies between 70 and 110 beats per minute. The breathing component was modeled as the absolute value of a sine function whose frequency is much lower than the assumed heart rate. Further details of the motion model and its parameters are provided in [18].

C. PERFORMANCE METRICS
For measuring the accuracy of radius estimation in those cases where a ground truth r t is available (data sets S and F), we use the root mean square error (RMSE). Assuming K simulation runs, each producing an estimated radius sequencê r t,k (t = 1, 2, . . . , T ; k = 1, 2, . . . , K ), the RMSE is defined as the square root of the average of the squared error (r t,k − r t ) 2 taken over all frames t and all simulation VOLUME 8, 2020 runs k, i.e., We also consider the inter-realization empirical standard deviation of RMSE k , which is denoted SD and defined as

D. SYNTHETIC DATA SET S
Data set S consists of a synthetic ultrasound sequence, referred to as S ∞ , as well as K = 100 realizations of four ''noisy'' synthetic ultrasound sequences, referred to as S + 25 , S + 15 , S * 25 , and S * 15 . Each sequence comprises T = 960 ultrasound frames. With a frame rate of 32 fps, this corresponds to 30 seconds of continuous recording. Sequence S ∞ was derived from a reference image, shown in Figure 4(a), through a time-varying deformation. The reference image was chosen as a real ultrasound image of the CCA  transverse section of a healthy subject, which was acquired as described in Section VI-F. The frames of S ∞ were then obtained by deforming the reference image in accordance with the motion model described in Section VI-B. More specifically, the CCA circle radius for each frame was calculated from the motion model, a displacement field (describing the tissue displacement at each pixel of the reference image) was calculated from the CCA circle radius, and the reference  image was transformed according to the displacement field. A more detailed description of this procedure can be found in [18]. Finally, the noisy sequences of data set S were obtained by corrupting S ∞ by additive Gaussian noise (in the case of S + 25 and S + 15 ) or by multiplicative Gaussian noise (in the case of S * 25 and S * 15 ). The subscript, 25 or 15, indicates the signal-to-noise ratio in decibels. Similarly to [17], we use noisy sequences to evaluate the robustness of the methods within a simple, controlled, and reproducible setting. In real B-mode (log-compressed) images, the noise characteristics  would be more complex than the additive and multiplicative Gaussian characteristics employed here, although additive noise would still be the most prominent component [32]. We created K = 100 realizations of each of the four sequence types S + 25 , S + 15 , S * 25 , and S * 15 by using 100 realizations of the respective noise process (i.e., the same ultrasound sequence S ∞ was corrupted by different noise realizations). Note that for S ∞ , formally, K = 1. For each of the noisy sequences S + 25 , S + 15 , S * 25 , and S * 15 , the first frame of one realization is shown in Figures 4(b)-(e), respectively.
For each sequence or sequence type, we calculated RMSE and SD for the radius estimates that we obtained with the proposed method and with the reference method [18]. Except for S ∞ , this calculation was based on the K = 100 realizations of the respective noise process. The results are presented in Table 3, along with the average RMSE and SD taken over all sequences. One can see that the proposed method always outperforms the reference method in terms of RMSE. In addition, the SD results show that the estimates produced by the proposed method are much more consistent across multiple realizations. In Figure 5, we show the individual estimated radius sequencesr t,k obtained with either method for ten realizations of S + 15 . The figure also shows the time-varying RMSEs, i.e., RMSE t 1 K K k=1 (r t,k − r t ) 2 . It can be observed that the estimates obtained with the proposed method are almost always closer to the ground truth. The large variability of the results of the reference method (across the individual realizations) is primarily caused by feature drift [18]. Figure 6 shows analogous results obtained for S * 15 . One can observe that the proposed method exhibits an even larger performance gain relative to the reference method. This is due to an even stronger feature drift, which is effectively compensated by the proposed method but less well by the reference method. We note that the increase in the RMSE of the proposed method that is observed in Figure 6(c) between 12 s and 14 s is caused by a temporary circle localization error; the RMSE starts decreasing again around 22 s, which is however not shown in Figure 6(c).

E. SYNTHETIC DATA SET F
Data set F consists of three synthetic ultrasound sequences referred to as F 0 , F 2 , and F 5 . Each sequence again comprises T = 960 ultrasound frames or 30 seconds of continuous recording (assuming a frame rate of 32 fps).
Sequence F 0 was created by means of the Field II simulation program [29], [30]. First, a three-dimensional tissue model (reference phantom) was generated. This reference phantom consists of infinitesimally small objects called scatterers, whose positions are uniformly distributed in the phantom with a density of about 50 mm −3 . With the resolution cell volume given by 0.2 mm 3 , about ten scatterers occupy one resolution cell. Each scatterer has an amplitude that is randomly drawn from a zero-mean Gaussian distribution. The standard deviation of that distribution was chosen equal to the intensity of the real ultrasound image that was used to create data set S (see Figure 4(a)), at the position corresponding to the scatterer's lateral and axial coordinates relative to a simulated ultrasound probe. The simulated ultrasound probe consisted of 192 elements, 72 of which were active during each single-line scan. The center frequency of the probe was set to 10 MHz. In each frame, 64 lines were scanned. With this set of parameters, we achieved a lateral resolution of 1.6 mm and an axial resolution of 0.32 mm.
An individual phantom was then created for each frame by displacing the scatterers of the reference phantom according to the motion model described in Section VI-B. The resulting sequence of phantoms was passed to the Field II simulation program, which produced the ultrasound sequence F 0 . The first frame of F 0 is shown in Figure 7.
To simulate speckle decorrelation, we then created sequences F 2 and F 5 by replacing in each frame of F 0 2% and 5% of the scatterers, respectively, with new scatterers. 2 The new scatterers were again uniformly distributed in the phantom of the respective frame, and their amplitudes were chosen as explained above. A more detailed description of this generation process is provided in [18]. Because of the high runtime of the Field II simulations, we generated only 2 Speckle decorrelation is more commonly associated with a nonuniform motion of the scatterers [9]. However, our approach of randomly replacing a proportion of the scatterers with new ones can be motivated by the fact that the complex ultrasound image (before envelope detection) can be modeled as the sum of complex phasors whose amplitudes and phases are determined by the scatterer's reflectivity and distance to the transducer, respectively [32]. Due to phase wrapping, it is irrelevant whether a scatterer travels, e.g., half or one and a half of the ultrasound wavelength. Thus, the effect of a nonuniform scatterer motion is similar to that of randomly drawing new scatterers. a single realization for each of the sequences F 0 , F 2 , and F 5 , i.e., the number of realizations for data set F is K = 1. The frames of sequences F 2 and F 5 are similar visually to those of F 0 (see Figure 7).
In Table 4, we list the RMSE obtained with the proposed method and the reference method for data set F. As in the case of data set S, the proposed method is seen to consistently outperform the reference method. Furthermore, Figure 8 shows the estimated radius sequences and the absolute estimation errors obtained for sequence F 2 . Again, the proposed method is seen to estimate the radius much more accurately than the reference method. Finally, still for sequence F 2 , Figure 9 shows the CCA circle estimates obtained with the proposed method and the reference method at four different frames times. It is seen that, while the estimates of the two methods are generally similar, the estimates of the reference method are slightly biased towards lower values of the radius. This confirms the results shown in Figure 8(a).

F. REAL ULTRASOUND SEQUENCE
To complement the results obtained for the synthetic data sets S and F, we consider a real ultrasound video sequence that shows the CCA transverse section of a healthy test subject (sex: male, age: 27, weight: 64 kg). This sequence was measured using an Ultrasonix OP device with linear probe L145/38 (Ultrasonix Medical, Richmond, BC, Canada). The sequence length is 787 frames, and the frame rate is 112 fps. The sequence was obtained with the subject's informed consent and with approval from the ethics committee. Figure 10 shows the estimated radius waveform, while Figure 11 shows the CCA circle estimates at four different frame times. From Figure 10, one can conclude that the estimated radius waveform obtained with either method exhibits the typical characteristics described for example in [26], where the waveform during a cardiac cycle consists of two consecutive peaks (local maxima), the first associated with the forward wave and the second associated with the reflected wave. One can furthermore see in Figure 10 that after an initial convergence period-which spans approximately two cardiac cycles-the waveform produced by the proposed method is smoother than that produced by the reference method. Both Figure 10 and Figure 11 show that the proposed method generally yields larger radius estimates than the reference method; however, no conclusion regarding estimation accuracy can be drawn from this observation as no ground truth is available. From the radius waveform obtained with the proposed method, we can derive the following approximate values of diagnostically relevant parameters: the heart rate is obtained as 80 beats per minute, the diastolic artery diameter D d as 6.00 mm, the systolic artery diameter D s as 6.48 mm, the pulse diameter D s − D d as 0.48 mm, and the circumferential strain (D s − D d )/D d as 0.08. All these values are in the typical range for a healthy male subject [26,Ch. 9]. If measurements of the pulse pressure were available, it would also be possible to calculate various arterial stiffness indices such as arterial distensibility and compliance [26,Ch. 9].

VII. CONCLUSION
The method we proposed in this paper is able to continuously and accurately track a circular approximation of the common carotid artery (CCA) wall, based on an observed B-mode ultrasound sequence of arbitrary length. We used an unscented Kalman filter (UKF) to track a composite state characterizing CCA wall motion. The operation of the UKF relies on a new state-space model that describes the dynamics of CCA wall motion and a related measurement process in a stochastic manner. The UKF fuses two different sets of measurements, of which one is produced by the pyramidal Lucas-Kanade optical flow algorithm and the other by a recently proposed CCA wall localization method. This fusion enables an effective, temporally consistent compensation of feature drift, which normally impairs the performance of optical flow based methods. The temporally consistent compensation of feature drift is a major advantage of our method over state-of-the-art methods.
We performed a quantitative evaluation of the accuracy of the proposed method for synthetic data sets containing additive and multiplicative noise and emulating speckle decorrelation. Our results demonstrate significant performance advantages relative to the state-of-the-art method presented in [18]. We also validated our method for a real ultrasound sequence.
The CCA circle radius waveform estimated by the proposed method can be used to calculate various indices that provide information about the health of the cardiovascular system [6], [7]. Moreover, because by its effective feature drift compensation the proposed method can be applied to arbitrarily long ultrasound sequences, it is possible to average clinically relevant indices over many cardiac cycles, thereby improving accuracy and reliability. In addition, the Fourier series coefficients calculated by our method can be used for a frequency analysis of the CCA radius waveform sequence, which is potentially useful for the diagnosis of cardiovascular diseases [33]. Developing an extension of our method that allows the estimation of other clinically important parameters of the CCA, such as radial and longitudinal strain [26], is an interesting direction for future research.
The proposed method can be extended in various other ways. For example, some or all of the parameters of the method can be included in the state and estimated online. The CCA wall can be characterized by a more sophisticated family of contours (instead of circles). Finally, an improved optical flow algorithm can be employed. Indeed, the pyramidal Lucas-Kanade algorithm used in our method is based on the assumption that tissue motion in the neighborhood of an FP can be approximated sufficiently well by a simple translation. The same assumption underlies our definition of the error metric ε (n) t in (1). To better account for nonuniform tissue motion, the translation assumption can be relaxed by using an affine model for the velocity field [34] or a global optical flow algorithm such as [35]. The image patch I (ỹ t ) t in (1) can then be replaced by an appropriately warped and translated version of the original image patch.