Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 May 4.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2009; 12(Pt 2): 331–338.
PMCID: PMC2863152

Incompressible Cardiac Motion Estimation of the Left Ventricle Using Tagged MR Images


Interpolation from sparse imaging data is typically required to achieve dense, three-dimensional quantification of left ventricular function. Although the heart muscle is known to be incompressible, this fact is ignored by most previous approaches that address this problem. In this paper, we present a method to reconstruct a dense representation of the three-dimensional, incompressible deformation of the left ventricle from tagged MR images acquired in both short-axis and long axis orientations. The approach applies a smoothing, divergence-free, vector spline to interpolate velocity fields at intermediate discrete times such that the collection of velocity fields integrate over time to match the observed displacement components. Through this process, the method yields a dense estimate of a displacement field that matches our observations and also corresponds to an incompressible motion.

1 Introduction

To measure regional function in the heart, magnetic resonance tagging [1] generates images that can be processed to find two-dimensional in-plane motion at every pixel location, and when short-axis and long-axis images are combined, three-dimensional (3D) motion of each point in the left ventricle (LV) can be inferred through interpolation of the sparse imaging data. Although dense methods for directly imaging 3D myocardial motion have been developed, they require too much time for routine acquisition of dense, three-dimensional myocardial motion in scientific research or clinical medicine. Therefore, interpolation methods are likely to be required in practice and may well be the critical element in promoting routine imaging of dense, 3D myocardial function in the heart. It is widely accepted that the volume change of myocardium during the cardiac cycle is no more than 4% [2]. Since materials that are incompressible undergo deformations that preserve volumes at all scales and have divergence-free velocity fields, it is natural to assume that one can improve interpolation by exploiting this constraint. Song et al. [3] first applied this property in building the 3D velocity of the heart from cine CT images. Denney et al. [4] directly applied the divergence-free constraint to reconstruct the 3D displacement field of the LV in an estimation theoretic approach. Recently, Bistoquet et al. [5] constructed nearly incompressible cardiac motion field from non-tagged MR images using a vector spline with a divergence-free matrix-valued function. There is a key problem with these approaches, however. Because the temporal resolution of the image sequences, the deformation between two neighboring time frames may be large. A velocity field that is approximated as the displacement field divided by the time interval is not theoretically predicted to be divergence-free. When this fact is ignored and the underlying field is interpolated in a divergence-free fashion this can lead to considerable errors when reconstructing motion fields in a time sequence since the errors in earlier time frames propagate to later time frames. In [5], this error was reduced by interpolating both forwards and backwards in time and then computing a weighted average of these solutions. However, solutions generated this way do not yield motions that have divergence-free velocity fields or correspond to incompressible motions. In this paper, we present a new approach to reconstruct a 3D, dense, incompressible deformation field in the LV of the heart from tagged MR images based on divergence-free vector spline with incomplete data samples. A key novelty of our approach is that, instead of computing divergence-free displacement field, we seek a sequence of divergence-free velocity fields from which the final displacement field is computed by integration. We also adopt a multi-resolution strategy and adaptive smoothing to reduce the computation and improve the accuracy. Our method was validated using both numerical simulation and in vivo cardiac experiments.

2 Background

2.1 Smoothing Divergence-Free Vector Spline

Given N points in space xn = [xn, yn, zn]T , n = 1,…,N, and vector-valued observations vn, n = 1,…, N, at these points, vector splines (VS) [6] interpolates a smooth vector field over the whole space. Specifically, the smoothing VS finds a vector field v(x) that minimizes

C(v)=ρJα,β(v(x))+1Nn=1Nv(xn)vn2, with


where ρ is the smoothing parameter, α and β are the weighting coefficients, div yields the divergence of a vector field, and rot yields the curl. It has been shown that (1) has the closed form solution [6] v(x)=n=1NK(xxn)·cn+p(x), where cn are the unknown vectorial coefficients and K(x) is the matrix-valued kernel function given by


where I is the identity matrix, h(x) = ‖x2k+1 is the solution to Δk+1 h(x) = δ(x) and Δ is the Laplace operator. p(x) is a polynomial of order k. The coefficients in the smoothing VS can be solved using known vector values at the sample points.

As a special case, the VS can be used to interpolate the divergence-free vector field by constraining the vector field to be divergence-free, i.e., divv(x) = 0. The divergence-free vector spline (DFVS) solution is similar to that of VS except that the kernel matrix becomes KDFVS(x) = [ΔI[nabla][nabla]T]h(x), and p(x) is also constrained to be divergence-free.

2.2 Smoothing VS from Incomplete Samples

In some applications—such as the displacement computed from tagged MR images—only selected components of the vector field are observed at the sample points. This incomplete sample data can be written as: {xn, ln, wn} for n = 1, 2,…,N, where ln is a normal vector representing a projection direction, and wn = ln · v(xn) is the projection of v(xn) on ln. The minimization problem of a smoothing VS given incomplete data can be expressed as


Arigovindan [7] showed that the solution to this problem is


where the coefficients cn are scalars, and K(x) and p(x) are the same as in VS. By replacing K with KDFVS, Eqn. (5) describes the solution to smoothing DFVS with incomplete samples.

3 Method

3.1 Experiments and Data Processing

We acquired CSPAMM cardiac image sequences using a breath-hold scenario on a Phillips 3T Achieva MRI scanner (Philips Medical Systems, Best, NL). An approved IRB protocol was used and informed consent was obtained. Both short axis (SA) and radial long axis (LA) images were acquired on a healthy subject. Two sets of images with orthogonal tag directions were acquired separately on each SA slice. The LA images were tagged in a direction perpendicular to the SA image planes. The imaging parameters were: tag spacing = 12 mm, pixel spacing = 1.25 mm, temporal resolution = 30 msec, time frames=20.We acquired twelve SA image slices with a 4 mm slice separation, and six LA image slices. The SA slices were divided into two interleaved groups (even and odd slice numbers) so that the slice separation within each group is 8 mm. The first group of six SA slices and all six LA slices were used to reconstruct a 3D, dense, incompressible displacement field of LV. The slices in the second group were used for validation. The relative locations of the slices used in the motion reconstruction are illustrated in Fig. 1.

Fig. 1
A tagged (a) SA and (b) LA image. For visualization purposes, the SA image shown is the product of the separately acquired horizontal and vertical tagged images. The overlaying lines depict the geometry of the (a) LA and (b) SA images.

All of the images were first processed using the harmonic phase (HARP) method [8] to yield sequences of HARP images. Let us define the time that the tags are just applied and have not deformed as the reference time t0. At t0, the tagging phase [var phi] is a linear function of the point’s coordinate x, and wrapped to the range [−π, π), i.e., [var phi](x, t0) = W(kx · l + [var phi]0) , where k is the known tagging frequency, l is a unit vector representing the tagging orientation, [var phi]0 is an unknown phase offset, and W is a phase wrapping operator. By assuming the tags do not deform much at the first time frame, [var phi]0 can be estimated from the HARP images at the first time frame [9]. Therefore, the tissue points at each time frame can be tracked back to t0 using HARP tracking [8]. For tagging direction l and a 3D spatial point xj imaged at t, if it comes from Xj = xj(t) − u(xj (t), t) at t0, then HARP tracking computes the projection of its displacement u(x, t) onto l as:


Because the SA images are acquired with two tag orientations, two projections of the displacement of each point are computed. For tissue points in the LA images, only one projection is computed. Therefore except for points at the intersections of LA and SA image planes, only partial knowledge of the displacement is available for any other pixel on the observed images. (Of course, no observations are available at 3D points that do not lie on an observed image plane.)

Intermediate image frames are not used to assist in tracking later frames because the observed tissues are not the same, primarily due to through-plane motion; thus the Lagrangian framework that is used to carry out incompressible interpolation (see next section) cannot take advantage of these observations.

3.2 3D Incompressible Displacement Field Reconstruction

The myocardium can be considered incompressible because it is composed mainly of water. For an incompressible elastic body subjecting to a deformation x = X+ u(X), the Jacobian determinant of the deformation satisfies det(I + [nabla]X u(X)) = 1 for any material point X, where [nabla]X is the material gradient operator. As well, the spatial velocity field v(x) giving rise to such a deformation must be divergence free, i.e., divv(x) = 0. Based on this physical property of the heart, we propose an approach based on DFVS to reconstruct the displacement field of the LV of the heart. HARP tracking provides N incomplete and non-uniform data samples {xn, ln, wn} at any time frame. Our goal is to reconstruct the 3D, incompressible displacement field u(x) such that ln · u(xn) = wn.

The deformation is calculated through the integration of the velocity. Let us define the integration variable as s, which takes on values in the interval [0, 1], and let the velocity of x at s be v(x(s), s). We define

w(x,s)=0sv(x(τ),τ)dτ, and v(x(τ),τ)=dw(x,τ)dτ,

with x(s) = x + w(x, s). Then the displacement field u(x) = w(x, 1). The problem can be reduced to a finite-dimensional problem by dividing the integration into discrete steps, i.e., sm = mδ for m = 0, 1,…,M with δ = 1/M. When M is reasonably large, the velocity is assumed constant within each interval, so u(x)=w(x,1)=δm=0M1v(x(sm),sm). Note v is not the true myocardial velocity, but rather a computational tool for the estimation of the displacement.

We use DFVS from incomplete data samples to interpolate the divergence-free velocity fields separately over each interval. The velocity fields are computed sequentially starting from s0 = 0 through sM = 1. Let us denote rn(sm) = ln · v(xn(sm), sm) for any step sm, and the data samples at sm are written as {xn(sm), ln, rn(sm)} for n = 1,…,N. The velocity at any sample point xn(sm) at sm is approximated by taking the first order expansion


so that


With the N data samples, the continuous velocity field v(x, sm) is interpolated with smoothing DFVS using Eqns. (4) and (5).

From Taylor’s expansion, the first order approximation of the velocity is accurate up to the order (1 − δm)2. Therefore, it is less accurate at smaller s and more smoothing is required at earlier steps. So the smoothing parameter ρ should be chosen to be large at small s and grow smaller as s approaches 1. At sM−1, ρ should be set to 0 so that the final displacement w(xn, 1) matches the original data samples exactly—i.e., ln · u(xn) = ln · w(xn, 1) = wn for n = 1,…,N. In practice, we choose the smoothing parameter as ρm=Mm1M1ρ0, where ρ0 is determined empirically.

To reduce computation time, a multi-resolution scheme is adopted. For smaller m, the sample points are downsampled so that only a subset of the samples is used in the interpolation. Since the computation time is dominated by solving the coefficients from Eqn. (5) with complexity O(N3), the multi-resolution scheme can greatly reduce the computation while not affecting the accuracy of the displacement field reconstruction.

The algorithm can be summarized as follows.

Algorithm 1: Dense 3D incompressible displacement field reconstruction

Given the data samples {xn, ln, wn}, and wn = ln · u(xn) for n = 1,…,N, and a dense 3D grid of data points yk for k = 1,…,K of which the 3D displacement vectors are to be computed, carry out the following steps:

  1. Initialize ρ0, M, yk (0) = yk, w(yk, 0) = 0, and w(xn, 0) = 0 for all k and n.
  2. for m = 0 to M − 1
    1. set sm = mδ, ρm = ρ0(Mm − 1)/(M − 1);
    2. downsample the data points x(sm) if needed;
    3. compute rn (sm) using Eqn. (9) for all sample points;
    4. compute the interpolating coefficients with samples {xn(sm), ln, rn(sm)} and ρ = ρm;
    5. compute the velocities v(yk (sm), sm) and v(xn (sm), sm) using Eqn. (5);
    6. set w(xn, sm+1) = w(xn, sm) + δv(xn (sm), sm), w(yk, sm+1) = w(yk, sm) + δv(yk (sm), sm), xn (sm+1) = xn(sm) + w(xn, sm+1), and yk(sm+1) = yk(sm) + w(yk, sm+1);
  3. Set u(xn) = w(xn, sM) and u(yk) = w(yk, sM), and the algorithm ends.

4 Results

4.1 2D Numerical Simulation

A 2D incompressible vector field was simulated using second order polynomials:


on an image with ranges x [set membership] [−50, 50] and y [set membership] [−50, 50]. This vector field is incompressible because det(I + Δu) = 1. The two displacement components are shown in Figs. 2(a) and (b). We picked a grid of sample points, shown in Fig. 2(c), and assumed the vectors on these sample points were known. To mimic the tagged image acquisition, these points were distributed densely in one dimension and sparsely in the other. Our method was then applied to reconstruct the 2D motion field in the whole image with M = 20 and ρ = 0.1.

Fig. 2
(a) The x and (b) y components of the simulated motion. The (c) sample points used in the interpolation. The Jacobian determinant from the motion field constructed using (d) our method and (e) one step divergence-free interpolation.

The reconstructed motion field was then compared with the closed-form solution. In this simulation the mean interpolation error was 0.064 pixel, and the standard deviation was 0.036 pixel. Fig. 2(d) shows the Jacobian determinant of the deformation. The average error of Jacobian determinant of the reconstructed deformation from 1 was 0.0046. For comparison, we also reconstructed the motion field using the direct divergence-free interpolation approach of Bistoquet et al. [5]. This approach yielded a mean interpolation error of 0.95 pixel and a mean error of Jacobian determinant (Fig. 2(e)) of 0.099.

4.2 Cardiac Motion Experiments

We applied Algorithm 1 to the cardiac images shown in Fig. 1 using a smoothing parameter ρ0 = 0.1 and M = 20 integration steps. The reconstructed 3D displacement field in the LV regions of three SA slices are shown in two views in Fig. 3 at different times. The reconstructed displacements of points on the LV in the validation slices were compared with the 2D displacement projection computed using HARP. Over all time frames, the mean displacement error was 0.446 mm, and the standard deviation was 0.310 mm. For comparison, we also computed the displacement field using Bistoquet’s approach, i.e., direct DFVS interpolation with backward-forward averaging [10]. The mean displacement error was 0.803 mm, and the standard deviation was 0.652 mm. Figs. 4(a) and (b) show the displacement error maps on the 5th validation slice at time frame 10 of our method and Bistoquet’s method, respectively. We also compared the incompressibility of the reconstructed motion fields. At time frame 10 when the heart deforms the most, the average absolute difference between the Jacobian determinant and unity of our approach was 0.030. The difference was mainly caused by both spatial and temporal discretization. The average absolute difference of Bistoquet’s method was 0.067. Figs. 4(c) and (d) show the Jacobian determinants at the 5th validation slice.

Fig. 3
The 3D displacement field illustrated using three SA slices. From top to bottom: two different views; From left to right: the displacement field at different time frames.
Fig. 4
The displacement error map on one slice using (a) our method and (b) Bistoquet’s method. The Jacobian determinant of the deformation field computed on the same slice from (c) our method and (d) Bistoquet’s method.

5 Conclusion

We presented an approach to reconstruct a 3D, dense, incompressible displacement field of the left ventricle of the heart using tagged MR images. Our method uses a divergence-free vector spline on incomplete and non-uniform sample data to interpolate the velocity fields at discrete integration steps, and the displacement field is achieved by integrating these velocity fields. Our method was validated with both numerical simulation and in vivo cardiac experiment.


*This work was supported by NIH grant R01HL047405.


1. Axel L, Dougherty L. MR imaging of motion with spatial modulation of magnetization. Radiology. 1989;171:841–845. [PubMed]
2. Yin F, Chan C, Judd R. Compressibility of perfused passive myocardium. Amer. J. Physiol.-Heart Circ. Physiol. 1996;8:1864–1870. [PubMed]
3. Song SM, Leahy RM. Computation of 3-D velocity fields from 3-d cine CT images of a human heart. IEEE Trans. Med. Imag. 1991;10(3):295–306. [PubMed]
4. Denney T, Prince JL. Reconstruction of 3D left ventricular motion from planar tagged cardiac MR images: an estimation theoretic approach. IEEE Trans. Med. Imag. 1995;14(4):625–635. [PubMed]
5. Bistoquet A, Oshinski J, Skrinjar O. Myocardial deformation recovery from cine MRI using a nearly incompressible biventricular model. Med. Imag. Anal. 2008;12:69–85. [PubMed]
6. Amodei L, Benbourhim MN. A vector spline approximation. J. Approx. Theo. 1991;67(1):51–79.
7. Arigovindan M. PhD thesis. École Polytechnique Federale de Lausanne; 2005. Variational reconstruction of vector and scalar images from non-uniform samples.
8. Osman NF, Kerwin WS, McVeigh ER, Prince JL. Cardiac motion tracking using CINE harmonic phase (HARP) magnetic resonance imaging. Magn. Reson. Med. 1999;42:1048–1060. [PMC free article] [PubMed]
9. Tecelao SR, Zwanenburg JJ, Kuijer JP, Marcus JT. Extended harmonic phase tracking of myocardial motion: improved coverage of myocardium and its effect on strain results. J. Magn. Reson. Imag. 2006;23(5):682–690. [PubMed]
10. Bistoquet A, Oshinski J, Skrinjar O. Left ventricle deformation recovery from cine MRI using an incompressible model. IEEE Trans. Med. Imag. 2007;26(9):1136–1153. [PubMed]