Myocardial strain computed at multiple spatial scales from tagged magnetic resonance imaging: Estimating cardiac biomarkers for CRT patients

HighlightsA framework to analyse cardiac strain at multiple spatial scales using a motion atlas is proposed.Optimal scales are found for the identification of clinical biomarkers for cardiac resynchronisation therapy patients.Combining strains from multiple scales improves the accuracy of CRT response prediction. Graphical abstract Figure. No Caption available. Abstract Abnormal cardiac motion can indicate different forms of disease, which can manifest at different spatial scales in the myocardium. Many studies have sought to characterise particular motion abnormalities associated with specific diseases, and to utilise motion information to improve diagnoses. However, the importance of spatial scale in the analysis of cardiac deformation has not been extensively investigated. We build on recent work on the analysis of myocardial strains at different spatial scales using a cardiac motion atlas to find the optimal scales for estimating different cardiac biomarkers. We apply a multi‐scale strain analysis to a 43 patient cohort of cardiac resynchronisation therapy (CRT) patients using tagged magnetic resonance imaging data for (1) predicting response to CRT, (2) identifying septal flash, (3) estimating QRS duration, and (4) identifying the presence of ischaemia. A repeated, stratified cross‐validation is used to demonstrate the importance of spatial scale in our analysis, revealing different optimal spatial scales for the estimation of different biomarkers.


Introduction
Cardiac deformation is driven by the electromechanics and perfusion of the myocardium, and is intrinsically tied to the disease state of the heart. As such, cardiac deformation has been widely analysed for the characterisation and detection of different diseases, both in more clinical ( Hsu et al., 2012;Jackson et al., 2014;Sohal et al., 2014 ) and methodological Peressutti et al., 2017 ) work.
An important determinant of cardiac function, and in-turn deformation, is the inter-play of multi-scale mechanisms controlling electrical conduction, perfusion and contraction of the heart . Multi-scale analysis spans fields of research including, but not limited to, computer vision ( Lindeberg, 1994;Starck et al., 1998 ), materials and structures ( Fish and * Corresponding author. E-mail address: m.sinclair@imperial.ac.uk (M. Sinclair).
Shek, 20 0 0 ) and biology ( West et al., 1997;Lesne, 2013 ). While the coronary arterial network and cardiac conduction system (CCS) display multi-scale properties ( Bassingthwaighte and van Beek, 2002;Sinclair et al., 2015;Sebastian et al., 2013 ), a multi-scale analysis of cardiac deformation, which is driven by these systems, has not been extensively investigated in the literature.
In this paper we build on our recent work on multi-scale analysis of cardiac strains  within the established framework of a cardiac motion atlas, and apply our analysis to the prediction and estimation of several cardiac biomarkers for cardiac resynchronisation therapy (CRT) patients. We start by reviewing aspects of multi-scale structure in the heart and its relation to disease, motivating our work for a multi-scale analysis of cardiac deformation. A description and overview of motion atlases in the literature is then presented, highlighting the advantages of highfidelity imaging data in the assessment of cardiac deformation. We then present the motivation of assessing CRT patients as a target cohort for our analysis. The left-most box shows primarily the large epicardial coronary arteries in red, and the right-most box includes many of the intramural arterioles. Image adapted with permission from Sinclair et al. (2015) . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Multi-scale cardiac structure
The branching structure of the coronary vasculature follows power law relationships ( Bassingthwaighte and van Beek, 2002 ), and vessel generation has been shown to follow a power law with downstream myocardial tissue volume ( Marxen et al., 2009;Sinclair et al., 2015 ). To help illustrate the relationship of vessel generation and the downstream tissue volume supplied by these vessels, Fig. 1 shows the increasing level of detail in the coronary arterial circulation in the same porcine heart as more vessel generations are considered. The mean tissue volume, V , supplied by the terminal vessels is indicated above each box, and is given in terms of the mean tissue volume supplied by the terminating vessels of the large coronary arteries, V 1 , in the left-most box. Given that the porcine coronary arterial network branches almost exclusively with bifurcations (at 98% of junctions) ( Kassab et al., 1993 ), the mean volume V approximately halves with the addition of each further generation of vessels. Thus, with the addition of every two generations, as in Fig. 1 between each box, the mean tissue volume supplied by each terminal vessel is approximately quartered.
Disease in the coronary circulation manifests at different generations of the coronary arteries, visualised in Fig. 1 , and also at the microvascular scale ( Camici and Crea, 2007 ). Similarly, block of the CCS can occur at any scale of the conduction network ( Park and Fishman, 2011 ). Depending on the vessel generation at which disease has affected the coronary circulation, the downstream myocardial tissue will be affected by perfusion abnormalities, in turn resulting in cardiac deformation abnormalities associated with the spatial scale represented by the supplied volume of tissue. This suggests that, depending on the disease, there may be a characteristic spatial scale at which the deformation is most altered from a healthy case, and may be most predictive for different clinical applications.
Spatial scale in the analysis of cardiac deformation has not been extensively investigated, despite the importance of scale in cardiac structure and function. Clinical analysis of regional characteristics in the LV is often performed using a standardised 17 segment AHA model of the LV ( Cerqueira et al., 2002 ). AHA segments are identified based on anatomical landmarks, and have coarse correspondence with underlying coronary perfusion territories. While analysis using AHA segments is easy to interpret, averaging cardiac deformation features into fixed, coarse segments may not be the best way of identifying complex deformation abnormalities.

Cardiac motion atlases
Earlier cardiac motion studies with a clinical focus often involved the estimation of global deformation parameters. In Park et al. (1996) , a finite element model of the LV was deformed using tagged magnetic resonance (T-MR) data and parameters corresponding to radial and longitudinal shortening as well as twisting were estimated. In Bardinet et al. (1996) , the evolution over time of LV wall thickness, cavity volume and twist were computed from a deformable model fitted to sequences of single photon emission computed tomography and computed tomography (CT) data.
In more recent works, motion atlases emerged as a versatile approach to analyse cardiac motion, allowing for more spatially and temporally localised cardiac deformation features to be assessed. Cardiac motion atlases have been used for characterising the motion of varied patient cohorts ( Hoogendoorn et al., 2013 ), for abnormal motion characterisation to identify disease ( Duchateau et al., 2011;De Craene et al., 2012 ) and for prospectively predicting treatment outcome ( Peressutti et al., 2015;. A cardiac motion atlas requires the normalisation and reorientation of subjects' cardiac geometry and motion both spatially and over time, thus removing biases due to heart orientation, size, shape and cardiac phase. This normalisation creates a coordinate system within which motions (i.e. displacements and velocities) and deformations (i.e. strains and strain rates) can be directly compared between subjects at corresponding anatomical points over time. This framework allows for a more flexible assessment of regional LV function compared to the consideration of global and AHA segment deformations.
The comparison of cardiac motion across cohorts has been facilitated in recent years by increased accessibility to cardiac imaging data from modalities including ultrasound (US), MR, and CT. The richness of the motion features encoded in a cardiac motion atlas depends on the quality and spatial/temporal resolution of the imaging data. For example, atlases built from CT data typically suffer from a limited temporal resolution but high spatial resolution due to the use of ionising radiation. An example of a CT-based motion atlas is Ardekani et al. (2009) , in which the end-diastolic (ED) and end-systolic (ES) phases alone were used to distinguish healthy subjects from patients with ischaemic heart disease. In Hoogendoorn et al. (2009) ; 2013 ) 15 phases were used to build a 4D (i.e. 3D+t) cardiac motion atlas, although the number of frames was still considered a limitation in terms of temporal resolution.
By contrast, US data typically has a high temporal resolution but low signal-to-noise ratio. Motion atlases built from US data have largely been constrained to 2D views due to the lower image quality of 3D views. Duchateau et al. (2011) proposed a motion atlas formed from 2D US image sequences to characterise motion abnormalities in CRT patients. This work was extended using a constrained manifold learning framework , with a further multi-scale function formulation introduced to address the issue of embedding new patients and healthy volunteers in a manifold created from a separate subgroup of patients . In Duchateau et al. (2014) 2D spatiotemporal maps based on displacements, velocities and strains derived from 2D speckle-tracking were used to identify septal flash in animals with strict left bundle branch block (LBBB). While producing important findings, their approach was limited to motion analysis from a 2D view, missing more complex 3D motion features which may be important for characterising different diseases.
Using short-axis (SA) cine-MR data, Perperidis et al. (2005) proposed a voxel-based 4D cardiac motion atlas built from 26 healthy volunteers' MR image sequences to classify hypertrophic cardiomyopathy patients. However, cine-MR images suffer from low through-plane resolution, although Perperidis sought to circumvent this issue using a shape-based interpolation. T-MR images, by contrast, have high resolution both temporally and spatially, and are designed specifically for tissue tracking and strain analysis ( Osman et al., 1999 used T-MR images to create a 4D motion atlas to characterise abnormal left ventricle (LV) velocity patterns in patients with heart failure compared to healthy subjects. Also using a 4D motion atlas built using T-MR, Peressutti et al. (2015) proposed an ensemble learning approach using random projections of LV displacements to predict CRT super-response in a cohort of 23 CRT patients. Additionally, a dictionary learning framework using local motion descriptors (based on LV displacements) was proposed to determine the location of LV infarct . Displacements, velocities and strains from a motion atlas built using T-MR were more recently combined with clinical information in a multi-kernel learning (MKL) framework to improve the prediction of CRT response .
In this study, we consider strains in a 4D motion atlas based on T-MR data. Strain represents the degree of local tissue deformation and therefore assumes an implicit spatial scale for its computation. This scale has been mostly overlooked so far but we introduce a method for its analysis. We explicitly investigate the issue of spatial scale in characterising cardiac strain with respect to different clinical biomarkers in a cohort of CRT patients.

Cardiac resynchronisation therapy
CRT is used to treat patients with electro-mechanical dyssynchrony, which diminishes systolic function and can result in heart failure ( Kirk and Kass, 2013 ). The current clinical selection criteria include a New York Heart Association functional class of II to IV, a QRS duration > 120 ms , and a LV ejection fraction (EF) ≤ 35%. However, these criteria result in approximately 30% of patients failing to improve after treatment ( Daubert et al., 2012 ). There is an ongoing need for better predictors of patient response, and as recent studies have demonstrated, CRT response is influenced by a range of different clinical indicators ( Parsai et al., 2009;Jackson et al., 2014;Sohal et al., 2014 ).
CRT patients can present with a range of symptoms which may manifest as mechanical abnormalities determined from imaging and electrical abnormalities determined from an electrocardiogram (ECG). Electrical abnormalities are the result of disease of the conduction system, while mechanical abnormalities may result from other underlying disease such as ischaemia or indeed conduction blocks. In the presence of LBBB, it is fairly common for a mechanical abnormality known as septal flash to be present, which is easily identifiable by 'eyeballing' in echocardiography ( Corteville et al., 2017 ). Septal flash is the result of what has been described as a 'true left bundle branch block' ( Corteville et al., 2017 ), and its presence has been shown to be a strong predictor of response to CRT ( Parsai et al., 2009;Doltra et al., 2014 ). Septal flash can be identified by a fast contraction/relaxation of the septal wall inward then outward from the LV bloodpool during isovolumetric contraction, i.e. within the QRS width. The septal wall is then pas-sively stretched as the lateral LV wall contracts against an increase in intraventricular pressure during systole, resulting in an inefficient, dyssnychronous pumping mechanism. CRT usually helps correct septal flash by circumventing the LBBB with pacing leads and restoring synchrony to LV contraction, thus improving the pump function and causing reverse remodeling . However, septal flash alone is not a sufficient criterion for predicting CRT response, and patients with other motion abnormalities have also been shown to benefit from CRT ( Parsai et al., 2009;Doltra et al., 2014 ). Patients treated for heart failure with CRT also have a more favourable outcome in the absence of ischaemic cardiomyopathy .
Identifying septal flash has been a subject of several of the aforementioned motion atlas studies ( Duchateau et al., 2011;2012;2013;. While septal flash is characterised by a well-defined large-scale deformation abnormality, it is yet unclear whether other symptoms such as ischaemia and other forms of conduction block may manifest in some consistent way which predisposes a patient to respond to CRT, and which may be detectable through the associated motion abnormalities. In this study we demonstrate that by considering strains at multiple spatial scales, we are able to identify different scales at which septal flash and ischaemia can be best detected. We also show that CRT response can be best predicted by combining strain information at multiple spatial scales, and with comparable performance to other recent studies.

Our contributions
We build on our recent work , in which we computed strain at different spatial scales in the LV, and subsequently used a motion atlas and dimensionality reduction to identify the spatial scales at which myocardial strain was most strongly predictive of CRT response and LBBB. In addition to a more indepth explanation of the proposed multi-scale strain, the novelties we present in this study include: (1) extending the CRT patient cohort from 34 to 43 patients, (2) the comparison of our multiscale strain computation to a standard method for computing myocardial strain, (3) the application of our framework to identifying patients with ischaemia and septal flash, and the estimation of QRS duration from strain, (4) the combination of different scales in the estimation of clinical biomarkers, and (5) the comparison to a moving average method for computing multi-scale strain.
This paper is structured as follows: Section 2 details the data acquisition protocol, Section 3 explains the creation of the motion atlas, the multi-scale strain computation, and the subsequent dimensionality reduction and classification. In Section 4 details of experiments performed and each clinical biomarker assessed are presented, with results given in Section 5 followed by discussion and conclusions in Section 6 .

Imaging data
Imaging and clinical data were collected for 43 CRT patients prior to treatment at St Thomas' Hospital, London. The study was approved by the institutional ethics committee and all patients gave written informed consent. 1 All patients underwent cardiac MR imaging using a 1.5T scanner (Achieva, Philips Healthcare, Best, Netherlands), with the acquisition of ECG-gated, breath-hold cine-MR, T-MR (3D-tagged) and DE-MR (delayed-enhancement) sequences. A single multi-slice short axis (SA) and three singleslice long axis (LA) cine-MR sequences were acquired (TR/TE =  ( ms ) 147.8 ± 21.6 -3.0/1.5ms, flip angle = 60 °). Slice thickness was 8mm and 10mm for SA and LA sequences respectively, with an in-plane resolution for both of ≈ 1.4 × 1.4 mm 2 . Three orthogonal T-MR sequences were acquired (TR/TE = 7.0/3.2ms, flip angle = 19 − 25 • , tag distance = 7mm, slice thickness = 7mm) and combined to produce a 3 D + t image with ≈ 1.0 mm isotropic resolution. The 3D T-MR had a reduced field-of-view (FOV) relative to the cine-MR, typically tightly enclosing the left ventricle. On average 22 T-MR frames were acquired per subject, with an average time of 50ms between each frame.
DE-MR images were acquired between 15 and 20 min after the administration of 0.1 to 0.2mmol/kg gadopentate dimeglumine (Magnevist, Bayer Healthcare, Dublin, Ireland) using conventional inversion recovery sequences. A multi-slice SA and three singleslice LA 2D images were acquired (TR/TE = 5.6/2.0ms, flip angle = 25 °). The same FOV and orientation as the cine-MR sequences was used.
The SA and LA cine-MR sequences were spatially aligned to the T-MR coordinate system, compensating for motion occurring between sequential breath-holds. The T-MR sequence was chosen as reference as it was free from respiratory motion artifacts. The ED frame of the cine-MR sequence was used to segment the myocardium given the high-resolution detail of myocardial boundaries. The T-MR images were used to estimate cardiac deformation, and the DE-MR images were used to quantify myocardial scar. Table 1 lists the patient biomarkers which were investigated. CRT response was determined from a 6-month follow-up after baseline examination. 2D echocardiography imaging was performed both at baseline and at 6-month follow-up to compute end-diastolic volume (EDV) and end-systolic volume (ESV), where the volumetric response of a patient to CRT is given by the improvement in ESV post-CRT relative to pre-CRT: (ESV pre − E SV post ) /E SV pre ≤ −15% . 13 of the 43 patients treated did not exhibit a volumetric response, and 25 patients presented with septal flash. The number of patients with an ischaemic aetiology was 17, based on a scar burden > 15% determined from the DE-MR images by a clinical expert using the cmr42 software (Circle Cardiovascular Imaging Inc.). The mean and standard deviation of the patient QRS durations ( QRS d ) are also reported in Table 1 .

Spatio-temporal motion atlas
A pipeline for constructing a 4D cardiac motion atlas has been detailed in previous studies Peressutti et al., 2017 ) and is briefly summarised here. The steps involved include (1) estimation of cardiac geometry, (2) motion-tracking, (3) computation of deformation features, and (4) transformation of geometry and deformation features to an unbiased atlas. Following step (4), dimensionality reduction is often used to transform the data to extract a suitable number of important features for further analysis, such as the use of linear embeddings like PCA ( Perperidis et al., 2005 ), non-linear manifold embeddings , or a more stochastic approach such as random projections ( Peressutti et al., 2015 ). The novelty of this study with regards to the above pipeline lies in (3), where detail is provided for the computation of cardiac strain at multiple spatial scales in Section 3.1.3 .

Estimation of LV geometry
In order to get an accurate representation of the LV myocardium, both the ED SA stack and 2 −, 3 − and 4 −chamber LA images were manually segmented. Since the SA and LA slices had been rigidly registered to the T-MR volume, artifacts due to image acquisition under different breath-holds were minimised. We took advantage of the high in-plane resolution of the SA and LA images in their respective orientations by fusing the segmentations and manually smoothing the resulting mask using ITK-SNAP ( Yushkevich et al., 2006 ), as illustrated in Fig. 2 , visualised with Eidolon 2 ( Kerfoot et al., 2016 ). A statistical shape model (SSM) surface mesh consisting of approximately 22,0 0 0 vertices was optimised to fit to the smoothed LV mask following an initial rigid registration using anatomical landmarks . The use of landmarks ensured that the LV mesh was aligned to the same anatomical features for each patient, and the endocardial and epicardial surfaces of the mesh were smoothly fitted to the myocardial mask . In order to reduce the number of vertices of the SSM surface mesh, a medial surface with regularly sampled vertices ( ≈ 3, 0 0 0) was generated using a combination of ray-casting and homogeneous downsampling followed by cell subdivision. Point-correspondence was retained using a consistent approach applied to the original SSM vertices for each patient. Note that for many patients, the 3D field-of-view (FOV) of the T-MR data did not entirely enclose the fitted SSM. Typically around the apical and basal regions a fraction of the medial surface mesh was not enclosed. Thus, these points were not considered for any patients in the analysis.

Motion tracking
The high resolution 3 D + t T-MR sequence was used for motion tracking. DICOM header information was used to determine the fraction of the cardiac cycle over which each 3 D + t T-MR sequence was acquired. Specifically, the trigger time of each image relative to the mean RR-duration recorded for each sequence was used to determine the temporal position of each frame in the cardiac cycle. A 3 D GPU-based B-spline free-form deformation (FFD) registration was used ( Rueckert et al., 1999 ) to estimate LV deformation between consecutive frames of the T-MR sequence. Subsequently, the inter-frame transformations were composed to estimate deformation between each time frame and the ED (first) time frame, producing a 3 D + t B-spline transformation, ψ. In order to compare cardiac phases between patients, the reference ED medial surface was warped using ψ over the cardiac cycle at 30 equally spaced cardiac phases for each patient, where missing frames towards the end of the cardiac cycle were interpolated between the last available frame and the ED state. Across the 43 patient cohort, the available fraction of the cardiac cycle from each T-MR sequence was typically at least 70%, so the first 21 frames of each transformation were used. The patient-specific LV deformations were therefore fully represented by their respective meshes deformed to each of the first 21 frames.

Multi-scale strain
In addition to producing a medial surface, the myocardial volume enclosed by the fitted SSM was sampled in a regular grid with half the resolution of the T-MR ED volume (i.e. with an isotropic  spacing of ≈ 2mm). This produced a LV point-cloud represented in Cartesian coordinates, P pc ∈ R N T ×3 , with the total number of points sampled in the LV point cloud N T ≈ 30, 0 0 0 for each patient, as illustrated in Fig. 3 . The 3 D + t motion transformation ψ was applied to the point-cloud to transform it to each of the first 21 frames of the cardiac cycle for each patient. We denote the first time frame (at ED) by t = 0 , and the time frames which make up the remaining 70% of the cardiac cycle by t = [1 , . . . , 20] .
Six spatial scales at which to analyse strain were then chosen following a power law, (1) corresponding to 1%, 2%, 4%, 8%, 16%, 32% of the total LV volume. At time t = 0, at each point, i , on the medial surface, P m i,t=0 , at each scale, s , a neighbourhood of K s nearest-neighbour points in the point-cloud were selected, Note that the subscript k is related to the spatial scale s via the number K s of nearest points at that scale. Since the point-cloud P pc is sampled in a uniform grid, the number of points K s in a given nearest-neighbour point-cloud, { P pc i,k,t=0 } , is proportional to the per-centage volume at each spatial scale in Eq.
(1) , that is: or between ≈ 300 at the V s = 1% scale and ≈ 10, 0 0 0 at the V s = 32% scale. Fig. 4 illustrates the shape and extent of the neighbourhoods at different spatial scales. The method for computing strain at each medial surface point P m i,t from its deforming point-cloud From large deformation mechanics, the deformation gradient tensor F maps the relative spatial position of two neighbouring particles before deformation ( d X ) to their relative spatial position after deformation ( d x ) ( Bonet and Wood, 2008 ). The mapping from the relative position at an undeformed state ( d X ) to that at a deformed state ( d x ) can be expressed as d x = F d X , from which the Green-Lagrange strain can be computed as E = 1 We cast this general representation into the context of computing strain at medial surface points, i , with respect to point-cloud neighbourhood points, k , at times, t . At ED, considering a point on the medial surface P m i,t=0 and its point-cloud neighbourhood at a given scale s , { P pc i,k,t=0 } , k = 1 , . . . , K s , the undeformed difference vector between P m i,t=0 and each point-cloud neighbour is expressed as, The deformation gradient F i,t,s ∈ R 3 ×3 at medial surface point, i , time, t , and scale, s , is then computed from the least-squares minimisation over all of the K s points in the neighbourhood, The least-squares minimisation gives more weighting to points in the neighbourhood which are further away from the medial surface point, since those points contribute larger values to d x i, k, t and d X i, k . This smoothly captures the intrinsic tissue deformation at each scale by including in the computation, but giving less weighting to, neighbourhood points from smaller scales. There are similarities in this approach to the use of a least-squares approximation to compute affine parameters from a local velocity field in McLeod et al. (2015) . Neighbourhood strain is computed using the Green-Lagrange strain tensor, For each patient, this computation was performed for every time frame, t , at every medial surface vertex, i , and for neighbourhoods at each spatial scale, s . Compute time for strain at all medial surface vertices for a given patient for a single frame ranged from approximately 3 s at V s = 1% ( ≈ 300 points) to 100 s at V s = 32% ( ≈ 10, 0 0 0 points) on 8 CPUs, i.e. computation time scaled in proportion to the size of the point-cloud neighbourhood, due to the least-squares computation in Eq. (6) .

Spatial normalisation
Spatial normalisation accounts for the differences in patientspecific LV geometries to enable an unbiased comparison of LV deformation. Strains were reoriented from the patient-specific to the atlas coordinate space, similarly to how displacements ( Perperidis et al., 2005 ) and velocities Duchateau et al., 2012 ) have previously been reoriented. Introducing the subscript, n , to denote each patient, the Green-Lagrange strain tensor in atlas space, E atlas i,t,s,n , was computed from the reoriented deformation tensor F atlas i,t,s,n = J i,φn F pat i,t,s,n J −1 i,φn , where J i,φn is the Jacobian at the location of the medial surface point, i , of the patient-to-atlas transformation φ n , following the derivation in Peressutti et al. (2017) .
Finally the reoriented strain E atlas i,t,s,n was projected into a local atlas coordinate system in radial, r , longitudinal, l , and circumferential, c , directions by the transformation E atlas i,t,s,n = P T i · E atlas i,t,s,n · P i , where P i = [ r i , l i , c i ] , a row-wise concatenation of the unit column vectors in the radial, longitudinal and circumferential directions at each medial surface point, i . The cylindrical components of strain were taken as the diagonal terms of E atlas i,t,s,n , giving e atlas i,t,s,n = [ e r i,t,s,n , e l i,t,s,n , e c i,t,s,n ] , for each patient, n , scale s , cardiac phase, t , and medial surface point, i .

Dimensionality reduction and classification
The local strains in atlas space e atlas i,t,s,n were concatenated into a single row vector for each spatial scale, s , separately, such that for patient n , ˆ e n,s ∈ R 1 ×M , where M = (3 × T × N m ) is a scalar with T the number of cardiac phases and N m the number of points in the atlas medial surface mesh. The row vectors ˆ e n,s for each patient were stacked to produce a matrix where N p is the number of patients. PCA was subsequently used to reduce the dimensionality of the data, and linear discriminant analysis (LDA) was used for classification tasks and linear regression for regression tasks of different clinical biomarkers. Details of the experiments performed are given in Section 4 .

Comparison of strain methodology
To provide a comparison to an existing method, we also computed strain as performed in a recent study , which we will refer to as Jacobian strain. Jacobian strain is computed directly from the Jacobian matrix obtained from the B-spline non-rigid deformation field ( Rueckert et al., 1999 ) at each point of interest, specifically at the vertices on the medial surface for each patient through time. Instead of computing strain from the deformation gradient tensor derived from a least-squares computation as in Eq. (7) , the Jacobian strain is computed from where ˜ J i,t is the Jacobian of the B-spline deformation field ψ at time t at medial surface point i . Spatial normalisation as described in Section 3.2 is then applied to each subject n , producing local Jacobian strains ˜ e atlas i,t,n . We demonstrate how Jacobian strain compares with our proposed approach at different spatial scales using (1) strain-time curves at two points in the LV wall for an individual patient and (2) computing the root mean squared difference and mean average difference between peak systolic strain magnitude of the two methods across the entire cohort, at every medial surface point.
As a further comparison to our proposed multi-scale strain computation, we also compute a multi-scale representation of strain using a moving average (MA) filter on the Jacobian strains over the medial surface. The filter kernels used have the same spatial extent as the neighbourhoods at each scale described in Section 3.1.3 . The points considered within a given kernel at each scale include a total of A s neighbouring medial surface vertices, denoted by j , for the vertex under consideration, i , This approach is more in-line with the Gaussian smoothing approach used to represent multi-scale features in computer vision ( Lindeberg, 1994 ), unlike the proposed least-squares approach which computes the relative deformation at a point compared to its neighbours at different scales. Dimensionality reduction as described in Section 3.3 is applied to the MA strains followed by PCA and LDA (or linear regression) for comparison to the proposed approach. Additionally, we plot the multi-scale MA strain-time curves at two points in the LV wall for the same patient as above as a comparison.
Mean accuracy, sensitivity and specificity, as well as their standard deviations (SD), were calculated for different classification tasks across different spatial scales V s with different numbers of PCA modes D j . Given the relatively small cohort, a typical crossvalidation method is to use a leave-one-out analysis. However, in order to quantify SDs and produce smoother grids of metrics (accuracy, sensitivity and specificity) across D j and V s , a repeated, stratified cross-validation (RSCV) was performed. This involved dividing X s into training X train s and test X test s data by randomly sampling from X s without replacement, while ensuring balanced classes in both training and test datasets ('stratified'). PCA was then performed on the training data X train s , into which space the test data X test s were subsequently embedded. Linear discriminant analysis (LDA) was trained on the embedded training data and used to classify patients from the embedded test data at different values of V s and D j . The RSCV was repeated at each D j and V s to obtain accuracy, sensitivity and specificity values for each repetition of the RSCV for each classification task. Similarly, a linear regression was fitted to the training data, and used to estimate the continuous variable QRS d of subjects in the test data in each repetition. In order to prevent over-fitting, the number of PCA dimensions, D j , was limited to 5, following the rule-of-thumb that the number of fitting parameters should not exceed more than ≈ 10% of the number of observations in the dataset. The average percentage of variance explained by the first 5 PCA modes in the training data across all RSCV folds were [59%, 57%, 55%, 55%, 61%, 70%] at V s = [1% , 2% , 4% , 8% , 16% , 32%] , respectively. Approximately 15 modes were required to explain > 90% of variance in the strain data, but the inclusion of more PCA modes (up to 10 modes) in the RSCV analysis did not improve results. Additionally, using a linear classifier (LDA) minimises the number of fitting parameters, and makes the results easier to interpret.
For each clinical biomarker, a RSCV with 100 repetitions was performed for every combination of (a) spatial scales, V s ∈ [1%, 2%, 4%, 8%, 16%, 32%], and (b) number of PCA dimensions, D j ∈ [1, 2, 3, 4, 5]. For the classification tasks, the results are reported in a grid across V s and D j for each of the metrics: the mean accuracy, sensitivity, specificity and their respective SDs over the 100 repetitions. For the regression task, results are reported in a grid across V s and D j for the metric of mean absolute error (MAE) and its SD. Note that since MAE is computed at every repetition from the test set (similarly for the classification task metrics), over 100 repetitions of the RSCV, the mean MAE is reported.
For the RSCV analysis, we found that a split of our cohort into train/test sets with a ratio of 7:1 (38 train samples and 5 test samples) produced the best results for all of the above tasks in terms of accuracy, sensitivity, specificity, and their respective SDs. This was due to (a) retaining a large proportion of the patients for the train data set (important for a relatively small cohort) and (b) ensuring a sufficient number of subjects from each class for each clinical biomarker were retained in the test sets (given the limited cohort size and imbalanced classes, shown in Table 1 ). In separate experiments, metrics generally decreased and SDs increased when a more balanced train/test split was used. This is to be expected as the number of the train samples is reduced, and the test samples contain more variance. Furthermore, increasing the train/test ratio above 7:1 meant that the class imbalance for certain tasks could not be adequately represented in the test data. Experiments also indicated that metrics stabilised after approximately 100 repetitions of the RSCV analysis for all experiments, so 100 repetitions were used for all reported results.

Predicting CRT response
The prediction of CRT response was performed using only strains across the first 70% of the cardiac cycle (at multiple spatial scales) derived from pre-treatment data. The 30 responders ( R ) were given a class label of 1 and the 13 non-responders ( NR ) a label of 0. By correctly identifying non-responders, we could avoid unnecessarily treating patients who do not benefit from CRT.

Identifying septal flash
Septal flash is a well-documented and visually identifiable motion abnormality . We chose to identify this motion abnormality as a test for the sensitivity of the classifier to strains computed at different scales. The 25 patients with septal flash were given a class label of 1, and the other 18 a class label of 0.

Identifying ischaemic aetiology
Patients treated with CRT who do not have an ischaemic aetiology have a more favourable outcome . Pa- tients with a ≥ 15% ischaemic burden will typically exhibit a significant cardiac motion abnormality. However, the abnormality depends largely on scar location and extent. Patients with an ischaemic aetiology were given a class label of 1, and 0 otherwise.

Estimating QRS duration
QRS d is a clinical determinant of cardiac health, which is used to identify LBBB and is a criterion for CRT patient selection as discussed in Section 2 . Intuitively, the contraction pattern of the LV is tied to the underlying electrical activation, for which QRS d is an indicator. In routine clinical practice, QRS d is determined from ECG readings. The estimation of QRS d from cardiac deformation in this work is included to demonstrate a correspondence between a key biomarker of electrical activation and the cardiac deformation which it influences. We seek to estimate QRS d based on LV strains at different spatial scales.  (bottom) strains are computed at a point in the septum (left) and in the lateral wall (right). Longitudinal and circumferential strains both decrease in magnitude at larger spatial scales in the septum, while longitudinal strain magnitude in the lateral wall decreases by a smaller margin at larger scales, and circumferential strain in the lateral wall changes very little. The radial strain computed in the septum and lateral wall switches sign at higher scales.

Combining multiple scales
As an extension to performing an RSCV at each spatial scale V s separately, we combined strains from different spatial scales to try to improve classification/regression accuracy. The intuition is that there may be features across different scales which together better characterise subjects for each task. We produced new matrices of strain values across all possible combinations of two or more scales, giving a total of 6 k =2 6 C k = 57 combinations ( 6 C 2 = 15 with 2 scales, 6 C 3 = 20 with 3 scales, 6 C 4 = 15 with 4 scales, 6 C 5 = 6 with 5 scales and 6 C 6 = 1 with all 6 scales). For each combination of scales, the strain matrices were simply concatenated, for example with the first two scales (corresponding to V s = 1% and V s = 2% ): All combined-scales strain data were analysed with a RSCV using the approach presented earlier in Section 4.2 for dimensionality reduction with PCA and classification with LDA (or regression), but considering only different numbers of dimensions D j ∈ [1, 2, 3, 4, 5]. The scale combinations which produced the highest accuracy for each task are reported.
In order to compare the performance of the best combined scales and the best single scale for each biomarker, statistical tests were performed on the accuracies obtained from the same 100 folds of data within each biomarker test. A two-tailed paired t -test was used to determine significant differences between normally distributed samples, and the Wilcoxon signed ranked test was used for non-normally distributed samples. A value of p < 0.01 was used to identify statistically significant differences, and normality was tested using the Shapiro-Wilk test.

Comparison of strain methodology
Fig. 5 compares Jacobian strains with strains computed using the proposed least-squares approach at V s = 1% . Strain curves are compared at two points on the LV of a patient with septal flash, demonstrating a close match of the two methods. Note how the Table 2 Results for the classification tasks (above) and regression task (below) for the clinical biomarkers of the CRT patient cohort. The best RSCV results are shown for single spatial scales of the proposed method as well as combined scales. Additionally, the best accuracies (and mean MAE) with a single scale and with combined scales are shown for the MA method beneath the dashed line for each task. The peak accuracy for each classification task is shown in bold in the second column, with the corresponding values of V s and D j used to achieve them. Similarly, the lowest mean MAE values are shown for the regression task with corresponding values of V s and D j . Statistically significant ( p < 0.01) differences between the best accuracies obtained from single versus multiple scales for the proposed method are indicated by ( * ), namely for CRT response and QRS d .  . 7. A comparison of mean strain magnitude at peak systole, across all vertices within the FOV for all patients, between strain at different spatial scales (dots) and Jacobian strain (dashed line). The mean absolute difference (blue line × ) and root mean squared difference (green line +) increases at higher scales. As we would expect, Jacobian strain is most similar to strain at the lowest scale V s = 1% . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) radial, longitudinal and circumferential strains are opposite in sign between the strain curves at the two points. Septal flash is characterised by a brief, early contraction of the septum, where during systole the septal wall does not contract. Thus a thinning of the septal wall (negative radial strain), and stretching in the circumferential and longitudinal directions (positive strains) are observed during systole, as the rest of the LV contracts and stretches the septal wall.

Description Class
To illustrate the differences between strains computed at different scales, Fig. 6 shows the multi-scale radial, longitudinal and circumferential strain curves computed in the septum and lateral wall for the same patient shown in Fig. 5 . The longitudinal and circumferential strains in the septum decrease significantly at larger spatial scales, whereas changes are less pronounced in the lateral wall. This difference is largely due to the difference in size of the septal region affected by the motion abnormality compared to the lateral wall, where the latter is considerably larger and hence longitudinal and circumferential strain are more constant across scales. The radial strain conversely exhibits a significant change in both the septum and the lateral wall at larger scales, switching signs in both cases. This effect and its implications are discussed further in the Discussion and Conclusions.
The mean absolute difference and root mean squared difference in peak-systolic strain magnitudes (i.e. the norm of the three local directions) between the Jacobian strain and strain at each spatial scale was computed across the 43 patient cohort and is displayed in Fig. 7 . The mean strain magnitude across all vertices for all patients is also shown, indicating that the Jacobian strain compares well with strain computed at V s = 1% , while strain magnitudes generally decrease with larger spatial scales. Given that the Jacobian is computed locally from the smooth B-spline interpolation of the deformation field from surrounding control points, it is unsurprising that strains computed at the 1% scale would most closely approximate the Jacobian strains.
In Fig. 8 , the strain-time curves of the MA method are presented as a comparison to the proposed approach shown in Fig. 6 . Longitudinal and circumferential strains follow a similar trend of decreasing magnitude at higher scales to the proposed approach in Fig. 6 . However, longitudinal strain in the lateral wall appears to decrease by more with the moving-average method. The greatest difference is that with the MA method radial strains have a monotonic decrease in magnitude at higher scales without any significant flip in sign, in contrast to Fig. 6 . Table 2 summarises the best results obtained with the RSCV for each classification task and the regression task along with the best spatial scale V s and number of dimensions D j used to obtain these results. Statistically significant differences between combined versus single scale results for the proposed method are also indicated. The MA method performs worse than the proposed method on all tasks and, for the purpose of conciseness, is left out for the remainder of the Results section, but addressed further in the Dis-cussion and Conclusions. Furthermore, the best results of combining scales for the MA method did not produce any significant differences to the peak single scale results for each task. Note that the label '( many )' for the combined scales result of the aetiology identification task using the MA method indicates that many different combinations of scales yielded the same result. The following sections summarise the results for each of the experiments from Section 4.2 . Fig. 9 shows annotated results from the prediction of CRT response across spatial scales and number of PCA dimensions. Peak accuracy at a single scale of 83.2% was achieved at V s = 8% using D j = 2 , with sensitivity and specificity of 92.8% and 64.0%, respectively. A peak sensitivity of 98%, which relates to the correct identification of responders, is achieved at V s = 1% using D j = 1 . A peak specificity of 64.5%, which relates to identification of nonresponders, is at V s = 8% using D j = 3 . Peak accuracy achieved with strains from combined spatial scales was found with scales [4,8,32]% and D j = 1 , with accuracy, sensitivity and specificity of 86.7%, 96.5% and 67.0%, respectively, outperforming peak accuracy from any single scale.  Fig. 10 shows annotated results from the identification of septal flash across spatial scales and number of PCA dimensions. The peak accuracy at a single scale of 94.0% was achieved with V s = 32% and D j = 1 , also producing the peak sensitivity and specificity of 93.0% and 95.0%, respectively. Best results achieved with a combination of strains from different scales was with scales [8, 32]% and D j = 1 with an accuracy, sensitivity and specificity of 92.2%, 93.3% and 91.0%, respectively. There was no statistical significance between peak accuracies from single versus combined scales. Fig. 11 shows annotated results from the identification of patients with an ischaemic aetiology across spatial scales and number of PCA dimensions. For the single scales, the highest accuracy of 70.5% was achieved at V s = 16% and D j = 1 . A peak sensitivity of 48%, related to the identification of patients with ischaemia, was achieved at V s = 32% and D j = 2 , while a peak specificity of 87.8% was achieved at V s = 4% using D j = 1 . The peak accuracy obtained from the RSCV on strains from combined scales was 71.0% with scales [8,16]% at D j = 1 , with a sensitivity and specificity of 50.0% and 81.5%, respectively. There was no statistical significance between peak accuracies from single versus combined scales. Fig. 12 shows annotated results from the estimation of QRS d across spatial scales and number of PCA dimensions. For single scales, the lowest mean MAE of 12.4ms was achieved at V s = 32% using D j = 2 . The best result using combined scales was achieved with [4, 16, 32]% and D j = 2 , producing the lowest mean MAE of 12.2 ms , which produced a statistically significant difference to the best single scale MAE.

Discussion and conclusions
We have presented a novel method for computing strain at multiple spatial scales, together with a framework for determining the optimal scale(s) with which to estimate different clinical biomarkers. We have demonstrated the utility of our method in predicting the outcome of CRT, in identifying a known motion abnormality (septal flash), in identifying a category of patients with varied motion abnormalities (ischaemic aetiology), and in estimating a continuous variable for patients ( QRS d ) via a RSCV regression.
The approach taken for all classification/regression tasks at separate spatial scales involved a subset of the first 5 PCA modes of strain at that scale. We see different patterns across the accuracy, sensitivity and specificity grids of V s and D j for each task, indicating that information from different scales and different PCA modes are being used for each classification. This is an interesting finding and shows that different clinical biomarkers are best encoded by different PCA modes and spatial scales.
For the prediction of CRT response, we demonstrate with our method that the highest sensitivity (identifying responders) of 98% is achieved with just the first PCA mode at a spatial scale of  [8,32]% denoted by C . Classification of subjects with septal flash is best at the larger scales, requiring only the first PCA mode at 32% to produce the peak classification performance. The yellow box shows the peak accuracy, and arrows indicate the general direction of increasing accuracy across V s and D j . A similar trend is seen for sensitivity and specificity. Bottom: The respective SDs for each RSCV, with a trend of decreasing SD following an increase in accuracy.. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 1%, with an accuracy and specificity of 80.2% and 44.5%, respectively. These results are comparable to using displacement, velocity, and strain kernels in an MKL framework, achieving an accuracy, sensitivity and specificity of 88.2%, 100% and 50%, respectively . The strain computed in that study used the Jacobian method to which our method was compared in Section 4 , and found to be closest to strains computed at the 1% spatial scale. However, we demonstrated here that an improved accuracy and specificity could be achieved (at the expense of sensitivity) at a higher spatial scale of 8% with 2 PCA modes, achieving an accuracy, sensitivity and specificity of 83.2%, 92.8% and 64.0%, respectively. Furthermore, our results showed that by combining strains at different scales, a peak accuracy of 86.7%, sensitivity of 97.0% and specificity of 67.0% were achieved. As demonstrated by Peressutti et al. (2017) , by combining multiple motion kernels with non-motion kernels, the specificity was improved to 62.5% while sensitivity was maintained at 100%. This suggests that by combining the strain features obtained across scales with other motion and/or non-motion features, a higher prediction accuracy could potentially be achieved.
For the identification of septal flash, we observe that the highest values of accuracy, sensitivity and specificity of 94%, 93% and 95%, respectively, are achieved at the largest scale 32% with only the first PCA mode. This suggests that PC1 at V s = 32% very strongly encodes the variation due to the motion abnormality across the cohort, and is less strongly encoded at smaller scales.
We also observe that combining strains from other scales worsens the accuracy, which suggests that the inclusion of additional scales does not provide any additional useful information for identifying septal flash, unlike for identifying CRT responders and, to a lesser extent, patients with an ischaemic aetiology. This is not surprising, as septal flash manifests in a consistent way with respect to its timing in the cardiac cycle and its anatomical location . The same cannot be said for motion abnormalities related to ischaemic scar or evidently to CRT response.
The more challenging class categories to identify -namely (i) non-responders in the CRT response prediction (specificity), and (ii) subjects with ischaemic aetiology (sensitivity) -were best identified at single scales with more components from the PCA embedding (3 and 2 respectively), and suffer from high SDs. Conversely, identifying the other class categories -(iii) CRT responders, (iv) with/without septal flash and (v) non-ischaemic aetiology -mostly performed best with just 1 PCA mode, and had relatively low SDs (see Table 2 ). The difficulty in correctly classifying (i) CRT non-responders and (ii) ischaemic aetiology is likely due to the diversity of underlying deformation patterns associated with those patients. Ischaemia, for example, can manifest in different regions of the heart with different extents, and may produce more subtle motion abnormalities depending on its severity and location. Therefore, it is unsurprising that these variations are not strongly encoded by the first few PCA modes in this analysis, especially given that ischaemic hearts are being com- pared to other diseased hearts which have other inherent motion abnormalities. Previous studies have demonstrated that diagnosis of infarct ( Suinesiaputra et al., 2017 ) and infarct localisation  can be achieved with high accuracy when comparing ischaemic patients to a healthy cohort. In future work, our framework could be extended to include a control group of healthy subjects for the improved identification of patients with ischaemic scar.
The regression task of estimating QRS d (with μ ± SD = 147 . 8 ±21.6ms across the cohort) produced the lowest error in terms of mean MAE of 12.4 ± 4.4ms at a single spatial scale of 32% with the first 2 PCA modes. A mean MAE of 12.4ms suggests an average error of < 10% for QRS d estimation across the cohort, where the mean QRS d is 147.8ms. The lowest value of mean MAE of 12.2ms was achieved by combining scales [4,32]% at D j = 3 and also [4,16,32]% at D j = 2 . While QRS d is routinely measured with ECG and does not require an analysis of cardiac deformation, it is interesting to see that QRS d can be determined best at a single scale of 32%, and then slightly better when including information from other scales.
Another finding which points to more complex deformation abnormalities being present in non-responder patients is the improvement in classification accuracy when combining strains from different scales. As mentioned above, CRT response prediction accuracy improves from 83.2% at V s = 8% and D j = 2 to 86.7% with scales [4,8,32]% and D j = 1 . Identification of scar sees a more marginal improvement in accuracy from 70.5% at V s = 16% and D j = 1 to 71.0% with scales [8, 16]% and D j = 1 . Similarly, QRS d estimation produced slightly lower errors when combining scales [4, 16, 32]%. Conversely, septal flash identification worsened when combining scales. These findings suggest that there are variations in cardiac deformation at different scales which more strongly encode differences between patients in the CRT response prediction, while septal flash is most detectable at a single scale. We observe that while QRS d estimation and ischaemic aetiology identification show only small improvements in peak performance using combined scales compared to single scales, this also means that combining strains across scales does not diminish performance (as with septal flash identification), pointing to potential limitations including the PCA embedding and linear classification/regression techniques.
One limitation of our analysis is that PCA does not take into account the relative importance of the strains at different scales, since the eigendecomposition of the covariance matrix is agnostic to the classification/regression tasks which we seek to perform using the PCA modes (unlike, for example, partial least-squares (PLS) analysis which would produce loadings on the data specific to the tasks). Despite this, classification/regression using the PCA modes from concatenated multi-scale strains out-performs classification/regression from any single scale strain (except for detecting septal flash), demonstrating a utility for combining strain across scales in the PCA analysis. Applying PLS in the future may allow for the strains from different scales to be weighted more appropriately, resulting in better performance on the different tasks. However, PCA also provides generic modes which are used in all classification/regression tasks for a given scale (or combination of scales), allowing observation of which (combinations of) features are useful for one task compared to another.
An observation from Section 4.1 was that the sign of the radial strain curves flips at the larger scales ( 8 − 32% ) for patients with septal flash, both in the septum and lateral wall. This can be explained as a function of the point-cloud neighbourhood extent in the radial direction, and the dyssynchrony of the contraction. Referring to Fig. 4 , at the 1% scale the extent of the neighbourhood in the radial direction is effectively constrained between the endocardial and epicardial boundaries. Thus for strain at a point in the septum at a small spatial scale (e.g. 1%), the radial strain is representative of a local myocardial thickening/thinning, where there is a negative radial strain (i.e. a thinning of the wall) during systole in the septum for patients with septal flash. This is because during systole the septum is passively stretched by the surrounding contracting myocardium while also being pushed outward by an elevated LV cavity pressure. Conversely in the lateral wall, radial strain at the 1% spatial scale is positive, indicating a thickening of the myocardium. However, at the larger spatial scales, neighbourhoods extend further circumferentially around the LV, and since myocardium further from the septum is moving away from the septum in the radial direction, this produces a positive radial strain. While this sort of strain no longer corresponds exactly with a local myocardial thickening, it provides a novel measure of dyssynchrony, which evidently better delineates septal flash patients from others at the larger scales.
Classification and regression tasks performed using the MA strain data produced consistently worse results than those produced by the proposed method. As demonstrated in Fig. 8 , the radial, longitudinal and circumferential strains generally decrease at higher scales. This is a natural result of the averaging process, which effectively smooths infinitesimal strain values over different kernel sizes. This is fundamentally different to the proposed approach which considers the relative deformation between a point and its neighbourhood, where the neighbourhood size is changed. Not only are the methods different in this way, but the MA method is susceptible to including more information from smaller scales at larger scales. Furthermore, using a Gaussian kernel as performed in image analysis ( Lindeberg, 1994 ) would only amplify this effect. This susceptibility produced (1) less pronounced patterns in the performance metric grids for the classification/regression tasks at different scales for the MA method (not shown for conciseness); (2) generally required more PCA modes to achieve peak accuracy at a given scale (see Table 2 ); (3) little difference between the best single scale versus combined scales results (see Table 2 ). Given this difference, the scales are not entirely comparable between the two methods, and as such we see a mismatch in the scale ( V s ) in Table 2 used to achieve peak accuracies with the two methods. The proposed method is able to generate more distinct scale-specific deformation features which contribute to better performance on the biomarker classification/regression tasks.
A limitation of our method is that strains are computed from neighbourhoods determined by nearest-neighbour distance, producing neighbourhood regions as shown in Fig. 4 . While computing strain at multiple scales with this approach provides novel deformation features which have allowed for different clinical biomarkers to be better identified in our framework, a future challenge would be to select or learn neighbourhood shapes as well as sizes which best reflect the underlying disease which we seek to identify. For example, Bai et al. (2016) proposed a method for parcellating the LV based on motion fields into 17 segments which showed significant differences for subjects of different age and gender. This method could potentially be extended to find LV regions of different scales associated with disease.
Another limitation is the relatively small cohort size, which is restricted by the difficulties of acquiring high-fidelity CRT patient data at a single site. For example, T-MR image quality can be negatively affected by patient movement in the scanner, and such cases had to be excluded from the analysis. Sample size is a common limitation for single-centre CRT patient studies Jackson et al., 2014;Peressutti et al., 2017;Ogano et al., 2014 ), and the research community would benefit from large open databases of CRT data, as is available for other cardiac pathologies through the Cardiac Atlas Project ( Fonseca et al., 2011 ). With a larger cohort, more advanced machine learning techniques could potentially be employed for patient classification and motion abnormality characterisation. We note that performing PCA on combined scales data still improved performance for identifying CRT responders and estimating QRS d , but generated little improvement for ischaemia identification. Experiments performed using a nonlinear Laplacian Eigenmap embedding and higher-order polynomial and RBF SVMs provided no improvement in results, and instead produced worse results in many cases, likely due to over-fitting given the limited number of samples in the training data. Furthermore, results obtained with the PCA embedding were more intuitive to explain. With a larger cohort, it is possible that manifold learning embeddings or neural network approaches for example could potentially represent the deformation patterns of the con-tracting LV across different scales as more compact non-linear features than PCA, and non-linear support vector classifiers could potentially more effectively utilise these features.
Finally, our method could be extended to other cardiac diseases to identify scales at which strain most strongly encodes clinical biomarkers, and in turn aim to improve clinical diagnosis.