Modified optical flow technique for cardiac motions analysis in echocardiography images.

The quantitative analysis of cardiac motions in echocardiography images is a noteworthy issue in processing of these images. Cardiac motions can be estimated by optical flow (OF) computation in different regions of image which is based on the assumption that the intensity of a moving pattern remains constant in consecutive frames. However, in echocardiographic sequences, this assumption may be violated because of unique specifications of ultrasound. There are some methodsapplying the brightness variation effect in OF. Almost all of them have presented a mathematical brightness variation model globally in the images. Nevertheless, there is not a brightness variation model for echocardiographic images in these methods. Therefore, we are looking for a method to apply brightness variations locally in different regions of the image. In this study, we proposed a method to modify ausual OF technique by considering intensity variation. To evaluate this method, we implement two other OF-based methods, one usual OF method and a modified OF method applying brightness variation as a multiplier and an offset (generalized dynamic imaging model [GDIM]) and compare them with ours. These algorithms and ours were implemented on real 2D echocardiograms. Our method resulted in more accurate estimations than two others. At last, we compared our method with expert'spoint of view and observed that three distance metrics between them was appropriately smaller than other methods. The Haussdorff distance between the estimated curve defined by the proposed method and the expert defined curve is 4.81 pixels less than this distance for Lucas-Kanade and 2.28 pixels less than GDIM.

INTRODUCTION

Cardiac motion analysis and tracking is surveyed in various researches. Several methods have been proposed to estimate heart motions from echocardiograms. One of the most popular and effective methods for motion tracking in cardiac images is optical flow (OF). This technique is utilized for the purpose of motion estimation of 2D cardiac images in many studies.[1-6]

OF calculation was first introduced by Horn and Schunck.[7] This technique is one of the most important and well-known methods in the context of machine vision. OF can be defined as the motion of brightness patterns in a sequence of images.[8] There are four groups of methods for OF calculation: differential methods, block-matching methods, energy-based methods, and phase-based methods. A complete survey of these four methods is presented in.[9] Most of these methods are based on the assumption that intensity of each region of a moving image remains constant in consecutive frames. In fact, in OF formulation, the main idea is that the intensity of a set of pixels is preserved between the times t and t+δt. However, this basic assumption is not satisfied in all cases. The intensity of an image may be inconstant for some reasons like motions of the object (translation or rotation), motions of the light source, or variations of brightness in the entire image.[10]

In echocardiography images, the intensity constancy assumption is violated in many regions because of complicated motions of beating heart, including translation, rotation, shear, contraction, and expansion.[1112] These complex motions, leading to displacement of the heart myocardial fibers in the direction of transducer, may result in entering or removing some pixels in some areas in the next frame. In fact, cyclic variations of the echocardiography images, as a function of the cardiac muscle behavior, depend on the angle of myocardial fibers respect to ultrasonic wave propagation. This fact leads to intensity variation of these areas in two consecutive frames. In addition, the intensity of image may alter for the effect of large displacements like the displacements of mitral valve. Also, speckle noise and its variations in each frame can influence the intensity. Intensity inconstancy, which leads to violation of the basic assumption of OF calculations, is an important source of error in motion estimation. Considering the intensity variations of moving patterns in OF computations could bring about more precise and less erroneous motion estimations.

In some of the previous studies, the effect of intensity variation in OF computations is surveyed.[1013-19] These studies resulted in more accurate estimations than common OF techniques. Nevertheless, none of them have presented a method for considering the intensity inconstancy in OF in echocardiography images. This matter motivated us to present a new method for considering intensity variations in echocardiography sequences.

A physical model is presented for infrared images, which models the intensity variationsas an exponential function.[19] In another approach,[17] namely generalized dynamic imaging model (GDIM), the hypothesis is that the intensity varies as product of a multiplier in previous amount of intensity, and an offset is added to it. In other words, the intensity variation is modeled with a multiplier and an offsetwhich is more accurate than the other works[2021] considering the intensity variations just as a multiplier or just an offset. In GDIM, initially, these quantities are unknown and then, will be computed associated with two more unknowns belonging to OF; the unknowns that are related to OF are velocity field components (u, v).

Also, there are some other researches accomplished in this context, which consider brightness inconstancy using different constraints and models.[22-24]

Finding a precise mathematical function as a model of intensity variation in echocardiography images is difficult due to the fact that these images have unique specifications and myocardial fibers have complicated motions in every cardiac cycle. These fibers have a non-uniform motion, in different directions and also different sizes in each direction.

As mentioned above, cyclic variations in echocardiography images are a function of the cardiac muscle motion and depend on the orientation of myocardial fibers toward ultrasonic wave.

None of the mentioned approaches has considered the specifications and variations of echocardiograms and they have not addressed the problem of intensity inconstancy in this category of images. Therefore, in this research, we motivated to propose a method which considers the variations in echocardiography images. To state the matter differently, in this approach, the intensity variation that occurs between two consecutive frames is calculated pixel-by-pixel in the image and afterwards, the effect of the variation is combined with OF calculations.

In this study, we investigate whether our proposed method improves the motion estimation of echocardiography images by OF techniques. In subsequent parts of this paper, in section 2, firstlythe OF method and its equations which are required for cardiac motion estimation are described; then our proposed method will be clarified. In section 3, validation results are expressed; these results first arise from a comparison between our method and two otherOF-based methods and then from a comparison between the proposed method and expert point of view as a golden standard. Finally, Section 4 is assigned to summarizing and conclusionof the paper and future works.

MATERIALS AND METHODS

The main idea of OF is based on the constancy of image intensity in consecutive frames. If I(x, y, t) is the intensity of pixelsat location expressed by (x, y) and time t, the basic assumption of OF can be declared by (1):[17]

I(x, y, t) = I(x + Δx, y + Δy, t + δt)      (1)

In this equation, Δx and Δy are determined as below:

(u, v) are the OF or velocity components in the horizontal and vertical directions at the point (x, y), and Δx, Δy are displacements at mentioned directions, respectively. Moreover, δt is the small time interval between two consecutive frames.

First-order Taylor expansion of the right hand of (1) will result in (2):[17]

In this equation, I are spatial derivatives of intensity in x, y directions, respectively, and I is the temporal derivative of intensity.

In (2), there are two unknowns (u, v) in one equation. Consequently, the solution of this equation is not possible with usual methods. Some methods are proposed to solve the (2) in an efficient way.[1725] In this paper, we use the method proposed by Lucas and Kanade[25] which estimates the motion locally assuming it to be invariant within a spatial window. They suggested that this equation could be solved by minimization methods like weighted least squares. The advantage of this method is that the OF calculations are not implemented in the entire image. In other words, calculations are implemented for small windows in the image. Therefore, the method will be computationally efficient.

The Weighted Least Squares equationis represented in (3):

In this equation, X = (x, y) and W(X) is a 2D window with dimension of L×L, including n pixels of X ∈ Ω that Ω is the domain containing the neighborhood of X. Minimizing this equation with respect to V =(u, v), (4) is obtained:[6]

A2AV=A2b      (4)

In this equation:

A = [▽I(x1),…,▽I(x)]      (5)

W = diag[W(x1),…,W(x)]      (6)

b = –(I(x1),…,I(xn))      (7)

The solution of (4) is expressed in (9):

V = [A]–1A

In all of these equations, the intensity variation influence is not considered.

In this study, a coarse-to-fine approach[1] is implemented. Utilizing this technique allows more accurate estimation especially for large motions. Number of the coarse-to-fine pyramid levels in our implementation is selected by comparison of efficiency between different choices. The optimum value obtained for the number of pyramid levels is 4. It means that a coarse-to-fine pyramid with four levels has the maximum accuracy and needs the minimum calculation time.

Each frame of the image is sub-sampled by the factor 1/2 at each level and a pyramid of images with four different resolutions is created. At the next step, the motion estimation algorithm, mentioned in previous equations, is accomplished starting at coarsest level (the level having less resolution) continuing to finest level (the level having the most resolution). The estimation of each level transfers to the next level as an initial estimation; the initial value in first level of the pyramid is set to zero.

The aim of this research is combining the intensity variation influence with mentioned equations to improve the motion estimation in echocardiography images.

In this study, the intensity variation is calculated locally for each pixel of the echocardiography image. Subsequently, the obtained values are replaced with the right hand of (2) that is zero. This means that in our calculations (2) is altered to (10), as below:

In this equation, ΔI is the intensity variation in the pixel located at (x, y) and at the time t. In fact, ΔI is calculated locally at each pixel of the image. The considerable point in (10) is that the actualvalues of Δx and Δy are unknown. Actually, (Δx,Δy) is the displacement vector or the motion field of each pixel which is in close relation with the principal unknown of the OF equation, namely the velocity vector (u, v). Consequently, we do not have (Δx,Δy) for computation of ΔI at the first step. As a result, an estimation of these unknowns is required. This estimation is done in the coarse-to-fine pyramid mentioned above. Firstly, for each frame of the sequence, the motion estimation algorithm based on the weighted least squares approach is applied to every level of the coarse-to-fine pyramid, from the coarsest to finest level, respectively. At each level, the frame is up-sampled with the factor 2 and the algorithm is applied again. This process is repeated while reaching to the pyramid apex, namely the finest level. The initial value for (Δx,Δy) at the first level is considered zero, and in the other levels, the initial value is obtained by doubling the estimated value in previous level.

In Figure 1, motion estimation steps in a four level coarse to fine pyramid is demonstrated. In this Figure, G00 and G10 are the initial sizes of two consecutive frames that in each level are sub sampled with the factor ½. G01, G02, G03 and G11, G12, G13 are the coarser levels of frame 1 and 2, respectively. In addition, Dx and Dy are the same (Δx,Δy) which are updated in each level of the pyramid.

Figure 1

Motion estimation steps in a four level coarse to fine pyramid

RESULTS

For validation of the proposed method, the algorithm was implemented on real echocardiography images. Data is acquired from this site:

“www.yale.edu/imaging/echo_atlas/contents/index.html” and the resolution of imaging is 3.5 MHz.

Our method was compared with two more OF-based methods, GDIM[17] and Lucas-Kanade.[25] All of the methods have comparable parameter settings. In the calculations related to our proposed method and Lucus-Kanade method, the optimum values obtained for the window W(X) dimensions were 17×17 pixels. It was obtained by calculation of the optimum value of mean square error (MSE) when a given displacement is applied to the image.

The criterion of comparison, MSE, is shown in (11):

In this equation, N is the number of points for which OF is calculated. t(k) is the set of true displacement values, c(k) is the set of calculated values, and e(k) is the error of each point.

We selected a number of landmarks on a frame of echocardiography sequence to draw a curve calculated by our method. These landmarks were selected on endocardial border of heart and were interpolated by B-spline curves. We compared it with expert defined curve by calculation of three distance metrics, mean, rms and Hausdorff, between two curves.

Motion components were estimated by our proposed method and two other methods, namely Lucas-Kanade[25] and GDIM.[17] The Lucas-Kanade method is one of the common OF techniques and the intensity variation is not considered in it. The OF calculations of the Lucas-Kanade method are implemented by (2). Nevertheless, in the GDIM method, the intensity variation is considered as a multiplier in first frame's intensity and an offset. In other words, the OF equationrelated to this method is written as (12), while in our proposed method, the intensity variation is considered locally and pixel-by-pixel:

In this equation, δm,δc respectively, are the multiplier and offset values. Moreover, u, v are the unknowns in this equation. Actually, these four parameters are calculated simultaneously using a least squares approach by (13):

The complete procedure of this calculation is expressed in the related reference.[17]

Validation Results for Applied Motions in Echocardiography Images

Practically, determination of a ground truth for heart motions is a challenging issue in validation of tracking methods. To evaluate the proposed method for echocardiography images, a normal heart echocardiography was used. Owing to the fact that the real value of displacements in different regions of an echocardiography is not available, three given displacements, 5,10, and 15 pixels, were applied on this image for validation purpose.

Our proposed method and two other methods (Lucas-Kanade and GDIM) were implemented to estimate the mentioned given displacements. Table 1 shows the MSEs related to displacements calculated by these algorithms. We observed that the MSE of our method was clearly less than two other methods.

Table 1

Mean square error for three different methods of optical flow calculation in an echocardiography image for three applied displacements

In addition of displacement error, we computed angular error for the proposed method and compared it with two other mentioned algorithms. We applied three different optional angles, 0°, 45°, 75°, on an echocardiography image and estimated them with three methods. Table 2 shows the MSE of estimated angles with three algorithms (proposed, Lucas-Kanade and GDIM). It is observed that the angular error related to our method is clearly less than two other methods. Therefore, the method will be accurate in estimation of the motion direction.

Table 2

Mean square error for three different methods of optical flow calculation in an echocardiography image for three applied angles

Finally, we implemented our proposed method on consecutive frames of a sequence of echocardiography images with a frame rate of 24 frames per second. This sequence was related to a normal heart. Figure 2 shows the OF calculated by our method between frames 12, 13 of this sequence which are related to expansion phase (diastole). In Figure 2, expansion of left ventricle can be observed by velocity arrows computed by our proposed OF estimation method. For more appropriate visualization of OF arrows we adjusted the background intensity in this Figure.

Figure 2

Frame 12 of an expanding heart motion estimation by the proposed method in an echocardiography image; optical flow is calculated between frames 12 and 13

To evaluate the proposed method for real heart motions in echocardiography images, we drew the endocardial border of the left ventricle calculated with the proposed method and compared it with the manual curve defined by an expert who was blind to automatic estimation result. For this purpose, at first, we defined a set of landmarks on the endocardial border and interpolated them by B-spline curves; secondly, we compared it with the curve defined by the expert. We calculated three distance metrics, namely mean, root mean square (rms) and Hausdorff distance, between automatic and manually defined curves. Mean distance is defined as (14):

In this equation, a and b are radial distances from a common central point and N is the number of points. Second metric distance, rms, is defined as (15):

The last and the most precise metric, Haussdorff distance, is defined as (16):

where d is the Euclidean distance between automatic and manual curves.

In Figure 3, panels (a) to (c) show the B-spline interpolated curve related to frame 13 of the echocardiographic sequence calculated by Lucas-Kanade, GDIM and the proposed method, respectively, and panel (d) shows the expert defined curve. Table 3 shows the results of the comparison of three automatically defined curves with the manually defined one. It is observed that the Haussdorff distance between the estimated curve defined by the proposed method and the manually defined curve, as a golden standard, is 4.81 pixels less than this distance for Lucas-Kanade and 2.28 pixels less than the GDIM method. Accordingly, we can say the proposed method is appropriately working when compared with the other methods.

Figure 3

B-spline interpolated curve related to frame 13 of the echocardiography sequence calculated by (a) Lucas-Kanade; (b) GDIM; (c) the proposed method; and (d) expert defined curve

Table 3

Distance metrics between three automatic and one expert defined curves

CONCLUSION

Cardiac motions in echocardiography images may lead to variation of intensity in different areas of the image. Intensity variations result in violation of the main assumption of the OF calculations. In this paper, we considered intensity variation in OF calculations and tried to compensate the effect of these variations to improve the estimations in echocardiography images. In our proposed method, the intensity variation was considered locally and pixel-by-pixel since heart has non uniform and complicated motions and definition of a mathematical function to model the intensity pattern is a problematic issue.

We evaluated our method by comparison with two other OF-based tracking methods, Lucas-Kanade and GDIM. Results show that precision of the proposed method was more than two other methods. In addition, we compared the proposed method with expert point of view and observed that the curve which was defined by our method was close to the curve which was defined by expert, when compared to two other methods.

This research can be a prelude to approach a specific model of intensity variation for echocardiography images. Investigating the way of intensity variation in different regions of echocardiography images may lead to obtaining such a model.

BIOGRAPHIES

Paria Torkashvand received the B.S degree in Electrical Engineering from Khaje Nasireddin Toosi University, Tehran, Iran, in 2007 and the M.S. degree in Medical Engineering from Iran University of Science and Technology, Tehran, Iran in 2010. Her research interests are Medical Image processing, specially, Echocardiography tracking.

E-mail: paria.torkashvand@gmail.com

Hamid Behnam received the B.S. degree in Electrical Engineering from Iran University of Science and Technology, Tehran, Iran, in 1988, the M.S. degree in Medical Engineering from Sharif University of Technology, Tehran, Iran, in 1992, the Ph.D. degree in Applied Electronics from Tokyo Institute of Technology, Tokyo, Japan in 1998. Since 1998 till 2004, he was a researcher at Iran Research Organization for Science and Technology and from 2004 he has been a faculty member at Iran University of Science and Technology, in Tehran, Iran. Currently, he is an Assistant Professor of Biomedical Engineering at the IUST. His research interests are Ultrasound in Medicine, Medical Image processing and Medical signal processing.

E-mail: behnam@iust.ac.ir

Zahra Alizadeh Sani earned a medical degree at Mashhad University of Medical Sciences, Mashhad, Iran in 1998. She was Resident of Cardio Vascular Diseases in Mashhad University of Medical Sciences, Mashhad, Iran in Oct. 2002-Sep. 2006 and completed fellowships in Echocardiography (Sep.2007-Mar.2009), at Iran University of Medical Sciences, and in Basic Cardiac MR (June. 2009), St. Marien Hospital Bonn Germany. Her research interests are Cardiovascular MRI, Physics, Pulse sequence, clinical application, Cardiovascular CT, Angiography, physics, technique, clinical application Echocardiography, physics, clinical application. Since 2009, she works as an assistant professor of Cardiovascular Imaging in Cardiovascular Imaging department, Shaheed Rajaei Cardiovascular and Medical center, Tehran University of Medical Sciences.

E-mail: d_zahra_alizadeh@yahoo.com