A method for direct estimation of left ventricular global longitudinal strain rate from echocardiograms

We present a new method for measuring global longitudinal strain and global longitudinal strain rate from 2D echocardiograms using a logarithmic-transform correlation (LTC) method. Traditional echocardiography strain analysis depends on user inputs and chamber segmentation, which yield increased measurement variability. In contrast, our approach is automated and does not require cardiac chamber segmentation and regularization, thus eliminating these issues. The algorithm was benchmarked against two conventional strain analysis methods using synthetic left ventricle ultrasound images. Measurement error was assessed as a function of contrast-to-noise ratio (CNR) using mean absolute error and root-mean-square error. LTC showed better agreement to the ground truth strain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\varvec{R}}}^{2}=0.91)$$\end{document}(R2=0.91) and ground truth strain rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\varvec{R}}}^{2}=0.85)$$\end{document}(R2=0.85) compared with agreement to ground truth for two block-matching speckle tracking algorithms (one based on sum of square difference and the other on Fourier transform correlation; strain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\varvec{R}}}^{2}=0.70)$$\end{document}(R2=0.70), strain rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\varvec{R}}}^{2}=0.70)$$\end{document}(R2=0.70)). A 200% increase in strain measurement accuracy was observed compared to the conventional algorithms. Subsequently, we tested the method using a 53-subject clinical cohort (20 subjects diseased with cardiomyopathy, 33 healthy controls). Our method distinguished between normal and abnormal left ventricular function with an AUC = 0.89, a 5% improvement over the conventional GLS algorithms.

Pairwise Cross-Correlation. Image cross-correlation provides a statistical estimate of the translation of an image pattern between two frames. This method is used in image registration 16 , speckle tracking 17 , particle image velocimetry 18 , and image correlation 19 . The 2D discrete spatial cross-correlation between two images, I n and I n+1 , is expressed as, where R x, y is the correlation plane, (i,j) are the summation indices of the correlation, and N and M are the image height and width, respectively. The cross-correlation can be performed in the spectral domain, as, where F is the 2D Fourier transformation (FT) and F is the complex conjugate of the FT. The expanded form of Eq. (6) is written as, where (u, v) are wavenumbers proportional to positions x, y and G is the image FT.
The Fourier transform affine theorem stipulates that rotation, stretch, and shear occurs on the FT magnitude and phase 20 . We use this to establish how the affine transform affects the rigid translation estimate by replacing G n+1 in Eq. (7) with the relationship for G n using Eq. (1), Here, u ′ = a 11 u + a 21 v , v ′ = a 12 u + a 22 v , and |F| = det(F) . Equation (8) provides a correlation plane with the peak shifted from the plane center by x, y , directly related to the translations t 1 and t 2 and the local deformation gradient tensor F. The correlation peak shifts are written as, The deformation gradients a ij produce the strain captured by the GLS measurement.
Translation-invariant FT magnitude correlation for GLS estimation. The components of F can be estimated separately of − → T using the magnitude, |G(u, v)| , of the FT 20,21 . The cross-correlation of the FT magnitudes is translation invariant yielding the terms of F with no contribution from − → T . The log-polar basis Fourier-Mellin transform 22,23 , popular in image registration, decouples terms of F to estimate image rotation and stretch. (1) x n+1 y n+1 z n+1 = F x n y n z n + − → T = a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 x n y n z n + t 1 t 2 t 3 .
(9) x ∼ = a 11 t 1 + a 12 t 2 , y ∼ = a 21 t 1 + a 22 t 2 . www.nature.com/scientificreports/ Contraction and relaxation of the LV result in deformation akin to anisotropic image rescaling. By changing the FT magnitude image coordinates from cartesian (u, v) to orthogonal logarithmic coordinates logu, logv , the resulting displacements from the correlation between two FT magnitude images now correspond to horizontal u and vertical u rescaling such that, This correlation is affected by rotation and shear based on the terms present in Eq. (8). However, if these terms are minimized beforehand, the correlation peak shift estimates (�u, �v) can be related to the terms a 11 and a 22 from F through, We substitute a 11 and a 22 from Eq. (11) into F in Eq. (2) and solve for ∇u such that, Since deformation of the LV in long axis apical (ALAX) scans occurs along its length, we assume GLS occurs predominantly along the vertical direction, providing the GLS estimator, Direct global longitudinal strain estimation algorithm. We now describe our algorithm using the translation-invariant FT magnitude correlation to estimate the GLS based on Eq. (13). A schematic of our algorithm is provided in Fig. 1. The algorithm comprises two stages-the first performs an image registration to minimize shear and rotation that corrupt the correlation accuracy, while the second performs the GLS estimation.
The algorithm begins by selecting frames to analyze (Fig. 1a). Next, the user selects three points from the first frame (Fig. 1b), corresponding to approximate locations of the LV apex, the annulus septal, and annulus lateral positions. These points are tracked between consecutive frames using standard pairwise cross-correlation (Fig. 1c). For each frame, the geometric center from the tracked points is computed along with the vertical axis orientation angle for a line formed from the annulus center to the apex (Fig. 1d). Frames are then aligned based on the geometric center, and the orientation angle is corrected. A circular ROI for each frame is defined from the tracked points. This ROI is applied to filter out most of the tissue signal present from the right ventricle (RV) free wall and left atrium (LA).
In the second stage, for each pair of sequential registered images, t and t + 1 (Fig. 1e), their FT and FT magnitude is computed. The FT magnitudes are interpolated from the image grid onto a logarithm-scale grid (Fig. 1f). Because the FT logarithm transformed images are symmetric about the image diagonals, they are separated into four quarters to improve measurement accuracy. Each sub-image is filtered 24 to further minimize the influence of tissue signal that may remain from the RV free wall and LA, preserving the signal content from the LV septal wall and LV free wall, and their FT is computed (Fig. 1g). The FT sub-images are then correlated using the spectral cross-correlation kernel and ensemble-averaged to provide the displacements x, y (Fig. 1h). A dynamic phase-filtered kernel is applied to the cross-correlation to improve estimate accuracy [25][26][27] . The displacements x, y are adjusted based on the logarithm-scale grid, becoming x ′ , y ′ . The GLSr between frames is computed using y ′ and Eq. (13). Finally, GLSr across each frame pair is integrated in time using 4th-Order Runge-Kutta to obtain GLS (Fig. 1i). The integral operator provides inherent smoothing, which suppresses noise in the GLS measurements. Drift correction is applied to ensure the measurement returns to an undeformed state. We will hereafter refer to this method as the Logarithm-Transform Correlation (LTC) method.
Speckle tracking strain. This study uses two standard STE algorithms against which we benchmark our method. One algorithm uses the spatial cross-correlation kernel introduced in Eq. (5) 17 , referred to herein as the Direct Cross-Correlation or DCC method. The second uses the spectral cross-correlation (Eq. 6), hereafter referred to as the Fourier Transform Correlation or FTC method. Boundary tracking is performed by propagating the segmented boundary of the initial frame through the estimated displacement fields using 4th-Order Runge-Kutta. GLS is estimated as the measured change in arc-length between the segmented and the tracked boundary from each frame. Image co-registration is not performed, as it is not required to obtain a consistent GLS measurement based on the arc-length change calculation. Drift correction is applied to ensure the measurement returns to an undeformed state.
Artificial echocardiograms. Error analysis was performed using synthetic LV ALAX echocardiograms generated by and made publicly available from the Laboratory on Cardiovascular Imaging and Dynamics at KU Leuven 14,28 . The synthetic echocardiograms mimic images acquired from Genera Electrics (GE) Vivid E9, Hitachi-Aloka Prosound α 7 CV, Philips iE 33, Siemens SC2000, and Toshiba Artida vendor machines. The datasets provide error sources from vendor realistic noise and tissue speckle signal loss due to out of plane motion 14 . Frame rate was varied to match the specific vendor machine. Two-chamber (A2C), three-chamber (A3C), and four-chamber (A4C) LV ALAX views were provided in the dataset for one normal and four ischemic conditions. The specific occlusions include distal and proximal left anterior descending artery, left circumflex artery, and www.nature.com/scientificreports/ right coronary artery. Ground truth boundaries, displacements, and strains for each dataset were included with the synthetic images. Mean absolute error (MAE) and root-mean-square error (RMSE) for GLS and GLSr were quantified as a function of contrast to noise ratio (CNR) for all LV ALAX views under all ischemic or normal conditions. A GLS and GLSr measurement was obtained for each of the frames in every echocardiogram analyzed. A total of 4,380 measurements were examined for each of the different measurement methods. The MAE and RMSE were calculated based on a ground truth GLS value for each frame which was quantified from the average segmental strain around the LV boundary using the first frame in each echocardiogram as the reference or zero-strain frame. (e) The LV is cropped from each frame, and these sub-images are Fourier-transformed. (f) The FT magnitude is calculated, interpolated onto a logarithm-basis, and separated into four sub-images. (g) Each sub-image is Fourier-transformed and convolved with a phase filter. (h) Ensemble phase correlation is performed, producing a correlation plane with a peak shifted from the plane center. This shift corresponds to a frame pair strain rate. (i) Strain is computed by temporally integrating the strain rate estimates. www.nature.com/scientificreports/ CNR was defined as the ratio between the difference of the means and the variance between the tissue signal and the signal inside the LV 29 , where µ V and µ T are the mean of the signal intensity inside the ventricle and throughout the myocardial tissue, respectively, and σ V and σ T are the standard deviation of the signal inside the ventricle and through the myocardial tissue, respectively. Boundaries for the ventricle and myocardium are included with the ground truth data and were used here to perform the CNR calculation for each frame in the artificial echocardiograms.
Error quantities were normalized by the peak GLS or peak GLSr.
Clinical imaging. The method's clinical capabilities were demonstrated using a cohort of pediatric patients with confirmed cardiomyopathy and age-matched controls collected from a study conducted at the University of Nebraska Medical Center in Omaha, Nebraska, USA. The Institutional Review Boards of Purdue, Nebraska, and Johns Hopkins Universities each approved the study protocol. All procedures were performed in accordance with relevant guidelines and regulations. Informed written consent was obtained from study subject guardians for those under age 18 or from the subject themself for those over age 18. Each patient underwent a routine echocardiogram study on an iE33 ultrasound system (Philips Healthcare, Andover, MA, USA). Studies were collected based on the American Society of Echocardiography recommendations 30 . The 53-subject cohort included 4 patients with confirmed dilated cardiomyopathy (DCM), 16 patients with confirmed hypertrophic cardiomyopathy (HCM), and 33 age-matched controls. Information on the cohort demographics is provided in Table 1 and heart function indices in Table 2.
Doppler measurements were collected in the ALAX A4C view. B-mode ALAX A2C, A3C, and A4C view scans were performed conventionally, not explicitly collected for strain measurements. B-mode scan frame rates varied between 29 frames per second (FPS) and 100 FPS, with a median of 50 FPS. Images were stored in Digital Imaging and Communications in Medicine (DICOM) format for post-processing. Ventricle dimension measurements were computed by the Simpson biplane method using the GE EchoPAC software. Data were classified into control (CTRL) and cardiomyopathy (CM) groups.
Peak absolute GLS (GLS P ) and peak absolute systolic GLSr (GLSrs) were computed for each patient across all ALAX views. When appropriate, repeated measurements within each ALAX view were averaged. GLS P and GLSrs are reported as the average of the three standard ALAX views 6 . Patient data was similarly evaluated using commercially available software (Image-Arena Version 4.6 Build 4.6.2.12, TomTec, Germany). Only one scan www.nature.com/scientificreports/ from each ALAX view was analyzed per Nebraska's standard clinical practices. Statistical significance was tested by one-way analysis of variance (ANOVA) and Tukey Honestly Significant Difference (HSD). Receiver operating characteristic (ROC) curves and area under the curve (AUC) were computed by repeated measurement random sampling from each view to report curves with 95% confidence intervals. The ROC performance was characterized by quantifying the Youden Index and distance to corner parameters. The Youden Index (YI) is a measure of informedness, or how well informed a classification test would be based on predictions, defined as with a range of 0 to 1, where 0 is no predictive ability and 1 is perfect predictive ability 31 . Distance to corner is a measure of the ROC curve at its optimal threshold to a perfect prediction at a sensitivity of 1, or perfectly predicting true positive class, and 1-specificity of 0, or perfectly predicting the true negative class.

Results
Error analysis results.  www.nature.com/scientificreports/ GLSr = 0.95s −1 ). LTC showed a 1.5 to 2-fold improvement in accuracy than DCC and a 1.5-fold improvement compared to FTC. Additionally, the LTC method is unaffected by CNR. In contrast, the FTC and DCC methods show a CNR dependence, albeit less for the FTC method. GLS estimates (Fig. 2b-1) yielded a linear regression fit with slope and bias of m DCC =0.60 and b DCC = 0.30% for the DCC method, m FTC =0.65 and b FTC = 0.33% for the FTC method, and m LTC = 0.92 , b LTC = 0.09% for the LTC method. Regression fit qualities of R 2 DCC = 0.71 , R 2 FTC = 0.74 , and R 2 LTC = 0.91 were measured. The LTC method shows more than 200% improvement in measurement accuracy compared against the DCC and FTC methods as a function of CNR (Fig. 2b-2), where values are normalized by GLS = 8.47%. The error analysis demonstrates that the LTC method is unaffected by CNR, while the DCC and FTC methods are affected by signal quality.

Clinical analysis results. Results comparing the LTC method against conventional GLS methods and
TomTec are presented in Fig. 3a-1,b-1. The LTC median GLS P are higher compared to the conventional methods (GLS P,LTC-CTRL = 15.84% , GLS P,LTC-CM = 10.00% ; GLS P,FTC-CTRL = 9.59% , GLS P,FTC-CM = 6.25% ; GLS P,DCC-CTRL = 8.07% , GLS P,DCC-CM = 5.47% ), but lower compared to TomTec (GLS P,TT-CTRL = 18.19% , GLS P,TT-CM = 13.86% ). Significance tests by ANOVA indicated the GLS P group means are significant (F-value = 90.02; p < 1 × 10 −5 ). Post-hoc analysis by Tukey HSD indicated the group means between health states were significant ( p < 0.05) for all methods. Furthermore, the LTC method was statically significant ( p < 0.05) compared to the conventional methods for both health states. Differences were not statistically significant between the FTC and DCC methods.
Similarly, for the GLSrs measurements, the LTC medians are higher than compared to the conventional methods (GLSrs LTC-CTRL = 0.75s −1 , GLSrs ,LTC-CM = 0.54s −1 ; GLSrs FTC-CTRL = 0.51s −1 , GLSrs ,FTC-CM = 0.38s −1 ; GLSrs LTC-CTRL = 0.50s −1 , GLSrs ,LTC-CM = 0.36s −1 ) but lower compared to TomTec (GLSrs LTC-CTRL = 0.86s −1 , Figure 3. Distribution of measurements and significance tests for each GLS measurement method on observing (a-1) peak absolute GLS (GLS P ) and (b-1) peak absolute systolic GLSr (GLSrs). (2) Receiver operating characteristic (ROC) curves displaying the ability for the LTC method estimated parameters to distinguish between normal and abnormal cardiac disease states based on (a-2) GLS P and (b-2) GLSrs. The analysis was performed on a set of control subjects (CTRL) and subjects with cardiomyopathies (CM). The LTC method was compared with the conventional STE methods, FTC and DCC, as well as with the commercially available TomTec software. The bounded regions about the ROC curve provide the 95% confidence interval over which sensitivity and specificity are measured, contingent on the beat records and scans analyzed within subjects. www.nature.com/scientificreports/ GLSrs ,LTC-CM = 0.69s −1 ). Significance tests by ANOVA indicated the GLSrs group means are significant (F-value = 44.35; p < 1 × 10 −5 ). Tukey HSD post-hoc analysis indicated the group means between health states were significant ( p < 0.05) for the LTC and TomTec methods, but were not significant for the conventional methods. Moreover, the LTC method was statically significant ( p < 0.05) compared to the conventional methods for both health states. Results were not significant between the FTC and DCC methods. ROC curves are presented in Fig. 3a-2, b-2. The GLS P ROC curves (Fig. 3a-2) show AUCROC values of GLS P,LTC-AUCROC = 0.89 ; GLS P,FTC-AUCROC = 0.83 ; GLS P,DCC-AUCROC = 0.84 ; and GLS P,TT-AUCROC = 0.86 . The LTC ROC confidence interval indicated this method at the low end performed as well as TomTec when classifying patients, but could improve classification by 10% based on the beats analysed. The FTC and DCC methods show increased confidence intervals, but consistently performs worse than the LTC method. The LTC ROC performance was marginally better than the other methods based on performance metrics (YI LTC = 0.79 , CD LTC = 0.22 ; YI LTC = 0.42 , CD LTC = 0.36 ; YI DCC = 0.54 , CD DCC = 033 ; YI TT = 0.72 , CD TT = 0.25).
The GLSrs ROC curves (Fig. 3 b-2) show the LTC AUCROC is consistently higher compared to the other methods (GLSrs LTC-AUCROC = 0.84 ; GLSrs FTC-AUCROC = 0.77 ; GLSrs DCC-AUCROC = 0.80 ; GLSrs TT-AUCROC = 0.76 ). The confidence interval for the LTC method was reduced compared to the conventional strain methods and shows at the low end that LTC is comparable to TomTec but can offer a 15% improvement overall. Additionally, the LTC ROC performance was better than the other methods based on performance metrics (YI LTC

Discussion
This study presents a new algorithm, the LTC method, for computing GLS and GLSr estimates from ultrasound scans. Error analysis using synthetic ultrasound images quantified the LTC method's improvement over two conventional STE methods. A clinical cohort was analyzed using the LTC method, the two conventional STE methods, and a commercial software method. The LTC method does not rely on LV shape assumptions and avoids the use of boundary segmentation. Furthermore, because LV segmentation is not required, regularization to preserve the segmentation shape is avoided. Moreover, the entire image of the LV is used to compute strain which minimizes out-of-plane motion correlation loss.
These claims are supported by basic principles derived in particle image velocimetry (PIV), another FFTbased cross correlation application. In PIV, fluid tracing particle motion is measured between frames 32 . Outof-plane motion in PIV images, which are speckle-like in nature, is a function of sample volume and out-ofplane velocity gradients 33 . PIV frame rates are high (greater than 100 FPS) and the illumination volume is thin (approximately 1 mm). In B-mode imaging, the sample volume is thicker (upwards of 10 mm) and frame rates are sufficiently high (more than 30 FPS) relative to cardiac tissue velocity gradients (on the order of 100 mm/s, or roughly 3 mm/frame at the lowest frame rate settings). Thus, speckle pattern losses in the LTC should not be significant, more so with additional signal content provided from the whole LV image.
The LTC method is novel by directly computing the GLSr between sequential frames, ensuring a reliable rate measurement. Computing GLS by integrating the GLSr provides a smoothing operation that reduces noise. Commercial GLS algorithms have constraints that enforce tracked boundaries that are smooth in space and time, providing measurements that appear physically consistent, but become a function of the regularization process 10,34,35 . Thus, the GLS measurements do not adhere to the underlying deformation that occurs between frames.
The error analysis presented in Fig. 2 demonstrated that the DCC and FTC methods are dependent on CNR. As CNR approaches 1, the mean pixel intensity of the LV wall and the mean pixel intensity of speckle noise are almost equal. This means the speckle noise cannot be differentiated from the LV wall, making the physical features ambiguous. Both methods underestimated the GLS, supporting the notion that it may be best to avoid such computation 6 . In contrast, the LTC method enables robust GLS computation even with noisy images, as supported by the MAE and RMSE plots in Fig. 2.
For demonstration purposes, example images from two test subjects' clinical data are provided to show the quality differences that affect CNR, provided in Fig. 4. These subjects are of the same gender and similar ages. In, one scan (Fig. 4a) the LV interior is clear of noise and the myocardium can be clearly delineated, providing a CNR of greater than 8 dB. The other scan (Fig. 4b) has increased noise within the LV interior and signal loss in the LV myocardium, resulting in a CNR of less than 2 dB. These differences for the latter case can greatly impact conventional STE method measurement accuracy.
Comparison of each method for all health states, presented in Fig. 3a-1,b-1, provides the measurement distributions and their statistical significance. For each method, the variance of the CTRL and CM distributions overlapped while statistical significance of the means was observed. Between methods, statistical significance was observed for the LTC method compared to the conventional STE methods. These results indicate larger GLS P and GLSrs values for normal function than for abnormal function, but establishing improvement of the LTC method through this analysis alone is not possible.
The ROC curves, presented in Fig. 3a-2,b-2, are used to determine if clinical separation is possible. The LTC method GLS P ROC curve showed marginal classification improvement compared to the conventional STE methods and the commercially available method. However, GLSrs ROC curve shows a nearly 10% improvement in classification. Both tested parameters show improved correct diagnosis rates.
The LTC method GLS P and GLSrs measurements were below nominal ranges reported from literature (peak GLS > 18%; peak GLSrs > 1 s −1 ) 36 . Commercial methods rely on regularization steps, which force the measurements to fit in LV shape models 6, 10 , thereby causing the GLS measurements to be a function of the regularization instead of the actual underlying deformation, possibly leading to overestimation 37,38 .
LTC algorithm limitations stem from the assumption that GLS and GLSr can be reliably measured from the septal and lateral walls, ignoring shortening near the apex and possible effects of signal loss, attenuation, www.nature.com/scientificreports/ and phase aberration along the myocardium. Furthermore, the algorithm assumes that reorientation has been performed correctly with no misalignment. If misalignment exists, the measurements will be a combination of GLS and off-axis strain. Out-of-plane motion in the artificial dataset is modelled, but may cause speckle signal dropout, requiring a correction step to mitigate gaps in the images 14 . Thus, a correction step is used which may minimize out-of-plane motion and reduce influence of this error component. Frame rate is also known to affect strain measurement accuracy, with a typical operating range of 50-80 FPS 9,10 . The LTC method should be evaluated in future work to determine if this optimal range is suitable for application or if lower or higher frame rates must be maintained for accurate measurements. Robustness of the LTC correlation to noise and out-ofplane correlation loss should be tested in greater detail in a future study, however this does require rigorous and controlled test conditions which is not trivial. Finally, this study was limited by its clinical validation, which was performed using retrospective pediatric data. Pediatric echocardiograms are typically of good quality, without significant attenuation, clutter noise, or phase aberration which are common in adult echocardiograms for patients with cardiovascular disease. A dataset optimized for GLS estimation would help further substantiate findings. Future work will apply the LTC algorithm to study patient populations with ischemic heart disease and cardiotoxicity-related heart failure. We presented a new correlation kernel, the logarithm transform correlation or LTC, for quantifying GLS and GLSr from echocardiography scans. Our LTC-based algorithm does not use LV shape assumptions, is machineagnostic, automated, and free of heuristic inputs. We compared the LTC against STE algorithms using artificial scans, analyzing error against ground truth GLS and GLSr values, and validated using clinical data from a study of pediatric cardiomyopathies. Results showed that the LTC method is unaffected by the image quality, providing improved measurement accuracy against the STE methods for both the synthetic data and clinical cohort. www.nature.com/scientificreports/