Vortical features for myocardial rotation assessment in hypertrophic cardiomyopathy using cardiac tagged magnetic resonance

HighlightsA novel curl‐based rotation measurement built from tensorial magnitudes is proposed.Proposed rotation descriptor makes no assumption about the cardiac topology.Locally increased vorticity values are present in hypertrophied myocardial segments.Extracted vortical features have proven useful in cardiomyopathy discrimination. Graphical abstract Figure. No caption available. ABSTRACT Left ventricular rotational motion is a feature of normal and diseased cardiac function. However, classical torsion and twist measures rely on the definition of a rotational axis which may not exist. This paper reviews global and local rotation descriptors of myocardial motion and introduces new curl‐based (vortical) features built from tensorial magnitudes, intended to provide better comprehension about fibrotic tissue characteristics mechanical properties. Fifty‐six cardiomyopathy patients and twenty‐two healthy volunteers have been studied using tagged magnetic resonance by means of harmonic phase analysis. Rotation descriptors are built, with no assumption about a regular geometrical model, from different approaches. The extracted vortical features have been tested by means of a sequential cardiomyopathy classification procedure; they have proven useful for the regional characterization of the left ventricular function by showing great separability not only between pathologic and healthy patients but also, and specifically, between heterogeneous phenotypes within cardiomyopathies.

Abstract Left ventricular rotational motion is a feature of normal and diseased cardiac function. However, classical torsion and twist measures rely on the definition of a rotational axis which may not exist. This paper reviews global and local rotation descriptors of myocardial motion and introduces new curlbased (vortical) features built from tensorial magnitudes, intended to provide better comprehension about fibrotic tissue characteristics mechanical properties. Fifty-six cardiomyopathy patients and twenty-two healthy volunteers have been studied using tagged magnetic resonance by means of harmonic phase analysis. Rotation descriptors are built, with no assumption about a regular geometrical model, from different approaches. The extracted vortical features have been tested by means of a sequential cardiomyopathy classifi-

Introduction
Hypertrophic cardiomyopathy (HCM) ) is a relatively common heart muscle disease with a heterogeneous phenotypic expression that occasionally overlaps with other pathologies that also present left ventricular hypertrophy. Differentiating the underlying etiology of the ventricular hypertrophy is a frequent clinical problem with relevant implications since each etiology needs a specific management and presents a different prognosis. HCM is characterized by a hypertrophied and nondilated left ventricle (LV) (Baron, 2008), often with an asymmetrical wall thickness distribution. HCM occurs in the presence of myocyte hypertrophy and interstitial and replacement fibrosis, which cause the walls of the ventricles to thicken (Maron et al., 1992) and a reduction on the cavity volume is usually observed. This thickening may block blood flow out of the ventricle. Therefore, the main features of a HCM heart summarize in increased LV mass and thickened walls, especially in the interventricular septum (Urbano-Moral et al., 2014). These abnormalities lead to altered forces revealing a significant reduction in the diagonal components of the strain (Saltijeral et al., 2010). Previous studies have shown that regional LV dysfunctions predate over the morphologic changes related with the phenotypic expression of hypertrophy and obstruction (Dhillon et al., 2014).
As previously stated, etiological factors are of great importance in the cardiovascular disease detection (Maron et al., 2006). Global indices, such as the global longitudinal strain (Shimon et al., 2000), have been employed for cardiovascular disease identification, reporting noticeable prognostic value; however, local measurements could provide more insight to the behavior of fibrotic tissue (Piella et al., 2010). In this direction, it has been hypothesized that the presence of greater myocardial twist may be associated with a greater degree of myocardial fibrosis in HCM patients. Consequently, assessment of LV rotation mechanics as a characteristic of cardiac function may 4 A C C E P T E D M A N U S C R I P T help differentiate the presence of fibrosis (Young and Cowan, 2012). Consistently with these studies, we adhere to the appropriateness of local analyses and their clinical value on the basis that most heart diseases typically affect localized regions of the myocardium. In addition, local studies can be used to improve cardiac analytics, which may help predict the effects of specific cardiovascular diseases on the tissue.
Rotation parameters have recently gained increasing attention due to their simplicity and ease of quantification; they constitute interesting measures of cardiac performance which provide additional information on myocardial mechanics as a complement of standard pump function indices (Rüssel et al., 2009a). However, most of the rotation parameters described in the literature implicitly require an accurate description of an axis of rotation. The center of mass given by myocardial boundaries is widely used as such; however, the heart can translate during the cardiac cycle, which commonly results in misalignments of the center along subsequent frames, incurring in estimation errors. Hence, non biased calculation methods, which compensate centroid motion, are mandatory for the use of LV torsion as a measure of myocardial dysfunction quantification (Sengupta et al., 2008). Still, additional drawbacks have been reported; Young et al. (1994) state that, for HCM, the axis of rotation is shifted from the LV center of mass towards the inferoseptal region. In addition, for HCM patients, due to their characteristic asymmetrical wall thickness distribution, accurate centroid estimation could become a extremely challenging task.
Imaging techniques provide essential information for the study of these pathologies; several modalities have been proposed in an effort to measure advanced cardiac mechanics in the LV: speckle tracking echocardiography (Helle-Valle et al., 2005;Bansala and Kasliwalb, 2013), Cine Displacement Encoded (DENSE) Magnetic Resonance Imaging (MRI) (Zhong et al., 2010) or traditional cine Steady State Free Precession (SSFP) MRI, combined with feature tracking techniques, (Heermann et al., 2014), to mention a few. Nevertheless, myocardial tissue tagging with cardiovascular magnetic resonance is currently the gold standard for assessing regional myocardial function (Shehata et al., 2009). If it is not widely used in the daily practice is because it is time consuming, but to date is an accurate method to measure regional contractility (Jeung et al., 2012). MR-Tagging (Ibrahim, 2011) is usually performed by spatial magnetization modulation (SPAMM) (Axel and Dougherty, 1989) or a variant of this technique. SPAMM is grounded on the ability of altering the magnetization of the tissue (within the limitations of 5 A C C E P T E D M A N U S C R I P T relaxation times in MR) even in the presence of motion. The tagging procedure is based on the superposition of a spatial modulation over the applied gradients which may be subsequently tracked throughout the cardiac cycle, from which the cardiac function can be assessed.
Harmonic Phase (HARP) based methods (Osman et al., 2000) are widely used as a motion estimation technique in MR-Tagging (MR-T). These methods are capable of reconstructing displacement fields accurately, grounded on the assumption of constant local phase, which turns out to be more reliable than the constant pixel brightness assumption. This approach is based on the use of SPAMM tag patterns, which modulate the underlying image, producing a set of spectral peaks in the Fourier domain. Each of these spectral peaks carry information about a particular component of tissue motion, and this information can be extracted using phase demodulation methods, obtaining tensorial descriptors of deformation and, for our case, rotation estimations.
Curl is a differential operator that describes the infinitesimal rotation of a vector field. Its direction determines the axis of rotation while its magnitude shows the amount of rotation. The term vortex is commonly associated to a localized increased value on the magnitude of the given curl vector (this property will be hereafter referred to as vorticity). The local rotation measured by the curl operator should not be confused with the bulk angular velocity vector observed within the myocardial tissue with respect to a fixed cardiac axis.
Numerous 4D phase-contrast MRI (Köhler et al., 2013) studies have made use of the curl operator. Flow vortical patterns in the heart chambers, the aorta, the carotid sinus and pulmonary circulation are physiological, but can also be related to certain pathologies including aortic aneurysms, pulmonary hypertension and congenital heart defects. Vortical patterns often occur because of morphological alterations, vessel widenings or after stenosis (von Spiczak et al., 2015). These structures may alter the pressure and shear forces on the walls and trigger processes leading to cell death.
It is our understanding that curl can also quantify the local rotation within the muscle. Consequently, in this paper we introduce a novel local rotation descriptor based on robust tensorial measurements that relates the presence of increased vorticity values with the hypertrophic tissue in the heart. Rotation is estimated without influence of global myocardial parameters, such as axis of rotation or cavity radius, allowing a regional comparative study in patients with LV hypertrophy of different etiologies; HCM and Secondary forms of LV Hypertrophy (SLVH), as well as healthy subjects. To 6 A C C E P T E D M A N U S C R I P T the best of our knowledge, this is the first study that relates local vortices in myocardial tissue with the presence of fibrosis.

Materials
For the validation of the proposed approach, our study is a retrospective analysis based on a database of patients who underwent the ordinary clinical protocol according to their symptoms; the database consisted in 78 individuals who were affected by either primary HCM or SLVH (hypertensive heart disease, aortic stenosis or athlete heart disease) or were healthy volunteers. The number of pathologic patients was 56; 39 of them, with ages from 30 to 86, were diagnosed as primary HCM. These patients showed hypertrophy, predominantly in the septal region of the LV. Following the same protocol, 17 patients were diagnosed of SLVH according to chronic pressure overload. The differential diagnosis between primary HCM and SLVH was based on previous echocardiopraphic studies and clinical and familial records. About the healthy volunteers, 22 were included in the study with ages between 16 and 84; these subjects underwent the MRI protocol because of a previous suspicion of cardiac pathology but all of them turned out to be healthy.
All subjects signed the ordinary informed consent for the MR session and agreed in writing to share the resulting images for research purposes. Personal data were treated according to current legislation. Demographic data of both controls and cases, the latter indexed by pathology type, are given in Table 1. A C C E P T E D M A N U S C R I P T We have acquired short axis (SA) and long axis (LA) MR-T datasets for each patient, from apex to base, using a MR Complementary SPAMM (C-SPAMM) SENSitivity Encoding (SENSE) Turbo Field Echo sequence on a Philips Achieva 3T scanner. Regarding the tagging parameters, we validate the method for a fixed tag spacing of k i = 1/λ, with λ = 7 mm using two different orientations U i = (cos(θ i ); sin(θ i )) with θ = [π/4; 3π/4] for the stripe directions.
Additionally, we have also acquired a balanced SSFP SA MR-Cine (MR-C) sequence at the same spatial location for each patient; snapshots of the acquired sequences are shown in Fig. 1. The myocardium has been segmented in the end-diastolic (ED) phase of the MR-C sequence by two cardiologists. Cine segmentations are used to align the tagging orientations to a common reference system to correct for patient motion. The ED segmentation is used to define a region of interest (ROI) in which to compute meaningful measurements. Resolution details about these sequences are included in Table 2.
We have implemented a preprocessing pipeline in order to (a) propagate the ROI in MR-C from ED to the end-systolic (ES) phase -in which subsequent calculations will be carried out-and (b) align the MR-C and MR-T sequences at ES. These two steps are: • Registration The MR-C sequence is processed by means of a groupwise elastic registration procedure (Cordero-Grande et al., 2013a) in order to propagate the ED segmentations towards ES phase. The transformation is achieved by B-spline based Free Form Deformations (FFD) (Rueckert et al., 2006). A gradient-descent optimization scheme is performed where the sum of the squared differences of image intensities is used as registration metric. A smoothness penalty term has also been introduced to constrain the spline-based FFD transformation to be smooth.

• Alignment
An affine registration method is performed to align MR-T and MR-C images at ES phase. The MR-T sequence has been detagged by means of a homomorphic filtering procedure (Makram et al., 2015) prior to the alignment process.

Motion estimation
3D HARP motion reconstruction using the C-SPAMM technique requires a minimum of 3 linearly independent wave vectors (Osman et al., 2000). We have extended the aforementioned HARP methodology for the computation of the deformation gradient tensor using SA and LA images on the intersection of the slices as shown in Fig. 2. For points on which LA axis images were not available, 2D motion has been reconstructed.
The motion estimation technique is based on the extraction of the local phase of the grid pattern according to the method presented in Cordero-Grande et al. (2011. A windowed Fourier Transform (WFT) is applied to the image at ES phase. The WFT provides a representation of the image spectrum in the surroundings of each pixel of the original image, so HARP band pass filtering techniques can be directly applied on the spatially localized spectrum of the image. To adequately retrieve the shape of the spectral peaks, we have resorted to an anisotropic filtering approach combining Gaussian band-pass and all-pass filters as proposed in Sanz-Estébanez et al. (2016a). Finally, each of the image phase ϕ i (x) (two for each plane) can be extracted in the spatial domain from the inverse WFT of the aforementioned filtered spectrum.
As mentioned before, we have extended the HARP methodology by allowing the estimation of motion under the application of a set of four wave vectors. Therefore, 3D deformation gradient tensor can be robustly recovered at the intersection points of both axes, by applying the methodology presented in Cordero-Grande et al. (2016). The material deformation gradient tensor F(x) can be estimated from the gradient of the phase image as stated in Osman et al. (2000) as: where K represents the two stripe orientations of four given wave vectors corresponding to each tagged image. Robust estimation of F(x) is achieved through Least Absolute Deviation (LAD) procedure. The reconstruction is performed via Iteratively Reweighted Least Squares (IRLS): with W l (x) a diagonal weighting matrix, which is updated at each iteration by considering fitting residuals (Cordero-Grande et al., 2016).

A C C E P T E D M A N U S C R I P T
From this estimated tensor, the main cardiac function characteristics can be obtained through the Lagrangian strain tensor, defined as: The spatial resolution of the reconstructed tensors depends on the width of the HARP band pass filter Prince, 2003, 2004); the HARP method is upper limited by half of the tag spacing (small deformation assumption). However, WHARP methods (Sanz-Estébanez et al., 2016a;Cordero-Grande et al., 2016) try to accommodate the band pass filter to the the local frequency of the signal in order to approach to the maximum achievable resolution. Therefore, effective HARP resolution will vary dynamically, allowing large deformations, as those observed at ES, being captured at a maximal scale of 1.5 times half the tag spacing. These tensors have been calculated at ES, where the greatest deformation along the cardiac cycle takes place.

Rotation parameters
In addition to thickening and shortening, the myocardium also undergoes a wringing motion during systolic phases due to the obliquely oriented subendocardial and subepicardial myofibers. Many descriptors have been proposed to measure this motion that rely on either global information derived from simplified anatomical models or on tensorial strain and deformation magnitudes built from local motion estimates.
Measures based on global information. It is well known that the LV apex globally rotates anticlockwise at a relatively constant rate throughout systole. On the contrary, the base, initially rotating anticlockwise, reverses direction providing a net clockwise rotation at ES phase. The resulting difference of these two motions is defined as twist, defined to be positive by convention (Young and Cowan, 2012).
There is currently a lack of standardization for methods used to characterize the global LV twisting motion. These descriptors rely on geometrical models of the heart for torsion and twist calculation. Consequently, both a well-defined fixed axis of rotation and regular myocardial radii over the whole heart are mandatory. For example, torsion has been traditionally calculated as relative rotation in degrees (Lorenz et al., 2000), rotation per length in degrees/mm, torsional shear angle, also in degrees (Buchalter et al., 1990), and longitudinal-circumferential shear strain (dimensionless) (Fung, 1965). Traditional rotation indices are obtained by vectorial product between position vectors at ES u ES and ED u ED phases as: As stated above, twist computation depends on the exact locations of the apical and basal slices and requires accurate motion compensation, specially for centroid motion correction. Twist per unit length is also widespread, since torsion is relatively constant in the longitudinal direction (Young and Cowan, 2012). Nonetheless, this measure does not scale appropriately between hearts of different sizes and we have not observed a significant complementary value with respect to the twist.
The torsional shear angle is a measure of the change in angle between line segments which are initially aligned with the anatomical axes of the LV. Many studies have used the formula given by Aelen et al. (1997). However, it has been demonstrated that it usually overestimates deformation (Rüssel et al., 2009b), so we have resorted to an unbiased alternative formula based on circumferential displacements: where D is the distance between selected segments. 1 However, HCM characteristic endocardial irregularities may hinder the accurate estimation of both the myocardial radius and the axis, which are crucial in this formulation.
Tensorial descriptors. Another important group of rotation descriptors focus on the properties of the tissue that provide localized characterization of the motion by tensorial analysis. As stated in solid mechanics (Fung, 1965), the 3D strain state at any point in a body can be fully represented by three diagonal strains and three shear strains. From them, the longitudinalcircumferential shear (E lc ) is a useful measure closely related to torsion. According to this analysis, local torsion measures can be defined, i.e., the 3D local torsion shear can be given by: where E cc and E ll represent the circumferential and longitudinal strains, respectively; these components can be obtained by straightforward operations on the Cartesian components in (3). Nevertheless, shear strains are several magnitudes lower than diagonal strains, so factors other than inotropic state may greatly affect its estimation (Petitjean et al., 2005). Angular variation between two states of stress in the plane in Cartesian coordinates can be expressed by a single angle φ as stated in Fung (1965): where ε represents the Cauchy's strain tensor (ε) directly related with the stress tensor by the Lamé parameters (Fung, 1965). Particular values of φ show angular variation of stress principal directions between both states (ES and ED phases). Additionally, in Cordero-Grande et al. (2013b) a novel rotation parameter has been proposed built from the (longitudinal) transformation that suffers a local coordinate at ED through time in the LR plane as given by the material deformation gradient tensor F.
Our work ellaborates on vortical patterns widely employed for the identification of abnormal flow patterns. Nonetheless, the term vortex bears different interpretations as defined in the literature. For most flow studies performed in clinical practice, the term vortex denotes rotating motion, where stream or pathlines tend to curl back on themselves (Markl et al., 2011). In fluid dynamics, a vortex is a region in a fluid in which the flow is rotating around an axis line, which may be straight or curved. More explanatory notes on the theoretical definition of the term vortex can be found in Stalder et al. (2010). In this paper, the term vortex will be used in the solid mechanics context as a local abnormally increased rotation component. In order to find evidence of perturbations within the myocardial tissue, material vortical patterns can then be associated to increased values of an dynamic rotation parameter extracted from the deformation gradient tensor.

A C C E P T E D M A N U S C R I P T
The curl of a given deformation field u describes local spinning vectors (see Fig. 3) and can be calculated as: where F ab represents a component of the material deformation gradient tensor in Cartesian coordinates. Hereinafter, vortical parameters will be expressed in Cartesian coordinates as opposed to the cylindrical coordinates from the strain tensor used in (6).  (9) in basal and apical slices, left and right respectively. The arrows show the extracted cardiac displacement field while the colour represents the intensity of vorticity (unitless). Some outliers are observed near the boundaries due to the difficulty of HARP methods in tracking material points in the presence of great intensity changes. Scales are set to best accomodate the range of values on the given cardiac plane, since myocardial rotation varies in modulus and direction along the cardiac axis.
In this paper we hypothesize that local increasing of vorticity values (in modulus) arises within myocardial segments with fibrosis-related perturbations. Nonetheless, high vorticity values, irrespective of the fibrosis degree, are prone to appear in myocardial boundaries, giving rise to multiple outliers in rotation estimation.
These vorticity measures are insensitive to the definition of the rotation axis, although 3D deformations are needed for its proper reconstruction.

A C C E P T E D M A N U S C R I P T
When LA information is not available, only the longitudinal component (ω z ) of the vorticity vector can be estimated. In addition, if the cardiac axis is not planned properly, vorticity parameters will be estimated with a systematic error related to the angular error of the axis. However, as it is common in clinical practice, we will assume that the main axis of rotation will lie on the LA planes (i.e., will be normal to the SA image planes). Hence, the ω z component of the vorticity vector provides clinical compliance as it is aligned with the wringing motion of myocardial fibers. Thus, a phase increment due to the aforementioned local rotation can be extracted by integration: This parameter will be referred to as local rotation ϑ. Therefore, twist motion will be approximated as For comparative purposes, we will also estimate ϑ and β rotation distributions with two other different methodologies. First, we will analyse the capabilities of the elastic registration algorithm described in Section 2.2.1 applied to MR-C to detect rotation from the estimated deformation fields. Second, we have employed an atlas-based approach that consists of a spheroidal model (Young and Axel, 1992) fitted at ED and that deforms due to the forces exerted from a stripe tracking procedure (Young et al., 1995) on the MR-T sequence throughout the cardiac cycle. Deformations are formulated from continuous parameter functions which, in addition, include parameterized twisting and rotation axis deformation as given in Park et al. (1996) and, therefore, can be applied to any shape.

Classification
We have resorted to a classification method (Sanz-Estébanez et al., 2016b) to assess the discriminating ability of the rotation features previously described. The procedure, sketched in Fig. 4, consists in an automated processing pipeline to classify heterogeneous groups of ventricular hypertrophy (and controls) from myocardial functional descriptors. The proposed classification method is grounded on the idea that populations will be highly overlapped overlap strongly irrespective of the specific features selected for classification if the problem is addressed through a single stage. Our purpose

A C C E P T E D M A N U S C R I P T
is to classify a sample into one of three classes, namely, control, primary HCM and secondary hypertrophy (SLVH). Since secondary hypertrophy patients have subtle differences with respect to the other two classes, we have resorted to a two-stage classification procedure. Thus, we have divided the classification process in three stages and performed a feature selection step for each stage independently, following a sequential methodology that adapts to the characteristics of the population at every stage. Different machine learning methods (both supervised and unsupervised) have been implemented for each stage and all their possible combinations have been tested. Mechanical descriptors extracted from the aforementioned tensors are an essential part of the classification procedure. We have considered different groups of features. First, the components of the strain tensor in the cylindrical coordinate system {R, C, L} are accounted for. We also use twist and

A C C E P T E D M A N U S C R I P T
torsion features (see Table 3) built from the aforementioned rotation parameters as extracted from the MR-T and MR-C sequences, as well as using the spheroidal model we have previously referred to. Additionally, since some rotation-related components have opposite directions in apex and base, we consider the location of the zero crossing for these components as well.
For the feature extraction step, we have previously selected for each feature the most representative cardiac segments (Cerqueira, 2002) within clinical practice. For twist and torsion descriptors, we have considered septal segments, whereas for tensorial parameters, mid-ventricular or basal segments have been chosen. In Table 3 we show the segments involved in the calculation of each of the features; then the overall feature is, for robustness purposes, the median of the distribution within those segments, which will be the input to the classifier. For the twist parameters the feature extracted is the difference between medians on apical and basal septal segments. On the other hand, for the zero crossing parameter, the feature represents the height of the cardiac axis at which the given rotation parameter, estimated slice by slice, changes its sign. Notice that the feature extraction step is grounded on clinical knowledge, i.e., it is not data-driven. Table 3 Cardiac segments involved during the feature extraction stage for each one of the motion descriptors employed in the classification procedure. For each row, only one feature will be extracted, summarizing all the segments indicated below. The extracted features will be the input to the feature selection step, from which the final (survivor) feature vector (from 2 to 5 components) will arise. The number within braces indicates the equation that defines the parameter.

A C C E P T E D M A N U S C R I P T
As reflected in Fig. 4, after feature extraction we carry out a normalization stage in order to diminish the influence of possible outliers. A sigmoidal function, with its scale factor set according to Theodoridis and Koutroumbas (1999), is used to this end. Data are mapped on the interval (0, 1) by imposing a generalized logistic function; outliers will tend to appear at either of the two extremes, while maintaining a linear relation for the rest of the data. Then, normalized features are arranged in vectors with different number of components (2 through 5). All possible combinations of features indicated in Table 3 have been tested.
Feature selection and classification performance assessment has been carried out in a similar but sequential manner. For both, data samples have been randomized and a Leave-10-out method has been applied; the proportions of control/primary/secondary have been kept unaltered along trials. Specifically, for feature selection in the first classification stage, we classify the samples in controls and primaries and calculate the accuracy; the feature set with the highest figure is selected for this stage. In parallel, and for the second stage, controls and secondaries are classified on one branch (left branch) and primaries and secondaries on the other branch. In this case, features are selected with the criterion of maximizing the sensitivity to secondaries so as to avoid bias towards the groups with larger sample size, specially between HCM and SLVH. This procedure has been labelled in Fig. 4 as K-fold Feature Selection. As for finding classification performance, a similar cross validation procedure has been carried out (labelled in Fig. 4 as Leave-10-out Classification) on new randomizations. The classifiers tested have been Fuzzy c-Means (FCM) (Bezdec, 1981) and Support Vector Machines (Cortes and Vapnik, 1995) both with quadratic (SVMq) and Gaussian (SVMg) kernels (Vert et al., 2004).
In order to assess the performance of a given feature in the classification procedure we have measured its surviving percent rate. This parameter shows the membership probability of the given feature to the feature vector extracted from the feature selection step; in other words, it is the frequency that the feature is employed within any of the stages of the classification along trials.

Rotation analysis
Torsion is known to be dependent on LV shape, with reduced twist in more spherically shaped hearts and increased torsion with concentric hypertrophy due to an increased lever arm for myocardial fibers. In HCM, torsion has been reported to increase despite reduced circumferential and longitudinal shortening (Young and Cowan, 2012). These findings, together with others described in Section 1, can be observed in the results included in Table 4, where the mean and standard deviation of the twist and torsion distributions derived from the aforementioned rotation features are shown. Table 4 Twist and torsion parameters extracted from the MR-T sequence (mean ± std) for segments in Table 3 for each population. The number within braces indicates the equation that defines the parameter.

Populations HCM SLVH Control
Twist ϑ (11) 10.46±2.11 9.23±2.28 7.33±2.08 Twist β (4) 13.53±2.50 13.80±3.37 9.87±2.81 Twist φ (7 In Fig. 5 we show snapshots of mid-ventricular slices of the local rotation extracted by means of the vortical approach as described in (10) for a HCM and a SLVH case as well as for a healthy volunteer. In general, septal segments for HCM have shown higher vorticity, specially when compared to lateral segments, whereas for secondary cases this behavior can be observed in any of the cardiac segments. In healthy volunteers the extracted values are lower compared to HCM patients independently of the cardiac segment.
In Fig. 6, we show color codes of the mean ± standard deviation (respectively, inner and outer rings within each segment) of the rotation parameters over the 17-segment model (Cerqueira, 2002), estimated from both the vortical approach given in (10) and the traditional approach by (4) for A C C E P T E D M A N U S C R I P T Figure 5. Snapshots on vortical-based rotation measurement ω z from (9) for HCM, SLVH and a healthy volunteer over mid-ventricular slices (from left to right).
the resulting distribution of the deformation vector field. Additionally, and for the sake of comparison, rotation parameters extracted with the elastic registration procedure over MR-C and from the deformable model have been included as well in the second and third rows, respectively. Student t-tests 2 have been performed to highlight differences on the mean of the vorticity modulus and the aforementioned rotation distributions on each of the 17 cardiac segments. Each population of the study has been compared with the other two, separately for each rotation parameter; the numerical results are shown in Table 5 and graphically in Fig. 7 bull's eye display. It is noticeable that the vortical approach seems to show larger differences between populations compared to the traditional approach. Septal segments seem to bear higher discriminating capability, specially when twist is measured by the vortical approach. Additionally, we have performed (twoway) ANOVA tests over the ϑ distributions extracted from MR-C and MR-T sequences as well as using the deformable model (over MR-T) for each population and cardiac segment. Significant differences (p < 10 −3 ) have been found in the vast majority of the comparisons. ϑ and β distributions were not compared since the measured parameters does do not represent the same component of the physiological rotation motion.

Classification analysis
We have assessed the survival rate of the feature selection stage of the classifier (recall Fig. 4). Results are shown in Table 6 for the features de-A C C E P T E D M A N U S C R I P T Figure 6. Regional study of the aforementioned vortical rotation parameters obtained from the MR-C and MR-T sequences, as well as using the deformable model over the latter, for the different populations. Traditional rotation is also depicted in the last row. For each cardiac segment two colors are depicted, the inner showing mean + std and the outer for mean − std. Table 3; features obtained from both the MR-C elastic registration as well as from the spheroidal deformable model have also been included. Additionally, we have also included conventional indices of cardiac motion used in clinical practice, such as wall thickening (WT) over mid-ventricular slices (see Dong et al. (1994); Prasad et al. (2010) for more details), ejection fraction (EF) and LV volume, with the latter both at ED (EDLVV) and ES A C C E P T E D M A N U S C R I P T Table 5 p-values for the comparisons between distributions (H, S and C stand for HCM, SLVH cases and controls, respectively) of the given rotation parameters indexed by number of segment. When unspecified, MR-T is the image source. Significance level after Bonferroni correction is 0.017. phases (ESLVV). The most repeated configuration of the classifier consisted of FCM for stages 1 and 2.2 and SVMg for stage 2.1. The best accuracy figures were obtained when using diagonal strain tensor components on stage 1, whereas for stages 2.1 and 2.2, the best feature vectors turned out to be [E ll , Twistϑ T AG , Twistφ] and [E ll , E cc , Twistϑ T AG , || ω||], respectively. If we take into account not only the best classifier but we also rank performance and analyze the first, say, ten results, the composition of the selected feature vectors shows some degree of variability which seems very much in accordance with the results in Table 6.

scribed in
In terms of end-to-end performance, the obtained accuracy (86%) seems comparable in classification figures with other procedures (see Gopalakrishnan et al. (2014);Puyol-Antón et al. (2017)) although, in these cases, no SLVH are analyzed, so comparisons have to be made cautiously. In disaggregated terms, the sequential classifier has obtained sensitivity figures higher than 70% for each group (specifically, 81% for control, 72% for secondary hypertrophy patients and 95% for primary HCM patients). It is worth mentioning that no primary is classified as control and viceversa; therefore, the pipeline proposed seems a proper screening tool. Secondary patients performance is clearly lower as compared with both controls and primary HCM patients, possibly due to a smaller sample size as well as the subtle differences they show.

ACCEPTED MANUSCRIPT
A C C E P T E D M A N U S C R I P T Figure 7. p-value bull's-eye plots from intra-segment comparisons between primary HCM and SLVH patients for ϑ distributions obtained with different methodologies, as well as traditional rotation β distribution. Scale is defined as −log 10 (p-value).
Finally, in order to assess the relative strength of the different measures in classification performance, we have run the classification pipeline with different features subsets. In particular, in Table 7 we have compared confusion matrices obtained with the full feature set (MR-T + MR-C + Deformable model + Conventional clinical indices) and with MR-T features only. From the latter, we have also shown classification performance obtained discarding traditional and vortical rotation features, respectively.

Discussion
The relationship between myocardial fibrosis and local mechanics is important for the diagnosis and treatment of cardiomyopathies (Karamitsos and Neubauer, 2011). This paper shows that LV rotation is essential for proper myocardial function. In our case, most of the measurements shown in Table 4 indicate that LV rotation can be considered as a marker for cardiac disease identification and might be helpful for cardiomyopathy understanding, thus providing complementary information to standard pump function indices.

23
A C C E P T E D M A N U S C R I P T Table 6 Surviving percent rate of the features employed in the classification procedure (see notation in Fig. 4). Stage 2.2 is devised to SLVH-to-HCM sensitivity, while stage 2.1 refers to SLVH-to-Control sensitivity.

Features
Stage 2 The diagonal components of the strain tensor (E cc , E ll and E rr ), defined in (3), have provided the highest separability between pathologic and healthy groups (cardiomyopathy screening) as shown in Table 6; our results are in accordance with this finding (Saltijeral et al., 2010). Twist parameters seem to be also valuable in this step, showing higher survival rate than shear strains and torsion parameters. For the refinement step for controls/SLVH (stage 2.1 in Fig. 4) the E cc component in mid-ventricular areas is the most discriminative. Amongst the rotation-based parameters, vorticity modulus || ω|| and twist ϑ T AG remain the most discriminative features.
In parallel, for the classifier stage 2.2, lower figures on the global performance are obtained. Besides, the selected feature vector presented one 24 ACCEPTED MANUSCRIPT A C C E P T E D M A N U S C R I P T Table 7 Confusion matrices for classification performance with different feature subsets. Matrices have been normalized with respect to the total number of patients. Each column represents the instances in a predicted class while rows represent the instances in an actual class. more component and two curl-derived entries; it consisted of a combination of two diagonal strain parameters, the vorticity modulus (9) and the twist extracted from the vortical approach (10). Extending the analysis to more (than one) high-ranked features vectors, we observe a greater degree of heterogeneity as well as the frequent presence of ω z (9). Consequently, the proposed curl-derived parameters (9)-(11) have revealed turned out to be particularly useful for the discrimination of primary and secondary cases as reflected by sensitivity figures extracted from Table 7.
To the best of our knowledge, this is the first study in HCM patients that relates vorticity in cardiac deformation fields with local myocardial mechanics and its abnormalities; our results suggest that vorticity may help deepen on the underlying characteristics behind primary and secondary cases of LV hypertrophy. For these parameters, neither the length nor the center of the heart are needed, so no bias is introduced in their estimation; as we have described above, vorticity is directly related to the deformation gradient tensor and, consequently, it can be estimated from the same information used in the strain tensor (3) analysis. In addition, they show higher survival rates than techniques that use fixed LV axis and center representations, giving rise to a more reliable parameter.
The color-coded results shown in Fig. 6 reveal that patients with both forms of ventricular hypertrophy present greater rotation distributions as compared with controls over most segments both for vortical (10) and traditional (4) rotations. As for primary HCM patients, there is a clear increased rotation (in modulus) for all segments, but high vorticity areas are mainly located on septal segments. SLVH patients showed a somewhat different pattern in mid-ventricular and basal regions of the heart, presenting higher A C C E P T E D M A N U S C R I P T values than controls; those values are not focused on septal segments but a slight bias to lateral segments may exist. Other regional analyses have also been performed by means of automatic LV segmentation (Bai et al., 2016;Liang et al., 2015). However, the presence of hypertrophic tissue and other pathologies may bias the final parcellation of the cardiac segments. For this reason, we have made use of the 17-segment model as a consistent and well-established model for motion analysis (Smiseth et al., 2016). The usefulness of the vortical parameters has also been reflected by the improvement shown in Table 7 with respect to traditional rotation parameters and the minor degradation with respect to the full-feature option when curl-derived parameters are used in isolation.
Our results also indicate a higher performance of MR-T for motion estimation with respect to MR-C; however, it is well-known that HARP procedures have some difficulties in correctly estimating the phase in the vicinity of boundaries. In those areas, elastic registration procedures over MR-C is usually more robust. For this reason, a coordinated procedure that weighs both information sources according to position may potentially provide better figures.
Additionally, vortical measures were compared when extracted both from the HARP method and by deforming a spheroidal model previously fit. Similar vorticity values were obtained although the extracted vortical patterns from the spheroidal model showed higher spatial smoothness due to the regularized functions used to define the deformation, thereby reducing its usefulness for classification (see Table 6).
Finally, conventional global indices in HCM diagnosis (see Table 1) have also been tested in the classifier. Neither of them presented a very representative survival rate (even in stage 1), hence, their influence on final performance does not seem relevant.

Conclusions and future lines
In this paper we have related anomalies on local vortical patterns with the presence of fibrotic tissue by means of an image processing pipeline and a two-stage sequential classification method. Local rotation parameters are estimated by means of a robust motion and tensor analysis so that potential biases of global analyses are avoided.
Local rotation was significantly increased in primary HCM patients, specially in the septum, compared to controls and secondary cases; in the latter,

A C C E P T E D M A N U S C R I P T
vortical abnormalities may show a slight trend to lateral segments and with values less pronounced than primaries. These findings may provide important information in hypertrophic diseases to establish a differential diagnostic between these two classes.
Classification figures, although collateral in the paper, are promising; clearly, discrimination between primary HCM and secondary cases is more challenging than between HCM and controls. Therefore, figures related to the former problem have been lower than those related to the latter. A larger cohort may let us increase this number in the near future. Classification of SLVH cases has proven to be a challenging task but figures, despite not being remarkable, are likely to improve when equalizing the number of subjects in the study or by introducing features that take into account the position of vortical peaks.
A C C E P T E D M A N U S C R I P T