Patient-specific Heart Geometry Modeling for Solid Biomechanics using Deep Learning

A. Problem Formulation

Let $M=(V,E)$ be a target mesh M with vertices V and edges E. Here, $V\in {ℝ}^{N_V\times 3}$ and $E\in {\mathbb{N}}^{N_E\times (2+1)}$, where ℝ is real numbers, ℕ is natural numbers, N_V is the number of vertices, and N_E is the number of edges. Each entry of V is a vertex in 3D space. Each entry of E is (1) a pair of vertex indices that indicates a connection between those two vertices and (2) a mesh component index that indicates the edge belonging to a specific component of M. For our task, we defined 5 mesh components: the left ventricular (LV) myocardium, the aortic wall of the ascending aorta, and the three aortic valve leaflets. This allowed us to design more component-specific losses and adjust various hyperparameters to fit our modeling needs.

The goal of template deformation strategies is to find the optimal displacement vectors at every vertex of a template mesh $M_0=(V_0,E_0)$. Here, $V_0\in {ℝ}^{N_{V_0}\times 3}$ and $E_0\in {\mathbb{N}}^{N_{E_0}\times }{}^{(2+1)}$. Given a set of displacement vectors $\delta \in {ℝ}^{N_{V_0}\times 3}$, the predicted vertices $\overline{V}$ are simply $\overline{V}=V_0+\delta$. We also use an alternative notation for vertex displacement $\overline{V}=\delta (V_0)$, and further extend the notation for the predicted mesh $\overline{M}=\delta (M_0)$ because $E_0$ is identical before and after deformation.

“Optimal” deformation is defined by a loss function $ℒ$, which in our task takes input parameters M, M₀, and $\delta$:

·* denotes the final optimized values. The exact details of $ℒ$ will be further explained in later subsections. Assuming we have a proper $ℒ$, Eq. 1 can be used to optimize the deformation for each patient mesh separately. However, that would complicate the meshing task during inference time due to the need for the ground-truth mesh label M and the need to perform the optimization for each new set of geometries.

In this work, we instead used a neural network $h_{\theta }$ as a function approximator that generates the displacements with a single forward pass with the input image $(h_{\theta }(I;M_0)=\delta )$. Here, $I\in {ℝ}^{H\times W\times D}$ is the input image of height H, width W, and depth D, and $\theta$ is the network parameters. Then, our problem boils down to optimizing $\theta$ with respect to the expected loss over our paired image-mesh training set $\mathrm{\Omega }$:

The overall training process is summarized in Fig. 3. After the model $h_{\theta }$ has been successfully trained, each test set image can be processed via a single forward pass to obtain $\mathrm{\delta }$, from which $\delta (M_0)$ is the final predicted patient-specific mesh with mesh correspondence.

B. Space-deforming Field via CNN

Broadly, there are two approaches for mesh deformation: (1) directly predicting node-specific displacement vectors, or (2) implicitly deforming the embedded object using a space-deforming field [36]. Within our problem formulation, this mainly affects our design choices for $h_{\theta }$. Our final proposed model predicts a regular-grid space-deforming field $\varphi \in {ℝ}^{H\times W\times D\times 3}$, from which we can trilinearly interpolate at V₀ to obtain $\delta$.

Following existing works [15], [37], we modeled $h_{\theta }$ as a 3D CNN, in the form of a modified U-net. The input was a 3D image, and the output was a set of control point displacements at specified isotropic intervals of the original image coordinates. The control point displacements were used to fill in the rest of the image coordinates using a 3rd order b-spline interpolation, which was approximated using iterative mean-filtering [38]. The resulting output was then treated as a stationary vector field of the final diffeormorphic flow, estimated using the scaling-and-squaring method [38], [39].

There are two main advantages to this approach. First, by predicting uniformly spaced displacements using a 3D voxel-grid input, we can leverage robust 3D CNN architectures without any major modifications. In contrast, directly predicting node-displacement vectors has generally involved geometric deep learning, which requires more complex network modifications and additional model parameters [33], [34]. Second and more importantly, we are able to control the field smoothness using the control point spacing instead of explicitly increasing the weighting on the field regularization terms. This exact mechanism helped significantly improve our mesh quality while maintaining spatial accuracy (Section V-E).

C. Loss Functions for Label-efficient Volumetric Meshing

Our proposed loss $ℒ$ is a linear combination of an accuracy term ${ℒ}_{acc}$ and a volumetric mesh quality term ${ℒ}_{mesh}$:

P(M) denotes randomly sampled points on the surface of M [40], and ${\tilde{M}}_0$ is a partial mesh of the template mesh M₀ that spatially matches the provided mesh labels. To better explain ${\tilde{M}}_0$ and M, we first discuss ${ℒ}_{acc}$, which is defined as the mean symmetric Chamfer distances averaged across all mesh components:

In Eq. 4, m_i and ${\tilde{m}}_{i,0}$ denote the i^th mesh component for M and ${\tilde{M}}_0$, respectively. C is the total number of components. In Eq. 5, A and B are arbitrary pointclouds with potentially different numbers of sampled points. In our task, we set the number of samples to be [150K, 40K, 3K, 3K, 3K] for [LV, aorta, leaflet1, leaflet2, leaflet3], which is roughly 15 points/mm² * SurfaceArea(m_i). We used the same number of samples between the target and predicted pointclouds.

To mitigate the difficulty of obtaining training labels, we aimed to reduce the dependency of our model to volumetric mesh labels by using minimally sufficient surface labels for M. The selection criteria of minimal sufficiency depends on each mesh component’s unique structural characteristics. Thus, the definition of m_i and the associated ${\tilde{m}}_{i,0}$ is different for each mesh component. The rationale and selection process for each mesh component is described below.

For the leaflets, segmentation alone is often not spatially accurate enough due to the limitations in the label resolution and quality. This was confirmed in our own experiments, where methods using only segmentation labels often suffered from poor spatial accuracy of the leaflets. (Table I). Additionally, the native leaflet thickness information is often obscured in CTA scans due to the high levels of calcification and limited image resolution. From these observations, we wanted leaflet labels that (1) use a node-based representation and (2) delineate the overall leaflet positioning and curvature instead of the leaflet volume. Based on these criteria, we chose the “base surface” representation from Fig. 2 because it was the simplest label to obtain that satisfies both of our constraints.

For the aorta, segmentation labels can still be problematic. First, CTA scans only allow for good segmentation of the blood pool, which only delineates the inner surface of the aortic wall, not the aortic wall volume. Second, the 0-genus ellipsoidal blood pool segmentation cannot be directly compared against the desired 1-genus tube-like mesh. Therefore, we also chose the “base surface” representation as our minimally sufficient surface label for the aorta.

For the LV myocardium, where the muscle mass is clearly visible in the CTA and is thick enough to be labeled accurately via conventional segmentation methods, the surface label is simply the marching cubes of the segmentation label.

After deciding on our minimally sufficient label representation, we can select parts of M₀ that spatially correspond to those labels. For the leaflets and the aorta, this is the “base surface” (Fig. 2), or the quadrilateral surface elements comprised of the bottom 4 nodes of each hex element (Fig. 4). For the LV myocardium, this is the extracted outer surface of the original volume elements. Note that M and M₀ do not need mesh correspondence because we used the sampled points for measuring surface accuracy.

While the use of surface labels greatly increases the flexibility of our method, it also fails to capture important mesh constraints, such as (1) the element thickness for the leaflets and the aortic wall and (2) the positioning of the inner volume elements for the LV myocardium. This leads to suboptimal performance of existing methods for our task (Table I). To address this, we directly incorporated a volumetric mesh quality term to minimize undue deformation from the initial mesh template.

The key starting ingredient for calculating various deformation energies is the deformation gradient $\mathbf{F}\in {ℝ}^{3\times 3}$ [41]:

where ${\overline{x}}_i$ and $x_i={\delta }_i({\overline{x}}_i)$ are the original and transformed node coordinates, respectively, of each tetrahedral element. Note that the computation of F for hexahedral elements requires the use of quadrature points, but we achieved similar performance by simply splitting all elements into tetrahedra via triangulation and applying Eq. 6.

Using F, we can calculate various isotropic and anisotropic energies to measure the amount of element distortion in the specified directions. For overall isotropic element stiffness, we used the as-rigid-as-possible (ARAP) energy, where the energy density for each element is defined as:

Here, I is the identity matrix, $\parallel \cdot {\parallel }_F$ is the Frobenius norm, and R and S are the rotation and stretch components of F obtained via polar decomposition F = RS. More precisely, $F=U\text{Σ}{V}^T$ via singular value decomposition, from which $R=U{V}^T$ and $S=V\text{Σ}{V}^T$. The first two equalities of Eq. 7 are from the definitions of ARAP energy and polar decomposition, and the last equality is from the rotational invariance of the Frobenius norm.

To ensure end-to-end differentiability, F needs to be full rank (i.e. no degenerate elements) and Σ needs to have distinct singular values (i.e. F ≠ I). In our experiments, both conditions were satisfied as long as we initialized the network output with uniformly sampled values within [−1e-3, 1e-3].

While the sole addition of isotropic energy showed some promise in a previous study [11], the prediction quality was hindered by the trade-off between surface accuracy and element quality due to the over-stiffening effect. Since this effect was observed most prominently for the leaflets, where the element conformation can vary significantly based on the degree of valve opening, the previously proposed solution was to use two mesh templates in open vs. closed valve states and measure the distance-based weighted average of the isotropic distortion energies from the two templates.

In this study, we instead incorporated an anisotropic energy term to maintain the desired overall mesh quality while applying less penalization of the isotropic energy. This allowed for greater flexibility of the leaflet and aortic wall elements, which lead to increased spatial accuracy of the predicted meshes for the desired level of mesh quality, while using a single template.

The anisotropic energy we used was a slight variation of the anisotropic St. Venant Kirchhoff (AStVK) energy, dubbed anisotropic square root (ASqrt) [41]. The energy density for each element is:

where $d\in {ℝ}^3$ is the normalized direction vector for measuring the anisotropic energy. In our task, d was calculated for each leaflet and aortic wall element independently in the thickness direction using the template mesh. More precisely, the 4 direction vectors going from the bottom 4 nodes to the corresponding top 4 nodes of each hex element were averaged and normalized (Fig. 4).

To combine these energy densities in a way that is minimally restrictive, we applied different weighting factors across the thick (LV myocardium) and thin (leaflet and aortic wall) elements:

for J LV myocardium elements and K leaflet and aortic wall elements. This assumes equal weighting across all elements of the same class of structures, but different overall weighting across different classes. $[{\lambda }_0,{\lambda }_1,{\lambda }_2]=[5,1,10]$ for our proposed model. The final loss is now fully defined by combining Eq. 3, 4, and 9.

D. Calcification Modeling

Calcification (Ca2) is crucial for accurate structural modeling and stress simulations. In this work, we first used our proposed deep learning method to obtain the patient-specific heart geometry, and subsequently used traditional techniques to generate and stitch Ca2 meshes onto the predicted structure. A few reasons for the sequential process were: (1) deformation strategies are inadequate for Ca2 due to its arbitrary amount and positioning, (2) Ca2 can be easily segmented using traditional techniques, and (3) most Ca2 attachment points were spatially consistent with our accurately predicted heart structures.

The exact implementation of the Ca2 meshing process described below is available in our published code. However, it should be noted that the Ca2 meshing technique was not the focus of this work; we merely aimed to demonstrate that our heart geometry is accurate enough for reliable incorporation of calcification meshes. We will be investigating faster and more robust calcification meshing techniques in future works.

Ca2 was first segmented using a fixed-window intensity thresholding, where the thresholds were determined empirically for each patient to account for the contrast agent-induced intensity differences. The selected threshold was often around [590, upper limit] HU. The final segmentation was then produced by trimming the portions outside of the boundary of predicted geometries and removing small islands with fewer than 5 connected voxels.

To convert the Ca2 segmentation to stitched meshes, we performed a series of automated post-processing steps. First, we eliminated mesh collisions by subtracting the heart geometry from the Ca2 segmentation. Since mesh boolean operations are notoriously unreliable, we first voxelized the predicted heart mesh, performed the subtraction in the voxel space, and converted the resulting Ca2 segmentation back to surface meshes using isosurface extraction [42].

Then, we extended the Ca2’s attachment surfaces closer to the heart surface. To calculate the direction of extension, we first used surface-normal ray-tracing to determine the first non-self intersection of all heart surface nodes with both the heart and Ca2 meshes. Valid ray-tracings were selected based on distance thresholding ([0, 3] mm for heart-to-Ca2, and heart-to-Ca2 must be closer than heart-to-heart). For each Ca2 connected component, the final extension direction was calculated as the average direction of all valid ray-tracings. We selected Ca2 nodes based on the angular difference between their surface normals and the extension direction, and extended them slightly past the heart surface based on the ray-tracing distances. The extended meshes were morphologically closed and subtracted in voxel space as described in the previous paragraph to prevent mesh collisions. We performed re-meshing via voronoi clustering to uniformly space the Ca2 nodes, and merged the surface-contacting Ca2 nodes onto the closest neighboring heart surface nodes. Surface contact was defined as nodes closer than 0.3 mm. After a few smoothing and mesh cleaning operations, such as removing non-stitched Ca2 meshes, we used TetGen [43] to generate the final Ca2 tetrahedra.

A. Data Acquisition and Preprocessing

B. Baseline experiments

For baseline experiments, we implemented various combinations of existing approaches from Section III.

C. Implementation Details

We used Pytorch ver. 1.13.1 [47] to implement a variation of a 3D U-net for our CNN [48], and Pytorch3d ver. 0.7.3 [40] to implement the GCN. The basic Conv unit was Conv3D-InstanceNorm-LeakyReLu, and the network had 4 encoding layers of ConvStride2-Conv with residual connections and dropout, and 4 decoding layers of Concatenation-Conv-Conv-Upsampling-Conv. The base number of filters was 16, and was doubled at each encoding layer and halved at each decoding layer. The GCN had 3 layers of graph convolution operations defined as $ReLU(w_0^T{f}_i+{\sum }_{j\in 𝒩(i)}{w}_1^T{f}_j)$ and a last layer without ReLU. The input to the initial GCN layer was the concatenation of vertex positions and point-sampled features from the last 3 U-net decoding layers. The GCN feature sizes were 227 for input, 128 for hidden, and 3 for output layers. We found ${\lambda }_i$ for every experiment with a grid search based on validation error, ranging 7 orders of magnitude. We used the Adam optimizer [49] with a fixed learning rate of 1e-4, batch size of 1, and 4000 training epochs. The models were trained with a B-spline deformation augmentation step, resulting in 140k total training samples (35 · 4000). All operations were performed on a single NVIDIA GTX 1080 Ti, with around ~2 days of training time and maximum GPU memory usage of ~6.5 GB. Inference takes 0.13 seconds/scan on the same GPU.

D. Spatial Accuracy and Volumetric Mesh Quality

For spatial accuracy, we computed the mean surface distance (CD), the worst-case surface distance (HD), and landmark correspondence error using the Euclidean distance of specific mesh nodes from 3 commissures, 3 hinges, and 3 leaflet tips (Corr). For volumetric mesh quality, we first evaluated two standard FE mesh metrics: |Jac| and skew [42], [50]. Then, since zero-volume elements and element inversions break FE simulations, we counted their number of occurrences using |Jac| ≤ 0 across all test-set patients. As one measure of deviation from the representative anatomical parameters, we also calculated the magnitude of thickness differences for the thin elements post-deformation (Thick diff). Refer to Table. I for abbreviations.

Our method outperformed all existing methods in every category except landmark correspondence error, where it came in second (Table I, II). It also outperformed the previous state-of-the-art ARAP-only condition in every category, demonstrating its superiority in the single-template deformation task.

Even compared against the double-template condition, our method had the best overall performance. This can be attributed to the fact that the distance-based deformation energy weighting strategy proposed in our MICCAI work is only truly effective when we have the input images aligned at the valve region. The image preprocessing for this work was a simpler and more realistic center crop with respect to the LV-aorta-valve complex, which does not necessarily guarantee valve alignment.

We also tested the viability of our approach for using half of the original image resolution, i.e. 2.5mm³ and 64³ voxels. The method generally performed as expected, with similar element quality and slightly lower spatial accuracy compared to the main result.

DeepCarve’s stitched calcification meshes were evaluated for all test-set patients. We compared the initial calcification segmentation and the final stitched meshes. To convert between the two representations, we used isosurface extraction and voxelization [42]. Our calcification meshing technique maintained the spatial accuracy of the initial input segmentation and generated nice quality elements for simulations (Table III). For all patients, there were 0 elements with |Jac| ≤ 0.

Qualitatively, our model produced patient-specific meshes with great spatial accuracy (Fig. 7). Our approach also significantly reduced the amount of thickness alterations (Fig. 8) and adverse element distortions (Fig. 9) compared to the ARAP-only condition. The final mesh with stitched calcification elements also illustrate the high quality of our final simulation-ready geometry (Fig. 7). We were only able to show one test-set patient’s final mesh due to space constraints, but all other test-set predictions exhibited similar mesh quality for both the heart and stitched calcification elements.

E. Hyperparameter Trade-off Analysis

From Eq. 9, different values of hyperparameters ${\lambda }_i$ can lead to vastly different meshes. For ${\lambda }_i\gg 1$, the predicted meshes would be nearly identical to the template mesh due to high deformation penalties, whereas for ${\lambda }_i\to 0$, the predicted meshes would be unusable for downstream tasks due to poor mesh quality. This is true for the two top baselines as well, i.e. the smooth-field and ARAP-only conditions. We used the idea of Pareto front approximation sets from multiobjective optimization [51] to demonstrate the differences in the degree of this trade-off between the 3 top performing methods (Fig. 10).

The first observation for Fig. 10 is the effect of the b-spline control point spacing in the final mesh viability. Between $\sigma =1$ and $\sigma =3$, where $\sigma$ is a unitless scaling factor ($\text{b-spline control point spacing}=\sigma *\text{image spacing}$), we can clearly see that only the $\sigma =3$ condition is able to generate valid models within the relevant range of ${\lambda }_i$. Although there is some observable difference in the lowest possible surface accuracy attainable by $\sigma =3$ presumably due to the coarser displacement calculations, this difference is minimal compared to the benefits of generating viable meshes. This phenomenon is further detailed in Table IV for 3 different values of $\sigma$ across 2 top performing methods. Following these observations, we used $\sigma =3$ for all subsequent model training and analyses. As briefly mentioned in Section IV-B, this ability to improve field smoothness without additional explicit regularization losses is one of the main reasons why a space-deforming field is preferred over node-specific displacements.

The second observation for Fig. 10 comes from the differences in hypervolume boundaries between the three $\sigma =3$ conditions. Since lower is better for all metrics, boundaries closer to the origin is better. Our method consistently outperformed the others for all metrics, where the improvement is especially large from the smooth-field condition. From the ARAP-only condition, there are consistent improvements throughout the metrics, and they especially stand out for the thickness difference vs. chamfer distance plot. This is especially beneficial for generating valid thin elements for the leaflets, where adjusting only the isotropic energy leads to over-stiffening and reduced spatial accuracy.

In terms of computational costs, the execution time for the most dense b-spline interpolation is in the order of milliseconds, so there is essentially no noticeable effect associated with different control point spacings. The diffeomorphic flow estimation also runs in the order of milliseconds, so it does not incur significant costs.

F. Post-TAVR Stress/Strain Analysis

To demonstrate the feasibility of our auto-generated meshes for biomechanical computations, we performed finite element analysis (FEA) to simulate the TAVR deployment procedure with Abaqus. We used a 26 mm, self-expandable, first-generation Medtronic CoreValve device, and selected 6 test-set patients that are suitable for this device based on each patient’s aortic valve annulus diameter and the manufacturer’s sizing recommendations [52]. We were not able to test all test-set patients due to our limited inventory of TAVR stent geometries. The superelastic Nitinol stent was modeled with a built-in material library in Abaqus [53]. Briefly, the TAVR simulations were performed in these three main steps:

1. The TAV stent was aligned co-axially within the aortic root and centered into the aortic annulus. The stent was positioned such that the lowest point of the crimped stent is 5mm below the aortic annulus [52].

2. The stent was then crimped with a cylindrical sheath to an exterior diameter of 6 mm [12].

3. The crimped stent was released by axially moving the cylindrical sheath away from the aortic root. The tissue-device contact interaction was then enabled during this procedure, with the stent-to-tissue friction coefficient set to 0.1 [54]. The finite element nodes at the distal end of the aorta and bottom of the left ventricle were constrained, which allows radial expansion only.

The deployed configuration of the stent with a patient’s left heart model is shown in Fig. 1 11, 12. The maximum principal stress and the maximum principal logarithmic strain distributions are shown in Fig. 11, 12. The stress was 0.1~3 MPa, and strain on the aortic root was 0~0.4. The stress concentrations and the maximum strain and stress on the leaflets occurred near the calcification. These results are in line with previous studies [12], [13], which indicates that our auto-generated meshes produce accurate simulation results. The approximate run-time of each simulation was around 10 hours using 3 Dual Intel Xeon Gold 6226 CPUs @ 2.7 GHz (24 cores/node). In future work, we will implement an FEA input file generator to fully automate the simulations.