Cardiovascular Disease (CVD) is the leading cause of death in the modern world. The mortality rate associated with CVD was estimated to be 17 million in 2005 and continues to be ranked as the top killer worldwide. CVD is the result of under-supply of the cardiac tissue and can lead to malfunction of the involved myocardial territories and manifest as hypokinesia or akinesia. Several imaging methods such as X-ray CT, MRI, and Ultrasound have been used for visualization of the heart function. MRI and X-ray CT provide excellent spatial resolution but the cost and lack of wide-spread availability cause challenges in the clinical settings. Echocardiography is a popular technique for cardiac imaging due to its availability, ease of use, and low cost. Echocardiography shows the motion and anatomy of the heart in real time, enabling physicians to detect different pathologies. However, analysis of motion of the myocardium in echocardiographic images is based on visual grading by an observer and suffers from inter and intra-observer variability. To overcome the inter- and intra-observer variability, computerized image analyses can help by quantitatively interpreting the data. To that end, cardiac image processing techniques, mainly categorized as segmentation and registration, have been widely used for assessing the regional function of the heart [1–3]. To perform such analysis, two independent techniques, have been utilized; these are TDI (Tissue Doppler Imaging) and speckle tracking. TDI computes the tissue motion based on the Doppler phenomenon and is dependent on the angle of insonification. Speckle tracking, on the other hand,is an image processing method based on the analysis of the ultrasound B-mode or RF images. B-mode based algorithms are robust to the variation of the transducer angle but rely entirely on the properties of echocardiographic images which may be noisy or inaccurate. The physical principle underlying B-mode and TDI are to a large degree independent and therefore for myocardial motion estimation carry complementary information [4, 5].

Many methods such as optical flow [6], feature tracking [7], level sets [8], block matching [9], and elastic registration [10] have been utilized for quantitative assessment of myocardial motion in B-mode images. Table 1 shows a description of some of the current methods used in motion estimation in echocardiographic images [6–21]. Suhling et al. [6] integrated rigid registration in an optical flow framework in order to detect myocardial motion from 2D echocardiographic images. B-spline moments invariants were applied to echo images to achieve invariance to the translation and rotation. The motion estimation algorithm was then applied to the B-spline moments of the image instead of the image intensity in a coarse to fine strategy and was validated using open chested dogs after ligation of a coronary artery. Additional validations were performed on simulation and phantom images. Ellen et al. [10] used elastic registration on 3D B-mode echocardiography images to extract myocardial motion and strain values. The method was validated using simulated and real ultrasoundimages. Esther-Leung et al. [11] proposed two different methods (1. model-driven, 2. edge-driven) for tracking the left-ventricular wall in echocardiographic images. Their approach was motivated by the fact that in echocardiography images, visibility of the myocardium depends on the imaging view; so the myocardium may be, partly, invisible to the beam. Their technique relied on a local data-driven tracker using optical flow applied to the visible parts of the myocardium and a global statistical model applied to the invisible parts. It was concluded that the shape model could render good results for both the visible and the invisible tissues in ultrasound images. Myronenco et al. [12] proposed the so-called Coherent Point Drift (CPD) technique for myocardial motion estimation, constraining the motion of the point set in the temporal direction for both rigid and nonrigid point set registration. A set of point distribution was computed based on endocardium and epicardium locations. The point set was modeled with a Gaussian mixture model (GMM). The GMM centroids were updated coherently in a global pattern using maximum likelihood to preserve the topological structure of the point sets. A motion coherence constraint was added based on regularization of the displacement fields. The purpose of regularization was to increase the motion smoothness.

Most of the motion estimation techniques developed thus far, are either based on TDI or B-mode. Recently, Garcia et al. [22] considered the combination of cardiac B-mode images and intra-cardiac blood flow data for computing the blood flow motion in the heart using continuity equation and mass conservation in polar coordinates. Their paper focused on the blood flow computation and did not consider the cardiac tissue displacements. Dalen et al. [23] and Amundsen et al. [24] previously combined TDI with speckle tracking by integrating TDI in the beam direction and speckle tracking in the direction lateral to the beam. However, this method discarded the speckle tracking data in the beam direction. The authors reported that they were unable to improve the motion estimation performance compared to speckle tracking techniques. In this paper, we propose integration of tissue Doppler and speckle tracking within a novel optical flow framework, we call TDIOF (Tissue Doppler Optical Flow). Our experimental results indicate that TDIOF outperforms both TDI and speckle tracking approaches.

The organization of the rest of the paper is as follows: in Sect. 2, we review the mathematical and algorithmic basis for the proposed method. Section 3 is a description of datasets used for validation of the proposed method. These include computational simulations, US data collected in a cardiac phantom, and in vivo data collected in patients recruited from the echocardiography laboratory at the Robley Rex VA Medical Center to our IRB-approved study. In Sect. 4, strain computations are discussed and, in Sect. 5, results from validation of the proposed method are described. Finally, in Sect. 6 we discuss observations related to TDIOF and our findings.

TDI and B-mode speckle tracking are different in both their physical underpinning and data type. In speckle tracking, tissue motion is determined from motion of speckles in Ultrasound images–typically, using a block matching approach applied to 2-D B-mode images [give references]. Although speckle tracking provides both components of the motion in the spatial domain, it is based on noisy B-mode images. TDI on the other hand only computes the velocity of tissue in one direction as moving towards (displayed as red) or away (displayed as blue) from the transducer. This means that the computed motion is the projection of the real motion in the direction of the transducer and therefore TDI is angle dependent. In this section, we first describe a novel energy minimization framework for estimation of myocardial motion from B-mode images which incorporates a velocity constraint from TDI.

The proposed method is based on optimization of three energy functions: (1) intensity constancy assumption, (2) velocity smoothness, and (3) similarity with Doppler data. The framework is minimized using an incremental algebra in method [25, 26] as described in Appendix A. In order to show the performance of TDIOF, it is compared to two popular motion detection techniques Horn-Schunck optical flow [27] and block matching [28] (Appendix B). Block matching is utilized in several commercial software platforms [29].

Strain is a measure of deformation of the cardiac tissue. With representing the identity matrix, the Lagrangian strain tensor at a given myocardial point and for a specific time point can be expressed as:
where the elements of the deformation gradient tensor, F, are:
while represents the spatial coordinates in the undeformed coordinates (typically taken to be the end-diastolic frame), and is the accumulated motion vector at the corresponding spatial location relative to the undeformed state . For the echocardiography data, the reference frame for the strain computation was considered to be the end diastolic frame and was selected based on ECG trigger. The deformation field was then computed between each two frames and was added to the motion field from the previous frame in order to measure the accumulated deformation and strain . Since the deformation field of the consecutive frames do not represent the motion of the same pixels, spline interpolation was used to align the deformation fields. For the tagged MRI data, the end-diastolic frame was always the first acquired image which is collected immediately after the R-wave trigger.

The normal strain in the direction of the unit vector n can be calculated from the Lagrangian strain tensor through the quadratic form , where n is a unit vector and can point to any direction on the unit sphere. Due to the geometry of the left ventricle, the normal strains are usually calculated in radial, circumferential, and longitudinal directions.

Regional analysis is performed on 17 American Heart Association (AHA) prescribed segments. Figure 7 shows the different segments. The acronyms stand for antero-septal (AS), anterior (A), lateral (L), posterior (P), inferior (I), and infero-septal (IS). For a review of topics related to determination of strain from cardiac images, the reader is referred to [35, 37, 38].

As noted in Sect. 3, TDIOF was applied to three different datasets: simulated images, data collected in a physical phantom, and in vivo data (both set A and set B). To further elucidate the effect of the Doppler term, results from TDIOF were compared to Horn-Schunck (HS) optical flow (beta = 0 in Eq. (A.11)) and block-matching (BM) (see Appendix) with the latter being the basis for most commercial speckle tracking methods [13]. Since the performance of each technique depends on the parameters of the method, it was necessary to optimize the parameters. Based on simulated images, an exhaustive search was performed over the parameters of TDIOF, HS optical flow, and BM speckle tracking method (a large range was considered for each parameter) and the best values were selected experimentally. To analyze the performance of the techniques with different parameter settings, simulated images were compared to the next simulated frame after being warped using the estimated motion field. Relative mean absolute error was used for the comparison. Relative mean absolute error was computed as ; where and are the first and subsequent warped images, and is the total number of points. Figure 3 shows the performance of the TDIOF technique using different parameters. The methods were then applied to all the datasets using the resulting parameters: number of scales for multiscale implementation: (smoothness weight): (TDI similarity weight): 0.001, and (penalizer parameter): 0.1. The parameters for the HS technique were set as follows: number of scales 5, and smoothness weight 2000.