High-Resolution Maps of Left Atrial Displacements and Strains Estimated with 3D CINE MRI and Unsupervised Neural Networks

Abstract

The functional analysis of the left atrium (LA) is important for evaluating cardiac health and understanding diseases like atrial fibrillation. Cine MRI is ideally placed for the detailed 3D characterisation of LA motion and deformation, but it is lacking appropriate acquisition and analysis tools.

In this paper, we present Analysis for Left Atrial Displacements and Deformations using unsupervIsed neural Networks, Aladdin, to automatically and reliably characterise regional LA deformations from high-resolution 3D Cine MRI. The tool includes: an online few-shot segmentation network (Aladdin-S), an online unsupervised image registration network (Aladdin-R), and a strain calculations pipeline tailored to the LA. We create maps of LA Displacement Vector Field (DVF) magnitude and LA principal strain values from images of 10 healthy volunteers and 8 patients with cardiovascular disease (CVD). We additionally create an atlas of these biomarkers using the data from the healthy volunteers.

Aladdin is able to accurately track the LA wall across the cardiac cycle and characterize its motion and deformation. The overall DVF magnitude and principal strain values are significantly higher in the healthy group vs CVD patients: $\mathbf{2.85\pm 1.59~{}mm}$ and $\mathbf{0.09\pm 0.05}$ vs $\mathbf{1.96\pm 0.74~{}mm}$ and $\mathbf{0.03\pm 0.04}$, respectively. The time course of these metrics is also different in the two groups, with a more marked active contraction phase observed in the healthy cohort. Finally, utilizing the LA atlas allows us to identify regional deviations from the population distribution that may indicate focal tissue abnormalities.

The proposed tool for the quantification of novel regional LA deformation biomarkers should have important clinical applications. The source code, anonymized images, generated maps and atlas are publicly available: https://github.com/cgalaz01/aladdin_cmr_la.

Atrial Cine MRI, Left Atrial Function, Segmentation Neural Network, Image Registration Neural Network, Online Few-Shot Online, Left Atrial Displacements, Left Atrial Strains, Atlas of the Left Atrium, Atrial Fibrillation, Left Atrial Mechanics.

1 Introduction

The motion and deformation of the left atrium (LA) are impaired in conditions such as atrial fibrillation (AF) and heart failure (HF) [1, 2]. So far, these assessments typically rely on global or one-dimensional changes in atrial wall length measured from single-plane cardiac views, quantified as LA strain. This is in contrast to the left ventricle (LV), for which several dedicated 3D regional metrics of motion and deformation exist [3], with known clinical correlates [4, 5, 6, 7].

1.1 LA Strains

The LA function comprises three distinct functional phases: 1) reservoir phase, in which the LA is filled with oxygenated blood from the pulmonary veins, reaching its maximum volume; 2) conduit phase, where blood is passively emptied into the LV, reducing the LA volume; and 3) booster pump, in which the LA actively contracts, further emptying the chamber to its minimum volume [8].

To assess LA motion and deformation, 2D images in the cardiac long-axis view are typically acquired across the cardiac cycle using either echocardiography or dynamic (Cine) Magnetic Resonance Images (MRI). Using these, the LA function across its three phases is usually assessed based on the ”global” or ”regional LA strains” [9, 10]. The global LA strain is a single value that describes the change in LA perimeter in a 4-chamber view [9]. It is the simplest and mostly reported strain value [9]. In 2D regional strain analysis, the change in lengthen of specific segments of the LA contour is measured [9]. These LA strains, and their temporal derivatives, ”LA strain rates”, have been correlated with AF burden, HF burden, presence of LA fibrosis, heart valve disease and cardiomyopathies [1, 11, 12].

While 2D strains are valuable, they provide an incomplete characterization of tissue deformations compared to the definition of strain in the physical sciences. In engineering, strains are tensors, which quantify the relative material deformation (stretch and shear) of a body in response to applied stresses. 3D full coverage strain measurements of the LA have the potential to capture vital information concerning tissue biomechanics (e.g., the presence of stiffer, fibrotic LA regions), thereby expanding clinical applications [2]. As in the LV, it may be used to identify atrial dyssynchrony, even before changes in functional biomarkers, such as ejection fraction, become abnormal [13]. Thus, high resolution spatio-temporal 3D maps of LA deformation can potentially provide more specific signatures of LA pathology, with greater diagnostic and prognostic value.

1.2 LA Strain Imaging

1.2.1 Speckle Tracking Echocardiography

LA deformations are most often assessed using 2D Speckle Tracking Echocardiography (STE) [9] in the clinic. STE analyzes cardiac motion by monitoring the movement of distinctive bright or dark speckles in grayscale echocardiographic images throughout the cardiac cycle. As the heart contracts and relaxes, the STE software tracks the motion of the speckles in the 2D image plane to quantify tissue displacements and deformations. 3D STE [14], which includes a detailed representation of the whole LA, enables better motion quantification by measuring deformations along all three direction. It has been evaluated for LV strain analysis [3], but its clinical value for the LA is still under investigation [15, 12].

1.2.2 Cine MRI

An alternative solution is to use Cine MRI to characterise the LA motion, as it offers better signal to noise ratio and higher spatial resolution than STE, and is operator independent [16, 17]. Cine MRI allows the reconstruction of MR images at various cardiac phases, enabling a dynamic evaluation of cardiac motion. Cine MRI is included in most clinical MRI examinations. The LA is visible, although incompletely, in the 2-chamber, 4-chamber and 3-chamber views, each acquired in a separate breath-hold as thick slices. These images have been used to study LA motion and deformation [10, 18, 19, 20], usually analyzed using Feature Tracking (FT) techniques [10]. FT is a post-processing step [16] which involves the identification and tracking across the cardiac cycle of key features (e.g., manually or automatically delineated LA wall regions).

These deformation analyses are inherently incomplete, as they are conducted as thick single slice images that lack full LA coverage. Moreover, at the typical spatial resolution used for ventricular 2D Cine MR imaging (slice thickness: $8-10~{}mm$, in-plane resolution: $\sim 2~{}mm$ [21]), partial volume effects preclude an accurate identification and motion characterisation of the LA wall. Different slices are also typically acquired at different breath-holds, further complicating 3D motion analysis. Nonetheless, 2D FT has been found clinically useful for the LA [22]. 3D FT has not yet been explored due to the lack of 3D CINE MRI protocols for the LA. Previous studies have nevertheless analysed LA function using stacks of thick slices covering the LA acquired in a short-axis view [23, 24].

1.3 Displacement Vector Field

Both FT and STE track the position of regions of interest (ROIs) of the image across the cardiac cycle, using optical flow techniques such as block matching or image registration (alignment) [16]. All methods ultimately estimate Displacement Vector Fields (DVFs), the vectors that link material points in the myocardium across different phases of the cardiac cycle. Strain tensors can be estimated through mathematical manipulation of the DVFs.

Neural network (NN)-based image registration approaches offer significant potential for medical imaging applications [25, 26]. Recently, various NN-based image registration methods have been proposed to analyze deformation of the LV across the cardiac cycle from Cine MR images. Common architectures are based on: U-Net [27, 28, 29, 30, 31], variational autoencoders [32, 33, 34], Siamese autoencoders [35, 36], or pyramid networks [37]. Usually, the displacement field is learned in an unsupervised way by using a spatial resampling module [38].

Segmentation maps (partitions of an image into anatomically meaningful regions) of cardiac structures have been used with cardiac motion tracking networks to improve their performance. Typically, this can be accomplished by incorporating the segmentation outputs as part of the input [32], constraining the DVFs so the transformed segmentations agree with the target segmentation [39, 27], or integrating them into a joint segmentation-motion tracking task [36, 29]. Additionally, there are efforts to extract information from multiple cardiac MRI views (3D and 2D) to improve the motion estimation [40, 41]. More recently, there is an interest in using biomechanics-informed modeling, which learn the underlying deformation properties of the heart, to regularize the DVFs towards more realistic motion [30, 34, 42, 43, 44].

Research in estimating the 3D DVFs of the LV from MRI data has been ongoing for many years, and has resulted in mature and specialized methodologies. However, the LA noticeably differs from the LV, bringing specific unresolved challenges. Compared to the LV, the LA has much thinner myocardial walls [45, 46, 47], a more irregular and complex morphology [45, 48], and structural diversity between individuals [45]. Out-of-the-box techniques designed for the LV are therefore unlikely to succeed in the LA.

1.4 Atlas

The clinical value of LA regional motion characterisation hinges on the ability to identify regional LA motion abnormalities. For this, LA maps of reference values of regional LA deformation biomarkers (such as DVF magnitude or principal strain values) are needed. These LA maps can be constructed by building an LA atlas of motion biomarkers. In this context, an atlas is a statistical model that characterizes the features of a specific population [49].

An atlas is typically created using registration techniques to align multiple representations (e.g., segmentation maps) of a given organ across a specific population [50, 29, 51]. Regional biomarkers (such as LA displacements and strains) can be mapped onto the atlas common reference space, to enable statistical and visual analyses across subjects [50, 51]. Atlases serve as standardized references for analyzing and comparing imaging data, and can provide a method for comprehensive characterization of normal LA motion and deformation patterns [49, 52, 53, 51]. Furthermore, they allow for the identification of regional differences in biomarkers between pathological cases and the healthy population [54].

1.5 Overview and Aims

We propose Aladdin, Analysis of Left Atrial Displacements and Deformations with unsupervIsed neural Networks, to automatically characterize LA displacements and strains from high resolution 3D Cine MRI. Using these, we characterise LA displacements and strains in healthy and CVD patient cohorts. We also create the first atlas of regional LA displacements and deformations across the cardiac cycle. We present:

1. A segmentation NN, Aladdin-S, for the automatic segmentation of the LA across the cardiac cycle.

2. An image registration NN, Aladdin-R, to automatically estimate LA DVFs across the cardiac cycle.

3. An algorithm to calculate regional LA strains across the cardiac cycle, which treats the LA myocardium as an infinitesimally thin 2D surface.

4. The construction of a DVF and strain atlas using data from 10 healthy volunteers.

5. Proof-of-principle characterization of regional LA strains of cardiovascular patients using the created atlas.

Overview of the workflow can be seen in Figure LABEL:fig:overview.

2 Methods

The source code for Aladdin, the LA biomarkers maps and the atlas, as well as the anonymized Cine MR images are publicly available: https://github.com/cgalaz01/aladdin_cmr_la.

2.1 LA Cine MR Images

We use 3D Cine MRI balanced Steady State Free Precession (bSSFP) scans of the LA acquired using a novel high-resolution acquisition protocol [55]. The images are electrocardiogram (ECG)-gated and acquired in a single 25-second breath-hold using a 1.5T Philips Ingenia MRI scanner with a 32-channel cardiac coil. The typical field of view used is $400\times 270\times 70~{}mm^{3}$ and the acquisition matrix is $256\times 256\times 36$ to give a reconstructed resolution of $1.72\times 1.72\times 2.00~{}mm^{3}$. We use SENSE with factors of 2.3 and 2.6 along each phase encode direction, and 55% view sharing for a total of 20 phases across the cardiac cycle. Cardiac phase 0 corresponds to the ECG R-peak, representing ventricular end diastole. This corresponds to atrial end systole at which the LA is expected to have the smallest volume. A representative MR image can be seen in Figure 1.

We analyse images acquired, after informed consent and under ethical approval, from 18 subjects: 10 healthy volunteers (age: $30.4\pm 4.7$ years; female: $40.0\%$), and 8 patients with cardiovascular disease (CVD) (age: $53.8\pm 15.7$ years; female $62.5\%$). The pathologies present in the CVD group are: myocarditis (1 patient), history of syncope (2), myocardial infarction (3), non-ischaemic cardiomyopathy (1), and hypertrophic cardiomyopathy (1).

High Resolution 3D Cine MRI

Figure 1: Top: Example cross section of the LA from a healthy volunteer imaged using the dedicated 3D high resolution acquisition protocol. Bottom: LA volume across the cardiac cycle from the same healthy volunteer. Phases 0, 8 and 15, in which the LA was manually segmented (see Section 2.2), correspond approximately to the transitions between the reservoir, conduit and booster pump phases.

2.2 Preprocessing

The LA images are preprocessed to standardize their characteristics before being inputted to the NNs: they are cropped to a size of $96\times 96\times 36$ voxels, centered at the LA, and their intensity is min-max normalized between $[0,1]$. We also perform a rigid-body translation of the LA to ensure its centroid is stationary across the cardiac cycle.

Our tool utilizes LA segmentation maps from cardiac phases 0, 8 and 15, which are manually segmented in all subjects. Moreover, 3 healthy volunteers and 3 CVD subjects have their LA manually segmented across the whole cardiac cycle by an expert. These additional segmentation maps across the cardiac cycle are used as ground truth data for evaluating the segmentation and image registration models.

Figure 2: Overview of the proposed common architecture of Aladdin-S (top) and Aladdin-R (bottom) and their respective data inputs (left). In total, we have 3 downsampling (decoding) and upsampling (encoding) steps. Each encoding block consists of a $3\times 3\times 3$ convolutional layer followed by a $2\times 2\times 2$ maxpooling layer. The decoding block consists of a concatenation layer of the output of the previous layer and the respective skip connection, followed by a $3\times 3\times 3$ transposed convolutional layer and a $3\times 3\times 3$ convolutional layer. After each convolutional layer, a Leaky Rectified Linear Unit (ReLU) activation [56] with an alpha value of $0.3$ is applied to all except the final convolution, which utilizes a sigmoid activation. Glorot uniform initialization [57] is used to set the initial weights. Both networks are trained independently on a subject-by-subject basis.

2.3 Aladdin-S: Segmentation Neural Network

To avoid manually segmenting the LA across all phases of the cardiac cycle, we propose a few-shot segmentation network, Aladdin-S, implemented in TensorFlow v2.10 [58]. Aladdin-S takes as input an LA image in a given cardiac phase and predicts a corresponding LA binary segmentation map. As we have limited available data, we perform online learning [35] (subject-by-subject basis training) to independently train the network to each subject. We use the manual segmentation maps at the three cardiac phases (0, 8 and 15) as ground truth data for training. We apply random augmentations to the input to approximate and generalize to the remaining phases. These include: anisotropic scaling with factors chosen by uniformly sampling the interval [0.75, 1.15]; rotations around the cardiac long axis with angle between $[-30,~{}30]^{\circ}$; Rician noise addition with a mean of 0 and a standard deviation chosen from [-0.05, 0.05]; and a bias field correction with a coefficient range of [0, 0.1] [59]. Each augmentation has a 95% probability of being applied. We use trilinear and nearest neighbor interpolation for the images and segmentation maps, respectively, when resampling.

Aladdin-S is a 3D convolutional NN that is based on the U-Net architecture [60]. As with the U-Net, Aladdin-S uses encoding and decoding blocks to capture features at different spatial resolutions, and skip connections to preserve higher-resolution features. Rather than the standard skip connection between each encoder and decoder block we use squeeze and excitation (SE) blocks [61], which offer dynamic feature map-wise re-calibration, with a squeeze ratio of 8. The complete architecture can be seen in Figure 2.

We update the weights of the model by minimizing the Dice loss of the target segmentation using the Adam optimizer algorithm, with parameters specified from [62], for $500$ epochs, which are empirically sufficient for model convergence. Experimentally, we found that an initial learning rate of $0.0005$, a learning decay rate of $0.8$, and a batch size of $1$ lead to the fastest convergence.

During inference, we perform test-time augmentations [63, 64] of 8 in-plane rotations to improve LA representation. To remove noise, we use connected-component labeling, which groups neighboring voxels together, and select the largest segment. We further smooth the segmentation maps by applying a median filter of size $3\times 3\times 3$ voxels. Then, we re-align the rotated segmentations and carry out a voxel consensus voting approach, considering a 75% voting agreement as the threshold for determining the accepted final segmentation.

We evaluate the performance of Aladdin-S using the Hausdorff Distance (HD, in $mm$) and Dice score (DS). With HD we identify the maximum distance between pairs of nearest points on the GT segmentation against the predicted segmentation, with a value closer to 0 indicating better performance. With DS we assess the level of overlap between the GT segmentation and predicted segmentation, with a Dice score of 1 indicating perfect overlap.

2.4 Aladdin-R: Image Registration Neural Network

As we are only interested in the motion of the LA wall, we use masked MR images of the LA myocardium (see Figure 2 and Figure Supplementary S1.3) as inputs to the registration neural network. This approach removes any irrelevant tissue to reduce tracking errors while leaving enough context for the network to distinguish between the LA wall and surrounding tissues. We first extract the contours of the LA segmentation maps yielded by Aladdin-S across the cardiac cycle. The LA contours are then dilated using a spherical structure with a radius of 2 voxels. Using this dilated LA contour segmentation, the original LA images are masked out (see Figure 2 and Figure S1.3) and used as inputs to the image registration network, Aladdin-R - see Figure 2.

Aladdin-R takes in a pair of contour-masked 3D LA images as predicted from Aladdin-S from the same patient: one at phase 0 (moving image), and the other at another cardiac phase (target image). The outputs of the model are the deformed moving image and the associated DVFs that encode this non-linear transformation. For simplicity, Aladdin-R uses a similar architecture to Aladdin-S as described in section 2.3. For Aladdin-R, the 3D input pairs are combined into a unified 3D image, each represented as a different channel. The Aladdin-R final convolutional layer uses a linear activation function and outputs each vector component of the DVF as a distinct feature map. After this final convolutional layer, we use a spatial resampling module [38] that linearly resamples the moving image to the target image given the predicted DVF.

As with Aladdin-S, we employ online learning to train Aladdin-R on all the cardiac phases from a single subject for $500$ epochs. The weights are updated by minimizing the mutual information loss [65] (the negative of the mutual information) between the target and the resampled moving image with an empirically-set bin size of 128. The learning rate is set to $0.001$, a learning decay rate of $0.8$, and a batch size of $1$.

The image registration task is ill-posed because the image similarity loss does not guarantee a unique optimal DVF solution. Moreover, the expected deformations of the LA across the cardiac cycle should be spatio-temporally smooth to preserve its physical integrity. For these reasons, we add a bending energy regularization term to the loss function [66], with an empirically selected regularization weight of $0.1$. This regularizer will guide Aladdin-R towards a solution that has the desired mathematical properties and is physiologically plausible.

We assess Aladdin-R’s ability to track the non-dilated LA contour across the cardiac cycle using once again the HD and DS metrics. These are evaluated between the ground truth LA contour (created from a manual segmentation - see Section 2.2) and LA contour for the same cardiac phase outputted by Aladdin-R. These analyses are performed in the 6 subjects with ground truth segmentation maps across all phases of the cardiac cycle.

We also compared the performance of Aladdin-R with state-of-the-art image registration solutions proposed for similar tasks: cardiac image registration packages Medical Image Tracking Toolbox (MITT) [67] and Advanced Normalization Tools (ANTs) [68]; and image registration NNs Voxelmorph ($Vxm$) [69] and Voxelmorph constrained on the image segmentation ($Vxm\text{-}seg$) [31]. Further details of these comparisons are provided in the Supplement.

2.5 LA Strain Calculation

Due to the LA’s thinness relative to its other dimensions, we treat it as an infinitesimally thin surface [70]. Strains on each LA segment can therefore be represented by a 2D symmetric tensor with two principal strain directions (embedded on the surface LA mesh) and two corresponding principal strain values.

We first obtain a closed surface triangular mesh of the LA by applying the marching cubes algorithm to the LA segmentation from the first cardiac phase using scikit-image [71]. We then use the PyVista library for the mesh post-processing operations [72] detailed in the Supplement.

For each triangular cell in the mesh, we calculate the finite Green-Lagrangian strain tensors [73] using cardiac phase 0 as the reference (undeformed) state. We take the following steps:

1. We calculate the metric tensor, $g$, of the vertex points in curvilinear coordinates. For each LA mesh triangle with vertex positions $x_{0}$, $x_{1}$ and $x_{2}$, defined in a Cartesian basis, the normalised covariant base $v=[g_{1},g_{2},g_{3}]$ can be calculated using:

   $$\begin{split}&g_{1}=\frac{x_{1}-x_{0}}{||x_{1}-x_{0}||}\\ &g_{2}=\frac{x_{2}-x_{0}}{||x_{2}-x_{0}||}\\ &g_{3}=g_{1}\times g_{2}\\ &g=v^{T}\cdot v\end{split}$$ (1)

2. The metric tensors at the reference configuration (cardiac phase 0), $G_{ij}$, and in each of the deformed configurations, $g_{ij}$, can be used to estimate the Green-Lagrange strain tensor $E_{ij}$ [74, Chapter 5]:

   $$E_{ij}=\frac{1}{2}(g_{ij}-G_{ij})$$ (2)

3. $E_{ij}$ are represented in local covariant coordinates, which are difficult to interpret. We now rewrite them in Cartesian coordinates for the reference configuration, $E_{kl}$, represented by the orthonormal basis $(e_{1},e_{2},e_{3})$:

   $$\begin{split}&e_{1}=g_{1}\\ &e_{2}=\frac{g_{2}-g_{2}\cdot e_{1}}{||g_{2}-g_{2}\cdot e_{1}||}\\ &e_{3}=\frac{e_{1}\times e_{2}}{||e_{1}\times e_{2}||}\\ \end{split}$$ (3)

4. The Green-Lagrangian strain tensor is calculated as:

   $$E_{kl}=E_{ij}(e_{k}\cdot G^{i})(G^{j}\cdot e_{l})$$ (4)

5. We compute the eigenvalues ($\lambda$) and eigenvectors ($v$) of $E$ defined in step 4. As the LA is represented by a 2D surface, the strain eigenvalue corresponding to the direction orthogonal to the local LA surface is always 0. The principal strain value and direction of $E$ at each triangular cell in the LA are treated as the regional main biomarkers of LA deformation.

2.6 LA Atlas Generation

Atlas Generation

Figure 3: Construction steps of the average atlas from the LA segmentation maps of the 10 healthy volunteers. It consists of a forward step (green full line) that applies the transformation $T_{i}$ that registers all cases to the initial reference space ($LA_{0}$, chosen randomly from the healthy volunteers pool). The forward transformations are then inverted and averaged, to give the average inverse transformed $\bar{T}^{-1}$. The atlas is created by applying $T_{i}\bar{T}^{-1}$ to each case to bring it to the final atlas space.

We create an atlas of the LA with the segmentation maps from the 10 healthy volunteers, to allow the direct comparison and statistical analyses of the calculated DVF and strain metrics. To construct the atlas, we follow a similar strategy to the brain atlas creation steps detailed in [50] using the SimpleITK [75] package. An overview of the atlas generation process can be seen in Figure 3.

The first step consists of registering the LA segmentation in cardiac phase 0 of each case ($i$) to an arbitrary selected reference case, in a forward transformation ($T_{i}$) step. This consists of successive affine and non-rigid diffeomorphic registrations, as detailed in the Supplement 5.2.1.

Then, to avoid biasing the atlas to the randomly selected reference case, we then calculate the average inverse transformation ($\overline{T}^{-1}$): $\overline{T}^{-1}=\frac{1}{n}\sum_{i=1}^{n}{T_{i}^{-1}}$. This has proven to be a robust strategy in minimizing potential biases of the atlas towards the initial reference case [50, 76, 77].

By applying $T\overline{T}^{-1}$ to each $i$, we transform each case into the atlas space. We overlay these transformed cases to obtain the atlas consensus segmentation map which requires at least 50% agreement between all cases. Then, we obtain the mesh of the atlas segmentation map using marching cubes.

The final step consists of the registration of the DVF and strain fields estimated in each subject to the atlas space. The rotation component of $T_{i}\overline{T}^{-1}$ needs to be applied to the vector fields (DVF and principal strain directions) to ensure that they are correctly oriented in the atlas space.

We construct the final atlas by obtaining the mean, standard deviation and coefficients of variation (CV) of the motion and deformation metrics: magnitude of DVFs, principal strain directions and principal strain values. These are calculated at each LA atlas vertex, over all healthy subjects and across all cardiac phases.

To evaluate how much the LA metrics in a given subject differ from the atlas, we register that case directly to the atlas using the affine and non-rigid methods employed for atlas creation (see Supplement 5.2.1). Then, to quantify individual deviations from the atlas we use the Mahalanobis Distance ($MD$) [78] for each metric:

$$MD_{vt}(\mathbf{s_{vt}})=\sqrt{(\mathbf{s_{vt}}-\boldsymbol{\mu_{vt}})^{T}\Sigma^{-1}_{vt}(\mathbf{s_{vt}}-\boldsymbol{\mu_{vt}})}$$ (5)

where $\mathbf{s}$ is the subject’s metric, $\boldsymbol{\mu}$ is the atlas mean of the metric, $\boldsymbol{\Sigma}$ the atlas covariance. $\mathbf{v}$ indexes the vertex, $\mathbf{t}$ the cardiac phase.

As a demonstration of the clinical application of Aladdin, we calculate the MD for an imaged CVD patient with a diagnosis of myocarditis (see Video S5).

3 Results

Our proposed method successfully estimated high-resolution LA displacement and strain maps across the cardiac cycle from all 18 subjects. We were additionally able to create an atlas of LA motion and deformation biomarkers and identify significant differences between healthy volunteers and CVD patients.

3.1 LA Segmentation

Aladdin-S Results

Figure 4: Results of Aladdin-S across the cardiac cycle. Left: Short-axis view showing a representative predicted segmentation map for phase 3 (ground truth: gray; prediction: orange). Right: Performance of Aladdin-S across the cardiac cycle based on HD (Plot A) and DS (Plot B), computed in 6 subjects.

Our proposed few-shot Aladdin-S model demonstrates high accuracy in segmenting the LA throughout the cardiac cycle, despite learning from augmented segmentation maps from 3 cardiac phases only. A representative example of a segmentation (from phase 3) can be seen in Figure 4. In the 3 healthy cases tested, we achieve a HD of $2.62\pm 1.12~{}mm$ and DS of $0.95\pm 0.01$. We observe similar results on the 3 CVD cases with a HD of $2.46\pm 0.71~{}mm$ and a DS of $0.94\pm 0.02$. Although Aladdin-S performs better on cardiac phases that are closer to the manual segmentations used for training, the data augmentation scheme used successfully prevents overfitting.

3.2 LA Registration

Aladdin-R is able to accurately track the LA wall across the cardiac cycle (see Video S1), particularly in cardiac phases with large displacements. The LA exhibits the expected rapid volume increase during the reservoir phase (0-8), followed by passive and active emptying during the conduit (9-15) and booster pump phases respectively (16-19), as shown in the last panel of Video S1.

To assess the accuracy of the registration algorithm, we compare the segmentation of the tracked LA with the ground truth LA segmentations for the same cardiac phase. We obtain an average HD of $2.87\pm 0.95~{}mm$ and a DS of $0.87\pm 0.05$, outperforming all the literature methods (see Table 1 and Video S1). Aladdin-R is able to more successfully track the LA wall during high deformation phases than literature methods MITT, ANTS and Voxelmorph (Figure 5, Video S1 and Table S1).

Moreover, Aladdin-R generates temporally and spatially smooth DVFs, surpassing the DVFs obtained using other methods, whose DVFs show a more erratic and zigzagging trajectory across the cardiac cycle (see Figure S2).

The two participant groups (healthy and CVD patients) display significantly different motion characteristics, as observed from the average DVF magnitude and principal strain value across the cardiac cycle (see Figure 6). The observed differences are significant during the first half of the reservoir phase and booster-pump phase (t-test $p<0.001$). The CVD cases have reduced DVFs and principal strain values overall, as well as less marked inflection points between the conduit and pumping phases (phases $\sim$11-13). Both groups have spatio-temporally smooth point displacements and larger DVF magnitude at the anterior wall compared to other LA regions (see Figure Video S2 and Video S3 first row). This spatial heterogeneity in DVF magnitude is more marked in healthy volunteers than CVD cases.

Image Registration Results

Figure 5: Performance of the compared registration methods in tracking the LA across the cardiac cycle on the Dice score (top) and Hausdorff distance (bottom). The proposed Aladdin-R model obtains the overall best performance across the whole cardiac cycle.

Healthy and CVD Comparison

Figure 6: Comparison between healthy (green) and CVD (red) groups across the cardiac cycle and their respective t-test $p$-value (black). Plot A shows the DVF magnitude and plot B the first Green-Lagrangian principal strain values.

Table 1: Average image registration results of the LA contour mask across the cardiac cycle. The average and standard deviation is shown for the Hausdorff distance (lower the better, with 0 the best) and Dice score (higher the better, with 1 the best).

Model	Hausdorff Distance (mm)	Dice Score
ANTs	$3.67\pm 1.08$	$0.62\pm 0.08$
Vxm	$4.41\pm 1.17$	$0.70\pm 0.09$
Vxm-seg	$3.72\pm 1.33$	$0.62\pm 0.12$
Aladdin-R	$\mathbf{2.87\pm 0.95}$	$\mathbf{0.87\pm 0.05}$

3.3 LA Strains

The calculated principal strain values are spatially and temporally smooth. Healthy cases (see representative case Video S2 second row) have higher first principal strain values on the anterior wall, consistent with the observed higher DVFs there. The smallest principal strain values can be usually be found on the LA roof and posterior wall. In contrast, CVD cases (see representative case Video S3 second row) show more spatially homogeneous principal strain values.

The principal strain directions predominantly align to the circumferential direction across the cardiac cycle (see Video S2 and Video S3 second row). The alignment with the short axis plane is largest at the end of the reservoir phase.

As with the DVF, we also observe differences in the average Green-Lagrangian principal strain value between the healthy ($0.09\pm 0.05$) and CVD ($0.02\pm 0.04$) groups (see Figure 6B). These are significant during most of the reservoir phase (phases 0-6) and the booster-pump phase(15-18). The healthy cases have marked peak strains at the end of the reservoir phase (phases $\sim$7-9) and at the start of the booster-pump phase (phases $\sim$11-13). On the other hand, CVD cases lack a marked strain peak at the booster-pump phase and have overall smaller principal strain values.

3.4 LA Atlas

We were able to successfully register the LAs of all 10 healthy volunteers to a common space and calculate average DVF and strain biomarkers in this atlas. The resulting DVF and principal strain atlas can be seen in Figure 7 for two cardiac phases and spatial orientations (see Video S4 for across the cardiac cycle).

Atlas

Figure 7: Atlas of the LA Displacements and Principal Strains: The generated DVF (first panel) and strain (second panel) atlas for the anterior (first row) and posterior (second row) views. For each group, the end of reservoir (phase 8) and start of booster-pump (phase 15) phases of the anterior and posterior views are shown.

We obtain an average (across space, time and healthy volunteers) magnitude of $1.38\pm 0.65$ for DVF and $0.09\pm 0.05$ for the principal strain value. The atlas fields preserve the observed properties in the DVFs and strains, showing temporal and spatial smoothness and higher displacements on the anterior wall. The three phases of the LA (reservoir, conduit and booster pump) remain clearly distinct (see Video S4).

The atlas captures the normal values and ranges observed in the healthy cohort and allows to perform a regional assessment on individual cases. We exemplify this application of the atlas by identifying abnormal LA motion signatures on a patient with myocarditis (Figure 8). The regions with the most atypical deformation patterns correspond to the largest MD with respect to the healthy volunteer atlas (see Video S5). In this case, both MDs of the DVF and first principal strain occur on the anterior wall during the end of the reservoir and start of the booster-pump phases. There is also larger MD for both on the posterior wall across the cardiac cycle predominantly due to smaller displacements. While the locations of large MD correspond between the DVF and principal strain, the MD of the DVF typically covers a larger region.

CVD Regional Analysis

Figure 8: The generated DVF (first panel) and strain (second panel) atlas for the anterior (first row) and posterior (second row) views. For each group, the end of reservoir (phase 8) and start of booster-pump (phase 15) phases of the anterior and posterior views are shown.

4 Discussion

In this study, we present a comprehensive technical workflow, Aladdin, for the automatic analysis of LA motion and strains, from high-resolution 3D Cine MR images. The proposed image analysis workflow consists of: a segmentation network, Aladdin-S, to accurately identify the LA contours; an image registration network, Aladdin-R, to extract LA DVFs; strain calculation from the DVFs; and, finally, the creation of a representative DVF and strain biomarkers LA atlas. Our results underscore the potential of this integrated methodology to provide clinicians and researchers with new regional biomarkers of LA motion and deformation likely to have important clinical applications.

As we use a novel high-resolution imaging protocol, we have a comparatively small dataset, unsuitable for traditional AI analysis techniques. Online learning, used in both Aladdin-S and Aladdin-R, is a natural choice to overcome issues of limited data and effectively address the large (compared to LV) structural and functional variability of the LA. As the networks are trained on a subject-by-subject basis, they are immune to biases from the training data and can continue to be employed in future studies to analyse images from subjects with various clinical presentations. The proposed methods are also well suited for future analyses of the motion and deformation of the right atrium or right ventricle (which are often also approximated as thin surfaces).

In this study, we focus on DVFs and principal Green-Lagrangian strain values (with ventricular end-diastole as the Lagrangian reference frame) to characterise LA motion. We could have used other related physical variables instead, representing, for example, strains in an Eulerian frame, or explicitly calculating the rate of change of the strain metrics [8]. In the LV literature, it is also common to compute projections of these variables into a cylindrical reference frame [27]. These metrics can all be easily computed from the metrics we present in this paper, and future studies should identify which biophysical variables are most useful in different clinical applications.

The proposed atlas allows rapid and reproducible assessments of LA DVF and/or strain irregularities, as manual comparison of 3D DVFs and strains is not only labor-intensive but also prone to errors. In the current study, we use the MD to identify regions of the LA with abnormal biomarkers. This identification is likely to be more robust when several deformation metrics (e.g., principal strain value and DVF magnitude) are considered simultaneously, as will be explored in further studies. This analysis can be performed using the proposed atlas or other techniques for analysing regional LA information, such as the Universal Atrial Coordinate (UAC) 3D LA surface mapping system [79]

Regional analyses of LA motion and deformation are likely to be clinically valuable for the identification of dense fibrosis or scar in the LA, as these regions are stiffer than healthy myocardium. Current MRI-based techniques for LA fibrosis identification rely on late-gadolinium enhancement images, whose analysis is highly subjective and has low reproducibility [80]. Future studies will determine whether 3D Cine MRI could be an alternative or complementary technique to identify LA fibrosis, with important clinical applications in AF treatment stratification and the personalization of catheter ablations in AF [81].

4.1 Comparison with the Literature

As full-coverage Cine MRI of the LA are a new and emerging research field, there is limited previous work using this modality [82, 67, 55] to compare to Aladdin. [83] obtained the 3D displacements from high-resolution 3D Cine using the MITT [67] on the LA segmentation maps. [82], on the other hand, focused on the identification of fibrotic tissue using full LA coverage MRI acquired using thick 2D slices. Displacements were estimated using ANTs [68]. In the current study, Aladdin-R led to more accurate DVF estimates than either MITT or ANTS (Figure 5, Video S1, and Figure S2). It also outperformed Voxelmorph [31], a popular registration NN.

Aladdin’s LA motion and deformation measurements are in good general agreement with the literature from 2D STE. Specifically, highest strain values are achieved during the end of the reservoir phase and a smaller peak during the start of the booster-pump phase [84, 85]. Additionally, the largest motion is observed on the anterior rather than the posterior side of the LA wall, and smallest motion on the roof [86]. The principal strain directions are predominantly circumferential throughout the cardiac cycle [84, 85]. We speculate whether Aladdin may additionally contribute to the estimation of LA myocardial fibre directions, which are important determinants of atrial electrophysiology and biomechanics. It would be interesting to compare our imaging-based estimation of the principal strain directions during the active contraction phase with the LA fiber orientation directions proposed using other methods [79, 87, 88, 89].

A comparison between LA deformation metrics derived from Cine MRI and other modalities will be performed to determine 3D Cine MRI’s clinical value. STE is the current gold-standard technique to assess LA motion. 3D STD enables full coverage scanning of the LA, but, when compared to 2D STE, suffers from lower spatial resolution, reduced image quality, and decreased reproducibility [17, 90, 91].

4.2 Limitations and Future Work

Further improvements to Aladdin can be implemented in future studies. We could join the segmentation and image registration tasks, as they are closely intertwined, to improve overall performance as has been shown for the LV [29]. We will explore memory-efficient models that can handle the 3D+t data as inputs to incorporate all cardiac frames into the model, rather than relying on the current pair-wise registration. This will likely improve the temporal consistency of the results. Lastly, estimates may be further improved by using biomechanics-informed regularisation [42], perhaps in a physics-informed NN setting [92, 43]. All these improvements should contribute towards a more reliable DVF and strain estimation.

We plan to update the atlas as more images are acquired. The current atlas of LA motion uses images from 10 subjects, acquired in a single scanner, from a narrow age range and limited ethnicities that are not representative of the general population. As is typical of Cine MRI-based motion analysis, we are not tracking the displacement of material points and use image registration methods with regularization to determine the displacements of the LA wall. This is expected to have little impact in the motion analysis in most circumstances. It may however lead to inaccurate strain estimates in the presence of atrial devices, severe focal atrial pathology or in other circumstances in which the smoothness assumptions of the regularisation techniques may be compromised.

5 Conclusions

We propose Aladdin, a tool to reliably provide 3D DVF and strain maps of the entire LA and map the information to a common atlas space. We expect this will allow the discovery of novel clinical prognostic and diagnostic biomarkers of LA function and improve our understanding of the atria’s involvement in cardiovascular disease processes.

Supplementary Material

Methods

5.1 Aladdin-R: Image Registration Neural Network

Figure S1: The different data input types that were experimented with for image registration methods: (1) 3D Cine MR LA image, (2) LA segmentation map predicted by Aladdin-S, and (3) the LA imaged masked using the contour from the Aladdin-S segmentation map.

Table S1: Summary of the image registration results across the cardiac cycle for input types 1, 2 and 3 (see Figure S1). The average and standard deviation is shown for the Hausdorff distance (lower the better, with 0 the best) and Dice score (higher the better, with 1 the best). (* indicates not suitable input type for specific method)

Input Type	Model	Hausdorff Distance (mm)	Dice Score
1	MITT	$4.76\pm 2.38$	$0.49\pm 0.21$
	ANTs	$3.41\pm 2.38$	$0.67\pm 0.12$
	Vxm	$4.07\pm 1.21$	$0.75\pm 0.07$
	Vxm-seg*	-	-
	Aladdin-R	$\mathbf{3.40\pm 1.28}$	$\mathbf{0.83\pm 0.07}$
2	MITT*	-	-
	ANTs	$3.97\pm 1.74$	$0.59\pm 0.15$
	Vxm	$7.41\pm 1.90$	$0.38\pm 0.17$
	Vxm-seg*	-	-
	Aladdin-R	$\mathbf{3.09\pm 0.90}$	$\mathbf{0.84\pm 0.05}$
3	MITT*	-	-
	ANTs	$3.67\pm 1.08$	$0.62\pm 0.08$
	Vxm	$4.41\pm 1.17$	$0.70\pm 0.09$
	Vxm-seg	$3.72\pm 1.33$	$0.62\pm 0.12$
	Aladdin-R	$\mathbf{2.87\pm 0.95}$	$\mathbf{0.87\pm 0.05}$

5.1.1 Input Types

We experimented with 3 different inputs for the registration methods (see Figure S1: 1) 3D Cine MR LA image; 2) LA segmentation map, and; 3) dilated contour-masked LA image. MITT [67] was evaluated only on the whole LA image as it is based on tracking voxel intensities and therefore not suitable for inputs 2 and 3. Overall, better performance was achieved using input 3, the contour-masked LA images (see Table S1). These were used as the inputs for all subsequent registration tasks.

5.2 LA Atlas Generation

5.2.1 Forward Transformation

The forward transformation ($T_{i}$) step for the atlas generation consists of:

1. An affine registration step: performed using a 2-level scale (downscale factors: ([2, 1]; smoothing sigma: [0, 0]) that minimizes the mean squared error between each case and the reference case, using a gradient descent optimizer;

2. A non-rigid registration step: performed using a 4-level multi-scale (downscale factors: [8, 4, 2, 1]; smoothing sigma: [4, 2, 1, 0]) symmetric forces demons registration algorithm [93] with a gradient descent method and a Gaussian kernel with a standard deviation of 2 voxels as a spatial regularizer.

5.3 LA Strain Calculation

5.3.1 Mesh Post-Processing

The obtained mesh used to calculate the surface principal strains if first post-processed. These steps include: Taubin smoothing for 50 iterations with a passband windowed sinc filter value of 0.01 [94], to homogenise the LA curvature without affecting its volume; linear cell subdivision by a factor of 2; and interpolation of the DVFs in image space across the cardiac cycle to each vertex of the mesh using a Gaussian kernel with a radius of 2 voxels.

Results

5.4 LA Registration

Figure Video S1: (video: https://github.com/cgalaz01/aladdin_cmr_la/figures/sup/la_reg.gif) The comparison of the image registration methods (from top to bottom: ANTs, Vxm-seg, Vxm and Aladdin-R) on registering the non-dilated contour-masked LA image across the cardiac cycle. A single slice of a representative case is shown. The last plot shows the approximate volume ($mm^{3}$) of the same case across the cardiac cycle for each of the methods. The volume is estimated based on the convex hull area of the transformed non-dilated contour-masked LA image.

Comparison of Point Tracking Between Registration Methods

Figure S2: Comparison of the displacement of randomly sampled surface points across the cardiac cycle of a representative healthy case between ANTs (A), Vxm-seg (B), Vxm (C) and Aladdin-R (D). Each line represents the trajectory of a single point across the cardiac cycle. The trajectory lines are colored based on the DVF magnitude.

Figure Video S2: (video: https://github.com/cgalaz01/aladdin_cmr_la/figures/sup/healthy_case.gif) The DVF (first row) and first principal strain (second row) of a representative healthy case across the cardiac cycle. The anterior (first column) and posterior views (second column) are shown. The vectors in the DVF panels show the magnitude (color and size) and direction of displacement. The black lines in the principal strain panels shown the principal strain directions. We observe larger DVF magnitude on the anterior wall rather than the posterior wall, and during the end of the reservoir phase and to a lesser extend at the start of the booster-pump phase. Higher first principal strain values are present during the reservoir phase and mostly on the anterior wall.

Figure Video S3: (video: https://github.com/cgalaz01/aladdin_cmr_la/figures/sup/cvd_case.gif) The DVF (first row) and first principal strain (second row) of a representative CVD case with myocarditis across the cardiac cycle. The anterior (first column) and posterior views (second column) are shown. The vectors in the DVF panels show the magnitude (color and size) and direction of displacement. The black lines in the principal strain panels shown the principal strain directions. We observe larger DVF magnitudes on the anterior wall rather than the posterior wall, and during the end of the reservoir phase. Higher first principal strain values are present during the reservoir phase and mostly on the anterior wall.

5.5 LA Atlas

Figure Video S4: (video: https://github.com/cgalaz01/aladdin_cmr_la/figures/sup/atlas.gif) In the first two columns, DVF and first principal strain of the atlas are shown in the first and second row respectively. Both DVF magnitude and principal strain values are higher on the anterior wall and during the end of the reservoir phase (phase 8) with a spike during the start of the boost-pump phase (phase 15). The last two columns show the Correlation of Variation (CV) of the DVF and first principal strain values respectively in the first and second row. Larger CV is present at regions with small displacement, such as the anterior wall during the start of the reservoir and the end of the boost-pump phases, and the posterior wall.

Figure Video S5: (video: https://github.com/cgalaz01/aladdin_cmr_la/figures/sup/cvd_md.gif) The regional analysis of a representative CVD case with myocarditis against the atlas distribution using Mahalanobis distance (MD) for the DVF (first row) and principal strain values (second row). The anterior (first column) and posterior (second column) view. For both DVF and principal strain, there is higher MD when there are smaller displacements (phases at the start of the reservoir and end of the boost-pump). Across the cardiac cycle, larger DVF MD are observed across the anterior wall and partially on the posterior wall with peaks between the end of the reservoir and start of the booster-pump phases (phases 8-15). The principal strain MD is larger on the central region of the anterior wall during up to the mid reservoir phases (phases 0-6). Additionally, higher strain value MD are observed at regions of the posterior wall between the start of the reservoir and start of the booster-pump phases (phases 0-15).