Analytical method to measure three-dimensional strain patterns in the left ventricle from single slice displacement data
- Abbas Nasiraei Moghaddam^{1, 2}Email author,
- Nikoo R Saber^{2},
- Han Wen^{3},
- J Paul Finn^{1},
- Daniel B Ennis^{1} and
- Morteza Gharib^{2}
https://doi.org/10.1186/1532-429X-12-33
© Moghaddam et al; licensee BioMed Central Ltd. 2010
Received: 12 January 2010
Accepted: 1 June 2010
Published: 1 June 2010
Abstract
Background
Displacement encoded Cardiovascular MR (CMR) can provide high spatial resolution measurements of three-dimensional (3D) Lagrangian displacement. Spatial gradients of the Lagrangian displacement field are used to measure regional myocardial strain. In general, adjacent parallel slices are needed in order to calculate the spatial gradient in the through-slice direction. This necessitates the acquisition of additional data and prolongs the scan time. The goal of this study is to define an analytic solution that supports the reconstruction of the out-of-plane components of the Lagrangian strain tensor in addition to the in-plane components from a single-slice displacement CMR dataset with high spatio-temporal resolution. The technique assumes incompressibility of the myocardium as a physical constraint.
Results
The feasibility of the method is demonstrated in a healthy human subject and the results are compared to those of other studies. The proposed method was validated with simulated data and strain estimates from experimentally measured DENSE data, which were compared to the strain calculation from a conventional two-slice acquisition.
Conclusion
This analytical method reduces the need to acquire data from adjacent slices when calculating regional Lagrangian strains and can effectively reduce the long scan time by a factor of two.
Keywords
Introduction
The left ventricular deformation associated with the pumping of blood from the heart is known to be complex [1]. For example, it is known that there are significant regional variations in the time course, magnitude, and spatial pattern of deformation within the ventricular wall [2]. Moreover, in the course of heart failure and other cardiac diseases, wall motion abnormalities become very pronounced and deviate from the normal, healthy function of the myocardium. Therefore, it is of significant clinical interest to quantify these regional wall motion abnormalities with high temporal and spatial resolution and therefore derive indices of the nature and extent of cardiac dysfunction to better inform clinical decision-making.
Cardiovascular magnetic resonance (CMR) motion encoding strategies have emerged in recent years as clinically useful tools for assessing derangements in regional cardiac function. In particular, CMR phase contrast (PC) [3, 4], tissue tagging [5, 6], HARP [7], and DENSE [8] have arisen as powerful techniques for encoding cardiac motion in a CMR exam and each has specific advantages and disadvantages. DENSE displacement encoding is similar in principle to the PC velocity-encoding method in that the motion components of small tissue volumes are encoded in the phase image at every image pixel location [9]; hence, the number of measures across the LV wall is typically higher than the number of tags across the wall in a tagging experiment. As an example, DENSE has been used successfully for detecting small focal regions of abnormal contraction in patients [10], in characterizing mouse ischemia models [11], and in acquiring detailed functional data in canine hearts in vivo [12].
The DENSE method has been implemented in different ways. In the meta-DENSE technique, which is used for data acquisition in this study, a series of phase-labeled images are acquired at the same time point in the cardiac cycle. The tissue can be phase-labeled at different points in the cardiac cycle prior to the image acquisition. The result is that the same pixel location in each image always corresponds to the same physical tissue. This ensures that the behavior (displacement and deformation) of the same "particles" within the tissue are followed within the frame of reference. This is oftentimes referred to as the "Lagrangian" framework for "tissue tracking." Kim et al. have proposed a variation of the DENSE sequence for two-dimensional (2D) breath-hold cine DENSE imaging [13] which was later extended to 3D by using the slice following technique [14]. Compared to meta-DENSE, cine DENSE acquires the data faster but with lower signal-to-noise ratio (SNR) and not in a Lagrangian framework.
The three-dimensional (3D) principal shortening strain is typically not oriented in any standard short or long axis imaging plane, therefore 2D measures typically underestimate the maximal shortening undergone at any point. From the 3D Lagrangian displacement vector for each pixel in the image, we are able to calculate the Lagrangian strain, which provides a complete description of myocardial deformation. The Lagrangian strain is calculated from the spatial gradients of the displacement vector field. In general, this requires that there exist two or more displacement vectors along each of the orthogonal directions along which the gradients will be computed. With typical 2D images, however, there are insufficient displacement data in the through-plane direction and calculation of the full 3D strain tensor is not possible unless we acquire displacement data in orthogonal directions or parallel planes. Acquiring multiple planes of data, however, increases the number breath holds, thereby extending the exam duration. The analysis can also be complicated in multiple plane acquisition because of image misregistration [15].
Several efforts have sought to make 3D strain calculation more efficient. Osman et al. acquired the out-of-plane strain using SENC [16], which derives the through-plane strain from the frequency of tag planes applied parallel to and within the image plane. Abd-Elmoniem et al. proposed zHARP [17], which is a modified version of slice following-CSPAMM [18, 19] and uses HARP processing to provide 3D displacement in a Lagrangian frame. They recently used "stacks of zHARP images" to derive the 3D strain tensor [20]. Hess et al. proposed a method that utilizes 2D DENSE in-plane displacement measurements in conjunction with a SENC through-plane strain measure [21]. These techniques are based on two or more slices and require at least 4 breathhold acquisitions. Therefore, they suffer from inter-breathhold myocardial position variability, which can cause a significant increase in the strain error (RMS error of 0.04 in through-plane normal strain) [21]. Sampath et al. presented a combination of HARP and SENC for 3D strain measurement [22]. The latter is the only technique that requires no more than one image slice for strain measurement; however, it does not calculate the complete 3D strain tensor because the shear components are not measured.
The goal of this work was to define the mathematics for calculating all components of the 3D strain tensor from a single slice of 3D displacement data by assuming tissue incompressibility. The importance of calculating out-of-plane strains from a single slice is emphasized by noting that in many strain imaging techniques (such as zHARP, slice following 3D cine DENSE and DENSE-SENC), the distance between adjacent slices may be too large to accurately calculate the through-plane strains due to the curvature of the heart in the through-plane direction. Herein, we derive a mathematical approach and the image processing tools to increase the derived 3D deformation information acquired from a DENSE CMR exam. Our methods are sufficiently general and can easily be extended to other 3D deformation imaging methods (e. g. zHARP). Specifically, we calculate the out-of-plane components of the Lagrangian strain tensor in addition to measuring the in-plane components from a single-slice DENSE CMR dataset covering a complete cardiac cycle.
The method has been validated using a simulated data set, which has an analytical deformation gradient tensor and mimics the incompressible deformation of the moving myocardium. It has also been compared to conventional methods that use two adjacent slices for strain measurements. Finally, we demonstrate the feasibility of the method in a human subject and conclude by comparing our results to those of other studies.
Theory
Deformation gradient tensor
where I is the identity tensor. Tensor F is calculated in the Cartesian coordinate system and is subsequently transformed to a locally polar coordinate system, which is more popular in cardiac strain measurement, accounting for the elongated ellipsoid shape of the left ventricle. The three unit vectors of this system lie in the local Radial, Circumferential, and Longitudinal (RCL) directions at any given point of the heart. Once the long axis of the heart and the local radial direction are known, these three directions can be defined uniquely for any point within the myocardium, including those on the epicardium and endocardium. These directions can be robustly calculated almost everywhere in the myocardium, apart from near the apex where they are not well defined (e. g. the radial and longitudinal directions become coaxial).
Displacement smoothing
The elements of the deformation gradient tensor are calculated using spatial derivatives (F_{ ij }= ∂x_{ i }/∂X_{ j }), which is a noise sensitive operator. Therefore, the spatial and temporal smoothing of the measured displacement should be applied prior to taking the derivative. Previously, it has been shown that the significant temporal frequencies of the heart motion are lower than five times the heart rate [23]. Hence, a fifth order Fourier basis function should be sufficient for temporal smoothing.
Noise can induce outlier displacements in the tissue tracking process, especially in regions of low signal intensity where phase measurements are especially sensitive to noise. These outliers can generate large errors in the calculation of the local strain. Median filters are quite popular and effective for suppressing this type of noise [24]. To decrease the possible adverse effects of filtering on the strain calculation, the smallest symmetric 2D median filter, which is the five-pixel plus sign-shaped filter, was applied to each component of the displacement vector to reduce the effect of the outliers. This is the minimum amount of smoothing that can be obtained with a median filter. For a displacement field with a higher noise level, however, we may need to increase the size of the applied median filter. Although components of Frelating to the out-of-plane direction are more noise-sensitive, no further filtering or smoothing was applied to the displacement in the out-of-plane direction. Instead, we invoked specific mathematical methods to calculate individual elements of the Ftensor, which determine the out-of-plane deformation.
Estimation of the deformation gradient tensor
It is important to choose a suitable region for the smoothing area. A large value for the smoothing area will over-smooth the physiologic spatial variations in elements of Frather than merely reducing the noise, whereas a small smoothing area does not reduce the noise sufficiently. This matter is discussed in more detail in the Discussion section.
In the conventional method for strain estimation, based on two or more adjacent slices, an equivalent procedure is followed, but points are chosen from both adjacent planes. With one slice, however, the small spatial variation of the data points in the direction normal to the surface nearly imposes a singularity condition, which makes the gradients in that direction highly sensitive to noise. It requires, therefore, exploiting some physical/mathematical constraints. This singularity condition does not exist when two adjacent slices are employed for strain measurement in the conventional methods.
To improve the accuracy of our calculation based on a single slice data, we compute the direction normal to the imaged surface for each neighborhood and redefine a local coordinate system using this normal vector as the out-of-plane unit vector. Calculation of Fin this new coordinate system ensures that the effect of the singularity condition only appears in the third column of the deformation matrix. It also implies that the expected error in the calculation of the third column of Fis greater than the expected error in the first two columns of F. Recall that Fis not symmetric and since the out-of-plane motion is measured by meta-DENSE with similar precision to the in-plane components, the calculation of the third row of Fis as well-posed as its first two rows. In order to compute the third column (out-of-plane elements) of Fwith reduced noise sensitivity, we employed the following mathematical methods and physical assumptions.
1) Incompressibility of the myocardium
Myocardial tissue volume may vary over the cardiac cycle due to blood volume changes. However, measurements in canine hearts in vivo show that the level of volume fluctuation is approximately 1% [26]. Therefore, tissue incompressibility is a reasonable assumption in the heart. The deformation gradient tensor for an incompressible body has a determinant equal to unity. This extra equation increases the robustness in calculation of the third column of Fand was enforced in all of our calculations.
2) Combining F with the inverse-deformation gradient tensor G
Fis the tensor that transforms the reference positions of the myocardial points within a neighborhood (smoothing area) to an arrangement, as close as possible, to their deformed shape. It is also possible to simultaneously define the tensor Gto perform the inverse operation. Gmaps the same set of points in the deformed shape to an arrangement, as close as possible, to their reference position. Ideally, the Gtensor should be the inverse of the Ftensor, but this is not actually the case, since we can only calculate a least-square estimate for each of these tensors. In fact, neither of the matrices performs the perfect matching of original and deformed points as a consequence of least-squares estimation.
E _{ 33 } Correction
where λ_{ iF }and λ_{ iG }are the stretch ratios derived from Fand Gwhen they are reduced to a 2 × 2 (in-plane) matrix.
3) Semi-symmetric rotation in the long axis view
In the short axis view, the effect of left ventricular rotation can be calculated with high precision, eliminating the need for a symmetry assumption. Left ventricular torsion, however, cannot be approximated and used for calculation of longitudinally associated shear strains.
Methods
DENSE CMR data acquisition
In this study, we used a specific version of DENSE CMR sequence called meta-DENSE [27]. DENSE CMR images were acquired on a three-Tesla (3T) whole body CMR scanner (Siemens Trio, Erlangen, Germany) at the Caltech Brain Imaging Center with an eight-channel cardiac array coil. Our study was approved by the institutional review board, and informed consent was obtained. A single-slice long axis dataset of high spatial and temporal resolution was acquired from the heart of a healthy 25-year-old female volunteer. The experiment was carried out under free-breathing conditions over a period of 40 minutes, which is directly proportional to the number of phases acquired (20 phases for this study). To synchronize the pulse sequence with the motion of the chest and heart, respiratory and heart monitoring were achieved using a pneumatic bellows and ECG monitoring.
In order to acquire directly the Lagrangian displacement vectors relative to a fixed time during the cardiac cycle, the image acquisition section of the DENSE sequence was positioned at a constant time in the late diastolic phase of the cardiac cycle when the myocardium is almost at rest. This fixed time was determined as 715 ms after the R-wave of the ECG trigger for this volunteer, while the encoding section was placed at a series of time steps in 30 ms increments after the ECG trigger.
The imaging parameters included FOV = 192 mm × 108 mm; acquisition matrix size = 128 × 72 pixels, pixel size = 1.5 mm × 1.5 mm; prescribed slice thickness = 5 mm; repetition time (TR) = 3.1 ms; mixing time (TM) = 600 ms; number of averages = 3; the displacement encoding sensitivity (k_{e}) = 0.80 and 0.93 cycle/cm for in-plane and out-of-plane directions respectively. 20 phases were acquired with increasing delays of 30 ms increments after the trigger, which together covered 600 ms of the cardiac cycle starting 30 ms after detection of the QRS complex and ending during diastasis (acquisition window = 850 ms). In this manner, the time-resolved deformations from the beginning of systole throughout most of diastole were covered. The heart is almost at rest during the part of the cycle not captured in the imaging process (almost 150 ms), such that negligible strain is expected within that period.
For validation purposes, two more data sets were acquired at UCLA (Department of Radiological Sciences) using a 1.5 Tesla scanner (Siemens Avanto, Erlangen, Germany). In one study, the subject was a normal, 35-year-old healthy male volunteer and the subject of the second study was a male patient with acute myocardial infarction (MI). Studies were approved by the institutional review board and conducted with the informed consent of the subjects. In the healthy volunteer, first a single 5 mm thick long-axis slice was acquired with meta-DENSE for processing with the proposed methods. Subsequently, two additional 5 mm slices, whose centers were separated by 8 mm, parallel to and centered about the first slice, were also acquired with meta-DENSE for processing with the conventional method. Imaging parameters, other than the number of slices, were similar for both acquisitions and only one phase was acquired for this validation study. We computed the strain from the single long axis (LA) slice using the proposed algorithm and compared it to the strain calculated from the conventional 3D method applied to the displacement field of these two latter slices.
For the patient study, which was a preliminary effort to evaluate the diagnostic value of the proposed approach, we applied our method to quantify the strain changes caused by myocardial infarction (MI). In this HIPAA compliant study, we first located the infarcted region in the myocardium of the patient using Late Gadolinium Enhancement (LGE) CMR. We next acquired DENSE images from the same single short axis slice to investigate myocardial contraction through strain mapping.
DENSE data pre-processing
After obtaining the DENSE images, segmentation was performed by semi-automatic masking of all parts of the anatomy except for the left ventricular myocardium. Since meta-DENSE always acquires the image at a specific time from the beginning of the cardiac cycle, each image pixel corresponds to one specific myocardium voxel in all cardiac phases. Therefore, performing segmentation for only one time frame is sufficient for masking all other time steps. Using the DENSEView software [28], the 3D displacement field was generated in a matrix format to be used with our algorithm for strain calculation. From the three components of the displacement for each image pixel, it is possible to generate an image of the displacement magnitude. If this image is created at a time point that has undergone a large displacement (end systole to end diastole), it becomes a robust indicator of myocardial anatomy and hence a good template for segmentation by manual masking. It is useful to have the displacement magnitude images of a few (3-4) phases, at maximum displacement, averaged together when performing the manual masking, as it helps to avoid inaccurate masking at noisy points.
Phase unwrapping was then performed on the segmented images by scanning the myocardial area through a region-growing algorithm, while searching for sudden changes in the displacement magnitude. These changes were fixed later by adding or subtracting the correction values, which are proportional to integer multiples of 2π changes in phase. This step was repeated separately for all three directions of displacement and for all phases. Spatial and temporal filtering of the displacement vectors was performed at this step, as explained in Displacement smoothing subsection.
Strain measurement
Using the algorithm described in the Theory section, the strain maps were then calculated from the three-dimensional displacement field for each of the myocardial material points at every time frame. To express the strain values in the local polar coordinate system, the RCL directions are required at all measurement points. In our method, these directions were calculated automatically for every single point of the myocardium. Since all deformations were calculated in comparison to the position of points at end diastole, the calculated directions at rest were used to express the strain tensor in the RCL coordinate system for all time steps.
The steps outlined in the Theory section are expected to lead to a tolerable strain calculation, such that in each local neighborhood, there exist a sufficient number of points with acceptable strains to estimate the elements of the strain tensor in that neighborhood. To void and recalculate the unreliable values within the strain tensor, strains whose magnitudes were different from the mean by more than 2 standard deviations were filtered out. Highly non-physiological strains with absolute values greater than one were eliminated before the aforementioned mean and standard deviation calculations. Next, the five pixel median filter, described in the Theory section, was applied to the strain map, in order to account for the outliers. To fill the voids created by the ignored points and to smooth the strain, a moving average filter acting over a 1 × 1 cm^{2} area was used.
Validation of the method
The proposed algorithm was validated first with a simulated displacement field that has its strain analytically calculated. Next, we validated our algorithm by comparison to strain calculation based on two adjacent slices. Finally, we compared our results to statistical estimations provided by Moore et al. [29] and Young et al. [30] over actual MR acquired data sets. Moore's measurement is based on 3D tagging dataset acquired for n = 31 hearts, while Young et al. has used a 3D finite element MR tagging-based model on 12 volunteers.
The simulated displacement field mimics the geometry and wall motion of LV excluding its apex. It is defined in a polar coordinate system where the Z direction coincides with the main long axis of the heart. The geometry of the LV before and after deformation is shown in Figure 2, with its parameters and analytical strain fully described in Appendix B. The simulated displacement field was fed into our algorithm and then the resulting values for the deformation gradient tensor, as well as the strain tensor, were compared to their analytical values. In the absence of noise, the previously described steps involving "displacement smoothing" and "filtering of the strain" were omitted to focus only on the effect of the singularity on the calculation of deformation gradient elements. To study the effect of noise on the algorithm, additive noises were simulated and independently added to each of the displacement elements. Each additive noise had a normal distribution, zero mean and varying standard deviation, which is normalized by the maximum displacement in the overall myocardial motion.
Results
Performance of the algorithm on simulated data
Accuracy of the elements of the deformation gradient matrix in the local coordinate system.
shear elements | normal elements | ||||||||
---|---|---|---|---|---|---|---|---|---|
F _{ 12 } | F _{ 13 } | F _{ 21 } | F _{ 23 } | F _{ 31 } | F _{ 32 } | F_{ 11 }-1 | F_{ 22 }-1. | F_{ 33 }-1 | |
Variation in Magnitude | 0.006 | 0.113 | 0.320 | 0.013 | 0.179 | 0.068 | 0.381 | 0.192 | 0.152 |
RMS Error | 0.000 | 0.026 | 0.107 | 0.013 | 0.008 | 0.002 | 0.001 | 0.027 | 0.026 |
Relative RMS | 3.7% | 23.1% | 33.5% | 98.2% | 4.6% | 2.3% | 0.3% | 13.9% | 17.1% |
Error in the strain element values, corresponding to the deformation gradient tensors reported in Table 1 for simulated data at maximum contraction.
Shear | Normal | ||||||
---|---|---|---|---|---|---|---|
RL | LC | CR | RR | LL | CC | ||
Analytical Magnitude | 0.136 | 0.033 | 0 | 0.523 | 0.163 | 0.133 | |
E | RMS Error | 0.047 | 0.005 | 0.017 | 0.036 | 0.023 | 0.025 |
Relative RMS | 34.6% | 13.8% | NA | 6.8% | 13.8% | 18.5% |
The first row of Table 1 shows the variations in the analytical values of the elements of the deformation gradient tensor from the identity matrix (F-I). The second and third rows respectively show the RMS of the error in our estimation and the relative RMS error (i.e. the RMS error normalized by the magnitude of the variation for the same element). These values were calculated for the elements of the deformation gradient tensor in the local coordinate system to avoid the propagation of the singularity condition, explained in the Theory Section, from the third column of the matrix Fto its other elements. The first and second direction base vectors are nearly equivalent to the Radial and Longitudinal directions, respectively. Nearly all element values are in agreement with analytical calculations and RMS errors are less than 0.03. The only exception is the in-plane element of F_{ 21 }which has a RMS error of 0.11, corresponding to 33.5% relative RMS error.
Table 2 shows that the RMS error for all six elements of strain is less than 0.05. In particular, the agreement of our algorithm with analytical values for those elements associated with the out-of-plane (circumferential) direction is notable. The highest error was for in-plane shear strain, E_{ RL }, which corresponds with the high relative error in F_{ 21 }.
Effect of noise and smoothing area on strain calculation for simulated data at maximum contraction.
Shear | Normal | |||||
---|---|---|---|---|---|---|
Noise (σ/max displacement) | RL | LC | RR | LL | CC | |
Relative RMS Error in E (smoothing area = 25 mm^{2}) | 0 | 34.6% | 13.8% | 6.8% | 13.8% | 18.5% |
0.01 | 34.6% | 15.7% | 8.1% | 14.1% | 21.6% | |
0.02 | 36.1% | 25.5% | 10.9% | 15.9% | 27.1% | |
0.03 | 36.0% | 32.6% | 13.3% | 19.0% | 37.0% | |
0.05 | 45.9% | 46.9% | 21.4% | 23.2% | 61.0% | |
Relative RMS Error in E (smoothing area = 15 mm2) | 0 | 27.4% | 14.4% | 5.7% | 9.4% | 12.2% |
0.01 | 28.0% | 21.2% | 8.0% | 10.3% | 17.0% | |
0.02 | 30.9% | 30.7% | 11.1% | 14.0% | 28.6% | |
0.03 | 34.2% | 42.0% | 15.3% | 21.1% | 46.8% | |
0.05 | 46.4% | 64.7% | 24.7% | 28.2% | 66.3% |
Comparison to the strain from two adjacent slices
High spatio-temporal strain measurement
Region | Sepal Wall | Lateral (free) wall | |||||
---|---|---|---|---|---|---|---|
Strain | Present method | Moore et al. | Young et al. | Present method | Moore et al. | Young et al. | |
E_{LL} | Basal | -0.12 | -0.14 ± 0.03 | -0.14 ± 0.03 | -0.22 | -0.15 ± 0.03 | -0.19 ± 0.04 |
Equatorial | -0.17 | -0.15 ± 0.03 | -0.15 ± 0.02 | -0.18 | -0.14 ± 0.04 | -0.16 ± 0.03 | |
Apical | -0.15 | -0.18 ± 0.04 | -0.19 ± 0.02 | -0.22 | -0.19 ± 0.03 | -0.18 ± 0.03 | |
E_{CC} | Basal | -0.17 | -0.17 ± 0.03 | -0.19 ± 0.03 | -0.18 | -0.21 ± 0.03 | -0.21 ± 0.03 |
Equatorial | -0.16 | -0.16 ± 0.03 | -0.20 ± 0.02 | -0.17 | -0.22 ± 0.03 | -0.21 ± 0.02 | |
Apical | -0.22 | -0.18 ± 0.03 | -0.20 ± 0.03 | -0.18 | -0.24 ± 0.04 | -0.22 ± 0.02 | |
E_{RR} | Basal | 0.40 | 0.45 ± 0.12 | 0.21 ± 0.10 | 0.58 | 0.52 ± 0.19 | 0.25 ± 0.09 |
Equatorial | 0.52 | 0.42 ± 0.19 | 0.16 ± 0.10 | 0.55 | 0.38 ± 0.18 | 0.21 ± 0.10 | |
Apical | 0.55 | 0.36 ± 0.22 | 0.07 ± 0.06 | 0.64 | 0.49 ± 0.29 | 0.10 ± 0.06 |
The peak values of shear strain averaged over six regions on the septum and LV free wall as shown in Figure 4, at the time of maximum deformation, compared to those reported by Moore et al.
Strain | Region | Sepal Wall | Lateral (free) wall | ||||
---|---|---|---|---|---|---|---|
Present method | Moore et al. | Variation coef. (%) | Present method | Moore et al. | Variation coef. (%) | ||
E_{LC} | Basal | 0.02 | 0.01 ± 0.05 | 500 | -0.03 | 0.03 ± 0.05 | 166 |
Equatorial | 0.03 | 0.03 ± 0.05 | 160 | 0.05 | 0.03 ± 0.04 | 133 | |
Apical | -0.03 | 0.04 ± 0.05 | 125 | 0.02 | 0.04 ± 0.03 | 75 | |
E_{LR} | Basal | -0.06 | 0.03 ± 0.09 | 300 | -0.11 | 0.05 ± 0.07 | 140 |
Equatorial | -0.05 | 0.08 ± 0.05 | 62.5 | 0.10 | 0.07 ± 0.06 | 85.7 | |
Apical | 0.09 | 0.06 ± 0.09 | 150 | -0.09 | 0.05 ± 0.08 | 160 | |
E_{RC} | Basal | 0.12 | 0.01 ± 0.06 | 600 | 0.04 | 0.05 ± 0.07 | 140 |
Equatorial | 0.03 | 0.02 ± 0.08 | 400 | 0.19 | 0.05 ± 0.06 | 120 | |
Apical | -0.20 | 0.11 ± 0.10 | 91 | 0.26 | 0.05 ± 0.08 | 160 |
Additional file 1: Strain mapping results from single slice DENSE data. This movie demonstrates the results of the proposed strain calculation method for an extended period of time. It displays the high spatio-temporal resolution maps for all six elements of the strain tensor and helps readers clearly observe how the local Radial, Circumferential, and Longitudinal (RCL) directions are defined for all myocardial points on a single slice of the LV in the long axis view. (WMV 9 MB)
Patient study guided by LGE
Discussion
Herein, we proposed an analytical method for strain calculation from a single slice of displacement data. The planar nature of the data points imposes a near-singularity condition when calculating gradients in the direction normal to the surface, which makes these gradients highly sensitive to noise. Although meta-DENSE measures the three-dimensional displacement field from a Lagrangian curved reference surface, rather than an Eulerian flat reference plane, the non-planar configuration of the data is generally insufficient to overcome the effect of singularity. We approached this problem by assuming some physical/mathematical constraints that are shown to lead to an improvement in the accuracy of the strain calculation, such that the points with acceptable strain are abundant enough to be used for the correction of strains in other points of the LV myocardium. In other words, in each local neighborhood, we only need a sufficient number of points with acceptable strains to estimate the elements of the strain tensor in that neighborhood.
The proposed method is a multi-step process comprised of median filtering of the displacement to application of the incompressibility and temporal filtering of the calculated strain. This results in different degrees of accuracy in the estimation of each element of the deformation and strain tensors. Some parameters, such as the smoothing area or size of median filter, may depend on the acquired data characteristics like SNR and resolution. In particular, the size of the median filter was chosen as minimum in the beginning. The overall outcome for the acquired data demonstrated that there is no need to increase the size of the median filter, since this filter along with the other applied methods, described above, provides sufficient smoothness to achieve reliable calculation for in-plane components of the deformation gradient tensor.
Simulated data
Verifying the algorithm's precision over analytically simulated data reveals the accuracy level for each step of the algorithm, as well as for the corresponding deformation elements. Moreover, testing the algorithm on the simulated data helps to investigate the effect of restraining factors, such as noise, on the accuracy of the calculated strain through the suggested procedure. Simulation is also important because there is no gold standard for strain estimation based on other computational methods. Although the methods that use multi-slice imaging can provide more accurate information about out-of-plane deformation, the large slice thickness and poor resolution in the out-of-plane direction, the curvature of the LV wall, and finally the inherent noise in the acquisition process substantially decrease the accuracy in the estimation of the shear and normal strains associated with that direction. Therefore, the validity of the overall algorithm has been verified in two ways: comparison of the results to the analytical strain available for a simulated displacement data and the calculated strain of a conventional two-slice acquisition.
The algorithm performed accurately for the simulated data with RMS errors of less than 0.03 for all elements of deformation gradient tensor except for F_{ 21 }. The RMS error of 0.11, corresponding to 33.5% relative RMS error, for this in-plane element of deformation gradient tensor is related to the filtering of its high spatial changes. This problem does not pertain to our method and remains the same for 3D methods that take advantage of more than one imaging slice depending on the smoothing required. Excluding the F_{ 21 }error, RMS errors for elements in the third column of matrix Fwere greater than the average error of other two columns, as expected. The error, however, was not as high in the third column as the singularity condition would cause in the absence of the applied constraints. In terms of relative error, F_{ 23 }had a large relative RMS error of almost 100% despite the small absolute error of 0.01. This was caused by the fact that F_{ 23 }was nearly zero in our simulated deformation. This problem was more pronounced in the calculation of E_{ RC }, which was exactly zero in our simulation. Therefore, the small RMS error of 0.017 in E_{ RC }was translated to an un-bounded relative error for this component of strain.
Smoothing area and noise
The main reason for the relatively large error in calculated F_{ 21 }was the sharp spatial changes of this element in the radial direction (Figure 2, lower right panel), prescribed by our simulation. Calculation of this element over an extended neighborhood smoothed these fast non-physiological changes and resulted in considerable relative RMS difference. Nevertheless, the global map of this quantity as well as its maximum and minimum was predicted well as shown in Figure 2, lower right panel. By decreasing the smoothing area from 25 mm^{2} to 15 mm^{2} in Equation 3, the RMS and relative RMS errors of this matrix element have dropped to 0.08 and 26.1%, respectively. This reduction of the smoothing area also decreased the 34.6% relative RMS error of E_{ RL }to 27.4%. In general, the smoothing in displacement and strain decreases the sensitivity of the calculation to sharp changes of the myocardial behavior. We cannot, however, decrease the smoothing area regardless of the noise level, since the extent of the smoothing area is proportional to the noise reduction, as shown in Table 3.
For the in-plane elements of strain, Table 3 shows that adding noise with σ = 0.01 times the maximum displacement does not considerably change the error for the larger smoothing area, since the effect of smoothing over a neighborhood is dominant. Although the relative RMS error is less for a smaller smoothing area at low noise, for noise with (σ/maximum displacement) equal or greater than 0.02, the error is less for the larger smoothing area. The optimum smoothing area depends on the noise level in the images as it was set to 20 mm^{2} for the in vivo acquired data.
The relative RMS error for in-plane strains increases almost linearly with the noise standard deviation, but in a quadratic fashion for strain elements associated with the out-of-plane direction. This illustrates the higher sensitivity of these elements to noise, as explained in the Theory section.
Comparison to the strain from two adjacent slices
In the second in-vivo study the proposed method was compared to the conventional two-slice method. The measured displacement was during mid-systole; therefore the calculated circumferential strain consisted of both positive and negative values in different regions. Figure 3 demonstrates that the outcome of the proposed method is in close agreement with the results of the standard 3D method that requires the acquisition of two adjacent slices. The pattern of the strain is similar between the two methods; however, the extreme values of the strain were slightly overestimated using our proposed method. This difference can be partially explained by the effect of curvature on the strain calculation in the common 3D method. That is, the shortening and lengthening are underestimated along the chord compared to their values on the curve itself. Furthermore, it should be noted that the wall curvature erodes the actual area for which we calculate the strain from adjacent slices. We compared the results of our proposed and standard 3D methods with analytical solutions for our simulated displacement field and verified that the discrepancies are partially derived from the aforementioned shortcomings of the 3D method. Nevertheless, discrepancies can be also related to the "singularity condition" of the proposed method, which was explained in the Theory section.
Overall, it seems that our methods, which have been validated through both the two-slice and the simulated displacement data, yield an acceptable estimation of the strain for in- and out-of-plane strains from the experimental meta-DENSE data.
Comparative study
To further validate our results, we compared them to those of Moore et al. [29] and Young et al. [30]. Moore et al. estimated wall thickening to be 0.36 to 0.52, which is more consistent with available physiological data (60% ejection fraction) and clinical evidences [31, 32].
Compared to Moore et al., 17 of our 18 strain calculations were within their range of [mean ± 2σ] and 14 were in the [mean ± σ] interval. One possible reason for the observed discrepancy is the averaging of the strains over regions of the myocardial wall. For example, the peak value of apical E_{CC} on the lateral free wall, calculated by our method is -0.18 in Table 4; however, if we only focus on the epicardium, this value becomes -0.20, which no longer falls out of [mean ± σ] range.
Lessick et al. [31] observed increasing gradients of thickening from base to apex. Our method shows this clearly on the entire sepal wall as well as on most of the free wall through the increase in E_{RR}, whereas Moore's results do not show this pattern on either wall.
For shear strains, the agreement of our results with Moore's is acceptable, although a little less so for normal strains. (5 were out of the [mean ± 2σ] interval.) Note, however, that the variation coefficients (σ/mean) of shear strains are in the range of 62.5% to 600% (average: 200%) in Moore's calculations. These large variations indicate either a high inter-subject variability or a large measurement error margin.
It is worth noting that three of the five largest discrepancies, i.e. calculations out of the [mean ± 2σ] interval, relate to E_{RC}. A closer look at the data set revealed a substantial asymmetric torsional movement from the equatorial to the apical level of the heart in this particular subject. This movement results in a significant positive E_{RC} in the lateral free wall and negative E_{RC} for the sepal wall close to apex. Thus, our method detected values which may not fall in the range given by Moore et al.
Besides the peak strain values, reported in the previous studies, the temporal changes of the strain components may also be important as they can reveal additional details of cardiac function. For instance, the mean value of E_{LL} as a function of time is shown in two regions of the free wall in Figure 4. It illustrates that although the contraction profiles are almost similar in the base and apex in the first 250 ms of the cycle, the peak strain at the apex happens at a later time (400 ms) and the strain decreases in magnitude rapidly afterwards. The importance of the temporal changes further justifies our proposed algorithm (strain calculation from a single slice), which saves time on the number of acquired slices. The saved time can be utilized to improve the temporal resolution.
Clinical study
The diagnostic potential of this non-invasive method in localization of the hypo-kinetic regions of the myocardium was shown by comparing it to LGE CMR. The comparison was performed in a short axis imaging plane without any assumption on the symmetry of torsional motion of the heart. The abnormal strain is seen in a wider section around the MI region as expected in a continuum, i.e. the lack of contraction in one region extends beyond its boundary through its mutual connection to the rest of the continuum [33]. For this reason, strain measurement is not supposed to replace the LGE. It is, however, important to note that the strain map is quite useful for infarct quantification in terms of its effect on the myocardial contractility and eventually the heart function. Figure 7 shows that the proposed method is able to reveal the impairment caused by infarction in the through-plane strain. There is another small area in the anterior part of the heart with unexpected strain values for radial and longitudinal strains. This is the location where the right and left ventricles are connected and thus we may expect some deviation from the normal strain in this particular region.
Limitations of the technique
The singularity condition, explained in the Theory section, poses an instability problem for the proposed method if the noise level is high. In general, the acquisition of the high spatio-temporal displacement fields with sufficient SNR is a non-trivial process and therefore data of sufficient quality for the proposed method may be difficult to obtain routinely. This highlights the value of the methods that help to extract as much information as possible from the acquired data.
A source of error in our results could be that we assume the heart function is completely periodic and hemodynamically stable during relatively long acquisition times. This assumption may not be sufficiently accurate and could result in errors.
Finally, it should be noted that the results in the region very close to the apex suffer from the fact that RCL directions are ill-defined there. This imposes some problems when the strain values are transformed to the polar coordinate system.
Conclusion
The proposed method is able to calculate all components of the 3D Lagrangian strain tensor for systolic and diastolic phases in the left ventricle from a single-slice of 3D displacement data by assuming incompressibility of the myocardium as a physical constraint. The feasibility and the validity of the calculation of all six elements of the myocardial strain tensor from a single-slice of 3D displacements has been demonstrated. The feasibility of obtaining very favorable spatio-temporal resolution for tissue tracking by DENSE CMR from human subjects has also been shown. Furthermore, the clinical potential of the method for the evaluation of myocardial damage was demonstrated through the comparison with LGE CMR.
The method developed in this study is general and can be used with other data acquisition techniques such as zHARP or slice following 3D cine DENSE. Our work provides a framework for the detailed analysis of cardiac function. For example, it is able to demonstrate the asymmetry of strains and non-uniform contraction patterns within the ventricular wall. It is also able to calculate 3D strains from a reduced dataset, using some assumptions, which helps reduce the total scan time relative to a normal examination.
Declarations
Acknowledgements
The authors acknowledge the assistance of Dr Juan Alvergue and Dr Alexander Sassani for helping with clinical part of the study. The authors acknowledge funding support from NIH-R00 HL087614 to DBE.
Authors’ Affiliations
References
- Parmley WW, Talbot L: Heart as a pump. Handbook of Physiology: Section 2: The Cardiovascular System. Volume 1: The Heart. Edited by: Berne RM. 1979, Bethesda, Maryland: American Physiological Society, 429-460.Google Scholar
- Humphrey JD: Cardiovascular Solid Mechanics. 2002, New York: SpringerView ArticleGoogle Scholar
- Bryant DJ, Payne JA, Firmin DN, Longmore DB: Measurement of Flow with NMR Imaging Using a Gradient Pulse and Phase Difference Technique. Journal of Computer Assisted Tomography. 1984, 8: 588-593. 10.1097/00004728-198408000-00002.View ArticlePubMedGoogle Scholar
- Van Dijk P: Direct Cardiac Nmr Imaging of Heart Wall and Blood-Flow Velocity. Journal of Computer Assisted Tomography. 1984, 8: 429-436. 10.1097/00004728-198406000-00012.View ArticlePubMedGoogle Scholar
- Axel L, Dougherty L: MR imaging of motion with spatial modulation of magnetization. Radiology. 1989, 171: 841-845.View ArticlePubMedGoogle Scholar
- Zerhouni EA, Parish DM, Rogers WJ, Shapiro EP: Tagging of the Human-Heart by Magnetic-Resonance - a Non-Invasive Method for Precise Analysis of Myocardial Motion. Clinical Research. 1988, 36: A330-A330.Google Scholar
- Osman NF, Kerwin WS, McVeigh ER, Prince JL: Cardiac motion tracking using CINE harmonic phase (HARP) magnetic resonance imaging. Magnetic Resonance in Medicine. 1999, 42: 1048-1060. 10.1002/(SICI)1522-2594(199912)42:6<1048::AID-MRM9>3.0.CO;2-M.PubMed CentralView ArticlePubMedGoogle Scholar
- Aletras AH, Ding SJ, Balaban RS, Wen H: DENSE: Displacement encoding with stimulated echoes in cardiac functional MRI. Journal of Magnetic Resonance. 1999, 137: 247-252. 10.1006/jmre.1998.1676.PubMed CentralView ArticlePubMedGoogle Scholar
- Aletras AH, Balaban RS, Wen H: High-resolution strain analysis of the human heart with fast-DENSE. Journal of Magnetic Resonance. 1999, 140: 41-57. 10.1006/jmre.1999.1821.View ArticlePubMedGoogle Scholar
- Aletras AH, Tilak GS, Natanzon A, Hsu LY, Gonzalez FM, Hoyt RF, Arai AE: Retrospective determination of the area at risk for reperfused acute myocardial infarction with T2-weighted cardiac magnetic resonance imaging - Histopathological and displacement encoding with stimulated echoes (DENSE) functional validations. Circulation. 2006, 113: 1865-1870. 10.1161/CIRCULATIONAHA.105.576025.View ArticlePubMedGoogle Scholar
- Gilson WD, Yang ZQ, French BA, Epstein FH: Gadolinium-enhanced displacement-encoded MRI can simultaneously image infarct area and myocardial function in mice. Circulation. 2002, 106: 436-436.Google Scholar
- Wen H, Bennett E, Epstein N, Plehn J: Magnetic resonance imaging assessment of myocardial elastic modulus and viscosity using displacement imaging and phase-contrast velocity mapping. Magnetic Resonance in Medicine. 2005, 54: 538-548. 10.1002/mrm.20589.PubMed CentralView ArticlePubMedGoogle Scholar
- Kim D, Gilson WD, Kramer CM, Epstein FH: Myocardial tissue tracking with two-dimensional cine displacement-encoded MR imaging: Development and initial evaluation. Radiology. 2004, 230: 862-871. 10.1148/radiol.2303021213.View ArticlePubMedGoogle Scholar
- Spottiswoode BSZX, Lorenz CH, Mayosi BM, Meintjes EM, Epstein FH: 3D myocardial tissue tracking with slice followed cine DENSE MRI. Journal of Magnetic Resonance Imaging. 2008, 27: 1019-1027. 10.1002/jmri.21317.View ArticlePubMedGoogle Scholar
- Sampath S, Prince JL: Automatic 3D tracking of cardiac material markers using slice-following and harmonic-phase MRI. Magnetic Resonance Imaging. 2007, 25: 197-208. 10.1016/j.mri.2006.09.033.View ArticlePubMedGoogle Scholar
- Osman NF, Sampath S, Atalar E, Prince JL: Imaging longitudinal cardiac strain on short-axis images using strain-encoded MRI. Magnetic Resonance in Medicine. 2001, 46: 324-334. 10.1002/mrm.1195.View ArticlePubMedGoogle Scholar
- Abd-Elmoniem KZ, Stuber M, Osman NF, Prince JL: ZHARP: Three-dimensional motion tracking from a single image plane. Information Processing in Medical Imaging, Proceedings. 2005, 3565: 639-651.View ArticleGoogle Scholar
- Rogers WJ, Shapiro EP, Weiss JL, Buchalter MB, Rademakers FE, Weisfeldt ML, Zerhouni EA: Quantification of and correction for left ventricular systolic long- axis shortening by magnetic resonance tissue tagging and slice isolation. Circulation. 1991, 84: 721-731.View ArticlePubMedGoogle Scholar
- Stuber M, Scheidegger MB, Fischer SE, Nagel E, Steinemann F, Hess OM, Boesiger P: Alterations in the local myocardial motion pattern in patients suffering from pressure overload due to aortic stenosis. Circulation. 1999, 100: 361-368.View ArticlePubMedGoogle Scholar
- Abd-Elmoniem KZ, Stuber M, Prince JL: Direct three-dimensional myocardial strain tensor quantification and tracking using zHARP. Medical Image Analysis. 2008, 12: 778-786. 10.1016/j.media.2008.03.008.PubMed CentralView ArticlePubMedGoogle Scholar
- Hess AT, Zhong XD, Spottiswoode BS, Epstein FH, Meintjes EM: Myocardial 3D Strain Calculation by Combining Cine Displacement Encoding With Stimulated Echoes (DENSE) and Cine Strain Encoding (SENC) Imaging. Magnetic Resonance in Medicine. 2009, 62: 77-84. 10.1002/mrm.21984.PubMed CentralView ArticlePubMedGoogle Scholar
- Sampath S, Osman NF, Prince JL: A combined harmonic phase and strain-encoded pulse sequence for measuring three-dimensional strain. Magnetic Resonance Imaging. 2009, 27: 55-61. 10.1016/j.mri.2008.05.020.PubMed CentralView ArticlePubMedGoogle Scholar
- Young AA: Epicardial deformations from coronary cine angiograms. Theory of heart: biomechanics, biophysics, and nonlinear dynamics of cardiac function. Edited by: Glass LL, Hunter P, McCulloch A. 1991, New York: Springer-Verlag, 175-207.View ArticleGoogle Scholar
- Gonzalez RC, Woods RE: Digital image processing. 2008, Upper Saddle River, NJ: Pearson/Prentice Hall, 3Google Scholar
- Lawson CL, Hanson RJ: Solving least squares problems. 1974, Englewood Cliffs, N.J.: Prentice-HallGoogle Scholar
- Rodriguez I, Ennis DB, Wen H:: Noninvasive measurement of myocardial tissue volume change during systolic contraction and diastolic relaxation in the canine left ventricle. Magnetic Resonance in Medicine. 2006, 55: 484-490. 10.1002/mrm.20786.PubMed CentralView ArticlePubMedGoogle Scholar
- Aletras AH, Aral AE: Meta-DENSE complex acquisition for reduced intravoxel dephasing. Journal of Magnetic Resonance. 2004, 169: 246-249. 10.1016/j.jmr.2004.05.008.View ArticlePubMedGoogle Scholar
- Wen H, Marsolo KA, Bennett EE, Kutten KS, Lewis RP, Lipps DB, Epstein ND, Plehn JF, Croisille P: Adaptive postprocessing techniques for myocardial tissue tracking with displacement-encoded MR imaging. Radiology. 2008, 246: 229-240.PubMed CentralView ArticlePubMedGoogle Scholar
- Moore CC, Lugo-Olivieri CH, McVeigh ER, Zerhouni EA: Three-dimensional systolic strain patterns in the normal human left ventricle: characterization with tagged MR imaging. Radiology. 2000, 214: 453-466.PubMed CentralView ArticlePubMedGoogle Scholar
- Young AA, Kramer CM, Ferrari VA, Axel L, Reichek N: 3-Dimensional Left-Ventricular Deformation in Hypertrophic Cardiomyopathy. Circulation. 1994, 90: 854-867.View ArticlePubMedGoogle Scholar
- Lessick J, Fisher Y, Beyar R, Sideman S, Marcus ML, Azhari H: Regional three-dimensional geometry of the normal human left ventricle using cine computed tomography. Annals of Biomedical Engineering. 1996, 24: 583-594. 10.1007/BF02684227.View ArticlePubMedGoogle Scholar
- Taratorin AM, Sideman S: 3d Functional Mapping of Left-Ventricular Dynamics. Computerized Medical Imaging and Graphics. 1995, 19: 113-129. 10.1016/0895-6111(94)00038-7.View ArticlePubMedGoogle Scholar
- Ashikaga H, Mickelsen SR, Ennis DB, Rodriguez I, Kellman P, Wen H, McVeigh ER: Electromechanical analysis of infarct border zone in chronic myocardial infarction. Am J Physiol Heart Circ Physiol. 2005, 289: H1099-1105. 10.1152/ajpheart.00423.2005.PubMed CentralView ArticlePubMedGoogle Scholar
- Maeda H: Quantification of synchronous contraction of left ventricle in normal subjects using ECG-gated-SPECT images. Physiological Measurement. 2004, 25: 71-84. 10.1088/0967-3334/25/1/007.View ArticlePubMedGoogle Scholar
- Kroeker CAG, Terkeurs HEDJ, Knudtson ML, Tyberg JV, Beyar R: An Optical-Device to Measure the Dynamics of Apex Rotation of the Left-Ventricle. American Journal of Physiology. 1993, 265: H1444-H1449.Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.