A three-dimensional left atrial motion estimation from retrospective gated computed tomography: application in heart failure patients with atrial fibrillation

Background A reduced left atrial (LA) strain correlates with the presence of atrial fibrillation (AF). Conventional atrial strain analysis uses two-dimensional (2D) imaging, which is, however, limited by atrial foreshortening and an underestimation of through-plane motion. Retrospective gated computed tomography (RGCT) produces high-fidelity three-dimensional (3D) images of the cardiac anatomy throughout the cardiac cycle that can be used for estimating 3D mechanics. Its feasibility for LA strain measurement, however, is understudied. Aim The aim of this study is to develop and apply a novel workflow to estimate 3D LA motion and calculate the strain from RGCT imaging. The utility of global and regional strains to separate heart failure in patients with reduced ejection fraction (HFrEF) with and without AF is investigated. Methods A cohort of 30 HFrEF patients with (n = 9) and without (n = 21) AF underwent RGCT prior to cardiac resynchronisation therapy. The temporal sparse free form deformation image registration method was optimised for LA feature tracking in RGCT images and used to estimate 3D LA endocardial motion. The area and fibre reservoir strains were calculated over the LA body. Universal atrial coordinates and a human atrial fibre atlas enabled the regional strain calculation and the fibre strain calculation along the local myofibre orientation, respectively. Results It was found that global reservoir strains were significantly reduced in the HFrEF + AF group patients compared with the HFrEF-only group patients (area strain: 11.2 ± 4.8% vs. 25.3 ± 12.6%, P = 0.001; fibre strain: 4.5 ± 2.0% vs. 15.2 ± 8.8%, P = 0.001), with HFrEF + AF patients having a greater regional reservoir strain dyssynchrony. All regional reservoir strains were reduced in the HFrEF + AF patient group, in whom the inferior wall strains exhibited the most significant differences. The global reservoir fibre strain and LA volume + posterior wall reservoir fibre strain exceeded LA volume alone and 2D global longitudinal strain (GLS) for AF classification (area-under-the-curve: global reservoir fibre strain: 0.94 ± 0.02, LA volume + posterior wall reservoir fibre strain: 0.95 ± 0.02, LA volume: 0.89 ± 0.03, 2D GLS: 0.90 ± 0.03). Conclusion RGCT enables 3D LA motion estimation and strain calculation that outperforms 2D strain metrics and LA enlargement for AF classification. Differences in regional LA strain could reflect regional myocardial properties such as atrial fibrosis burden.

LA endocardial surface meshes were generated from the LA segmentations using the Medical Image Registration Toolkit (MIRTK), as shown in Figure 1B.Surface meshes were post-processed using CemrgApp which gave an output depicted in Figure 2A.The sleeve extent of the PVs and LAA were labelled semi-automatically using clippers that were generated using a Voronoi diagram representation of the LA.The PVs and mitral valve (MV) were clipped using spherical clippers defined by the user.
Supplementary Figure 1.Segmentation and surface meshing of the left atrium (LA) from enddiastolic CT image.A: The output segmentation of the LA body, pulmonary veins (PVs) and left atrial appendage (LAA) from a CNN.B: Surface mesh of the LA endocardium extracted from the LA body segmentation.

Universal Atrial Coordinates and Fiber Fields
Universal atrial coordinates (UACs) were registered to each patient's anatomical model by manually selecting 4 landmarks on the endocardium that were used to define the orthogonal axes of the UACs.The landmarks were selected as outlined in Roney et al [28] and were defined as: the LSPV/posterior wall boundary, RSPV/posterior boundary, fossa ovalis and left atrial appendage.
Average endocardial fiber fields from a dataset of seven ex vivo human hearts imaged using diffusion tensor magnetic resonance imaging [29] (DT-MRI) were appended to each model using the registered UACs, shown in Figure 2B.This provided a myofiber vector direction for each element on the model.2. Refinement of the surface mesh of the LA endocardium and registration of a human atrial fiber atlas.A: The LA endocardium surface was refined in CemrgApp: the LAA and all four PVs and their extents were labelled semi-automatically.The mitral valve (MV) and PVs were manually clipped using spheres.B: An atrial fiber atlas from a human ex vivo dataset was appended to each model using the Universal Atrial Coordinates (UACs).The UAC registration enabled each surface elements to be assigned an estimated endocardial myofiber orientation, denoted by the grey arrowed glyphs.

Regional Separation of the LA Body using UACs
The LA body was separated into 5 different regions: posterior wall, septum, lateral wall, anterior wall and inferior wall using the UACs.The pulmonary vein positions were used as landmarks to create the borders between the 5 regions, as depicted in Figure 3. 3. Separation of the LA body in 5 regions (posterior, septum, lateral, anterior and inferior walls) using the UACs.A: The pulmonary veins (PVs) were used as landmarks to delineate the LA body into 5 regions, denoted by different colours.The PVs and LAA area shaded light gray to show these regions were omitted from the strain calculations.This ensured only strains from the LA body contributed to strain measurement.B & C: The 5 regions viewed on the LA anatomy.Abbreviations: LAAleft atrial appendage, LSPVleft superior pulmonary vein, LIPVleft inferior pulmonary vein, RSPVright superior pulmonary vein, RIPVright inferior pulmonary vein, MVmitral valve.

Imaging Data
Retrospective electrocardiogram (ECG)-gated computed tomography (CT) imaging was indicated for each patient prior to cardiac resynchronization therapy (CRT).Ten or twenty frame CT images, at 10% or 5% intervals between the R-R interval and starting from the R-wave, were created for each patient.In this study, ten CT image frames were used to measure LA motion and strain.The dimension of images varied between patients, however the majority of cases had in-plane shape of 512 x 512 voxels with a range of 0.32-0.49mm isotropic resolution, and a range of 94-399 slices with a through-plane thickness of 0.4-1.2mm.

Feature Tracking
This section outlines details on feature tracking usage, optimization and verification steps.

Left Atrial Feature Tracking in Retrospective Gated CT
The task of LA feature tracking in retrospective gated CT images can be re-framed as non-rigid image registration.To track LA features, all later CT image frames need to be registered to the initial end-diastolic (ED) CT image.This provides displacement fields to deform the LA model from ED to all the later time frames in the cardiac cycle.Strains can then be calculated on the deformed model.
In this study, we firstly optimized and applied the temporal sparse free form deformation (TSFFD) method [17] for LA feature tracking in retrospective gated CT.
Images were cropped to cover the LA body, PV ostia and LAA prior to feature tracking.This ensured that the images fit within computational memory constraints, and avoided spatial resolution downsampling which ensured sufficient image resolution of the LA was achieved.

Temporal Sparse Free Form Deformation Method
The TSFFD method was used as it is based upon the widely used free form deformation method, with additional sparsity and cyclicity constraints which encourage temporally and spatially smooth displacement fields that are appropriate for the application of cardiac motion measurement.It has previously been applied to LV feature tracking in retrospective gated CT.The TSFFD method performs nonrigid registration of all later time frames to the initial end-diastolic CT frame simultaneously using multiple levels of control points with increasing spatial resolution.Control point (CP) displacements are parametrized using cubic B-spline functions in both space and time which encourage spatially and temporally smooth displacement fields.The sparsity constraint encourages sparsely activated CPs.

TSFFD Optimisation
The workflow for the optimisation of the TSFFD image registration method is outlined in Figure 4.The sparsity weight (SW) and bending energy (BE) parameters, which parametrize the sparsity constraint and stiffness of deformation respectively, were optimized using a grid search through 63 different SW and BE combinations.The SW and BE values that were iterated through are displayed in Figure 5.The grid search was performed at a coarser and finer set of multi-level control point resolutions, which are outlined in Table 1.
Supplementary Table 1.Two sets of control point (CP) resolutions at which hyperparameters for the TSFFD method were searched through.A coarser and finer set of CP resolutions were tested.The LA feature tracking accuracy was quantified using three error metrics based upon tracking LA features from end-diastole (ED) to end-systole (ES).The ED-to-ES registration was chosen since it involves the largest change of shape of the LA and will therefore generate the largest errors in LA feature tracking.Surface meshes and segmentations of the LA body at ED were deformed to ES using the TSFFD deformation field and compared with the ground truth surface meshes and segmentations created from the target image at ES.A combination of the ASD, DHD and DSC errors with equal weighting was used to identify the optimal TSFFD hyperparameter combination for LA feature tracking.To combine the error metrics, errors were first normalized by calculating the Z-score, which represents the number of standard deviations away from the mean across the cohort:

𝑍 = 𝑥 − 𝜇 𝜎
For each hyperparameter combination, the Z-scores of each error metric were combined to give the unified score, Lower values of   represent more accurate LA feature tracking across the cohort.Five-fold cross validation was used with 25 patient cases to select the optimal hyperparameter combination.This ensured the optimal hyperparameter combination was robust to which cases were used to calculate the LA feature tracking errors using   .Results of the 5-fold cross validation used for the feature tracking optimization are shown in Figure 5.The hyperparameter combination which gave the lowest   score was SW = 0.0, BE = 1e-9.This combination gave an   score of -1.47.A hold-out test set of 5 patient cases was used to evaluate the accuracy of the feature tracking on an unseen dataset.Error metrics for the test set are shown in Figure 6.The mean ± standard deviation of the test set errors were as follows, ASD: 0.52 ± 0.16 mm, DHD: 4.22 ± 2.32 mm, DSC: 95 ± 2 %.

Feature Tracking Verification
Synthetically generated gated CT images with known LA deformation were used to verify the optimized feature tracking method.Calculated displacement fields using TSFFD for each patient were applied to the ED image to create a set of synthetic gated CT images.The synthetic CT images were then re-registered to yield inferred LA motion which was compared with the ground truth LA motion.This enabled the calculation of errors for material point tracking, and global and regional strain calculation.Synthetic CT images were created across the patient cohort, therefore our verification method captured inter-individual variability in atrial size and motion.
The Euclidean distance between inferred and ground truth positions for the LA endocardial surface vertex points was calculated using the ED to ES registration to give a measure of material point tracking error.The median and mean ± standard deviation of the root mean square error (RMSE) of the Euclidean distance between the inferred and ground truth positions were 0.4 mm and 0.6 ± 0.5 mm, respectively across the cohort.This suggests that our feature tracking method was able to predict vertex coordinates on the LA endocardial surface within roughly one voxel resolution.
The synthetic gated CT images enabled the comparison between inferred and ground truth strain measurements, as depicted in Figure 7. Ground truth (solid lines) and inferred (dotted lines) global and regional strain transients were compared for each case.
Supplementary Figure 7. Synthetic CT images with known LA deformation enabled the comparison of inferred and ground truth strain measurements.Shown above are inferred (dotted lines) and ground truth (solid line) strain transients for a sample HFrEF only case.The shape of the strain transients are similardifferences in strain transient shape can be quantified using the root mean square error (RMSE).Absolute differences in reservoir strain measurements can also be quantified.
Differences between the ground truth and inferred strain transients are described using the root mean square error (RMSE).The RMSE values across the cohort are outlined in Table 2 for global and regional strain transients.Absolute differences in reservoir strain measurements across the cohort are outlined in Table 3.

Two-dimensional Feature Tracking
Two-dimensional (2D) global longitudinal strain (GLS) was calculated by extracting 2-and 4chamber views from the RGCT images using manual landmarks, as shown in Figure 8.The coupled cross-hair rotation tool in CemrgApp was used to align imaging planes with the 2-and 4-chamber views.The TSFFD method was used to track LA features in the 2D images.
Where () represents the element area at time t, and ( = 0) represents the initial element area at end-diastole.Positive area strain indicates endocardial area stretch and negative area strain indicates endocardial area contraction.
Endocardial fiber strain denotes deformation parallel to the local endocardial fiber direction.The calculation for fiber strain,   , is described as follows.First, the green strain tensor,  , for each element was calculated from the deformation gradient,  , using Fiber strain is calculated by first finding the green strain,   , by projecting the green strain tensor onto the direction parallel to the element myofiber orientation,  : Green strain was then converted to linear strain to calculate fibre strain,   :

Difference between Strain Metrics
Area strain does not account for the direction of deformation, unlike fiber strain which measures deformation with respect to the local myofiber direction.Figure 9 depicts the differences between area and fiber strains using a sample of a LA endocardial surface.
Supplementary Figure 9. Depiction of the different information that can be extracted from the area strain and fiber strain metrics.A subset of the left atrial endocardium surface is shown at its end-diastolic (left) and end-systolic (right) configurations.Individual elements are outlined in pink.Area strain and fiber strain were both calculated, with red and blue colors representing positive and negative strain, respectively.The local myofiber orientation for each element is shown by the arrowed glyphs.Top row, right: Area strain is colored red for all elements which reflects the uniform increase in area.Bottom row, right: Fiber strain is shaded light red for two elements and light blue for the remaining elements since deformation along the fiber directions are increasing and decreasing, respectively.

Reservoir Strain Measurement
To avoid artefacts in strain measurements due to pacing leads in the gated CT imaging, the top 0.1% of elemental strains were omitted from the strain analysis when taking mean averages over elements for regional and global reservoir strains.This circumvented extreme unphysiological strains from using the feature tracking with imaging that included lead artefacts from CRT patients.

Two-dimensional Global Longitudinal Strain Measurement
2D GLS was measured in both the 2-and 4-chamber views using the percentage length changes of the LA body contour, excluding the MV, PVs and LAA.Final values for 2D GLS were taken using the mean strain value from both the 2-and 4-chamber views.Visualisations of LA body contours in the 4-chamber view is displayed in Figure 10.
Supplementary Figure 10.2D GLS calculated from the 4-chamber view for a sample case.A: Enddiastolic (t=0%) frame and overlaid LA body contour, exlcuding the PVs, MV and LAA.B: Endsystolic (t=40%) frame and overlaid tracked LA body contour at t=40% shown in green, which is larger than the original LA body contour, shown in white, from the end-diastolic frame.
Surface mesh-based and segmentation-based errors were calculated: the average surface distance (ASD) and directed Hausdorff distance (DHD) and the dice score coefficient (DSC).ASD is the mean normal distance between the tracked and ground truth meshes at ES, and DHD is the maximal minimum distance between points in the tracked and ground truth ES meshes.Lower values for ASD and DHD indicate more accurate LA feature tracking.DSC represents the overlap in LA body segmentation between the tracked and ground truth LA segmentations at ES. Higher values for DSC indicate more accurate LA feature tracking.Supplementary Figure4.Schematic that depicts quantification of LA feature tracking accuracy for optimization.A: Anatomical information of the LA body at the end-diastolic CT frame was captured using an endocardial surface mesh and segmentation.B: Surfaces and segmentations were deformed from end-diastole (t=0%) to end-systole (t=40%) using the deformation calculated from the feature tracking from the retrospective gated CT images.C: The tracked surface (in blue) was compared with the ground truth surface (in transparent grey) and the tracked segmentation (solid) was compared with the ground truth segmentation (transparent) end-systole.D: Differences in surfaces and segmentations at end-systole were quantified using surface-based errors (average surface distance and directed Hausdorff distance) and volume-based errors (dice score coefficient), respectively.

Figure 5 .
Results of the feature tracking optimization using 5-fold cross validation for the Bending Energy and Sparsity Weight hyperparameters.Lower values of   , indicated by darker blue shaded squares, represent more accurate LA feature tracking in the CT images.

Table 2 .
Statistics of root mean square error (RMSE) values between the ground truth and inferred strain transients.Statistics are given per patient group.

Table 3 .
Statistics of absolute error in reservoir strain measurements between ground truth and inferred reservoir strains.Statistics are given per patient group.