Keywords

These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

1 Introduction

Characterization of left ventricular myocardial deformation is useful for the detection and diagnosis of cardiovascular diseases. Conditions such as ischemia and infarction undermine the contractile property of the LV and analyzing Lagrangian strains is one way of identifying such abnormalities.

Numerous methods calculate dense myocardial motion fields and then compute strains using echo images. Despite the substantial interest in Lagrangian motion and strains, and some recent contributions in spatio-temporal tracking [1, 2], most methods typically calculate frame-to-frame or Eulerian displacements first and then obtain Lagrangian trajectories.

Therefore, we propose a novel dynamic shape tracking (DST) method that first provides Lagrangian motion trajectories of points in the LV surfaces and then computes dense motion fields to obtain Lagrangian strains. This approach aims to reduce the drift problem which is prevalent in many frame-to-frame tracking methods.

We first segment our images and obtain point clouds of LV surfaces, which are then set up as nodes in a directed acyclic graph. Nearest neighbor relationships define the edges, which have weights based on the Euclidean distance and difference in shape properties between neighboring and starting points. Finding the trajectory is then posed as a shortest path problem and solved using Dijkstra’s algorithm and dynamic programming.

Once we obtain trajectories, we calculate dense displacement fields for each frame using radial basis functions (RBFs), which are regularized using sparsity and incompressibility constraints. Since the trajectories account for motion primarily in the LV boundaries, we also include mid-wall motion vectors from frame-to-frame RF speckle tracking. This fusion strategy was originally proposed in [3] and expanded in [4]. Ultimately, we calculate Lagrangian strains and validate our results by comparing with sonomicrometry based strains.

2 Methods

2.1 Dynamic Shape Tracking (DST)

We first segment our images to obtain endocardial and epicardial surfaces using an automated (except the first frame) level-set segmentation method [5]. There are surfaces from N frames, each with K points, and each point has M neighbors (Fig. 1a, with \(M = 3\)). \(x \in \mathbb {R}^{N\times K \times 3}\) is the point matrix, \(F \in \mathbb {R}^{N\times K \times S}\) is the shape descriptor matrix (where S is the descriptor length) and \(\eta \in \mathbb {R}^{N\times K \times M}\) is the neighborhood matrix (shape context feature is described in [6]). The \(j^{th}\) point of the \(i^{th}\) frame is indexed as \(x_{i, j}\) (\(i \in [1:N]\) and \(j \in [1:K]\)).

Fig. 1.
figure 1

Points, shape descriptors and trajectories.

Full size image

Let \(A_j(i): \mathbb {N}^N \mapsto \mathbb {N}^K\) be the set that indexes the points in the trajectory starting at point \(x_{1, j}\). For any point \(x_{i, A_j(i)}\) in the trajectory, we assume that:

  1. 1.

    It will not move too far away from the previous point \(x_{i-1, A_j(i-1)}\) and the starting point \(x_{1, j}\). Same applies to its shape descriptor \(F_{i, A_j(i)}\).

  2. 2.

    Its displacement will not differ substantially from that of the previous point. Same applies to its shape descriptors. We call these the \(2^{nd} \textit{ order weights}\).

Trajectories satisfying the above conditions will have closely spaced consecutive points with similar shape descriptors (Fig. 1b). They will also have points staying close to the starting point and their shape descriptor remaining similar to the starting shape descriptor, causing them to be more closed. The \(2^{nd} \textit{ order weights}\) enforce smoothness and shape consistency.

Let \(\mathbb {A}\) represent the set of all possible trajectories originating from \(x_{1, j}\). Hence, \(A_j \in \mathbb {A}\) minimizes, the following energy:

$$\begin{aligned} \begin{aligned} \hat{A}_j =&\ \mathop {{{\mathrm{argmin}}}}\limits _{{A_{j} \in \mathbb {A}}} \sum _{i=2}^{N} \lambda _1||x_{i, A_j(i)} - x_{i-1, A_j(i-1)}|| + \lambda _2||x_{i, A_j(i)} - x_{1, j}|| \\&+ \lambda _3||F_{i, A_j(i)} - F_{i-1, A_j(i-1)}|| + \lambda _4||F_{i, A_j(i)} - F_{1, j}|| + 2^{nd} \textit{ order weights} \end{aligned} \end{aligned}$$
(1)

Graph based techniques have been used in particle tracking applications such as [7]. We set each point \(x_{i, j}\) as a node in our graph. Directed edges exist between a point \(x_{i, j}\) and its neighbors \(\eta _{i, j}\) in frame \(i+1\). Each edge has an associated cost of traversal (defined by Eq. 1). The optimal trajectory \(\hat{A}_j\) is the one that accrues the smallest cost in traversing from point \(x_{i, j}\) to the last frame. This can be solved using Dijkstra’s shortest path algorithm.

Algorithm. Because our search path is causal, we don’t do all edge weight computations in advance. We start at \(x_{1, j}\) and proceed frame-by-frame doing edge cost calculations between points and their neighbors and dynamically updating a cost matrix \(E \in \mathbb {R}^{N \times K}\) and a correspondence matrix \(P \in \mathbb {R}^{N \times K}\). The search for the trajectory \(A_J\) stemming from a point \(j = J\) in frame 1 is described in Algorithm 1.

figure a

2.2 Speckle Tracking

We use speckle information, which is consistent in a small temporal window (given sufficient frequency), from the raw radio frequency (RF) data from our echo acquisitions, and correlate them from frame-to-frame to provide mid-wall displacement values.

A kernel of one speckle length, around a pixel in the complex signal, is correlated with neighboring kernels in the next frame [8]. The peak correlation value is used to determine the matching kernel in the next frame and calculate displacement and confidence measure. The velocity in the echo beam direction is further refined using zero-crossings of the phase of the complex correlation function.

2.3 Integrated Dense Displacement Field

We use a RBF based representation to solve for a dense displacement field U that adheres to the dynamic shape (\(U_{sh}\)) and speckle tracking (\(U_{sp}\)) results and regularize the field in the same manner as [4]:

$$\begin{aligned} \begin{aligned}&\hat{w} = \mathop {{{\mathrm{argmin}}}}\limits _{w} {f_{\mathrm {adh}}(Hw, U^{\mathrm {sh}}, U^{\mathrm {sp}})} + \lambda _1 ||w||_1 + \lambda _2 f_{\mathrm {biom}}(U)\\&f_{\mathrm {biom}}(U) = ||\nabla \cdot U||^2_2 + \alpha ||\nabla U||^2_2. \end{aligned} \end{aligned}$$
(2)

Here, U is parametrized by Hw, H represents the RBF matrix, w represents the weights on the bases, \(f_{adh}\) is the squared loss function, \(f_{biom}\) is the penalty on divergences and derivatives, which along with the \(l_1\) norm penalty results in smooth and biomechanically consistent displacement fields. This is a convex optimization problem and can be solved efficiently. The \(\lambda \)’s here and in Eq. 1 are chosen heuristically and scaled based on the number of frames in the images.

3 Experiment and Results

4DE Data (Acute Canine Studies). We acquired 4DE images from 7 acute canine studies (open chested, transducer suspended in a water bath) in the following conditions: baseline (BL), severe occlusion of the left anterior descending (LAD) coronary artery (SO), and severe LAD occlusion with low dose dobutamine stress (SODOB, dosage: \(5\,\mu g \backslash kg \backslash min\)). All BL images were used while 2 SO and SODOB images were not used due to lack of good crystal data. Philips iE33 ultrasound system, with the X7-2 probe and a hardware attachment that provided RF data were used for acquisition. All experiments were conducted in compliance with the Institutional Animal Care and Use Committee policies.

Sonomicrometry Data. We utilized an array of sonomicrometry crystals (Sonometrics Corp) to validate strains calculated via echo. 16 crystals were implanted in the anterior LV wall. They were positioned with respect to the LAD occlusion and perfusion characteristics within the crystal area, which are defined accordingly: ischemic (ISC), border (BOR) and remote (REM) (see Fig. 2). Three crystals were implanted in the apical (1) and basal (2) regions (similar to [9]). We adapted the 2D sonomicrometry based strain calculation method outlined in [10] to our 3D case. Sonomicrometric strains were calculated for each cube and used for the validation of echo based strains.

Fig. 2.
figure 2

Crystals and their relative position in the LV

Full size image

Agreement with Crystal Data. In Fig. 3a, we show, for one baseline image, strains calculated using our method (echo) and using sonomicrometry (crys). We can see that the strain values correlate well and drift error is minimal. In Fig. 3b, we present bar graphs to explicitly quantify the final frame drift as the distance between actual crystal position and results from tracking. We compare the results from this method (DST) against that of GRPM (described in [4, 11]), which is a frame-to-frame tracking method, in BL and SO conditions for 5 canine datasets.

The last frame errors for crystal position were lower and statistically significant for both BL an SO conditions (\(p < .01\)).

Strain Correlation. Pearson correlations of echo based strain curves (calculated using the shape based tracking (SHP) only and the shape and speckle tracking combined (COMB) with corresponding crystal based strain curves, across the cubic array regions (ISC, BOR, REM) for all conditions are summarized in Table 1. We see slightly improved correlations from SHP to COMB method. Correlation values were generally lower for ischemic region and longitudinal strains for both methods.

Table 1. Mean correlation values across regions for SHP and COMB methods.
Full size table
Fig. 3.
figure 3

Peak strain bar graphs (with mean and standard deviations) for radial, circumferential and longitudinal strains - shown across ISC, BOR and REM regions for echo and crystal based strains.

Full size image

Since we only had a few data points to carry out statistical analysis in this format, we also calculated overall correlations (with strain values included for all time points and conditions together, \(n > 500\)) and computed statistical significance using Fisher’s transformation. Change in radial strains (SHP \(r = .72\) to COMB \(r = .82\)) was statistically significant (\(p < .01\)), while circumferential (SHP \(r = .73\) to COMB \(r = .75\)) and longitudinal (SHP \(r = .44\) to COMB \(r = .41\)) were not.

Physiological Response. Changes in the crystal and echo based (using the combined method) strain magnitudes, across the physiological testing conditions - BL, SO and SODOB, is shown in Fig. 3.

Both echo and crystal strain magnitudes generally decreased with severe occlusion and increased beyond baseline levels with low dose dobutamine stress. The fact that functional recovery was observed with dobutamine stress indicates that, at the dose given ischemia was not enhanced. Rather, it appears that the vasodilatory and inotropic effects of dobutamine were able to overcome the effects of the occlusion.

However, in average, the strain magnitude recovery is less in the ISC region compared to BOR and REM regions for both echo and crystals. For echo, the overall physiological response was more pronounced for radial strains.

4 Conclusion

The DST method has provided improved temporal regularization and therefore drift errors have been reduced, specially in the diastolic phase. A combined dense field calculation method that integrates the DST results with RF speckle tracking results provided good strains, which is validated by comparing with sonomicrometry based strains. The correlation values were specifically good for radial and circumferential strains.

We also studied how strains vary across the ISC, BOR and REM regions (defined by the cuboidal array of crystals in the anterior LV wall) during the BL, SO and SCODOB conditions. Strain magnitudes (particularly radial) varied in keeping with the physiological conditions, and also in good agreement with the crystal based strains.

We seek to improve our methods as we notice that the longitudinal strains and strains in the ischemic region were not very good. Also, the DST algorithm occasionally resulted in higher error at end systole. Therefore, in the future, we will enforce spatial regularization directly by solving for neighboring trajectories together, where the edge weights will be influenced by the neighboring trajectories. We would also like to extend the method to work with point sets generated from other feature generation processes than segmentation.