Improve myocardial strain estimation based on deformable groupwise registration with a locally low-rank dissimilarity metric

Abstract

Background

Current mainstream cardiovascular magnetic resonance-feature tracking (CMR-FT) methods, including optical flow and pairwise registration, often suffer from the drift effect caused by accumulative tracking errors. Here, we developed a CMR-FT method based on deformable groupwise registration with a locally low-rank (LLR) dissimilarity metric to improve myocardial tracking and strain estimation accuracy.

Methods

The proposed method, Groupwise-LLR, performs feature tracking by iteratively updating the entire displacement field across all cardiac phases to minimize the sum of the patchwise signal ranks of the deformed movie. The method was compared with alternative CMR-FT methods including the Farneback optical flow, a sequentially pairwise registration method, and a global low rankness-based groupwise registration method via a simulated dataset (n = 20), a public cine data set (n = 100), and an in-house tagging-MRI patient dataset (n = 16). The proposed method was also compared with two general groupwise registration methods, nD + t B-Splines and pTVreg, in simulations and in vivo tracking.

Results

On the simulated dataset, Groupwise-LLR achieved the lowest point tracking errors (p = 0.13 against pTVreg for the temporally averaged point tracking errors in the long-axis view, and p < 0.05 for all other cases), voxelwise strain errors (all p < 0.05), and global strain errors (p = 0.05 against pTVreg for the longitudinal global strain errors, and p < 0.05 for all other cases). On the public dataset, Groupwise-LLR achieved the lowest contour tracking errors (all p < 0.05), reduced the drift effect in late-diastole, and preserved similar inter-observer reproducibility as the alternative methods. On the patient dataset, Groupwise-LLR correlated better with tagging-MRI for radial strains than the other CMR-FT methods in multiple myocardial segments and levels.

Conclusions

The proposed Groupwise-LLR reduces the drift effect and provides more accurate myocardial tracking and strain estimation than the alternative methods. The method may thus facilitate a more accurate estimation of myocardial strains for clinical assessments of cardiac function.

Background

Myocardial strain is a sensitive measure of the left ventricular mechanical function, which can be useful to aid the diagnosis and prognosis of various diseases [1,2,3]. Cardiovascular magnetic resonance-feature tracking (CMR-FT) estimates myocardial strains from routinely performed cine MRI images using feature tracking algorithms [4]. Compared with the reference standard tagging-MRI, CMR-FT does not need additional time-consuming sequences and facilitates retrospective strain analysis [5]. However, it is well known that CMR-FT based on regular cine does not estimate regional strain accurately [6, 7]. Since regional strain is often hypothesized to be more sensitive to the early progresses of cardiac diseases, a more accurate estimation of the regional strain with advanced feature tracking algorithms is highly desirable.

Current commercial software tools for CMR-FT often use optical flow or pairwise registration for tracking of the myocardium across adjacent cardiac phases [8]. However, studies have shown that these frame-by-frame tracking methods often cause an accumulation of tracking errors through the early frames. These errors become greater for later frames, which eventually cause inaccurate feature tracking and strain estimation, often referred to as the drift effect [9]. Several methods have been proposed to address this challenge. One method is based on drift compensation, which approximates the drift by a linear function and subtracts the trend from the estimated strain [10]. However, since the linearity of the drift is assumed only for simplicity, the compensation is highly empirical and does not truly eliminate the accumulative tracking errors. Another type of pairwise CMR-FT methods seeks to directly align every frame to the end-diastolic frame [9, 11]. Although such methods can avoid the drift effect, large deformations between distant frames such as end-diastole and end-systole can cause even more severe tracking errors [12]. In addition, there are also attempts to mitigate the errors by combining forward and backward registration passes [13,14,15]. However, the tracking errors caused by pairwise alignment still exist in both passes and are not guaranteed to be offset when combined. Recently, groupwise registration has shown the potential to reduce the drift effect by simultaneously aligning all cine images to a common reference frame. For example, Metz et al. proposed a deformable groupwise registration method for tissue tracking, which used a sum of the voxelwise signal variances as the dissimilarity metric since a well-aligned set of images should have a low signal variation across time [16]. Qiao et al. implemented a similar groupwise registration method for myocardial tracking via cine MRI images and demonstrated an improved tracking accuracy over pairwise registration [17]. In addition to signal variance, signal low-rankness across the spatiotemporal domain can also be used as a dissimilarity metric for groupwise registration. If the registration is perfect and the temporal signal variation is spatially consistent, the aligned images should encompass a low-rank property, which can be used to guide the registration. For example, Peng et al. decomposed the spatiotemporal images into an error term, which should be sparse, and a set of well-aligned images, which should have the low-rank property [18]. They then combined the sparsity and low-rankness into a single dissimilarity metric to guide the groupwise registration, and found improved accuracy and efficiency when compared with other registration methods. Haase et al. extended this method by introducing nonparametric deformation and total variation regularization, and found improved accuracy compared with alternative groupwise registration methods for several motion tracking applications [19]. However, none of these studies on groupwise registration evaluated the accuracy and reproducibility for strain estimation.

In this work, we sought to develop a strain estimation method based on deformable groupwise registration with a locally low-rank (LLR) dissimilarity metric and evaluate its accuracy and reproducibility. Previous myocardial tracking studies have shown that global low-rankness can be used to guide groupwise registration [18, 19]; however, it is well-known that signals in cine MRI are often locally low-rank rather than globally low-rank [20]. Different tissues in cine MRI may encompass signal variations of different patterns, and thus the global rank can be relatively high even after motion compensation. For example, the myocardial intensity may present a periodic variation across the cardiac cycle due to through-plane motion and the fresh spin effect. The blood signals may also change irregularly due to turbulent flows. The banding and flow artifacts commonly observed in cine MRI can disturb regional signals in the image. Therefore, Locally Low Rankness (LLR) may be a better metric to measure the alignment of cine images over the cardiac cycle. The baseline methods for evaluation included three alternative CMR-FT methods (optical flow, pairwise registration, and global low-rankness (GLR)-based groupwise registration) and two general groupwise registration methods (nD + t B-Splines [16] and pTVreg [21]). We compared our methods with all the baseline methods in terms of myocardial tracking accuracy via a simulated dataset (XCAT) [22] and a public cine dataset (ACDC) [23]. We further compared our method with the alternative CMR-FT methods in terms of strain estimation accuracy and reproducibility via the simulated dataset, public cine dataset, and an in-house tagging-MRI patient dataset.

Methods

Groupwise registration formulation

The goal of groupwise registration is to simultaneously align all cine images to a common reference frame. Let $\:\text{f}\left(\mathbf{x},\text{t}\right)$ represent the cine image sequence, where the spatial coordinate $\:\mathbf{x}\in\:{\Omega\:}=\left\{\left(\text{x},\text{y}\right)\right|1\le\:\text{x}\le\:{\text{N}}_{\text{x}},1\le\:\text{y}\le\:{\text{N}}_{\text{y}}\}$ and the temporal coordinate $\:1\le\:\text{t}\le\:{\text{N}}_{\text{t}}$. $\:{\text{N}}_{\text{x}}$, $\:{\text{N}}_{\text{y}}$, and $\:{\text{N}}_{\text{t}}$ represent the pixel number in the x direction, pixel number in the y direction, and frame number, respectively. The deformed cine image sequence can be expressed as

$$\:\stackrel{\sim}{\text{f}}\left(\mathbf{x},\text{t}\right)=\text{f}\left(\mathbf{T}\left(\mathbf{x},\text{t}\right),\text{t}\right)=\text{f}\left(\mathbf{x}+\mathbf{d}\left(\mathbf{x},\text{t}\right),\text{t}\right),$$

(1)

where $\:\mathbf{T}\left(\mathbf{x},\text{t}\right)=\mathbf{x}+\mathbf{d}\left(\mathbf{x},\text{t}\right)$ is the transformation field and $\:\mathbf{d}\left(\mathbf{x},\text{t}\right)$ is the displacement field. We use the free-form deformation modeling of the displacement field [24]. Specifically, given a control point mesh $\bf{\Phi} ({\rm{i}},{\rm{j}},{\rm{t}})$ with uniform control point spacing $\:{\updelta\:}$, $\:\mathbf{d}\left(\mathbf{x},\text{t}\right)$ is modeled by a linear combination of control points with cubic B-spline coefficients as

$${\bf{d}}\left( {{\bf{x}},{\rm{t}}} \right) = \sum\limits_{{\rm{k}} = 0}^3 {} \sum\limits_{{\rm{l}} = 0}^3 {{{\rm{B}}_{\rm{k}}}\left( {\rm{u}} \right){{\rm{B}}_{\rm{l}}}\left( {\rm{v}} \right)\bf{\Phi} \left( {{\rm{i}} + {\rm{k}},{\rm{j}} + {\rm{l}},{\rm{t}}} \right)} ,$$

(2)

where $\:\text{i}=\lfloor\:(\text{x}-1)/{\updelta\:}\rfloor,\:$ $\:\text{j}=\lfloor\:(\text{y}-1)/{\updelta\:}\rfloor\:$, $\:\text{u}=(\text{x}-1)/{\updelta\:}-\lfloor\:(\text{x}-1)/{\updelta\:}\rfloor\:$, $\:\text{v}=(\text{y}-1)/{\updelta\:}-\lfloor\:(\text{y}-1)/{\updelta\:}\rfloor\:$, and $\:{\text{B}}_{\text{k}}\left(\text{u}\right)$ is the $\:\text{k}$-th basis function of the B-splines. By constraining the temporal average of the transformation or displacement field to be an identity map or zero, we register all images to an average common reference frame:

$${1 \over {{{\rm{N}}_{\rm{t}}}}}\sum\limits_{{\rm{t}} = 1}^{{{\rm{N}}_{\rm{t}}}} {{\bf{T}}({\bf{x}},{\rm{t}}) = {\bf{x}}\>\> \Leftrightarrow } \>\>{1 \over {{{\rm{N}}_{\rm{t}}}}}\sum\limits_{{\rm{t}} = 1}^{{{\rm{N}}_{\rm{t}}}} {{\bf{d}}({\bf{x}},{\rm{t}}) = 0,}$$

(3)

Which based on Eq. (2) is equivalent to

$${1 \over {{{\rm{N}}_{\rm{t}}}}}\sum\limits_{{\rm{t}} = 1}^{{{\rm{N}}_{\rm{t}}}} {\bf{\Phi} ({\rm{i}},{\rm{j}},{\rm{t}}) = 0.}$$

(4)

GLR and LLR dissimilarity metrics

To register all images in $\:\stackrel{\sim}{\text{f}}\left(\mathbf{x},\text{t}\right)$, we can use either the previously described GLR metric [18, 19] or the proposed LLR metric. To use the GLR metric, one firstly reformulates the entire series of images into a Casorati matrix:

$$\bf {\rm{\bf{C}}} = \left[ {\matrix{{{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{1}}}{\rm{,1}}} \right)} & {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{1}}}{\rm{,2}}} \right)} & \cdots & {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{1}}}{\rm{,}}{{\rm{N}}_{\rm{t}}}} \right)} \cr {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{2}}}{\rm{,1}}} \right)} & {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{2}}}{\rm{,2}}} \right)} & \cdots & {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{2}}}{\rm{,}}{{\rm{N}}_{\rm{t}}}} \right)} \cr \vdots & \vdots & \ddots & \vdots \cr {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{L}}}{\rm{,1}}} \right)} & {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{L}}}{\rm{,2}}} \right)} & \cdots & {{\rm{\tilde f}}\left( {{{\rm{\bf{x}}}_{\rm{L}}}{\rm{,}}{{\rm{N}}_{\rm{t}}}} \right)} \cr } } \right] \in \>{\mathbb{R}^{{\rm{L \times }}{{\rm{N}}_{\rm{t}}}}},$$

(5)

where $\:\text{L}={\text{N}}_{\text{x}}{\text{N}}_{\text{y}}$ is the total number of voxels, and the $\:\text{l}$-th row holds the signal in the $\:\text{l}$-th voxel. GLR assumes the Casorati matrix is rank-deficient if all images are well aligned. The GLR dissimilarity is thus determined by the nuclear norm of $\:\mathbf{C}$:

$${{\rm{D}}_{{\rm{GLR}}}} = \parallel{\bf{C}}{\parallel _{\rm{*}}} = \sum\limits_{{\rm{t}} = 1}^{{{\rm{N}}_{\rm{t}}}} {{{\rm{\sigma }}_{\rm{t}}},}$$

(6)

where $\:{{\upsigma\:}}_{\text{t}}$ is the $\:\text{t}$-th singular value of $\:\mathbf{C}$.

In cardiac cine imaging, however, signals may be spatially inconsistent even after registration, due to a set of confounders such as motion, blood flow, and artifacts. As a result, the post-registration global rank may be still high, rendering GLR less effective for the guidance of the registration. Figure 1 illustrates this problem using a realistic cine movie example. From the illustration, it is clear that the local signal rank is more sensitive to the elimination of motion than the global signal rank. For example, the rank for voxels surrounding the myocardium border is reduced significantly by the registration. Therefore, we propose to use the LLR metric, which uniformly partitions the spatial domain into $\:{\text{N}}_{\text{c}}$ overlapped square patches with a prespecified patch size and spacing, reformulates the patches into local Casorati matrices, and computes the sum of the nuclear norm of these local Casorati matrices:

$${{\rm{D}}_{{\rm{LLR}}}} = \sum\limits_{{\rm{k = 1}}}^{{{\rm{N}}_{\rm{c}}}} {\parallel {{\bf{C}}_{\rm{k}}}{\parallel _{\rm{*}}} = } \sum\limits_{{\rm{k = 1}}}^{{{\rm{N}}_{\rm{c}}}} {} \sum\limits_{{\rm{t = 1}}}^{{{\rm{N}}_{\rm{t}}}} {{{\rm{\sigma }}_{{\rm{k,t}}}}} ,$$

(7)

where $\:{\bf{C}}_{\rm{k}}$ is the local Casorati matrix associated with the $\:\text{k}$-th patch, and $\:{\rm{\sigma }}_{{\rm{k,t}}}$ is the $\:\text{t}$-th singular value of $\:{\bf{C}}_{\rm{k}}$.

Fig. 1

The maps of LLR cost—a surrogate of the patchwise signal rank—of a cine image sequence in the ACDC dataset before (Row 1) and after (Row 2) registration. Row 3 shows the relative difference between Row 1 and Row 2. Note that the rank reduction is highly region-dependent, with high reductions (~ 40%) in regions with large motion, such as the myocardial borders. Some regions have a high rank even after image registration, such as the blood pool. GLR, which is equivalent to LLR with a patch size of 128 (the rightmost column), assumes the rank is globally low after image registration. However, as the rightmost column shows, the post-registration global rank is still high, resulting in only a ~ 10% rank reduction

Regularization

To enforce the spatiotemporal regularity of the deformation field, we use a combination of spatial and temporal regularization terms. For the spatial smoothness, we introduce a bending-energy-based regularization term over the spatial domain [24]:

$$\eqalign{{{\rm{R}}_{{\rm{spatial}}}} & = \sum\limits_{{\rm{t}} = 1}^{{{\rm{N}}_{\rm{t}}}} {} \sum\limits_{{\bf{x}} \in \>\Omega } {} \left\| {{{\partial {\>^2}{\bf{T}}} \over {\partial \>{{\rm{x}}^2}}}({\bf{x}},{\rm{t}})} \right\|_2^2 \cr & + 2\left\| {{{\partial {\>^2}{\bf{T}}} \over {\partial \>{\rm{x}}\partial \>{\rm{y}}}}({\bf{x}},{\rm{t}})} \right\|_2^2 + \left\| {{{\partial {\>^2}{\bf{T}}} \over {\partial \>{{\rm{y}}^2}}}({\bf{x}},{\rm{t}})} \right\|_2^2. \cr}$$

(8)

For the temporal smoothness, we introduce a second-order regularization term over the temporal domain

$${{\rm{R}}_{{\rm{temporal}}}} = \sum\limits_{{\rm{t}} = 1}^{{{\rm{N}}_{\rm{t}}}} {} \sum\limits_{{\bf{x}} \in \>\Omega } {} \left\| {{{\partial {\>^2}{\bf{T}}} \over {\partial \>{{\rm{t}}^2}}}({\bf{x}},{\rm{t}})} \right\|_2^2.$$

(9)

We approximate the partial derivates in $\:{\text{R}}_{\text{s}\text{p}\text{a}\text{t}\text{i}\text{a}\text{l}}$ and $\:{\text{R}}_{\text{t}\text{e}\text{m}\text{p}\text{o}\text{r}\text{a}\text{l}}$ by the finite differences. Furthermore, given the cyclic motion of the heart, we adopt cyclic finite differences in $\:{\text{R}}_{\text{t}\text{e}\text{m}\text{p}\text{o}\text{r}\text{a}\text{l}}$ to enforce the consistency of the motion between the first frame and the last frame.

As a result, the proposed LLR-based groupwise registration (Groupwise-LLR) is formulated as

$$\begin{aligned}\:\underset{\varvec{\upvarphi\:}}{\text{min}}{\text{D}}_{\text{L}\text{L}\text{R}} & +{\uplambda\:}{\text{R}}_{\text{s}\text{p}\text{a}\text{t}\text{i}\text{a}\text{l}}+{\upmu\:}{\text{R}}_{\text{t}\text{e}\text{m}\text{p}\text{o}\text{r}\text{a}\text{l}} \\& \quad \text{s}.\text{t}.\frac{1}{{\text{N}}_{\text{t}}}\sum\:_{\text{t}=1}^{{\text{N}}_{\text{t}}}\varvec{\upvarphi\:}(\text{i},\text{j},\text{t})=0,\end{aligned}$$

(10)

where $\:{\uplambda\:}$ and $\:{\upmu\:}$ are the regularization coefficients, and the constraint is from Eq. 4. Inspired by a previous study on groupwise registration [25], we implemented the projected gradient descent algorithm to solve this problem in a coarse-to-fine multi-resolution framework. In this framework, the proposed method firstly estimates the deformation at the lowest resolution level to roughly align the images. Then, the method gradually refines the deformation at higher resolution levels to achieve more accurate alignment. We initialize the control point mesh $\bf{\Phi} ({\rm{i}},{\rm{j}},{\rm{t}})$ as zero, leading to zero displacements across all cardiac phases.

Strain estimation

After obtaining the temporally resolved transformation field, we estimate the transformation field from the first image (end-diastole) to each later image by a two-step process, in which the coordinate $\:\mathbf{x}$ in the end-diastolic image is firstly mapped to the common reference frame, and then to the $\:\text{t}$-th image. Mathematically, this process is represented by

$${{\bf{T}}_{1 \to {\rm{t}}}}\left( {\bf{x}} \right) = {\bf{T}}\left( {{{\bf{T}}^{ - 1}}\left( {{\bf{x}},1} \right),{\rm{t}}} \right),$$

(11)

where ${{\bf{T}}_{1 \to {\rm{t}}}}\left( {\bf{x}} \right)$ represents the transformation field from the first image to the $\:\text{t}$-th image. Based on the estimated ${{\bf{T}}_{1 \to {\rm{t}}}}\left( {\bf{x}} \right)$, the Green-Lagrange strain tensor [26] is evaluated by

$${\bf{E}}\left( {{\bf{x}},{\rm{t}}} \right) = {1 \over 2}\left[ {{{\left( {\nabla {\>_{\bf{x}}}{{\bf{T}}_{1 \to {\rm{t}}}}\left( {\bf{x}} \right)} \right)}^{\rm{T}}}\nabla {\>_{\bf{x}}}{{\bf{T}}_{1 \to {\rm{t}}}}\left( {\bf{x}} \right) - {\bf{I}}} \right],$$

(12)

where $\:{\nabla\:}_{\mathbf{x}}$ represents the Jacobian operator about $\:\mathbf{x}$ (implemented using finite differences) and $\:\mathbf{I}$ is an identity matrix of size two. The strain along a certain direction is then computed by

$$\:{\text{e}}_{\mathbf{u}}\left(\mathbf{x},\text{t}\right)={\left(\mathbf{u}\left(\mathbf{x}\right)\right)}^{\text{T}}\mathbf{E}\left(\mathbf{x},\text{t}\right)\mathbf{u}\left(\mathbf{x}\right),$$

(13)

where $\:\mathbf{u}\left(\mathbf{x}\right)$ is the strain direction field determined by segmentation of the left ventricular myocardium in the end-diastolic image. For those long-axis slices, $\:\mathbf{u}\left(\mathbf{x}\right)$ represents the longitudinal direction; for those short-axis slices, $\:\mathbf{u}\left(\mathbf{x}\right)$ represents the circumferential or radial direction. The global longitudinal strain (GLS) is estimated by averaging the strain values over the left ventricular myocardium within a single long-axis slice. The global circumferential strain (GCS) and global radial strain (GRS) are estimated by averaging the strain values over the left ventricular myocardium within a specified slice (e.g., base, mid-ventricle, apex) or across multiple slices covering the entire left ventricle.

Baseline methods

The proposed Groupwise-LLR was compared with the Farneback optical flow [27], a commercially available pairwise registration method [28], and the GLR-based groupwise registration (Groupwise-GLR) based on both simulated and in-vivo data. Among these alternative CMR-FT methods, the optical flow and pairwise registration estimate the transformation field between each pair of the neighboring images, i.e., ${{\bf{T}}_{{\rm{t}} - 1 \to {\rm{t}}}}\left( {\bf{x}} \right)$, $\:2\le\:\text{t}\le\:{\text{N}}_{\text{t}}$. The transformation field ${{\bf{T}}_{1 \to {\rm{t}}}}\left( {\bf{x}} \right)$ is then a composition of the neighboring transformation fields. Groupwise-GLR was implemented in the same way as Groupwise-LLR except for the use of a GLR dissimilarity metric. The proposed method was also compared with two general groupwise registration methods, nD+t B-Splines [16] and pTVreg [21], in simulations and in vivo tracking. Among them, nD + t B-Splines uses an nD + t B-spline deformation model and a sum of the voxelwise signal variances as the dissimilarity metric. pTVreg uses a similar LLR metric as the proposed one. However, pTVreg extracts nonoverlapped patches to compute the LLR metric and fixes the patch size for all resolution levels. In comparison, Groupwise-LLR partitions the spatial domain into overlapped patches to enhance the continuity of the estimated motion among neighboring patches, and adjusts the patch size along with the image resolution to improve the consistency of the multi-resolution optimization. Additionally, pTVreg is characterized by the total variation regularization, which is different from the second-order Tikhonov regularization used in this work. Furthermore, it should be emphasized that the proposed Groupwise-LLR is the first strain estimation method based on groupwise registration while the two baseline groupwise registration methods have not been applied for strain estimation before.

Implementation details

All methods were implemented with MATLAB (R2022a, MathWorks, Natick, MA, USA). The parameters of each method were manually tuned for the simulated and in vivo experiments. For the proposed Groupwise-LLR, the number of resolution levels was 3 and the downsampling factor between two adjacent resolution levels was 2. The regularization coefficients were not reweighted between different resolution levels. The cubic interpolation was used for image downsampling and the linear interpolation was used for image warping. The patch size and spacing used to sample the local Casorati matrices were 5 and 3 pixels at the lowest resolution level and doubled each time moving to the next level. The control point spacing, spatial regularization coefficient $\:{\uplambda\:}$, and temporal regularization coefficient $\:\mu\:$ were 6 and 7 pixels, 6х10^− 4 and 1х10^− 3, 0.06 and 0.1 for the simulated dataset and two in vivo datasets, respectively. We determined the regularization coefficients using a semi-quantitative approach based on small separate trial datasets. For the simulated experiments, we firstly set $\:{\uplambda\:}$ = 5х10^− 4 based on a manual inspection of the spatial smoothness of the displacement fields. We then varied the ratio $\:\mu\:/{\uplambda\:}$ from 0.1 to 1000, and found that $\:\mu\:/{\uplambda\:}$ = 100 minimized the EPE over the simulated trial dataset. Thus, we fixed $\:\mu\:/{\uplambda\:}$ = 100 hereafter. We varied $\:{\uplambda\:}$ from 2х10^− 4 to 1х10^− 3 with a step of 2х10^− 4, and found that $\:{\uplambda\:}$ = 6х10^− 4 minimized the EPE over the simulated trial dataset. For the in-vivo experiments, we tried different values of $\:{\uplambda\:}$ from 5х10^− 4 to 2.5х10^− 3 with a step of 5х10^− 4, and found that $\:{\uplambda\:}$ = 1х10^− 3 minimized the mean contour distances over the in-vivo trial dataset. The optimization at each resolution level was terminated when the relative value change of the cost function was less than a tolerance of 1х10^− 5.

Simulations

The simulated cine MRI dataset was generated by a digital phantom (XCAT) [22] and its MRI extension MRXCAT [29]. Specifically, XCAT was used to generate the cardiac-phase-resolved 3D tissue masks of 20 subjects, each with a different gender, cardiac anatomy, and cardiac motion pattern. The corresponding transformation field was also generated for each subject. Interpolation was then used to generate 2D cine movies in a 2-chamber slice and 7 or 8 short-axis slices, with a 1.5 mmх1.5 mm resolution and 24 frames. MRXCAT was then used to generate the signal intensity of each tissue in cine MRI. The studied feature tracking methods were performed in each slice to estimate the strains. The tracking quality was evaluated by end-point error (EPE), which is the difference between the estimated transformation field and its ground truth, averaged over the entire myocardium at end-systole or over the entire cardiac cycle. Strain error was computed at end-systole for both voxelwise (VSE) and global (GSE) levels, where the former was an average of absolute strain difference between the estimated strain and ground truth over the entire myocardium, and the latter was the absolute difference between the global estimated strain and global ground truth strain.

Experiments on the ACDC dataset

The ACDC dataset [23] contains 100 subjects evenly distributed over five study groups: dilated cardiomyopathy (DCM), myocardial infarction (MI), abnormal right ventricle (ARV), hypertrophic cardiomyopathy (HCM), and normal cardiac anatomy and function (NOR). For each subject, a series of short-axis cine images were acquired with a balanced steady-state free-precession (bSSFP) sequence on a 1.5T (Siemens Area, Siemens Medical Solutions, Germany) or 3T (Siemens Trio Tim, Siemens Medical Solutions, Germany) clinical scanner. Manual segmentation of the end-diastolic and end-systolic myocardium is provided by the dataset. To evaluate the tracking accuracy of a method, the end-diastolic epicardial and endocardial contours were first tracked to the end-systolic phase. Then, the mean distance between the tracked and annotated end-systolic contours was computed in the basal, mid-ventricular, and apical slices. To evaluate the inter-observer reproducibility of the measured strains, the end-diastolic endocardium and epicardium were independently contoured by a reader (3 years of cardiac MR experience) in a subset (n = 30) of the ACDC dataset, and the intraclass correlation coefficients (ICC) of the end-systolic global strains (radial and circumferential) based on the annotations by the dataset and the reader were calculated. To evaluate the intra-observer reproducibility, the same set of data was annotated by the same reader over one year later, and the ICC based on the two annotations by this reader were calculated.

Experiments on the patient dataset

The patient study was approved by the institutional review board. All patients provided informed written consent. The clinical dataset contains 16 patients (9 male, age 51 ± 14 years) scanned between February 2016 and February 2017 on a 3T scanner (Philips Ingenia Elition, Philips Healthcare, Nederland). Clinical indications included normal (n = 3), coronary heart diseases (n = 4), connective tissue diseases (n = 4), nonischemic heart diseases (n = 3), and other (n = 2). All patients were scanned with a bSSFP cine sequence and a tagging-MRI sequence. The bSSFP cine sequence used the following parameters: TR = 2.8 ms, TE = 1.4 ms, bandwidth = 1890 Hz/pixel, flip angle = 45°, FOV = 268х280 mm², matrix size = 134х140, slice thickness = 7 mm, number of phases = 30. The tagging-MRI sequence had the following parameters: TR = 5.9 ms, TE = 3.6 ms, bandwidth = 432 Hz/pixel, flip angle = 10°, tag spacing = 7 mm, FOV = 229х280 mm², matrix size = 154х188, slice thickness = 8 mm, number of phases = 9 to16. For each subject, only the matched cine and tagging images at basal, mid-ventricular, and apical short-axis slices were analyzed. The left ventricular myocardium in each end-diastolic cine or tagging image was annotated according to the AHA 16-segment model [30]. The tagging images were analyzed by pairwise registration following a previous work [31]. The end-systolic global and segmental strains estimated by the CMR-FT methods were compared with those by tagging using ICC.

Statistics

Continuous variables were reported as mean ± standard deviation. The Wilcoxon signed rank test was performed to determine whether two continuous variables were equivalent. The effect size of the Wilcoxon signed rank test was computed by dividing the absolute standardized test statistic z by the square root of the number of sample pairs. Following Cohen’s classification [32], the threshold of the effect size for a small, medium, and large effect are 0.1, 0.3, and 0.5, respectively. For evaluation of the inter/intra-observer reproducibility, ICC was calculated using a single rating, absolute agreement, and 2-way mixed-effects model; for evaluation of the cine-tagging correlation, ICC was calculated using a single rating, consistency, and 2-way mixed-effects model [33]. Bland-Altman analysis was also performed to evaluate the inter/intra-observer reproducibility of each method. All statistical analyses were performed with MATLAB and SPSS (Version 26, IBM, New York, USA). P values less than 0.05 were considered statistically significant.

Results

Simulations

Figure 2 shows changes of the EPE, VSE, and global strains over the entire cardiac cycle for one simulated subject. Groupwise-LLR achieved the lowest errors for tracking and strain estimation over the entire cardiac cycle. The largest improvement of Groupwise-LLR was found at end-systole (orange arrows). In late-diastole (blue arrows), the four groupwise methods were more accurate than the pairwise registration and Farneback optical flow, both of which suffered from the drift effect. Groupwise-LLR led to improved strain accuracies compared with the other methods, with a substantial improvement in the radial strain estimation.

Fig. 2

Plots of the end-point errors (EPE), voxelwise strain errors (VSE), and global strains (GS) over a cardiac cycle for one simulated subject. The curves were plotted as a function of normalized time. Groupwise-LLR achieved the lowest errors among the studied methods for tracking, voxelwise strain estimation, and global strain estimation. The largest improvement occurred at end-systole (orange arrows). In late-diastole, the four groupwise methods were more accurate than Farneback optical flow and pairwise registration (blue arrows) due to the elimination of the drift effect. LAX = long-axis; SAX = short-axis; LS = longitudinal strain; CS = circumferential strain; RS = radial strain

Table 1 compares the values of EPE_ES, EPE_all, VSE_ES, and GSE_ES from the studied methods over the entire simulated dataset. Groupwise-LLR achieved lower EPE_ES and EPE_all compared with the other methods for both long-axis and short-axis views (p = 0.13 against pTVreg for EPE_all in the long-axis view, and p < 0.05 for all other cases). Furthermore, Groupwise-LLR achieved lower VSE_ES and GSE_ES than the other methods for all strain types (p = 0.05 against pTVreg for the longitudinal GSE_ES, and p < 0.05 for all other cases). In particular, the improvement of radial strain accuracy by Groupwise-LLR was relatively substantial.

Table 1 Comparison of the end-point errors, voxelwise strain errors, and global strain errors between different methods on the simulated dataset

ACDC dataset.

Tracking quality

Figure 3 compares the tracking quality of all the studied methods, measured by the mean distance between each predicted contour and manually drawn contour at end-systole. For both epicardial and endocardial contours, Groupwise-LLR achieved significantly lower mean distances (1.72 ± 0.57 mm and 2.11 ± 1.33 mm) than Farneback optical flow (1.88 ± 0.72 mm and 2.50 ± 1.61 mm, both p<0.01), pairwise registration (1.80 ± 0.72 mm and 2.30 ± 1.47 mm, both p<0.01), nD + t B-Splines (1.87 ± 0.77 mm and 2.21 ± 1.51 mm, both p<0.01), pTVreg (1.80 ± 0.69 mm and 2.29 ± 1.56 mm, both p<0.01), and Groupwise-GLR (1.85 ± 0.73 mm, p<0.01 and 2.17 ± 1.48 mm, p<0.05). These results demonstrate the superiority of Groupwise-LLR in terms of tracking quality.

Fig. 3

Comparison of the tracking quality of different methods on the ACDC dataset. The height of each bar represents the mean value of the mean distance between the predicted contour by each method and the manually drawn contour at end-systole. The error bar indicates the standard deviation. For each subject, only the basal, mid-ventricular, and apical slices were used for the contour tracking. The p value and effect size (es) from the one-sided Wilcoxon sign rank test are labeled over each bar for the three alternative methods. Groupwise-LLR showed a significant reduction of mean distance compared with the other methods for both the epicardial and endocardial contours

Strain accuracy

Figure 4 shows the global strain curves of three subjects, each from a different study group. For every CMR-FT method, the MI and DCM subjects had lower strains than the control. However, whereas the two groupwise methods generated strains closer to zero at late-diastole, the two pairwise methods generated nonzero strains in late-diastole due to the drift effect (arrows). Figure 5 shows the mid-ventricular strain maps of an MI subject, who has wall thinning and akinesia in the septum. All four methods captured the reduction of radial strain in the septum; however, Farneback showed more inaccurate strains compared with the other methods. The late-diastolic displacement and strain estimates also appeared more accurate with the groupwise registration methods, due to the elimination of the drift effect. An additional movie file shows the estimated strain maps across the whole cardiac cycle [see Additional file 1].

Fig. 4

Global circumferential and radial strain curves of three subjects from different study groups in the ACDC dataset. The curves were plotted as a function of normalized time. The end-diastolic and end-systolic mid-ventricular images of each subject are shown at the top of each column. Normally, the strain at the end of the cardiac cycle should be nearly zero due to the periodicity of the cardiac motion. However, in the presence of the drift effect, the estimated late-diastolic strain can be significantly different from zero, which is observed for the Farneback optical flow and pairwise registration but not for the groupwise methods. ED = end-diastolic; ES = end-systolic; DCM = dilated cardiomyopathy; MI = myocardial infarction; NOR = normal cardiac anatomy and function; GCS = global circumferential strain; GRS = global radial strain

Fig. 5

Mid-ventricular cine images of a subject with myocardial infarction in the ACDC dataset and the corresponding end-systolic and late-diastolic displacement fields, circumferential strain maps, and radial strain maps. The septal wall (arrows) exhibited abnormal thinning and akinesia due to the presence of infarction. Farneback caused nonsmooth estimates of the circumferential strain and underestimates of the radial strain in the anterior wall. Both Farneback and pairwise registration suffered from the drift effect and manifested nonnegligible displacements and strains in the late-diastolic phase. The groupwise registration methods demonstrated a smoother circumferential strain map, more accurate radial strain estimates, and a reduction of drift effect. ED = end-diastolic; ES = end-systolic; LD = late-diastolic (corresponding to the last frame); CS = circumferential strain; RS = radial strain

Figure 6 shows the mean and 95% confidence interval of the global strain curves over each study group obtained by different methods. One can notice the systematic occurrence of drift effects in late-diastole with the two pairwise methods (arrows), which was largely reduced by the two groupwise methods. The mean absolute value of the late-diastolic GCS/GRS were 1.51%/4.11% for Farneback optical flow and 1.33%/4.09% for the pairwise registration. In comparison, the mean absolute value of the late-diastolic GCS/GRS were 0.74%/1.60% for Groupwise-GLR, and 1.00%/1.98% for Groupwise-LLR.

Fig. 6

The mean and confidence interval of the global circumferential and radial strain curves over each study group (n = 20) in the ACDC dataset. Each curve represents an average global strain, with different colors representing different study groups. The shade around each curve represents the 95% confidence interval at each time point. The arrows highlight the systematic drift effect of the Farneback optical flow and pairwise registration. GCS = global circumferential strain; GRS = global radial strain; DCM = dilated cardiomyopathy; MI = myocardial infarction; ARV = abnormal right ventricle; HCM = hypertrophic cardiomyopathy; NOR = normal cardiac anatomy and function

Table 2 shows the end-systolic global strains obtained by different methods for each study group. Groupwise-LLR generated a significantly greater GRS for every study group than the other methods, while the GCS was more similar. These trends were consistent with the findings from the simulated dataset.

Table 2 Comparison of the End-Systolic Global Circumferential and radial strains between different methods on the ACDC dataset

Reproducibility

Figure 7 shows the inter-observer reproducibility of each CMR-FT method using ICC [33]. For all methods, the ICC of GCS was “moderate” to “excellent” while the ICC of GRS was “excellent”. Figure 8 shows the results of the intra-observer reproducibility evaluation. For all methods, the ICC of both GCS and GRS were “excellent”. These results suggest that the four methods have similar inter/intra-observer reproducibility.

Fig. 7

Bland-Altman analysis for the inter-observer reproducibility evaluation of the end-systolic global strains based on the ACDC dataset. The intraclass correlation coefficient (ICC) and its 95% confidence interval were labeled at the top of each plot. The four methods had similar inter-observer reproducibility. The disagreement of global circumferential strain (GCS) was greater than that of global radial strain (GRS) for every method. ES = end-systolic; DCM = dilated cardiomyopathy; MI = myocardial infarction; ARV = abnormal right ventricle; HCM = hypertrophic cardiomyopathy; NOR = normal cardiac anatomy and function

Fig. 8

Bland-Altman analysis for the intra-observer reproducibility evaluation of the end-systolic global strains based on the ACDC dataset. The intraclass correlation coefficient (ICC) and its 95% confidence interval were labeled at the top of each plot. The four methods manifested similarly high intra-observer reproducibility. GCS = global circumferential strain; GRS = global radial strain; ES = end-systolic; DCM = dilated cardiomyopathy; MI = myocardial infarction; ARV = abnormal right ventricle; HCM = hypertrophic cardiomyopathy; NOR = normal cardiac anatomy and function

Patient dataset

Table 3 compares the ICC of different methods with tagging for end-systolic global strains. For end-systolic GCS, pairwise registration, Groupwise-GLR, and Groupwise-LLR all showed a stronger correlation with tagging than Farneback optical flow. For end-systolic GRS, both groupwise methods showed a stronger correlation with tagging than the optical flow and pairwise registration in the mid-ventricular slice (r = 0.80 and 0.78 vs. r = 0.55 and 0.72), apical slice (r = 0.50 and 0.57 vs. r = 0.46 and 0.39), and mean over the three slices (r = 0.65 and 0.69 vs. r = 0.49 and 0.56). All feature tracking-based global strains had a significant correlation with those based on tagging-MRI, except the apical radial strain measured by the pairwise registration (p = 0.06).

Table 3 Comparison of the correlation with tagging between different methods for end-systolic global strains on the patient dataset

Figure 9 shows the bull’s eye plots of the per-segment ICC of each method relative to tagging-MRI. For the circumferential strain, pairwise registration, Groupwise-GLR, and Groupwise-LLR achieved similar performances and outperformed Farneback optical flow. For the radial strain, Groupwise-LLR yielded higher correlation coefficients than Farneback optical flow, pairwise registration, and Groupwise-GLR in 13, 10, and 13 of 16 segments, respectively, suggesting an improved accuracy in regional radial strain estimates.

Fig. 9

Bull’s eye plots of per-segment intraclass correlation coefficients with tagging for end-systolic strains on the patient dataset. Groupwise-LLR demonstrated stronger correlations with tagging compared with the other methods, especially for the radial strain. *: Significantly correlated with the estimates by tagging-MRI. ES = end-systolic; CS = circumferential strain; RS = radial strain

Discussion

In this work, we propose a novel CMR-FT method based on deformable groupwise registration with an LLR dissimilarity metric. We evaluated the myocardial tracking quality, strain accuracy, and strain inter/intra-observer reproducibility of the proposed method via simulated and in vivo data. The simulations show that the proposed Groupwise-LLR method was more accurate in both tracking and strain estimation than the other methods. Results on the ACDC dataset show that the proposed method improved tracking quality and suppressed the drift effect common for optical flow and pairwise registration. However, the studied CMR-FT methods were similarly reproducible. The patient data show that Groupwise-LLR achieved a stronger correlation with tagging for the radial strains in the median slice, apical slice, and most myocardial segments.

Feature tracking based on cardiac cine MRI has drawn increasing attention due to its reproducibility, ease of use, and excellent value for risk stratification and prognosis of cardiac diseases [4, 34, 35]. Optical flow and pairwise registration are currently the major algorithms used in existing commercial software [8, 36]. However, since both of them track the myocardium frame-by-frame, small errors within each tracking step can accumulate to a larger error, resulting in a less accurate evaluation of myocardial contractility in the later cardiac phases [9]. On the other hand, although groupwise registration has been postulated to provide more accurate myocardial tracking [16, 17], little data exist to evaluate the accuracy of strains measured by the groupwise registration methods. Our findings filled this niche by showing that groupwise registration methods indeed generate more accurate strain estimates than those pairwise methods, especially for the radial strain which has been proven to be highly related to amyloidosis [37], myocarditis [38], aortic stenosis [39], etc. Furthermore, by minimizing the drift effect, the proposed method provides a more accurate estimation of the strains in the diastolic phase. The diastolic function of the heart can be affected by various diseases, such as heart failure and cardiomyopathies [40, 41]. Therefore, the improved accuracy in the estimation of diastolic strains and strain rates may translate to more accurate diagnosis and prognosis of the relevant diseases. In addition to the evaluation of the left ventricular function, the proposed method may also be used to analyze the function of other chambers, such as the right ventricle [42] and left atrium [43]. For these chambers, the assessment of diastolic function is also very important. For example, the strains for the left atrium are usually divided into reservoir strain, conduit strain, and booster strain, where the latter two occur in the diastolic phase [44]. A more accurate estimation of these diastolic strains may provide more accurate evidence for diagnosis and prognosis of various diseases.”

Recently, due to the rapid advances of machine learning, many CMR-FT methods based on deep learning have been proposed [17, 45,46,47,48,49,50]. However, these learning-based feature tracking methods were almost all based on pairwise registration, and thus still suffer from the same problem as their optimization-based counterparts. On the other hand, although some learning-based groupwise registration methods have also been proposed recently for cardiac cine MRI [51,52,53], none of them have been applied to strain estimation. The improvement of our method over traditional pairwise feature tracking indicates that the application of learning-based groupwise registration methods could potentially yield a more accurate assessment of the strains. Additionally, learning-based strain estimation methods usually require a relatively diverse training dataset to ensure generalizability for different diseases, scanners, and image qualities. However, for some anatomies, such as the right atrium and pulmonary artery, such a diverse dataset has not been publicly available. As our method is completely unsupervised, its application to strain estimation for other chambers and vessels can be more easily available.

The study and technique have limitations. Firstly, we did not evaluate the longitudinal strain over the in vivo datasets. However, the results on the simulated dataset preliminarily validated the potential of the proposed method for longitudinal strain estimation. Secondly, we only validated our method using a relatively small patient cohort. Further evaluation based on a larger and more diverse patient cohort is highly warranted. Thirdly, we only applied our method to 2D cine images and thus could not avoid the errors caused by the out-of-plane motion. Application of the method while exploiting the correlation between adjacent image planes may partially account for the influence of the through-plane motion [54]. Furthermore, the development of 3D cine imaging methods [55, 56] and strain estimation methods may provide an ultimate solution to this issue. Fourthly, since we did not assess the proposed method for different field strengths, slice thicknesses, flip angles, resolutions, sequences (e.g., GRE), and artifacts (e.g., severe banding and flow artifacts [57]), it is not yet known how robust the method is to images of different qualities. Although the proposed groupwise registration method may be relatively robust to low-quality images due to the exploitation of the context information of the whole image sequence, this awaits systematic investigation in large-cohort clinical studies. Finally, the computational time of the proposed method is relatively long. The typical computational time of the four CMR-FT methods for processing a single cine movie was 0.8s (Farneback optical flow), 33s (pairwise registration), 24s (Groupwise-GLR), and 38s (Groupwise-LLR). However, the computational time of the proposed method is still within a reasonable range. Further accelerations can be obtained by adopting patchwise parallel computation for the LLR dissimilarity and a more efficient optimization algorithm. Moreover, learning the complex groupwise optimization process with an end-to-end neural network may also reduce the computational time to a large extent.

Conclusions

In conclusion, we developed a CMR-FT method based on deformable groupwise registration with an LLR dissimilarity metric. Our results based on simulated and in vivo datasets suggest that the method improves the accuracy of myocardial tracking and strain estimation compared with previous methods. Application of the proposed method may thus facilitate a more accurate assessment of the myocardial function in the clinical setting, especially along the radial direction and in the diastolic phase.