Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2863152

Formats

Article sections

Authors

Related links

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

NIHMSID: NIHMS158348

See other articles in PMC that cite the published article.

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.

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.

Given *N* points in space x* _{n}* = [

$$C(\mathbf{v})=\rho {J}_{\alpha ,\beta}(\mathbf{v}(\mathbf{x}))+\frac{1}{N}{\displaystyle \sum _{n=1}^{N}\Vert \mathbf{v}({\mathbf{x}}_{n})-{\mathbf{v}}_{n}{\Vert}^{2},\text{\hspace{1em}with}}$$

(1)

$${J}_{\alpha ,\beta}(\mathbf{v})={\displaystyle \int [\alpha \Vert {\nabla}^{k}(\text{div}\mathbf{v}(\mathbf{x})){\Vert}^{2}+\beta {\displaystyle \sum _{i=1}^{3}\Vert {\nabla}^{k}{(\text{rot}\mathbf{v}(\mathbf{x}))}_{i}{\Vert}^{2}}]d\mathbf{x}}$$

(2)

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]
$\mathbf{v}(\mathbf{x})={\displaystyle {\sum}_{n=1}^{N}\mathbf{K}(\mathbf{x}-{\mathbf{x}}_{n})\xb7{\mathbf{c}}_{n}+\mathbf{p}(\mathbf{x})}$, where **c*** _{n}* are the unknown vectorial coefficients and

$${\mathbf{K}}_{\text{vs}}(\mathbf{x})=[\frac{1}{\beta}\mathrm{\Delta}\mathbf{I}+(\frac{1}{\alpha}-\frac{1}{\beta})\nabla {\nabla}^{T}]h(\mathbf{x})$$

(3)

where *I* is the identity matrix, *h*(**x**) = ‖**x**‖^{2k+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., div**v**(**x**) = 0. The *divergence-free vector spline* (DFVS) solution is similar to that of VS except that the kernel matrix becomes **K**_{DFVS}(**x**) = [Δ**I** − ^{T}]*h*(**x**), and **p**(**x**) is also constrained to be divergence-free.

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: {**x*** _{n}*,

$$\text{arg}\phantom{\rule{thinmathspace}{0ex}}\underset{\mathbf{v}}{\text{min}}\phantom{\rule{thinmathspace}{0ex}}C(\mathbf{v})=\rho J(\mathbf{v}(\mathbf{x}))+\frac{1}{N}{\displaystyle \sum _{n=1}^{N}{({\mathbf{l}}_{n}^{T}\mathbf{v}({\mathbf{x}}_{n})-{w}_{n})}^{2}.}$$

(4)

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

$$\mathbf{v}(\mathbf{x})={\displaystyle \sum _{n=1}^{N}\mathbf{K}(\mathbf{x}-{\mathbf{x}}_{n}){\mathbf{l}}_{n}{c}_{n}+\mathbf{p}(\mathbf{x}),}$$

(5)

where the coefficients *c _{n}* are scalars, and

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.

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 *t*_{0}. At *t*_{0}, the tagging phase is a linear function of the point’s coordinate **x**, and wrapped to the range [−π, π), i.e., (**x**, *t*_{0}) = *W*(*k***x** · **l** + _{0}) , where *k* is the known tagging frequency, **l** is a unit vector representing the tagging orientation, _{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, _{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 *t*_{0} using HARP tracking [8]. For tagging direction **l** and a 3D spatial point **x _{j}** imaged at

$${w}_{j}=\mathbf{l}\xb7\mathbf{u}({\mathbf{x}}_{\mathbf{j}},t)=\mathbf{l}\xb7({\mathbf{x}}_{j}-{\mathbf{X}}_{j}).$$

(6)

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.

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** + _{X} **u**(**X**)) = 1 for any material point **X**, where ** _{X}** is the material gradient operator. As well, the spatial velocity field

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

$$\mathbf{w}(\mathbf{x},s)={\displaystyle {\int}_{0}^{s}\mathbf{v}(\mathbf{x}(\tau ),\tau )d\tau ,\text{\hspace{1em}and\hspace{1em}}\mathbf{v}(\mathbf{x}(\tau ),\tau )=\frac{d\mathbf{w}(\mathbf{x},\tau )}{d\tau}},$$

(7)

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., *s _{m}* =

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 *s*_{0} = 0 through *s _{M}* = 1. Let us denote

$$\mathbf{u}({\mathbf{x}}_{n})-\mathbf{w}({\mathbf{x}}_{n}({s}_{m}),{s}_{m})=\mathbf{v}({\mathbf{x}}_{n}({s}_{m}),{s}_{m})(1-\delta m),$$

(8)

so that

$${r}_{n}({s}_{m})={\mathbf{l}}_{n}\xb7\mathbf{v}({\mathbf{x}}_{n}({s}_{m}),{s}_{m})=\frac{{w}_{n}-{\mathbf{l}}_{n}\xb7\mathbf{w}({\mathbf{x}}_{n}({s}_{m}),{s}_{m})}{1-\delta m}.$$

(9)

With the *N* data samples, the continuous velocity field **v**(**x**, *s _{m}*) 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 *s*_{M−1}, ρ should be set to 0 so that the final displacement **w**(**x*** _{n}*, 1) matches the original data samples exactly—i.e.,

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*(*N*^{3}), 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 {**x*** _{n}*,

- Initialize ρ
_{0},*M*,**y**(0) =_{k}**y**,_{k}**w**(**y**, 0) = 0, and_{k}**w**(**x**, 0) = 0 for all_{n}*k*and*n*. - for
*m*= 0 to*M*− 1- set
*s*=_{m}*m*δ, ρ= ρ_{m}_{0}(*M*−*m*− 1)/(*M*− 1); - downsample the data points
**x**(*s*) if needed;_{m} - compute the interpolating coefficients with samples {
**x**(_{n}*s*),_{m}**l**,_{n}*r*(_{n}*s*)} and ρ = ρ_{m};_{m} - set
**w**(**x**,_{n}*s*_{m+1}) =**w**(**x**,_{n}*s*) + δ_{m}**v**(**x**(_{n}*s*),_{m}*s*),_{m}**w**(**y**,_{k}*s*_{m+1}) =**w**(**y**,_{k}*s*) + δ_{m}**v**(**y**(_{k}*s*),_{m}*s*),_{m}**x**(_{n}*s*_{m+1}) =**x**(_{n}*s*) +_{m}**w**(**x**,_{n}*s*_{m+1}), and**y**(_{k}*s*_{m+1}) =**y**(_{k}*s*) +_{m}**w**(**y**,_{k}*s*_{m+1});

- Set
**u**(**x**) =_{n}**w**(**x**,_{n}*s*) and_{M}**u**(**y**) =_{k}**w**(**y**,_{k}*s*), and the algorithm ends._{M}

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

$$\begin{array}{c}{u}_{x}=\frac{1}{90}{(\frac{8}{281}x+y)}^{2}+\frac{1}{15}(x-0.5y)-15\hfill \\ {u}_{y}=\frac{1}{900}{(\frac{8}{281}x+y)}^{2}+\frac{1}{15}(2x-y)-3\hfill \end{array}$$

on an image with ranges *x* [−50, 50] and *y* [−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.

(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.

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.

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.

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |