Strainger things: discrete differential geometry for transporting right ventricular deformation across meshes

. Cardiac dynamics have been a focus of image analysis, and their statistical models are used in a wide range of applications: generating synthetic datasets, derivation of specific biomarkers of pathologies, or atlas-based motion estimation. Current representations of dynamics, mainly based on displacements, often overlook the physiological basis of cardiac contraction. We propose to use local strain as a more accurate representation, and demonstrate this on 3D echocardiography surface meshes of the right ventricle. Our methodology, based on a differential geometry algorithm, deforms the surface mesh according to a pre-imposed strain field. This approach allows for a clearer disentan-glement between cardiac geometry and dynamics, better differentiating deformation changes than those due to changes in cardiac morphology. The methodology is demonstrated in two toy examples: transporting deformation from one individual to another; and simulating the effects of a pathology on a healthy patient, namely an akinetic right ventricular outflow tract.


Introduction
The heart ventricles are made up of myocytes, which are muscle cells that contract during systole to generate force, thereby raising the pressure inside the chamber and ejecting blood. The shape and function of the heart are closely linked, and any change in shape can affect how the heart works, even if each individual myocyte is still working properly [14]. This is crucial to understand because some illnesses primarily alter the myocardium shape (such as hypertrophy) and, subsequently, affect the cardiac function, while others primarily alter This a pre-print version. The final document will be available at http://www.springerlink.com the conditions under which the myocardium works (such as infarction), leading to altered shape as well as impaired function.
Finding ways to represent cardiac function that are independent to its shape is a significant area of research with various applications. These include identifying functional abnormalities in populations with diverse shape variations, creating synthetic populations of shapes with different functional patterns, and transferring function between different individuals or imaging modalities (for example, using an atlas to estimate the function of a patient for whom only computed tomography (CT) images are available).
Researchers in medical image analysis have proposed various methods for better considering independent components of shape and deformation when analyzing the cardiac function, such as spatio-temporal statistical shape models [9,1,11]. These approaches describe the pointwise displacement of every point of the mesh. More modern approaches use parallel transport to transport displacement fields among two meshes while accounting for the non-linear geometry of the mesh spaces [7,10]. However, this strategy does not consider the physiological basis of cardiac mechanics: it ignores that the resulting pointwise displacement is the aggregation of the local contraction and relaxation of individual myocytes, and thus movement at a certain point is affected by the contraction of many myocytes, which will pull the shape. Therefore, strain-based indices are often used in the clinical community as they are size-independent and better comparable between individuals, allowing for better identification of abnormalities [13]. Strain-based analysis has been performed by authors, but these representations cannot be mapped back to the original space, generating synthetic representative elements for visualisation [5].
While there are standardized strain representations for the left ventricle (LV), which is divided into segments [3], analyzing strain in the right ventricle (RV) is more challenging due to its thin walls and difficulty in accurately imaging it with magnetic resonance imaging (MRI) or ultrasound [12]. This complicates the calculation of strain as the RV is represented as a surface and requires the use of differential geometry, which is the field of mathematics that studies smooth surfaces. Additionally, the RV has higher anatomical variability and irregularity [8] compared to the regular ellipsoid shape of the LV, making standardization of segments more difficult.
In the present work, we propose a strain-based methodology to represent deformation across the cardiac cycle and analyse the cardiac function across a population. In particular, the developed strain representation can be imposed to a given surface mesh at end-diastole (ED), to generate realistic dynamics. This methodology is based on differential geometry, where we solve the inverse problem of strain computation in a local anatomical system of reference, to construct a mesh satisfying the desired strain. To allow a direct transportation and comparison of strain between different individuals, strain is expressed in a local anatomical system of coordinates, consisting of the longitudinal and circumferential coordinates. We demonstrate the utility of this representation on two toy examples: to transport deformation between different subjects, and how the strain of a healthy individual can be modified to model a pathology.

Data
To illustrate the methodology, we used patient-specific triangular surface meshes of the RV obtained from a commercial software (4D RV-Function 2.0, TomTec Imaging Systems GmbH, Germany) for segmentation and tracking along 3D echocardiographic sequences. The meshes of different individuals are in approximate point-to-point correspondence thanks to the model-based segmentation.
We tested the methodology on the data from two male controls. They were selected from a small control cohort for being the subjects with the largest global longitudinal and circumferential strain, respectively. We focused on the meshes at end-systole (ES) and ED.

Discrete differential geometry
Discrete differential geometry studies 2D surface meshes embedded in the 3D space. In discrete differential geometry, a triangular surface is characterized (modulo translations and rotations) by its first and second fundamental forms [4], as in classical differential geometry. The first fundamental form is defined at every triangle and consists of the length of each of its sides (this uniquely defines a triangle, modulo translations and rotations). The second fundamental form is defined at each edge, and corresponds to the dihedral angle, which is the angle between the normals of the incident triangles. It is worth noting that not all combinations of dihedral angles and lengths correspond to a valid mesh, but they require to satisfy integrability conditions, which are the discrete equivalent to the Codazzi-Mainardi equations (also known as Gauss-Codazzi) [15]. Informally, they correspond to ensuring that for every possible loop traversing triangles, the sum of the traversed dihedral angles must be a multiple of 2π and is that the common edge of two triangles must have the same length in both triangles.
At each cell of the mesh, we define a local system of coordinates that corresponds to the longitudinal and circumferential directions. These directions are estimated using a heat diffusion method, where the longitudinal direction is defined at each triangle as the gradient of the stationary solution to the heat equation from the apex to the pulmonary and tricuspid valves, as explained in [6]. These directions are orthogonal and defined in the tangent space of the mesh, which in a discrete mesh corresponds to the plane in which the cell is contained.

Deformation representation
Given an ED mesh and its corresponding ES mesh in point-to-point correspondence, we can compute all parameters needed to represent its ED-ES deformation. The different steps of these computations can be seen in Figure 1. Basically, for each i-th cell, we need to compute the deformation tensor A i , which is the linear mapping from the undeformed cell to the deformed cell. We express this tensor in the local anatomical system of coordinates, leading to a 2 × 2 matrix. Since we use the Lagrangian formulation, this is defined in the system of coordinates of the ED. The second part of the deformation consists of the changes of dihedral angles, which are represented for each edge, as the difference of the angles at ED minus the ones at ES.

Strain transport algorithm
Here we describe the method to impose a strain (expressed as defined in Section 2.3) from a given source individual to a new target individual for which only the ED mesh is required. This strain can come from an acquisition where both ED or ES are available (that might correspond to another individual, but also to the same individual in a different modality), but it may also come from an atlas, or be a purely synthetic strain.
A linear mesh reconstruction has been proposed [15], which given the desired edge lengths and dihedral angles, finds a mesh (with a fixed topology) that satisfies as much as possible these constraints in the L 2 manner. This approach works by optimizing the 3D coordinates of each node, and defining a dummy matrix variable F, defined at each triangle, that represents the rotation defining how each triangle is positioned in space. In the original approach, F was not forced exactly to be a rotation, which allowed efficient solver algorithms, but produced shrinking artefacts. In an extension [2], F was constrained to be a rotation by using a parametrization in the Lie Algebra of rotations SO(3), although this increased the computational cost. The steps of the algorithm to transport a strain deformation from a mesh to another are defined as follows, and a schematic view of this process is provided in Figure 2. : 1. Compute the longitudinal and circumferential directions of the ED target mesh at each cell. 2. Express each point of each triangle of the target mesh using its anatomical system of coordinates and deform it using the anatomical deformation tensor of the prescribed strain A i . Compute the length of each triangle. 3. Compute the new dihedral angles, by adding the desired variation to the ones of the target ED mesh. 4. Find the ES mesh by using the mesh reconstruction algorithm to recover a mesh satisfying, as much as possible, the prescribed strain.
Note that a full match to the prescribed strain is not possible, since the desired lengths and dihedral angles computed in steps (1)-(3) above might not exactly satisfy the Codazzi-Mainardi equations (note that they depend on both the target mesh, and the desired strain). This is because, even if we can locally deform every triangle/edge to satisfy the desired new lengths and dihedral angles, imposed strain needs to be globally consistent, as expressed in the Codazzi-Mainardi equations.This implies that the reconstruction step performed in step (4) above can be seen as a projection to the space of dihedral angles and edge lengths that exactly satisfy the equations.

Experiment design
We assessed the utility of the proposed methodology in two toy examples: 1. We first examined how to synthetically simulate the impact of pathology on strains. A specific region was selected to become akinetic and undergo no deformation. We selected the right ventricular outflow tract, which can be impacted by diseases such as pulmonary hypertension or arhythmogenic right ventricular dysplasia. The procedure was as follows: we calculated the deformation (consisting of deformation tensors and dihedral angle variations) using the original ED and ES data, and manually selected a region near the outflow tract. We then removed the kinetics from that region by setting the deformation tensors to the identity matrix and the dihedral angle variation to 0. Finally, this modified deformation was imposed to the ED mesh. 2. We tested the capacity of our method to impose a real deformation to a new geometry. We selected the mesh with the largest longitudinal strain as the target, to which we imposed the strain from the mesh with the largest circumferential strain. The evaluation metrics consisted of the residual of the longitudinal and circumferential strains.

Code availability
The Python code is publicly available as a library in the github (https://github. com/gbernardino/rvmep) for non-commercial academic usage.

Synthetic deformation
As described above, we first estimated the deformation of a real sequence, synthetically removed deformation at the outflow tract, and applied it to the original ED mesh. This can be seen in Figure 3. We can also observe that the outflow tract has low deformation, both quantitatively (Figure 3.c, blue points), and visually, since it has the bloated form characteristic of ED, instead of a more squeezed shape which would be expected at non-pathological ES. As expected, this region has high displacements, as the rest of the ventricle is contracting and pulling it.

Inter-individual transport
We also assessed the ability of our method to correctly transport the real ED-ES deformation of a given individual to a new target ED mesh, both qualitatively and quantitatively, by comparing the original target ES mesh to the one reconstructed after transporting the source deformation. This can be seen in Figure 4. We observe that the resulting mesh undergoes large circumferential motion, especially near the apical region (marked with arrows), as for the reference. This is visible in the global circumferential strain (GCS), which increased from 16% in the original sequence to 22% in the transported sequence, with the GCS of the target strain being 20%. The mean absolute difference between the desired and obtained longitudinal and circumferential strain (averaged over all cells) was 4% (± 6%) and 3% (±3%), respectively. In comparison, the standard deviation of the strain across all triangles was 19%. This mismatch is mainly due to the fact that for a deformation to be exactly applicable to a surface, its elements need to jointly satisfy the Codazzi-Mainardi equations.

Discussion
Our strain-based methodology was able to transport a given deformation (either synthetic, or from another individual) to the mesh of another individual. We obtained very promising results in our preliminary experiments. Perfect reconstruction was not obtained, as it is not possible to exactly impose a strain field and recover a valid mesh: the strain field must satisfy the discrete Codazzi-Mainardi conditions. Contrary to previous methods [1,7,9,11], we use strain to model cardiac dynamics. While this is very well established in the clinical community, and has strong physiological reasons, the engineering community has focused on modelling displacement, and afterwards computing strain. This approach has issues, as displacement has a high dependence on the shape as well as the strain. Our strain-based method allows modelling easily how some pathologies would affect the deformation of an individual, such as a region becoming akinetic. This would be more difficult to model in a displacement-based approach, as the akinetic region would still experience displacement due to pulling from the other contracting areas. The strain representation of the cardiac dynamics is more decoupled from the cardiac shape than displacement, since it represents the actual process of myocyte contraction.
Our method has assumptions: it relies on point-to-point correspondences and assumes that the computed anatomical directions are consistent among different individuals. For the former, we think its effect is lower than for displacementbased methods, as strain is less dependent to shape. Regarding the latter, while it can be a problem for pathologies such as some congenital heart diseases, where anatomy suffers large abnormalities, we believe it represents a more appropriate assumption compared to displacement-based strategies.
Limitations of our methods is that local strain data is noisier than displacement, being its spatial derivative, and the fact that our method is not able to capture any rigid-body motion. The latter is not a big issue, since in most cases, the RV apex is fixed due to pericardial constraints, and rigid motion of the heart is mostly due to breathing artefacts.

Conclusion
We proposed a discrete differential geometry method to deform a surface mesh using an imposed deformation field, based on a mesh reconstruction algorithm.
The method defines strain in a local anatomical system of coordinates, allowing transportation of strain across different individuals. In the future, we would like to apply the method to the left atria and show its usefulness for multimodality: there are modalities with higher temporal resolution (such as 2D echocardiography), from which local strain can be estimated and transported to a higher spatial resolution modality, such as MRI or CT.