A Reconstruction Method of Blood Flow Velocity in Left Ventricle Using Color Flow Ultrasound.

Vortex flow imaging is a relatively new medical imaging method for the dynamic visualization of intracardiac blood flow, a potentially useful index of cardiac dysfunction. A reconstruction method is proposed here to quantify the distribution of blood flow velocity fields inside the left ventricle from color flow images compiled from ultrasound measurements. In this paper, a 2D incompressible Navier-Stokes equation with a mass source term is proposed to utilize the measurable color flow ultrasound data in a plane along with the moving boundary condition. The proposed model reflects out-of-plane blood flows on the imaging plane through the mass source term. The boundary conditions to solve the system of equations are derived from the dimensions of the ventricle extracted from 2D echocardiography data. The performance of the proposed method is evaluated numerically using synthetic flow data acquired from simulating left ventricle flows. The numerical simulations show the feasibility and potential usefulness of the proposed method of reconstructing the intracardiac flow fields. Of particular note is the finding that the mass source term in the proposed model improves the reconstruction performance.

1. Introduction

Vortex flow imaging has recently attracted much attention in the field of clinical cardiac assessment owing to reports of its feasibility for analyzing intraventricular vortex flows [1-3]. The vorticity of intraventricular blood flow describes a rotational flow pattern that offers possible clinical indices of cardiac functions such as sphericity, vortex depth, vortex length, and vortex pulsation correlation.

There are several methods to compute and visualize the velocity fields of blood flow inside the left ventricle (LV), with echo particle image velocimetry (E-PIV) being representative of the commonly used noninvasive methods [4]. It tracks the speckle patterns of blood flow to estimate blood motion within the imaging plane. Although it is generally unable to measure out-of-plane particle motion from 2D echocardiography data (called B-mode images), a recent study extending E-PIV to 3D volume data demonstrated the possibility of out-of-plane assessment [5]. However, E-PIV is not completely noninvasive because it requires the intravenous injection of a contrast agent to obtain images suitable for the speckle-tracking algorithm.

To develop less invasive techniques, methods to reconstruct blood flows from color flow images (also called C-mode images, color Doppler images, color Doppler data, or Doppler echocardiography) have been proposed. The color flow images reflect the projected velocity components in the direction of ultrasound beam propagation [6]. To compute the flow velocity from color flow images, Garcia et al. [7] assumed a 2D divergence-free condition on the velocity fields; they decomposed each 2D velocity vector into a radial component obtained from the color flow data and an unknown angular component, which was computed using their assumption of the 2D flow. Ohtsuki and Tanaka [8] also assumed 2D flows and recovered the 2D velocity fields from the color Doppler data using the concepts of stream function and streamline in 2D fluid flow. However, the assumption of a 2D divergence-free condition is an oversimplification that ignores out-of-plane flows.

Arigovindan et al. [9] proposed a velocity reconstruction method using color Doppler data acquired from beams in two different directions. To cope with the nonuniformly sampled data of multiple imaging planes, they used 2D B-spline on each of the velocity components to be estimated, and the unknown coefficients of the 2D B-spline were calculated from the measured color Doppler data using least squares. Similar to the 2D reconstruction, Gomez et al. [10] recovered 3D velocity fields from multiple registered color Doppler images using 3D B-spline and least squares. The registration of multiple imaging planes for the above two methods remains very challenging in a practical environment.

Recently, a new imaging modality (Doppler vortography) based on 2D color Doppler data was introduced by Mehregan et al. [11], who assumed that a vortex flow pattern has axisymmetric features in the neighborhood of its center. Their method employs a simple kernel filter designed to find the positions of axisymmetry in the 2D color Doppler images. The vortex flow was recovered using a color Doppler-variable vorticity function that directly computes vorticities from color Doppler values. However, the assumption of axisymmetry does not reflect detailed flow patterns, and it may lead to inaccurate vortex positions and vorticity values in patients with severe dysfunction where axisymmetry cannot be assumed at all.

In this paper, we propose a 2D Navier-Stokes model to reconstruct intraventricular flows using color flow images and LV boundaries extracted from echocardiography data. Although the use of the full 2D Navier-Stokes equations in this setting has already been proposed and evaluated for 2D flow field regularization [12], the originality of this work is the inclusion of a source-term to deal with the out-of-plane flow component. The proposed model considers both in-plane and out-of-plane blood flows for an imaging plane in apical long-axis three-chamber (A3CH) view. Particular attention is given to the appearance and disappearance of the out-of-plane components in the imaging plane, which is modeled as a mass source term of a source-sink distribution. Blood flows in the imaging domain are reconstructed through solving a system of equations, which include a 2D incompressible Navier-Stokes equation for the mass source term and the color flow data measurement equation describing the projected velocity component for the color flow data. The boundary conditions required to solve the system of equations are given by the LV borders extracted from echocardiography data.

The performance of the proposed method is evaluated numerically using synthetic flow data with LV motion. The proposed method is shown to be feasible and potentially valuable for reconstructing intracardiac flow fields.

2. Materials and Methods

Commonly used ultrasound systems can provide not only 2D echo images but also color flow images, which represent the scanline directional components of the velocity fields using the phase-shift estimated by a standard autocorrelation algorithm [13]. Our flow reconstruction method is to reconstruct the intraventricular flows using color flow images and the LV boundaries extracted from the echo images. In this section, we describe the overall outline of our flow reconstruction method using those ultrasound measurements based on two assumptions as follows:

the time difference between sequential color flow imaging frames is very small;

the echo and color flow images are acquired simultaneously and separately for the entire heart cycle.

2.1. Mathematical Model on 2D Imaging Plane

We focus on the dominant vortex flow appearing in the A3CH view, which passes through the apex and the mitral and aortic valves as shown in Figure 1(b), and mathematical model for blood flows inside the LV on the imaging plane of the A3CH view.

Figure 1

A sector scanning probe and apical long-axis views. The vector a in (a) means ultrasound beam propagation directions in sector-type scanning for cardiac application. (b) shows three (4-chamber, 3-chamber, and 2-chamber) apical imaging views for the cardiac scanning. They have the relationship of rotation by approximately 60 degrees toward each other. The imaging plane of the apical long-axis 3-chamber view passes through both of mitral and aortic valves.

Let D be a 2D imaging domain and Ω(t) the cross-section of the LV region in the A3CH view so that they satisfy Ω(t)⊆D⊆ℝ 2. Color flow data are practically measured in the imaging plane D. Figure 1(a) describes the scanline directional vectors a : = (a 1, a 2) for the 2D imaging plane D as an example of sector scanning. As shown in Figure 2, Γ(t) and Γ(t) denote the mitral and aortic valves, respectively. The parameter t denotes the LV region and valves as time-varying boundaries. The superscripts I and O of Γ(t) and Γ(t) stand for inlet and outlet valves, respectively. The 3D coordinate system takes the xy-plane to contain Ω(t) and the z-axis to be normal to this plane. Figure 2 describes a basic LV structure and the diastolic/systolic motion of its wall observed in the A3CH view.

Figure 2

Description on the LV region Ω(t) in the A3CH view. They show the diastolic and systolic motion of  Ω(t), respectively. Here, Γ(t) and Γ(t) refer to the mitral and aortic valves (as inlet and outlet valves), respectively.

Let c(x, t) be the measured color flow data and (u(x, t), v(x, t)) the velocity fields of flow at the position x ∈ D and time t, respectively. Then the color flow data c(x, t) can be expressed as the inner product of the scanline vector and the velocity vector:To recover (u, v) from knowledge of c on the imaging plane D, we propose the following 2D Navier-Stokes model:where ρ = 1050 kg/m3 and μ = 0.00316 Pa·s are the density and viscosity of the blood flow, respectively [14]. This model (2) is equivalent to a 2D incompressible flow having a source-sink distribution s(x, t) [15]. In fact, the 3D blood flows appear or disappear in the imaging plane. Therefore, we reconstruct the 2D velocity fields in the imaging plane by solving (1) and (2) with LV boundary conditions. We impose the following boundary conditions:where u wall and v wall are velocity components computed by the motion of  ∂Ω(t) and n is the unit outward normal vector to ∂Ω(t). Here, ∂Ω(t) is the boundary of  Ω(t), and its motion can be extracted from echocardiography data as described in the above subsection. Further mathematical explanation of (1) and (2) is given in the Appendix A.

2.2. Reconstruction Algorithm

For numerical implementation, we write the system of   (1) and (2) as the following linear second-order system:We discretize the 2D imaging region Ω into the mesh grid elements with M × N nodes and apply the standard finite difference method for the linear equation (4). Then we obtain the discretized linear system of the following form:where D , D , and L are the x-derivative, y-derivative, and Laplace operators of finite difference, respectively; the superscript (n) denotes the nth time-step. Note that the measured data c on the right-hand side are the values at the (n + 1)th step, not the (n)th step.

Let I and O be the N × N identity and zero matrices, respectively. Also set D and L to and let 𝔸 1, 𝔸 2, 𝔻 , 𝔻 , and 𝕃 be the (M × N)×(M × N) matrices defined by The above discretized linear system (5) can then be written in the form of 𝕂 U = F, where

For notational simplicity, we drop the superscripts (n) and (n + 1). Then, in the right-hand side, u, v, p, s, c, k 1, and k 2 are the (M × N) × 1 vectors defined as follows:

However, the coefficient matrix in (8) is nonsymmetric, sparse, and large scale. Moreover, it has a bad condition number and is almost singular. Hence, we use the least-squares approach with Tikhonov regularization to solve the minimization problem with a regularization term of the form:By setting Γ = λ𝕀 (4×, the minimization problem is written asand its solution is given byWe perform the reconstruction algorithm using the heuristically selected Tikhonov regularization parameter λ shown in the reconstruction algorithm (11).

2.3. Numerical Simulations with Moving LV Boundaries

We obtain synthetic flow data inside a virtual moving simplified LV wall. To obtain the synthetic intraventricular flows, we construct a 3D moving LV region of the stroke volume of about 20 mL. As shown in Figure 3(a), the time-dependent LV volume ranges from 123 mL to 143 mL during the entire cycle of 0.5 s. Figure 3(b) shows a LV shape model corresponding to the volume of 133 mL at the preset pressure state. We then perform a numerical simulation of the forward problem of the Navier-Stokes equation inside the 3D moving LV for the beat cycle. For the forward simulations, the “fluid-structure interaction model” of the COMSOL software is used, and 3D intraventricular velocity fields are computed as a solution of the forward problem. A no-slip boundary condition is used for computing the velocity fields. We assume that the blood flow during the filling of the LV is ejected in the normal direction to the inlet valve surface Γ(t), that the flow velocity is uniform on the entire inlet boundary, and that the LV volume change is equal to the total amount of its net inflow: where n is the normal vector to the mitral valve (inlet valve) boundary denoted by Γ(t). The boundary portion containing the mitral and aortic valves was fixed not to be moved. Neumann boundary conditions of velocity fields were given by setting the zero normal derivatives at the outlet Γ and inlet Γ valves for the LV diastole and systole, respectively. For simplicity, pressure is assumed to be a constant at the inlet, while the Neumann condition for the pressure is used at the outlet. (The work of [16, 17] will also help readers create synthetic flow data.)

Figure 3

LV model for numerical experiments. (a) LV volume curve according to diastolic and systolic motions, (b) 3D LV shape model at the end systole, and (c) 2D imaging plane of A3CH view positioned in the constructed 3D LV model. The dashed arrows means the scanning lines transmitted from a linear array transducer for 2D imaging.

We project the 3D synthetic velocity fields on the imaging plane of A3CH view, which is set to be located in the constructed 3D LV model (see Figure 3(c)). The projected 2D velocity fields are used as reference data to evaluate the proposed 2D velocity reconstruction algorithm. Representative cases of these synthetic 2D flows are depicted in Figure 4, which show the time-varying dominant vortex patterns inside the virtual moving LV wall. While the LV model is relaxing, flows enter into the LV through the inlet valve. The main stream impinges against the LV wall, and a large vortex is simultaneously formed near the center. The large vortex moves upward during the contraction process and weakens at the end of the systole cycle.

Figure 4

Synthetic intraventricular velocity fields projected on the imaging plane. (a) and (b) show the velocity fields for the diastole and systole process of LV, respectively.

Figure 5 shows the vorticity fields obtained by taking the curl operator to the synthetic flow fields. At the early stage of relaxation, small vortices are formed on the both sides of the main incoming stream. The higher vorticity near the wall is related to the fact that the high momentum incoming flow results in increasing the wall shear stress when the incoming flow impinges to the wall; we think. While the vorticity gradually grows during the expanding process, it eventually shrinks and weakens during the process of contraction.

Figure 5

Vorticities corresponding to the synthetic 2D intraventricular velocity fields. (a) and (b) show the vorticity changes according to the diastole and systole process of LV, respectively.

The inner product of the synthetic 2D velocity fields and the scanline directional unit vectors gives the scanline directional velocity components, which can be regarded as the color flow data of the intraventricular flow on the A3CH view measured by a real ultrasound system. The one-directional velocity component data are used as the input data for our proposed method. Figure 6 illustrates the scanline directional velocity components. The figures in the second column corresponding to the end-diastole and end-systole reflect the overall similar flow patterns that are comparable to real color patterns (see Figures 12(c) and 12(d) in Appendix B).

Figure 6

One-directional velocity component data. (a) and (b) show the scanline directional components of the generated velocity fields during the diastole and systole process of LV, respectively.

Figure 12

Real data acquisition using LV phantom and border extraction. (a) and (b) are echo images in ventricular shape at end-diastole and end-systole phases, respectively. (c) and (d) are color flow images corresponding to (a) and (b), respectively.

2.4. Choice of the Parameter

In this study, all the experiments were performed by setting the parameter . To optimize the value of this parameter, we investigated the errors of the reconstructed velocity fields with respect to the change of λ 2. We adjusted the value of λ 2 to each power of 10 from 10−3 to 101. Figure 7 shows that the L -norm errors between the reconstructed velocity fields and the reference data have formed a convex plot with the minimized L -norm error attained at λ 2 = 0.1.

Figure 7

Error change depending on the parameter λ. (a) L -norm errors between u recon and u exact, (b) between v recon and v exact, and (c) between (u recon, v recon) and (u exact, v exact) with respect to λ 2. Here, (u recon, v recon) and (u exact, v exact) stand for the reconstructed and the reference velocity fields, respectively.

3. Results

We demonstrated the feasibility of the proposed method by showing the L 2-norm errors between the synthetic flow and the reconstructed flow fields. We reconstructed the velocity fields from the synthetic scanline directional velocity components by repeatedly solving the minimization problem (11) five times for the entire beat cycle. Figure 8 shows the reconstruction results of velocity fields obtained from the input data of the one-directional velocity components. As in the forward numerical simulations, the incoming stream impinged to the lower right wall. Overall, the reconstructed flows showed very similar patterns to the reference flows shown in Figure 4. However, the main vortex near the center was shifted upward relative to the reference vector field. This shifting phenomenon was also continued in the remaining steps.

Figure 8

Reconstruction results. The reference vector field (a) and the reconstructed 2D vector fields by the proposed method (b).

Figure 9 illustrates the performance of the proposed model. Part (a) compares the pointwise error magnitude images of the reconstructed velocity fields with the reference data. While the reconstructed u-component errors were distributed in the range of  −0.39~0.25 m/s, the errors of the v-components showed larger differences of  −0.24~0.27 m/s for the entire cycle. Figure 9(b) shows the distribution of the mass source term s and reflects the out-of-plane components of the flows. The ratio s/ρ ranged from −151.1 s−1 to 82.5 s−1. From these contributions of s, the averaged pointwise errors for u- and v-components were 0.06 m/s and 0.02 m/s, respectively. The averaged pointwise error-magnitude was 0.065 m/s. We also compared the reconstruction results between the 2D Navier-Stokes models with and without the mass source term s. When s was used, the largest pointwise errors for the u- and v-components were 0.39 m/s and 0.27 m/s, respectively; reconstruction results not considering s showed larger pointwise errors of 0.71 m/s and 0.78 m/s for these components, respectively, for the entire cycle. These results support that the proposed model improves reconstruction performance.

Figure 9

Errors of the reconstruction results. (a) 2D plots for error-magnitude of reconstructed (u, v) relative to the reference velocity fields and (b) the mass source term s.

For further quantitative comparison, we computed L 2-norm errors of each velocity field on the given whole region at each 1/10 time step of the whole cycle. Table 1 shows the relatively large velocity error (≈47%) for the end of the contraction process compared with the expanding process. While the pointwise errors of velocity may be larger than its L 2-norm errors at some local regions, those pointwise errors at the local regions did not affect the major features of the vortex flows. This implies that the proposed model is reasonably accurate for estimating the unsteady features of the blood flows when the reconstruction algorithm of flow is performed repeatedly over several cycles of diastole and systole processes.

Table 1

Comparison of L 2-norm errors for the velocity fields recovered from the synthetic and real data. The third and fourth columns represent normalized L 2-norm errors for pointwise and global energy estimates of velocity. Here, , Δv = v recon − v exact, and Δk/k exact = |‖v exact‖2 − ‖v recon‖2|/‖v exact‖2.

t/T	‖v exact‖2	‖v recon‖2	‖Δv‖2/‖v exact‖2	Δk/k exact
0.1	0.012	0.011	28.2%	11.1%
0.2	0.026	0.023	17.6%	8.9%
0.3	0.033	0.029	21.2%	9.9%
0.4	0.028	0.024	28.2%	12.0%
0.5	0.018	0.015	39.6%	17.4%
0.6	0.015	0.013	39.8%	16.8%
0.7	0.015	0.013	36.6%	13.4%
0.8	0.014	0.012	36.3%	12.0%
0.9	0.011	0.010	39.3%	13.0%
1.0	0.010	0.007	47.5%	18.6%

Although the proposed reconstruction method produced nonsmooth distributions of vorticity inside the LV region (as shown in Figure 10), the evolutions of the large vortical structure inside the LV region were clearly observed, and the vorticities overall showed very similar patterns to the reference results described in Figure 5. The pointwise and global errors for the vortex fields are listed in Table 2. These nonsmooth vorticity distributions may be induced by uncertainties of the boundary geometry and the regularization term in the minimization problem (11).

Figure 10

Reconstruction results of vorticities. These were obtained by taking the 2D curl operator to the reconstructed velocity fields.

Table 2

Comparison of L 2-norm errors for the corresponding vortex fields. The third and fourth columns represent normalized L 2-norm errors for pointwise and global energy estimates of vorticity. Here, Δω = ω recon − ω exact and Δk /k exact = |‖ω exact‖2 − ‖ω recon‖2|/‖ω exact‖2.

t/T	‖ω exact‖2	‖ω recon‖2	‖Δω‖2/‖ω exact‖2	Δk ω/k ω exact
0.1	1.09	1.19	36.8%	8.7%
0.2	2.35	2.48	24.7%	5.3%
0.3	2.82	2.90	21.7%	2.8%
0.4	2.31	2.34	21.2%	1.4%
0.5	1.67	1.67	32.1%	0.4%
0.6	1.38	1.40	40.6%	1.7%
0.7	1.33	1.40	46.7%	5.2%
0.8	1.20	1.30	51.4%	7.8%
0.9	0.95	1.02	52.0%	6.7%
1.0	0.79	0.82	47.8%	4.2%

4. Discussion and Conclusions

We propose here a new 2D Navier-Stokes model that reconstructs intraventricular flows using color flow data and LV boundaries extracted from echocardiography data. The proposed model considers both in-plane and out-of-plane blood flows on the imaging plane of an apical long-axis three-chamber view. The out-of-plane components moving out of the imaging plane were modeled as the mass source term of a source-sink distribution. We reconstructed blood flows in the imaging domain by solving a system of equations, which includes the 2D incompressible Navier-Stokes equation of the mass source term and a color flow measurement equation describing the one-directional velocity component of the color flow data. The boundary conditions that are required to solve the system of equations are given by the LV borders extracted from echocardiography data. To evaluate the proposed method, numerical experiments were performed on synthetic flow data following a virtual LV motion. The results showed that the proposed method is feasible and potentially valuable for reconstructing intracardiac flow fields.

The numerical experiments used an imaging domain of size 64 × 64 pixels, because the computational costs of the proposed reconstruction method were large. We divided the whole cycle into 1000 time-steps for stable numerical implementation and repeatedly performed the proposed algorithm five times for the entire heartbeat cycle. We observed that the solutions obtained after two repeats were very similar to each other. A detailed mathematical description of the convergence of the solution will be given in our next work.

All the experiments were conducted using MATLAB 7.10.0, and the computational time for each step was about 12 s using an Intel i7-4702MQ Quadcore CPU running at 2.0 GHz with 8 GB of RAM. The speed of the proposed method needs to be improved for its practical application, and a study to reduce the processing time is under way.

To reconstruct 2D flow patterns using the color flow images requires an experimental ultrasound system that can acquire echo and color data independently; however, we had only a commercial ultrasound system that cannot support data acquisition for the experiments. The frame-rate limitation of the color flow data acquired by the ultrasound imaging system should also be overcome. The color flow images from the ultrasound imaging system are obtained by applying the autocorrelation algorithm [13] after repeatedly transmitting and receiving the ultrasound beam along the same scanline; this then causes the low frame-rate of color flow data acquisition. The proposed method assumed that the time difference between the sequential color flow images is very small and that each color value reflects the exact velocity components in the direction of ultrasound beam propagation. For the proposed method to be practically applicable, the frame-rate needs to be improved. The development of an ultrafast imaging system is ongoing. Such a system is expected to have a frame-rate higher than 1000 fps using plane wave beam-forming; it would acquire echo and color flow image data separately and support the performance of the proposed method in real application. Some recent works by other researchers have examined ultrafast flow imaging [18-20]. A review of these references will help us to perform our future work on an ultrafast imaging system. In Appendix B, we describe an in vitro phantom experimental setup for empirical data (color flow ultrasound and LV boundary data) to demonstrate the feasibility of our reconstruction method.

To advance the proposed algorithm, we are studying on mitral valve tracking. Based on reports [21, 22] that the motion of the mitral valve affects the vorticities of the intraventricular flows, studies on the boundary conditions containing the valve motion are under way.

This work is the first that describes overall our proposed method. More detailed mathematical analysis and further validation tests may be required to verify its scientific validity and practical applicability. We will perform several follow-up works to achieve this.