Robust Myocardial Motion Tracking for Echocardiography: Variational Framework Integrating Local-to-Global Deformation

Abstract

This paper proposes a robust real-time myocardial border tracking algorithm for echocardiography. Commonly, after an initial contour of LV border is traced at one or two frames from the entire cardiac cycle, LV contour tracking is performed over the remaining frames. Among a variety of tracking techniques, optical flow method is the most widely used for motion estimation of moving objects. However, when echocardiography data is heavily corrupted in some local regions, the errors bring the tracking point out of the endocardial border, resulting in distorted LV contours. This shape distortion often occurs in practice since the data acquisition is affected by ultrasound artifacts, dropouts, or shadowing phenomena of cardiac walls. The proposed method is designed to deal with this shape distortion problem by integrating local optical flow motion and global deformation into a variational framework. The proposed descent method controls the individual tracking points to follow the local motions of a specific speckle pattern, while their overall motions are confined to the global motion constraint being approximately an affine transform of the initial tracking points. Many real experiments show that the proposed method achieves better overall performance than conventional methods.

1. Introduction

In company with the development of real-time three-dimensional echocardiography (RT3DE), the demands for automated analysis methods of left ventricle (LV) assessment such as ejection fraction, motion analysis, and strain analysis are rapidly increasing. Nevertheless, most of the analysis methods are still based on the measurements in a few two-dimensional (2D) slices, because they are available in clinical practice [1, 2]. In general, the quantitative assessment for heart function is performed by manually tracing endocardial border in some 2D slices of different view at frames (such as end-systole (ES) or end-diastole (ED) frames) selected from the entire cardiac cycle and automatically tracking the traced LV contour over the remaining frames [3, 4]. The motion tracking of LV is carried out by observing the speckle pattern associated with deforming tissue. Speckle pattern is an inherent appearance in ultrasound imaging and its local brightness reflects the local echogeneity of the underlying scatterers. Since it is a difficult task to automatically track the motion of endocardial border in ultrasound images due to ultrasound artifacts, dropouts or shadowing phenomena, low contrast, and so on, user intervention is somewhat required for stable and successful tracking of endocardial border.

In the last decades, there have been numerous studies for tracking of LV wall motion such as the tracking methods using deformable models [5–8], active shape models [9–11], and optical flow methods [2, 12–15]. Those methods have some limitations to practical application of endocardial border motion tracking. In deformable models, their methods are relatively time consuming due to iterative contour evolution with stopping criteria and often need preprocessing for speckle reduction before wall motion tracking. Active shape models are the statistical methods based on the dataset of trained images so that they require additional effort to train on many images. Both deformable models and active shape models provide the motion information of LV border and enable user to measure the volume inside LV, whereas they are somewhat inadequate for strain analysis related to the motion and deformation of heart, because they are not speckle tracking-based methods providing motion information of local region on the myocardium but shape-based tracking methods.

On the other hand, optical flow methods, which use the assumption that the intensity of a moving object is constant over time, provide the local motion information of myocardium. They are capable of measuring the LV volume as well as the myocardial wall motion analysis or strain analysis to detect LV abnormalities. After an initial contour of endocardial border is traced, each point on the contour tracks the specific intensity and speckle pattern in sequential images. However, it is problematic to track the endocardial border in ultrasound images with unclear speckle pattern or weak signals. In practical environment, there often exist some incorrectly tracked points due to ultrasound artifacts, dropouts, or shadowing phenomena of cardiac wall [16]. When edge dropout or indistinguishable speckle pattern is present in a local neighborhood of a tracking point, the errors bring the tracking point out of the endocardial border, resulting in distorted LV contours throughout the entire cardiac cycle as shown in Figure 1 or irregular distances between the tracked points in Figure 2. These distorted results affect LV volume measurement or strain analysis.

(a)

(b)

(c)

(d)

Figure 1 

The estimation of endocardial border by the Lucas-Kanade optical flow method. Case  1: a tracking point getting out from the real LV shape distorts the whole shape near the border with weak edges; (a) initially traced endocardial border and its tracking points at an ED frame, (b) the tracked result at the ES frame, (c) at the frame between ES and ED, and (d) at the next ED frame.

(a)

(b)

(c)

(d)

Figure 2 

The estimation of endocardial border by the Lucas-Kanade optical flow method. Case  2: the tracked points are irregularly spaced by indistinguishable speckle patterns; (a) initially traced endocardial border and its tracking points at an ED frame, (b) the tracked result at the ES frame, (c) at the frame between ES and ED, and (d) at the next ED frame.

In order to cope with these problems, we develop a new optical flow method equipped with a global motion constraint that is designed to prevent each tracking point from getting out of the endocardial border. In the proposed model, the Lucas-Kanade (LK) optical flow method [17] and a global motion constraint being approximately an affine transformation of the initial tracking points are incorporated into a variational framework. So the individual tracking points follow speckle patterns (corresponding to each tracking point) and their overall motions are confined to the global motion constraint. The global motion constraint is based on the results [18, 19] that heart motion is regarded as the nonrigid motion by rotation, contraction/expansion, and shear. Typically, nonrigid motion consists of global deformation and local deformation. The global deformation is modeled by an affine transformation while the local deformation is described by a free-form deformation.

The proposed algorithm is capable of tracking LV border in real-time since its movement is directly computed from the difference between two sequential images via a simple matrix multiplication. For performance evaluation, we carry out various real experiments with Samsung Medison R&D Center (http://www.samsungmedison.com/). Numerous experiments show better performance of the proposed tracking methods compared to the conventional tracking methods.

2. Methods

2.1. Conventional Optical Flow Tracking Methods

Let 𝐼(𝐫,𝑡) represent the intensity of echocardiography at the location 𝐫=(𝑥,𝑦) and the time 𝑡. Optical flow tracking methods are based on the assumption that the intensity of a moving object is constant over time, so that the noisy time-varying images 𝐼(𝐫,𝑡) approximately satisfy 𝐮()()+𝜕𝐫,𝑡⋅∇𝐼𝐫,𝑡𝐼()𝜕𝑡𝐫,𝑡≈0, (1) where 𝐮(𝐫,𝑡) is the velocity vector to be estimated. Based on (1), numerous approaches for estimating the velocity vector 𝐮(𝐫,𝑡) have been proposed and those were applied to LV border tracking in echocardiography [2, 13–15].

Horn and Schunk [20] proposed the optical flow technique incorporating the smoothness of the motion vector in the entire image as a global constraint. In their model, the velocity 𝐮(𝐫,𝑡) at each time 𝑡 is determined by minimizing the energy functional: 𝐸_𝑡(𝐮)∶=_Ω𝐮(𝐫)()+𝜕⋅∇𝐼𝐫,𝑡𝐼()𝜕𝑡𝐫,𝑡²+𝜆|∇𝐮(𝐫)|²𝑑𝐫, (2) where Ω is the image domain and 𝜆 a regularization parameter which controls the balance between the optical flow term and the smoothness on 𝐮. The velocity 𝐮(𝐫,𝑡) at each time 𝑡 can be computed by solving the corresponding Euler-Lagrange equation that is a reaction-diffusion equation. In [21], it has been observed that this global method with the global smoothness constraint is significantly more sensitive to noise than the local method used by Lucas and Kanade [17].

Lucas and Kanade [17] used the assumption of locally constant motion to compute the velocity 𝐮(𝐫₀,𝑡) at a target location 𝐫₀=(𝑥₀,𝑦₀) and time 𝑡 by forcing constant velocity in a local neighborhood of a point 𝐫₀=(𝑥₀,𝑦₀), denoted by 𝒩(𝐫₀). Following Lucas and Kanade, Barron et al. [21] estimated the velocity 𝐮(𝐫₀,𝑡) by minimizing the weighted least square criterion in the neighborhood 𝒩(𝐫₀): 𝐮𝐫₀,𝑡∶=argmin_{𝐮𝒩(𝐫0)}𝑤𝐫−𝐫₀()+𝜕𝐮⋅∇𝐼𝐫,𝑡𝐼(𝜕𝑡𝐫,𝑡)²𝑑𝐫, (3) where 𝑤 is a weight function that enables to give more relevance to central terms rather than the ones in the periphery. Here, “arg min” stands for the argument of the minimum, that is, the vector 𝐮 for which the right integral attains its minimum value. Since this method determines 𝐮(𝐫₀,𝑡) at each location 𝐫₀ by combining information from all pixels in the neighborhood of 𝐫₀, it is reasonably robust against image noise. We used (3) as the Lucas-Kanade method, because this weighted window LK method is essentially same as the LK method. When the weight function 𝑤 is uniform, the form is the same as the Lucas and Kanade one, in fact.

As we mentioned in Section 1, there often exist some incorrectly tracked points due to weak signal on cardiac wall since echocardiography data is acquired through transmitting and receiving ultrasound signals between the ribs, causing considerable shadowing of cardiac wall [16]. Due to these incorrectly tracked points, LK method may produce significantly distorted LV shape.

Recently, Sühling et al. [13] improved the weighted window LK method (3) by introducing a linear model for the velocity along the time direction, and the displacement 𝐮(𝐫₀,𝑡) is obtained by evaluating 𝐮 such that 𝐮,𝐛∈ℝ² and 2×2 matrix 𝐴 minimize the following energy functional: 𝐸_𝑡(𝐮,𝐴,𝐛)∶=^∞_−∞ℝ2×𝑤𝐫−𝐫₀,𝑠𝐮+𝐴𝐫−𝐫₀()+𝜕+𝑠𝐛⋅∇𝐼𝐫,𝑡+𝑠𝐼()𝜕𝑡𝐫,𝑡+𝑠²𝑑𝐫𝑑𝑠, (4) where 𝑤 is the symmetric window function, which gives more weight to constraints at the center of the local spatiotemporal region than to those at the periphery. Since this method uses multiple frames centering around the time 𝑡, it is more robust than the LK method (3) using the single frame at 𝑡. However, the same problem of LV shape distortion as in LK method still remains.

Compared with the approaches based on the LK method, Duan et al. [15] used the region-based tracking method (also known as the block matching or pattern matching method) with the cross-correlation coefficients as a similarity measure. For given two consecutive images 𝐼 at time 𝑡 and 𝑡+Δ𝑡, the velocity vector 𝐮=(𝑢,𝑣) for each pixel 𝐫=(𝑥,𝑦)∈Ω is estimated by maximizing the cross-correlation coefficients: 𝐮𝐫₀,𝑡∶=argmax_𝐮×∫_𝒩(𝐫0)[𝐼(𝐫,𝑡)𝐼(𝐫+𝐮,𝑡+Δ𝑡)]𝑑𝐫∫_𝒩(𝐫0)[𝐼(𝐫,𝑡)]²𝑑𝐫∫_𝒩(𝐫0)[𝐼(𝐫+𝐮,𝑡+Δ𝑡)]².𝑑𝐫 (5)

Instead of maximizing the cross-correlation coefficients, the velocity vector can be estimated by minimizing the sum-of-squared difference (SSD) [21] as follows: 𝐮𝐫₀,𝑡∶=argmin_{𝐮𝒩(𝐫0)}×𝑤𝐫−𝐫₀[𝐼(𝐫,𝑡)−𝐼(𝐫+𝐮,𝑡+Δ𝑡)]²𝑑𝐫. (6)

The block matching method uses similarity measures that are less sensitive to noise, of fast motion, and of potential occlusions and discontinuities [15].

The above three local methods have drawback in dealing with the problem of the contour shape distortion in the presence of locally weak signal corrupted by rib shadowing and other factors. Hence, we need to develop a method alleviating shape distortion.

2.2. Proposed Method

The proposed method uses an affine transformation to describe a global motion that is synthesized by integrating local deformations. We denote the endocardial border traced at initially selected frame (e.g., end-systole or end-diastole frame) by a parametric contour 𝒞^∗={𝐫^∗(𝑠)=(𝑥^∗(𝑠),𝑦^∗(𝑠))∣0≤𝑠≤1} that can be identified as its 𝑛 tracking points 𝐫^∗₁=𝐫^∗(𝑠₁),…,𝐫^∗_𝑛=𝐫^∗(𝑠_𝑛). Here, 0=𝑠₁<𝑠₂<⋯<𝑠_𝑛=1. Let 𝒞(𝑡)={𝐫(𝑠,𝑡)=(𝑥(𝑠,𝑡),𝑦(𝑠,𝑡))∣0≤𝑠≤1} be the contour deformed from 𝒞(0)=𝒞^∗ at time 𝑡. The motion of the contour 𝒞(𝑡) will be determined by an appropriately chosen velocity 𝐔(𝑡) indicating a time change of tracking points (𝐫₁(𝑡),…,𝐫_𝑛(𝑡)): 𝐔(𝑡)𝐮∶=¹(𝑡)⋮𝐮_𝑛(𝑡)=𝑑𝐫𝑑𝑡¹(𝑡)⋮𝐫_𝑛(𝑡)𝐫with¹(0)⋮𝐫_𝑛(0)=𝐫^∗₁⋮𝐫^∗_𝑛. (7)

Here, we identify the contour 𝒞(𝑡) with tracking points (𝐫₁(𝑡),…,𝐫_𝑛(𝑡)).

In our method, 𝐔(𝑡) for each time 𝑡 is a minimizer of the following energy functional reflecting local-to-global deformation: ℰ_𝑡(𝐔)1∶=2^𝑛_𝑖=1×∫_{𝒩(𝐫𝑖(𝑡))}𝑤𝐫^′−𝐫_𝑖(𝑡)𝐮_𝑖⋅∇𝐼𝐫^′+𝜕,𝑡𝐼𝐫𝜕𝑡^′,𝑡²𝑑𝐫^′+𝜆𝐫_𝑖(𝑡)+𝐮_𝑖−𝑎₁(𝐔)𝑎₂(𝐔)𝑎₃(𝐔)𝑎₄(𝐔)𝐫^∗_𝑖−𝑎₅(𝐔)𝑎₆(𝐔)², (8) where 𝜆 is a nonnegative parameter, 𝑤 is the weight function as used in the LK method, and the affine coefficients 𝑎₁(𝐔),…,𝑎₆(𝐔) at time 𝑡 are given by 𝑎₁(𝐔)𝑎₃(𝐔)𝑎₂(𝐔)𝑎₄(𝐔)𝑎₅(𝐔)𝑎₆(𝐔)=Φ𝒞^∗𝑇Φ𝒞^∗−1Φ𝒞^∗𝑇×𝐫₁(𝑡)+𝐮₁^𝑇⋮𝐫_𝑛(𝑡)+𝐮_𝑛^𝑇, (9) where Φ𝒞^∗𝐫∶=^∗1𝑇1𝐫⋮⋮_∗𝑛𝑇1. (10)

The first term in (8) controls the individual tracking points to follow the local motions of a specific speckle pattern, while the second term controls their overall motions to be confined to the global motion constraint being approximately an affine transform of the initial tracking points.

The first term in (8) reflects the well-known LK optical flow (3) that probes local motions using blood-to-tissue intensity ratio.

The second term concerns a misfit between the estimated tracking points and their projection onto the space 𝒲, the space of affine transforms of the initial tracking points, given by 𝐫𝒲=^∗1𝑇1𝐫⋮⋮_∗𝑛𝑇1𝑎₁𝑎₃𝑎₂𝑎₄𝑎₅𝑎₆∶𝑎₁,…,𝑎₆∈ℝ⊂ℝ^𝑛×2. (11)

To be precise, a careful computation yields 𝐫theprojectionof¹(𝑡)+𝐮₁^𝑇⋮𝐫_𝑛(𝑡)+𝐮_𝑛^𝑇𝒞onto𝒲=Φ^∗Φ𝒞^∗𝑇Φ𝒞^∗−1Φ𝒞^∗𝑇𝐫₁(𝑡)+𝐮₁^𝑇⋮𝐫_𝑛(𝑡)+𝐮_𝑛^𝑇=𝐫^∗₁^𝑇1𝐫⋮⋮_∗𝑛𝑇1𝑎₁(𝐔)𝑎₃(𝐔)𝑎₂(𝐔)𝑎₄(𝐔)𝑎₅(𝐔)𝑎₆(𝐔). (12)

Hence, the second term in (8) with the above identity reflects a global motion involving contraction, expansion, translation, and rotation.

To compute the minimizer 𝐔 of the energy functional (8), we need to derive the Euler-Lagrange equation which can be obtained by taking partial derivative of ℰ_𝑡 with respect to each 𝐮_𝑗: 𝟎=𝜕ℰ_𝑡𝜕𝐮_𝑗=_{𝒩(𝐫𝑗(𝑡))}𝑤𝐫−𝐫_𝑗(𝑡)()×𝐮∇𝐼𝐫,𝑡_𝑗()+𝜕⋅∇𝐼𝐫,𝑡𝐼(𝐫𝜕𝑡𝐫,𝑡)𝑑𝐫+𝜆_𝑗(𝑡)+𝐮_𝑗−^𝑛_𝑖=1𝑑𝑖,𝑗𝐫_𝑖(𝑡)+𝐮_𝑖,for𝑗=1,…,𝑛, (13) where 𝑑(𝑖,𝑗) is the (𝑖,𝑗)-component of the 𝑛×𝑛 matrix 𝒫𝒞^∗∶=Φ𝒞^∗Φ𝒞^∗𝑇Φ𝒞^∗−1Φ𝒞^∗𝑇. (14)

The derivation of the Euler-Lagrange equation is given in the appendix.

For numerical algorithm, we replace the integral over 𝒩(𝐫_𝑗(𝑡)) in (13) by summation over pixels around 𝐫_𝑗(𝑡). Assuming that the neighborhood 𝒩(𝐫_𝑗(𝑡)) consists of 𝑚 pixels 𝐫_𝑗1,…,𝐫_𝑗𝑚, (13) becomes 𝟎=𝐴^𝑇_𝑗𝑊_𝑗𝐴_𝑗𝐮_𝑗+𝐴^𝑇_𝑗𝑊_𝑗𝐛_𝑗×𝐫+𝜆_𝑗(𝑡)+𝐮_𝑗−^𝑛_𝑖=1𝑑𝑖,𝑗𝐫_𝑖(𝑡)+𝐮_𝑖, (15) where 𝐴_𝑗=[∇𝐼(𝐫_𝑗1,𝑡),…,∇𝐼(𝐫_𝑗𝑚,𝑡)]^𝑇, 𝑊_𝑗=diag(𝑤(𝐫_𝑗1−𝐫_𝑗(𝑡)),…,𝑤(𝐫_𝑗𝑚−𝐫_𝑗(𝑡))), and 𝐛_𝑗=[(𝜕/𝜕𝑡)𝐼(𝐫_𝑗1,𝑡),…,(𝜕/𝜕𝑡)𝐼(𝐫_𝑗𝑚,𝑡)]^𝑇.

For notational simplicity, let the time 𝑡 be fixed and let 𝐮_𝑗𝑢∶=^𝑗𝑣_𝑗,𝐫_𝑗(𝑡)=𝑥_𝑗𝑦_𝑗,𝛼_𝑗𝛽_𝑗𝛽_𝑗𝛾_𝑗∶=𝐴^𝑇_𝑗𝑊_𝑗𝐴_𝑗𝜉+𝜆𝐼,^𝑗𝜂_𝑗=𝐴^𝑇_𝑗𝑊_𝑗𝐛_𝑗+𝜆𝐫_𝑗(𝑡). (16)

Then, the system (15) can also be represented by 0=𝛼_𝑗𝑢_𝑗+𝛽_𝑗𝑣_𝑗−𝜆^𝑛_𝑖=1𝑑𝑖,𝑗𝑢_𝑖+𝜉_𝑗−𝜆^𝑛_𝑖=1𝑑𝑖,𝑗𝑥_𝑖,0=𝛽_𝑗𝑢_𝑗+𝛾_𝑗𝑣_𝑗−𝜆^𝑛_𝑖=1𝑑𝑖,𝑗𝑣_𝑖+𝜂_𝑗−𝜆^𝑛_𝑖=1𝑑𝑖,𝑗𝑦_𝑖. (17)

This can be concisely written by Λ−𝜆𝒫𝒞^∗𝑈+𝐵𝑉=−Ξ+𝜆𝒫𝒞^∗𝒞𝑋,𝐵𝑈+Γ−𝜆𝒫^∗𝑉=−Π+𝜆𝒫𝒞^∗𝑌, (18) where Λ=diag(𝛼₁,…,𝛼_𝑛), 𝐵=diag(𝛽₁,…,𝛽_𝑛), Γ=diag(𝛾₁,…,𝛾_𝑛), 𝑈=[𝑢₁,…,𝑢_𝑛]^𝑇, 𝑉=[𝑣₁,…,𝑣_𝑛]^𝑇, Ξ=[𝜉₁,…,𝜉_𝑛]^𝑇, Π=[𝜂₁,…,𝜂_𝑛]^𝑇, 𝑋=[𝑥₁,…,𝑥_𝑛]^𝑇, and 𝑌=[𝑦₁,…,𝑦_𝑛]^𝑇. Using the block matrix form, we can rewrite it as the system of linear equations: Λ−𝜆𝒫𝒞^∗𝐵𝐵Γ−𝜆𝒫𝒞^∗𝑈𝑉=−Ξ+𝜆𝒫𝒞^∗𝑋−Π+𝜆𝒫𝒞^∗𝑌. (19)

Therefore, we can directly compute the movement 𝐔=[^𝑈𝑉] of size 𝑛×2 from the formula: 𝑈𝑉=Λ−𝜆𝒫𝒞^∗𝐵𝐵Γ−𝜆𝒫𝒞^∗−1−Ξ+𝜆𝒫𝒞^∗𝑋−Π+𝜆𝒫𝒞^∗𝑌, (20) because the column vectors of the block matrix (Λ−𝜆𝒫(𝒞^∗))𝐵𝐵(Γ−𝜆𝒫(𝒞^∗)) of size 2𝑛×2𝑛 are linearly independent.

For the parameter 𝜆=0, the displacements 𝐮_𝑗(𝑗∈{1,…,𝑛}) by (15) are exactly the same as those by the LK optical flow. However, (20) has a distinction to be capable of controlling the global shape in that the bigger the parameter 𝜆 is, the stronger the shape constraint is imposed. The LK optical flow performs a role as the local deformation subject to the global shape constraint, which is represented by the relationship of all 𝑛 tracking points. Therefore, each point efficiently tracks maintaining the global deformation of initial LV contour.

2.3. Heuristic Choice of Parameter 𝜆

For heuristic choice of parameter 𝜆, we use various datasets of manually delineated LV borders by clinical experts. With manually defined data 𝐫_𝑗,𝒞^∗,and𝐔(𝑡) in a given image 𝐼, we define the parameter 𝜆 as a function of quantity 𝐫_𝑗,𝒞^∗,𝐔(𝑡),𝐼, and time 𝑡: 𝜆∶=∑^𝑛_𝑗=1𝐴^𝑇_𝑗𝑊_𝑗𝐴_𝑗𝐮_𝑗(𝑡)+𝐴^𝑇_𝑗𝑊_𝑗𝐛_𝑗²₂∑^𝑛_𝑗=1𝐫_𝑗(𝑡)+𝐮_𝑗(𝑡)−∑^𝑛_𝑖=1𝑑𝑖,𝑗𝐫_𝑖(𝑡)+𝐮_𝑖(𝑡)²₂. (21)

We should note that if 𝐮_𝑗(𝑡) satisfies (15) for all 𝑡 and 𝑗=1,…,𝑛, then 𝜆=𝜆, the constant independent of time 𝑡.

From numerous experiments, we observed that 𝜆 tends to depend mainly on the contrast of the image 𝐼, and its dependency on time 𝑡 is relatively small. We found a linear relationship between log(𝐼_tissue/𝐼_blood) and 𝜆, where 𝐼_tissue/𝐼_blood is an overall tissue/blood intensity ratio.

To investigate behavior of the parameter 𝜆, we generate synthetic speckle images consisting of tissue and blood regions and test them by changing conditions including tissue/blood contrast as shown in Figure 3. We use an apical long-axis view template shown in Figure 4(a). When the synthetic images are generated, it is assumed that speckle is fully developed so that the statistics of echo envelope follow the Rayleigh distribution ([22, 23]) and, by log-compression, the distribution of the intensities is changed into the Fisher-Tippett distribution ([24, 25]) 𝑓_𝜎(𝐼)=2𝛼₁2exp𝛼₁𝐼−𝛼₂−ln2𝜎²×2−exp𝛼₁𝐼−𝛼₂−ln2𝜎², (22) where 𝜎 is the distribution parameter represented in Rayleigh distribution and 𝛼₁, 𝛼₂ are the predetermined system parameters for log-compression of echo envelope. Finally, the synthetic images are smoothed by low-pass filter (in Figures 4(b), 4(c), and 4(d), resp.).

Figure 3 

Image frames by varying tissue/blood intensity ratio. We use echographic texture modeling and heart motion modeling to generate image frames with various contrasts.

(a)

(b)

(c)

(d)

Figure 4 

Synthetic images: (a) original LV template, (b) speckle image with Rayleigh distribution, (c) with Fisher-Tippet distribution, and (d) its smoothed image by a Gaussian filter.

For modeling of heart motion, we simulate a heart with the nonrigid motion integrating global and local deformations. Figure 5(a) illustrates the deformation of LV in four simulated images. LV contours are represented by 13 tracking points and a natural cubic spline connecting them, which are denoted by 𝐫_𝑖 at each time. Their 𝑥𝑦-coordinates and displacements 𝐮_𝑖 from previous tracking points to next tracking points are listed in Table 1. For the sake of convenience in computation, it is assumed that the global deformation is modeled by an affine transformation of coefficients 𝑎₁=0.92, 𝑎₂=−0.03, 𝑎₃=−0.01, 𝑎₄=0.89, 𝑎₅=20, and 𝑎₆=6, which illustrate a contraction, and the local deformation is modeled by a free-form deformation of 0.1% variants with respect to the global deformation. In the first row of Figure 5(a), the blue solid lines and the red asterisks are showed as LV contour and tracking points by the defined heart motion, respectively. The green lines and asterisks mean LV contour at the previous frame. To generate the sequential images indicating the heart deformation, we also generate the tracking points of the epicardial contours so that the wall thickness between two contours is changed from 20 to 25 pixels in the sequential images. Using node points containing the endocardial and epicardial tracking points, the Delaunay triangulation meshes are generated and the sequential images are filled using linear spatial transformation from each mesh at previous image to the corresponding mesh at next image (second row).

Table 1 

The tracking points and displacements used in sequential synthetic images.

(a)

(b)

Figure 5 

Result showing the independency of 𝜆 with time 𝑡: (a) sequential synthetic images for myocardial motion. The first row shows the synthetic initial image and the tracking points representing the sequential motion of heart, the second row the sequential images corresponding to the motion of the tracking points, and the last row the tracking points (yellow “o” marks) and LV contour (yellow dotted line) by the displacements computed using 𝜆=300 and the mean value of 𝜆 computed using the sequential images. (b) The change of 𝜆 according to 𝑡. The 𝜆 value varies within the range of [250,350].

We first test the dependency of 𝜆 on the time 𝑡. For the given sequential synthetic images and tracking points at each time step, we compute 𝜆 and plot the change of 𝜆 with time 𝑡. The parameter 𝜆 varies within the range of 250 to 350 as shown in Figure 5(b). Using 𝜆=300, the mean value of 𝜆, we again compute (20) and get the displacements having the errors within 1 pixel compared to the reference displacements in Table 1(b). In this test, we use the 2-dimensional Gaussian function of variance 𝜎²_𝑥=𝜎²_𝑦=5² (pixel size) for the weight function 𝑤 over the square neighborhood with side length 21 pixels. From this test, we observe that the dependency of 𝜆 on the time 𝑡 is negligibly small.

Next, we test the dependency of 𝜆 on the tissue/blood intensity ratio. We generate the two consecutive images by varying the intensity of tissues as mentioned in Figure 3 and evaluate the change of 𝜆 with respect to the tissue/blood intensity ratio. Figure 6 shows that the relationship between the image intensity contrast and log₁₀𝜆 is approximately linear. This linear relationship enables us to provide a way of choosing the parameter 𝜆 depending on the tissue/blood intensity ratio effectively.

Figure 6 

Graphs showing the relationship between 𝜆 and the tissue/blood intensity contrast. The stars are the points obtained from the simulation; the straight line fits these points.

3. Experimental Results

We test the proposed algorithm in clinical setting using many real data. We compare the performance of the proposed algorithm with some widely used tracking algorithms including the block matching tracking methods using sum-of-squared difference (SSD) and cross-correlation coefficient, and the LK optical flow. For experiments, we use the 35 cases of 240×320 size 2D echocardiography data acquired using a Samsung Medison V10 ultrasound system (Seoul, Republic of Korea) and a phased array transducer P2-4BA (2–4 MHz). We use 19 tracking points to track the endocardial border and make the LV contour connecting the points using the natural cubic spline. All the experiments were conducted using MATLAB 7.5 and laptop computer (Inter processor U7300 at 1.3 GHz and 1 GB RAM), and the computational time was about 40 milliseconds at each frame.

3.1. Assessment of LV Border Tracing

A quantitative evaluation on the performance of the proposed tracking algorithm is done on real 2D image sequences. For computation of 𝐮⋅∇𝐼+(𝜕/𝜕𝑡)𝐼, we use the standard finite difference method. We use the Hausdorff distance 𝜀_𝐻 [26, 27] to compare the automated LV contours produced by algorithms with manually traced contours by a clinical expert. Here, the Hausdorff distance between the contour 𝒞₁ and 𝒞₂ is given by 𝜀_𝐻𝒞₁,𝒞₂=maxsup_{𝐫1∈𝐶1}min_{𝐫2∈𝐶2}𝐫₁−𝐫₂,sup_{𝐫2∈𝐶2}min_{𝐫1∈𝐶1}𝐫₁−𝐫₂. (23)

For two representative cases among the 35 cases of 2D echocardiography data, the LV tracking results of the proposed method and the conventional methods are shown in Figures 7 and 8, respectively. For the sake of a name, we call them Cases I and II.

Figure 7 

Case I: real images with weak signals in endocardial border. The second and third rows are the results by region-based tracking methods using sum-of-squared difference (SSD) and cross-correlation, respectively. The fourth row is the result by the LK optical flow and the final row is the result by the proposed method (𝜆=120).

Figure 8 

Case II: real images with indistinguishable speckle patterns in endocardial border. The second and third rows are the results by region-based tracking methods using sum-of-squared difference (SSD) and cross-correlation, respectively. The fourth row is the result by the LK optical flow and the fifth row by the proposed method (𝜆=100).

In Figure 7, the first row is manually traced LV contours by a clinical expert for images at ED, ES, and ED frames in the entire cycle. The next two rows are results by two block matching tracking methods using the different similarity measures of sum-of-squared difference (SSD) and cross-correlation. The fourth and fifth rows are obtained by the LK optical flow and the proposed method, respectively. The three conventional methods produce distorted LV contours due to a few incorrect tracking points alienated from the real LV border. On the other hand, the proposed method successfully follows local speckle patterns without distorting the whole LV shape.

For initial 10 sequential images, we compute 𝜆 by manually identifying each tracking point to the corresponding position on each image. From the computed parameters 𝜆, 𝜆 is set to 120 according to the 𝜆-choice method described in Section 2.3.

In Figure 8, we test for real images having indistinguishable speckle patterns near endocardial border. Due to the presence of indistinguishable speckle patterns, the three conventional methods produce irregular distribution of tracking points as shown in second, third, and fourth rows in Figure 8. The proposed method keeps regular distribution of tracking points and successfully track local speckle patterns. For Case II, 𝜆 is set to 100.

Figure 9 shows the comparison results of four different methods using Hausdorff distance between contours drawn manually and contours generated automatically for the entire cycle from an ED frame to the next ED frame. The proposed method provides the smallest errors in final tracking results of both Cases I and II.

(a) Case I 

(b) Case II

Figure 9 

Comparison results of LV contours using the Hausdorff distance 𝜀𝐻 for the entire images of (a) Case I and (b) Case II, which are shown in Figures 7 and 8, respectively.

3.2. Assessment of Individual Tracking Point Errors

For performance evaluation of the proposed algorithm, we propose an additional assessment regarding the repeatability of local point along the forward and backward entire cardiac cycle. Let {𝐫^initial₁,…,𝐫^initial_𝑛} be the set of initial tracking points on a manually delineated contour (see the images of the left column in Figure 7). Let 𝑡_𝑅 be a time interval of a one cycle image between ED frame and the next ED frame. Using one cycle image 𝐼(𝐫,𝑡),0≤𝑡≤𝑡_𝑅, we generate a forward-backward image defined by 𝐼(̃)=𝐼()𝐫,𝑡𝐫,𝑡if0≤𝑡≤𝑡_𝑅𝐼𝐫,2𝑡_𝑅−𝑡if𝑡_𝑅≤𝑡≤2𝑡_𝑅. (24)

Using this forward-backward image 𝐼(̃𝐫,𝑡),0≤𝑡≤2𝑡_𝑅, we apply an automated tracking algorithm to get the returning tracking position 𝐫^returning_𝑗 at time 𝑡=2𝑡_𝑅. The local tracking point assessment is obtained by estimating the distance between the initial position 𝐫^initial_𝑗and the corresponding returning position 𝐫^returning_𝑗: ()=Forward-backwardpointtrackingerrorFBTE1𝑛^𝑛_𝑖=1𝐫^initial_𝑖−𝐫^returning_𝑖². (25)

For the previous two representative cases, Cases I and II, Table 2 shows the comparison results of the proposed method with the conventional methods using the FBTE.

Table 2 

The comparison results of the proposed method with the conventional methods using FBTE, for Case I and Case II (in pixels).

Table 3 shows the mean and standard deviation of the forward-backward point tracking errors of the results obtained by the three conventional tracking methods and the proposed method. Tables 2 and 3 reveal that the proposed method provides improved performance compared with the conventional tracking methods.

Table 3 

The comparison results of the tracking algorithms for the total experimental dataset of 35 cases. The errors are measured using the FBTE (in pixels).

4. Discussion and Conclusion

The proposed method controls the individual tracking points following optical flow by confining their overall motions by penalizing the misfit between the estimated tracking points and their projection onto the affine transform space 𝒲 in (11) of the initial tracking points.

We have experimentally demonstrated that the proposed method is capable of performing robust real-time LV border tracking even in the presence of indistinguishable portions of the LV walls in echocardiography data. In practice, echocardiography data often contains edge dropout or indistinguishable speckle patterns in a local neighborhood of a tracking point which may bring the tracking point out of the endocardial border, resulting in distorted LV contours. The proposed method effectively deals with these problems by taking advantage of an LV shape space describing a global motion that is synthesized by integrating local deformations governed by the LK optical flow model. Various experiments show that the proposed method achieves better overall performance than the widely used conventional methods including the block matching tracking methods using sum-of-squared difference (SSD) and cross-correlation, and the LK optical flow.

The proposed method performs the LV border tracking by directly computing the displacements between two sequential images via a simple matrix multiplication. The computational time is affected by the size of the matrix, depending on the number of tracking points.

We also proposed a new performance evaluation method for LV tracking that is based on the forward-backward tracking error estimation as shown in Section 3.2. The conventional evaluation of global tracking performance using the delineated LV contours has some limitations in estimating errors of individual tracking points; in the case when tracking points erroneously move along LV border, the LV contour connecting the tracking points cannot reveal those individual tracking errors. The forward-backward point tracking error estimation provides a better local tracking performance assessment in the whole cycle.

The proposed technique can be extended to three dimensions by using 3D affine transformation as a global deformation.

Appendix

Derivation of the Euler-Lagrange Equation (13)

In this appendix, we derive the Euler-Lagrange equation (13) from (8). From (9), we have 𝜕𝑎₁𝜕𝑢_𝑗𝜕𝑎₂𝜕𝑢_𝑗𝜕𝑎₅𝜕𝑢_𝑗=Φ𝒞^∗𝑇Φ𝒞^∗−1𝑥^∗_𝑗𝑦^∗_𝑗1,𝜕𝑎₃𝜕𝑢_𝑗𝜕𝑎₄𝜕𝑢_𝑗𝜕𝑎₆𝜕𝑢_𝑗=000,𝜕𝑎₁𝜕𝑣_𝑗𝜕𝑎₂𝜕𝑣_𝑗𝜕𝑎₅𝜕𝑣_𝑗=000,𝜕𝑎₃𝜕𝑣_𝑗𝜕𝑎₄𝜕𝑣_𝑗𝜕𝑎₆𝜕𝑣_𝑗=Φ𝒞^∗𝑇Φ𝒞^∗−1𝑥^∗_𝑗𝑦^∗_𝑗1, (A.1) for 𝑗=1,2,…,𝑛. From the first and fourth identities of (A.1), we have 𝑑𝑖,𝑗=𝜕𝑎₁𝜕𝑢_𝑗𝑥^∗_𝑖+𝜕𝑎₂𝜕𝑢_𝑗𝑦^∗_𝑖+𝜕𝑎₅𝜕𝑢_𝑗=𝑥^∗_𝑖𝑦^∗_𝑖1𝜕𝑎₁𝜕𝑢_𝑗𝜕𝑎₂𝜕𝑢_𝑗𝜕𝑎₅𝜕𝑢_𝑗=𝑥^∗_𝑖𝑦^∗_𝑖1Φ𝒞^∗𝑇Φ𝒞^∗−1𝑥^∗_𝑗𝑦^∗_𝑗1=𝑥^∗_𝑖𝑦^∗_𝑖1𝜕𝑎₃𝜕𝑣_𝑗𝜕𝑎₄𝜕𝑣_𝑗𝜕𝑎₆𝜕𝑣_𝑗=𝜕𝑎₃𝜕𝑣_𝑗𝑥^∗_𝑖+𝜕𝑎₄𝜕𝑣_𝑗𝑦^∗_𝑖+𝜕𝑎₆𝜕𝑣_𝑗. (A.2)

Using (A.1) and (A.2), the partial derivatives 𝜕ℰ_𝑡/𝜕𝑢_𝑗 and 𝜕ℰ_𝑡/𝜕𝑣_𝑗 of ℰ_𝑡 with respect to 𝐮_𝑗=(𝑢_𝑗,𝑣_𝑗) are computed by 𝜕ℰ_𝑡𝜕𝑢_𝑗=_{𝒩(𝐫𝑗(𝑡))}𝑤𝐫^′−𝐫_𝑗(𝑡)×𝐮_𝑗⋅∇𝐼𝐫^′+𝜕,𝑡𝐼𝐫𝜕𝑡^′𝜕,𝑡𝐼𝐫𝜕𝑥^′,𝑡𝑑𝐫^′𝑥+𝜆^𝑗+𝑢_𝑗−𝑎₁𝑥^∗_𝑗+𝑎₂𝑦^∗_𝑗+𝑎₅−∑^𝑛_𝑖=1𝑥_𝑖+𝑢_𝑖−𝑎₁𝑥^∗_𝑖+𝑎₂𝑦^∗_𝑖+𝑎₅𝜕𝑎₁𝜕𝑢_𝑗𝑥^∗_𝑖+𝜕𝑎₂𝜕𝑢_𝑗𝑦^∗_𝑖+𝜕𝑎₅𝜕𝑢_𝑗=_{𝒩(𝐫𝑗(𝑡))}𝑤𝐫^′−𝐫_𝑗(𝑡)×𝐮_𝑗⋅∇𝐼𝐫^′+𝜕,𝑡𝐼𝐫𝜕𝑡^′𝜕,𝑡𝐼𝐫𝜕𝑥^′,𝑡𝑑𝐫^′+𝜆𝑥_𝑗+𝑢_𝑗−𝑎₁𝑥^∗_𝑗+𝑎₂𝑦^∗_𝑗+𝑎₅−^𝑛_𝑖=1𝑑𝑖,𝑗𝑥_𝑖+𝑢_𝑖−𝑎₁𝑥^∗_𝑖+𝑎₂𝑦^∗_𝑖+𝑎₅,𝜕ℰ_𝑡𝜕𝑣_𝑗=_{𝒩(𝐫𝑗(𝑡))}𝑤𝐫^′−𝐫_𝑗(𝑡)×𝐮_𝑗⋅∇𝐼𝐫^′+𝜕,𝑡𝐼𝐫𝜕𝑡^′𝜕,𝑡𝐼𝐫𝜕𝑦^′,𝑡𝑑𝐫^′𝑦+𝜆^𝑗+𝑣_𝑗−𝑎₃𝑥^∗_𝑗+𝑎₄𝑦^∗_𝑗+𝑎₆−∑^𝑛_𝑖=1𝑦_𝑖+𝑣_𝑖−𝑎₃𝑥^∗_𝑖+𝑎₄𝑦^∗_𝑖+𝑎₆𝜕𝑎₃𝜕𝑣_𝑗𝑥^∗_𝑖+𝜕𝑎₄𝜕𝑣_𝑗𝑦^∗_𝑖+𝜕𝑎₆𝜕𝑣_𝑗=_{𝒩(𝐫𝑗(𝑡))}𝑤𝐫^′−𝐫_𝑗(𝑡)×𝐮_𝑗⋅∇𝐼𝐫^′+𝜕,𝑡𝐼𝐫𝜕𝑡^′𝜕,𝑡𝐼𝐫𝜕𝑦^′,𝑡𝑑𝐫^′+𝜆𝑦_j+𝑣_𝑗−𝑎₃𝑥^∗_𝑗+𝑎₄𝑦^∗_𝑗+𝑎₆−^𝑛_𝑖=1𝑑𝑖,𝑗𝑦_𝑖+𝑣_𝑖−𝑎₃𝑥^∗_𝑖+𝑎₄𝑦^∗_𝑖+𝑎₆. (A.3)(A.4)

Simple arrangement of the above identities leads to the Euler-Lagrange equation: 𝟎=𝜕ℰ_𝑡𝜕𝐮_𝑗=_{𝒩(𝐫𝑗(𝑡))}𝑤𝐫−𝐫_𝑗(𝑡)()×𝐮∇𝐼𝐫,𝑡_𝑗()+𝜕⋅∇𝐼𝐫,𝑡𝐼(𝐫𝜕𝑡𝐫,𝑡)𝑑𝐫+𝜆_𝑗(𝑡)+𝐮_𝑗−𝐠_𝑗(𝐔)−^𝑛_𝑖=1𝑑𝑖,𝑗𝐫_𝑖(𝑡)+𝐮_𝑖−𝐠_𝑖(𝐔), (A.5) where 𝐠_𝑖(𝐔)=^{𝑎1(𝐔)𝑎2}_𝑎^(𝐔)_{3(𝐔)𝑎4(𝐔)}𝐫^∗_𝑖+^𝑎5_𝑎^(𝐔)_6(𝐔). The Euler-Lagrange equation (13) can be obtained by a careful rearrangement of (13) using the identity 𝐠_𝑗(𝐔)=^𝑛_𝑖=1𝑑𝑖,𝑗𝐠_𝑖(𝐔). (A.6)

Hence, it remains to prove the above identity. It suffices to prove ^𝑛_𝑖=1𝑥^∗_𝑖𝑑𝑖,𝑗=𝑥^∗_𝑗,^𝑛_𝑖=1𝑦^∗_𝑖𝑑𝑖,𝑗=𝑦^∗_𝑗,^𝑛_𝑖=1𝑑𝑖,𝑗=1, (A.7) for 𝑗=1,2,…,𝑛. To compute ∑^𝑛_𝑖=1𝑥^∗_𝑖𝑑(𝑖,𝑗) simply, denote 𝜙_𝑗𝜑_𝑗𝜓_𝑗∶=^𝑛_𝑘∑₌₁𝑥_𝑘^𝑛∗2∑_𝑘=1𝑥^∗_𝑘𝑦∑^∗_𝑘^𝑛_𝑘=1𝑥^∗_𝑘𝑛𝑘∑₌₁𝑥^∗_𝑘𝑦∑^∗_𝑘^𝑛_𝑘=1𝑦_𝑘^𝑛∗2∑_𝑘=1𝑦^∗_𝑘𝑛𝑘∑₌₁𝑥∑^∗_𝑘^𝑛_𝑘=1𝑦^∗_𝑘𝑛⁻¹𝑥^∗_𝑗𝑦^∗_𝑗1. (A.8)

Then ∑^𝑛_𝑖=1𝑥^∗_𝑖𝑑(𝑖,𝑗) can be simplified by ^𝑛_𝑖=1𝑥^∗_𝑖𝑑=𝑖,𝑗_𝑛𝑖=1𝑥^∗_𝑖𝑥^∗_𝑖𝑦^∗_𝑖1×𝑥^∗₁𝑦^∗₁1⋮𝑥^∗_𝑛𝑦^∗_𝑛1^𝑇𝑥^∗₁𝑦^∗₁1⋮𝑥^∗_𝑛𝑦^∗_𝑛1⁻¹𝑥^∗_𝑗𝑦^∗_𝑗1=∑^𝑛_𝑖=1𝑥_𝑖^𝑛∗2∑_𝑖=1𝑥^∗_𝑖𝑦∑^∗_𝑖^𝑛_𝑖=1𝑥^∗_𝑖×^𝑛_𝑘∑₌₁𝑥_𝑘^𝑛∗2∑_𝑘=1𝑥^∗_𝑘𝑦∑^∗_𝑘^𝑛_𝑘=1𝑥^∗_𝑘𝑛𝑘∑₌₁𝑥^∗_𝑘𝑦∑^∗_𝑘^𝑛_𝑘=1𝑦_𝑘^𝑛∗2∑_𝑘=1𝑦^∗_𝑘𝑛𝑘∑₌₁𝑥∑^∗_𝑘^𝑛_𝑘=1𝑦^∗_𝑘𝑛⁻¹𝑥^∗_𝑗𝑦^∗_𝑗1=∑^𝑛_𝑖=1𝑥_𝑖^𝑛∗2∑_𝑖=1𝑥^∗_𝑖𝑦∑^∗_𝑖^𝑛_𝑖=1𝑥^∗_𝑖𝜙_𝑗𝜑_𝑗𝜓_𝑗=𝜙_𝑗^𝑛_𝑖=1𝑥_𝑖^∗2+𝜑_𝑗^𝑛_𝑖=1𝑥^∗_𝑖𝑦^∗_𝑖+𝜓_𝑗^𝑛_𝑖=1𝑥^∗_𝑖=𝑥^∗_𝑗. (A.9)

This completes the first identity of (A.7). Similarly, we can get the remaining two identities of (A.7).

Acknowledgments

The author would like to thank Samsung Medison R&D Center for supporting many ultrasound data and valuable suggestions. This work was supported by the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (no. 2012-8-1798 and no. 2011-8-1782).