Walter G. O'Dell, Christopher C. Moore, William C. Hunter,
Elias A. Zerhouni, and Elliot R. McVeigh
Walter G. O'Dell,
Christopher C. Moore, William C. Hunter,
Elias A. Zerhouni, and Elliot R. McVeigh
The Departments of Biomedical Engineering and Radiology
JHU Medical Imaging Lab
Johns Hopkins University School of Medicine
Baltimore, Maryland, 21205
Previously presented material point tracking schemes for MRI tag data require identifiable points within the images, such as intersections between tags [Axel89improved, McVeigh91] , between tags and myocardial contours [Moore92] , or points along striped tag lines [Moore93SMRI]. However, analyses based on such sparse tag intersection data neglect valuable information contained in the intervening portions of the tag lines. Accordingly, the resulting descriptions of deformation suffer from poor spatial resolution. Intersection methods also require images with high spatial resolution both in the frequency and phase directions, leading to acquisition times longer than one breath-hold and thus limiting their clinical applicability.
In this paper we present a novel method for reconstructing the 3D deformation field of the left ventricle (LV) from tagged MR images which uses position and displacement information along the entire length of each tag and relies on accurately defined tag profiles [Atalar94] rather than poorly defined heart contours [Bazille94] . This approach eliminates the need for accurate simultaneous measurement of displacement in two dimensions and permits the computation of the 3D deformation gradient tensor at any point in the heart wall.
To test the performance of the method, a computer-generated model of a prolate spheroid undergoing deformations simulating those measured in the beating heart were analyzed. Finally, the method was applied to parallel tagged, breath-hold cine data sets from ten human hearts. Noise propagation properties using a Monte Carlo analysis were then tested using human cardiac geometry and deformation.
|
Fig. 1 Images of an in vivo human heart from a normal volunteer
using the conventional tagging and imaging protocol. The progression (from left to
right) through three phases in the cardiac cycle are shown: early, mid, and
late systole. Two cardiac short axis and one long axis image are shown, each
displaying tag lines from a different set of mutually orthogonal tag planes.
|
The images were then processed by a semi-automated software package [Guttman94] to define positions along the tag lines and around the endocardial and epicardial heart contours. A typical tag and contour data set for a normal human after image processing is shown in Figure 2.
Fig. 2 Left: A typical 3D tag data set for x-displacement in
an in vivo human heart from a normal volunteer, showing the appearance
of the tags and the contours near end-systole on 7 short axis image planes.
Right: A typical 3D contour data set and the estimated prolate spheroidal
centroid (upper grey circle) and apical focal point (lower black circle).
For a tag separation in the reference state of 6.0 mm, approximately twelve tags were produced in the myocardium in each set. At a point separation along each tag of 1 mm, this gave over 4600 raw data points, of which every other point was used for the fitting.
|
Fig. 3 Depiction of a short image at some time after initial tagging.
The inset depicts an enlarged view of one deformed tag line after detection
of tag points at 1 mm intervals.
|
|
Fig. 4 Deformation of a tag plane into a tag surface in 3D. The
bold curves are the intersections of the tag surface within the heart
wall with four vertical image planes. The enlarged dots on these curves are
the tag points generated by an automated tag detection algorithm.
 xn represents the 1D displacement associated with
the nth tag point. The rectangle in the figure at the upper left shows the
region of the heart depicted in the lower figure.
|
The key concept of the method is shown in Figure 4.
In three dimensions, each observed tag line represents the intersection
of a deformed tag surface with the 2D image plane. The tag surface
results from the deformation of an initially flat tag plane. The figure
depicts a tag surface within the heart wall that is cut by four
parallel image planes. In each plane, tag points were marked along tag
lines between the endocardial and epicardial heart contours. For each tag point
on the deformed tag surface, the component of its displacement vector that
was normal to the tag plane (x in Figure 4) and
its position were measured. Reconstruction of the displacement field was
performed using the tag point displacement and position data to solve for the
coefficients of a series expression describing the field. This is analogous to
using digital sampling and Fourier analysis to reconstruct a continuous 1D
signal. Here, this expression described the displacement over the entire LV.
In the simplest case, each unidirectional displacement data set was used to
solve for an interpolation function for the displacement field in that one
direction, independent of the motion in the other Cartesian directions.
Once the three independent displacement field expressions were determined, they
could then be used to compute the 3D displacement of any point in the deformed
heart wall. At higher complexity, displacement data from multiple directions
can be used simultaneously to solve for the full displacement vector field in
an any coordinate system. The displacement field at each time frame was
computed independently, using the time at which the tag planes were generated
as the universal reference state.
x(x,y,z) (refer
to Figure 4), took the form:
x = a0
+ a1x + a2y + a3z
In vector notation this can be written [in this paper vectors are denoted with a lower case bold italic font, tensors are upper case bold italic, and scalars are plain italic.]
x = a p
where a = [a0,a1,a2,a3] and p = [1,x,y,z].
For each sampled tag point (xn,yn,zn), each term in the vector
p was evaluated and the x displacement value,
xn,
measured. By considering the x-displacements and p vectors of
all tag points, a series of simultaneous equations was created and the unknown
coefficients (a) were fit using singular value
decomposition (a least-squares fitting method)
[Press88]. Series expressions for the y and z
displacement fields were generated and fit analogously.
|
Fig. 5 The prolate spheroidal coordinate system. Surfaces of constant
 are ellipsoids and surfaces of constant longitudinal
coordinate,  , are hyperboloids.
|
Figure 5 shows the prolate
spheroidal coordinate system where is
the radial component, is the longitudinal
angle, 
is the circumferential angle, and
f is
the focal length. These are related to Cartesian coordinates via
x = f sinh()
sin()
cos()
y = f sinh()
sin()
sin()
z = f cosh()
cos()
For each prolate spheroidal coordinate displacement field, a series expansion
in prolate spheroidal coordinates (,
, )
was created using
a generating function that is analogous to the spherical harmonic series:
where Pl m is the associated Legendre polynomial function, and the unknown coefficients, ai, were numbered sequentially. For the typical tagging and imaging geometry, and typical image signal-to-noise characteristics, we found that a fit using N=1, L=4 (where N is the radial fitting order and L is the angular fitting order) was most appropriate [ODell94ffSMR], generating 50 free fitting parameters per coordinate direction.
To fit for the coefficients of the prolate spheroidal series expansions from measures of displacement in the Cartesian coordinate directions, the measured displacements were projected onto the local prolate spheroidal axes. The infinitesimal displacements in Cartesian and prolate spheroidal coordinates are related by the following directional derivatives:
where
x =
[x,
y,
z],
=[,
,
]
c
and the J is the Jacobian matrix for the this coordinate
transformation. For a, b, and c defined
as the vectors of unknown coefficients for the
,
,
and
expressions, respectively, and p the vector of
50 series terms derived from the generating function, then
=
a*p,
=
b*p,
and
=
c*p.
For the case where the Cartesian axes (as defined by the three orthogonal
sets of tag normals) are centered and aligned with the prolate spheroidal
coordinate system as described in
equations 1-3, the measured
Cartesian displacements can be expressed:
where now
=
[a0,...,a49, b0,...,b49, c0,...,c49],
x =
[ J00 p0, ... , J00 p49,
J01p0, ... , J01p49,
J02p0, ... , J02p49 ].
y and
z
were defined analogously. Similar to the method of solution of the
Cartesian power series coefficients above, each measured displacement value
was collected into a single row vector,
x, and the
corresponding
vectors were inserted as columns in a matrix P.
The unknown coefficients , were then
solved using the equation
x =
* P
In this way, all the unknown coefficients can be solved simultaneously
using all the available tag data and singular value decomposition.
Although for finite deformations the relations between the displacements
in x, y, and z and those in
, ,
and are not
exact, the simultaneous least-squares fitting will always find
the optimal fit of the given displacement data to the series expansion.
The fitting procedure was performed independently for each time frame in
the data set. Since the heart geometry changes throughout the cardiac
cycle, a new prolate spheroidal coordinate system was calculated
for every time frame based upon the measured epicardial contour point
data. This was done by first estimating the 3 coordinates of the centroid
and the 3 coordinates of the apical focus (point on the long axis, toward
the apex, at a distance of a focal length from the centroid) and computing
a least-squares optimal
value at this centroid and apical focus
for the given set of epicardial LV contour points. Within an iteration loop,
the 6 coordinates were
then independently adjusted until the combination of centroid and focal point
coordinates and radial coordinate, ,
was found that produced a prolate spheroid with the smallest fitting error to
the measured LV contour points.
Variation of the prolate spheroidal coordinate centroid and apical focus required a modification to Equation 6 to account for the coordinate transformation necessary to align the axes of the new spheroidal coordinate system with that defined by the normals to the undeformed tag planes. For the transformation to the new coordinate system (x') defined by the rotation matrix R, such that x' = R * x, Equation 4 becomes:
x = R-1 *
J *
and this modification was carried through the derivation. In addition, this formalism permits 3D reconstruction from sets of images with arbitrary tag normals, so long as three dimensions are spanned.
u, which was computed
numerically by taking partial
derivatives of the displacement expressions with respect to the three
coordinate directions. Because the terms in the series expansion were computed
based on the coordinates of the tag points in their deformed state,
u described the motion from the
deformed state back to the undeformed state. Thus, to compute F from
u required an
inverse operation:
F = (I - u)-1, where
I is the identity tensor. The Lagrangian finite strain tensor,
E, can then be computed from
E = 1/2 (FT F - I)
where the `T' superscript denotes matrix transpose.
For both the mathematical and experimental models, a set of material points
was defined within the heart wall in the first time frame and 3D displacements
and deformations were evaluated over time at these points as they moved
through space. A mesh of 96 material points regularly spaced in
,
, and
was selected, based on the imaging and heart
geometries. The material point locations at the first time frame were
constrained longitudinally by the most basal and most apical short axis
image slices and radially by the endocardial and epicardial contours. The
interpolation of the contour boundaries between image slice locations was
achieved by fitting the LV epicardial and endocardial contour points of the
first time frame to a 4th order prolate spheroidal coordinate expansion
for as a function of the two angular
coordinates.
,,
) as
a function of the initial coordinates
(
,
,
) and the
deformation parameters, were as follows:
where f is the focal length of the prolate spheroidal coordinate system
and Vol() is the volume within a
prolate spheroidal shell of constant .
Deformation in the radial component, , was
limited to a volume conserving (incompressible) contraction mode
by setting = 1.
The initial and deformed short axis endocardial
equatorial radii, f*sinh(endo)
and f*sinh(endo), were chosen to
approximate realistic heart geometry and resulted in an ejection fraction of ~50
percent. The prolate spheroidal deformation parameters are given in
Table 1.
Table 1 Axisymmetric deformation model parameters on a
prolate spheroid.
Table 2 Cartesian deformation parameters.
Bulk motions and linear deformations incorporating ellipticalization and shears in x, y, and z were then superimposed upon the axisymmetric modes in the manner of Arts, et al. [Arts92]. The magnitudes of these additional modes were chosen to give the overall model approximately physiological rigid body motions and to reproduce those deformations observed by Arts, et al. These deformation values are given in Table 2. Image prescription error, where the central axis of the imaging volume is not coincident with the long axis of the left ventricle, was also simulated with representative values of 5 degrees offset about the X axis and 8 degree about the Y axis.
Seven image slices and ten tag planes were simulated for each each image set. Separation between adjacent data points along a given tag was fixed at 2 mm. Over the entire 3D data set, this produced approximately 2400 tag data points. The prolate spheroid deformation model is shown in Figure 6. The two short axis image sets are spatially coincident, giving a grid tag appearance.
| Fig. 6 Simulated short axis (Left) and long axis tag data sets for the prolate spheroid model. Physiologic deformation modes with added Cartesian motions were modeled for 7 image slices and 10 tag planes in each view. The angles of tilt are due to simulated image prescription error. |
x,
y, and
z data sets, and the resulting 3D
material point
trajectories were computed for the 96 material points over 100 trials.
A single root-mean-square deviation (RMSD) value was computed for
each point from the 3D distribution of the cloud of estimated point
trajectories around the expected location.
Two measures of the displacement field-fitting performance were evaluated with the prolate spheroid deformation model: tracking and fitting. The 3D tracking performance was evaluated by comparing the estimated material point trajectories derived from the field-fitting reconstruction to the exact solution computed from the deformation model Equations 6-8. The mean absolute value and the standard deviation of the tracking error were computed for the set of 96 material points. The fitting performance was evaluated by computing the standard deviation of the error (SDE) (around a mean error of zero) between the estimated and the measured 1D displacements at each tag point. Thus, the fitting performance reflects the ability of the reconstruction to account for the variations in the displacement field at the locations where the measurements were made, while the tracking performance reflects both the quality of the fitting and the interpolation of the displacement field between the measurement locations.
Fig. 7 Short and long axis images of a human heart at 228ms into systole,
overlaid with tag point locations predicted from the 3D field-fit.
Figure 7 shows representative short and long axis human cardiac images 228 ms into systole and demonstrates the agreement between the predicted locations of the tag points and the actual tags in the MR images.
| Fig. 8 The root-mean-square deviation (RMSD) of the computed material point position as a function of input noise level on a human heart geometry and deformation field. The distribution of tracked points formed an isotropic cloud around the noise-free reference trajectory. The RMSD of the cloud at the endocardial and epicardial surface points was the same, but of a larger value than at midwall point locations. The error bars correspond to the standard deviation in the RMSD over the 100 Monte Carlo trials. |
The scatter of points around the expected value was found to be isotropic. The root-mean-square deviation (RMSD) of the computed material point position increased linearly as a function of input noise level, as shown in Figure 5. The scatter of points around the expected value was found to be isotropic. The RMSD at the endocardial and epicardial material points was found to be equal, and was greater than that found at the midwall. At an input noise standard deviation of 0.5mm, the midwall RMSD tracking error was 0.077 +/- 0.015 and that at the endocardium and epicardium was 0.126 +/- 0.030. The RMSD computation using a subset consisting of 50 trials gave results that were not significantly different, suggesting that 100 trials were sufficient to obtain a stable convergence of the estimate in tracking precision.
The combination of parallel-tag MR imaging and displacement field-fitting is an accurate and robust method for reconstructing the 3D deformation field over the entire left ventricle. Because of its reliance on parallel tag data, the displacement field-fitting method is inherently less susceptible than previous MRI reconstruction schemes to tag and contour detection error since only the centerline of any tag need be determined during image analysis. The incorporation of a greater number of points along each tag line increases the sampling density over the heart compared to techniques that use only those points located at intersections between tags or between tags and myocardial contours.
The fitting of displacement fields can also be realized using other methods, including energy minimization techniques [Young92, Denx94WNAM]. In this approach, the fits can be computed numerically for a predefined set of material points (nodes) and interpolation between these points accomplished via linear or finite element interpolation functions. Although field-fitting by solving for series coefficients has some mathematical similarities to energy minimization and other methods, the series method has several advantages. First, it is independent of material strain energy models. Second, the global fit provides the ability to express the 3D deformation at any myocardial point with a relatively small number of parameters, independent of a predefined set of material points.
The field-fitting approach is immediately suitable to other methods of 3D strain analysis where motions in multiple directions are independently measured. For example, three sets of velocity encoded images [Pelc91] could be used to solve for a set of expressions describing the velocity field over a prescribed region of interest, in a manner analogous to that described for tag-based displacement fields. This would have the advantage of smoothly interpolating the velocity field between pixels and of filtering out uncertainty in the velocity measurements at each pixel. Grid tag data sets can also be analyzed by interpreting the grid as two independent sets of parallel tags.