The role of computational continuum physics in cardiac research has expanded beyond an engineering analysis of devices or even idealized organs, to become a platform for discovery and integration. This thrust has been paralleled by a major research focus on biomedical image processing (www.itk.org
). Not surprisingly, biomedical geometries are increasingly derived from imaging modalities such as magnetic resonance or computed tomography. Not only can a connection to imaging data provide critical geometric information — atherosclerosis for example is arguably a disease with a strong geometric component — but also critical physiologic data, such as diffusion anisotropy, genetic expression or even metabolism. Biomedical geometries, however, are complex and hierarchical. Despite this, little attention in the biomedical literature has been paid to the enormous challenges posed by geometric computation of biomedical geometries.
Efficient unstructured mesh generation of these geometries in a way that is tuned to the physics of the problem is challenging. Finite computational resources dictate that the grid must be adapted to the geometry in order to minimize both computational cost and discretization error. At the same time, many biophysics problems require that the grid be organized both by scale and by intrinsic properties. The wall of the heart, for example, is a laminated muscle consisting of three separate layers, each with a separate family of fiber and myocyte orientations. Thus, a grid of the heart muscle must be similarly layered. That same layering persists at the smaller scale of the heart geometry, where, for example the musculature transitions into the connective tissue of the cardiac valves; or at an even smaller scale where the heart muscle involutes to become a network of coronary arteries. These transitions of scale are of course mirrored in the blood. The lumen of a coronary artery, for example, may be several orders of magnitude smaller than the span across a ventricle, and thus there is a need to equilibrate error over a range of meaningful scales. These challenges of course are not limited to the heart.
However, the heart does pose some unique challenges. First, the effects of blood flow on cardiac mechanics and physiology are profound. From a biological viewpoint the interface is the zone of greatest activity and communication; from a medical point of view it is the site of the greatest vulnerability to disease. Thus, biomedical fluid-structure problems must aim to resolve continuum fields accurately at the interface. To address these challenges, we have presented the numerical framework to accomplish a Lagrangian treatment of tissue-blood interfaces [12
]. Work in that area continues. Such an approach is greatly simplified if exact nodal compatibility exists between each region of the grid. Secondly, unlike arteries or veins in isolation, most of the heart structures cannot be approximated by a plane-stress membrane assumption. This implies the explicit discretization with volume elements such as tetrahedra, hexahedra and prisms. Quality discretization of complex geometries with hexahedra remains forbidding.
We have presented a method for computing high-quality scale-invariant grids from such imaging data, and have exemplified the approach with a demanding, perfusion fixed micro MRI dataset of the in situ mouse heart and its structure. The mesh quality achieved for this geometry is excellent by any measure. This quality owes to three factors: a careful and topologically correct segmentation to produce initial surfaces, the adoption of the gradient-limited feature size field, and variational optimization of marched prisms based on the face-offsetting method.
The concept of scale is directly embedded into all aspects of our approach by efficient computation of the gradient-limited feature size field. Were edge lengths based solely on some geometry-independent values, the number of elements required would increase by two or three orders of magnitude to capture the finest structures. In contrast, adoption of the GLFS field enables the efficient placement of mesh only where it is needed. Further savings are achieved by decoupling the mesh densities in the surface tangential and normal directions, though tangential density is locally isotropic. Adoption of the GLFS also enables adequate circumferential discretization of tubular regions of the geometry (Equation 4
), and a finite boundary layer to resolve the stronger gradients that are found at the wall.
To explicate our method, we have presented detailed pseudo-code for novel aspects of our methods. Further details may be found in earlier references [11
]. Our methods are applicable to a wide variety of multimaterial data sets and are particularly well-suited to segmentation of perfused tissues. The objective here is to enable others to reproduce our results on related data sets.
Despite the success of our approach, there remain at least four areas of potentially substantial improvement. First, the grid density tangential to the surface is isotropic. On structures such as the cordae or coronary arteries, where fluid flow is largely axial, it would be considerably more efficient to introduce anisotropy into the grid, such that the greater density were arranged in the direction with the highest curvature. Such an approach has been successfully applied to vascular trees in the past [30
]. A promising approach that is compatible with our current methods is described in [18
]. A second area of improvement would be to relax the requirement of constrained Delaunay tetrahedralization for the blood core. Discretization of this zone can be considered preliminary since adaption in our fluid solver is driven by an error metric, and thus discretization will evolve. Nevertheless, the initial quality of this interior grid is important. Although the current approach will typically generate high-quality hybrid meshes for accurate and efficient numerical simulations, the use of constrained Delaunay tetrahedralization can limit the quality of the tetrahedral core for some complex biological geometries. Thirdly, it must be noted that while the overall method does not depend on “good” surfaces, it does depend on “correct” surfaces. In most cases, surface adaptation and volume conserving smoothing improve the input surfaces sufficiently. Application of the latter is critical for mesh improvement, and has the advantage over conventional smoothing approaches that the volume conservation constraint limits drift. However, our smoothing algorithm does not guarantee against self-intersection, and therefore needs to be applied with care. With complex geometries, a robust volume-conserving smoothing with theoretical guarantees against self-intersection would be an important asset. Finally, prisms have inherent advantages with respect to tetrahedra in terms of orthogonality and their avoidance of numerical pitfalls, such as locking1
. Additionally, prisms can potentially reduce the element count in a given zone by a factor of 2–3. Thus, the face-offsetting method has additional potential as a discretization tool, if select nodes can be merged, collapsing prisms into a small number of pyramids and tetrahedra. Addressing these four areas of improvement is part of our future work.
Though computational geometry has received little attention in the biomedical literature, there are some notable exceptions. A number of references to discretization of biomedical geometries can be found in the numerical methods literature. We cite only a few here that have taken a comprehensive view of biomedical geometries that spans imaging and geometric analysis. For example, a method described by Labelle for “isosurface stuffing” [26
], fills an isosurface with tetrahedral mesh with excellent quality statistics. Though this method has not been applied to biomedical geometries, that we are aware of, it is certainly applicable. While this approach may be appropriate for some applications, it does not organize the grid by layer or scale and is generally unsuitable for fluid-flow simulations. A somewhat related approach that produces multiple material lattice-like tetrahedra on biological domains has been presented in [37
]. This approach is similarly suited to a range of applications that fall outside cardiac fluid-structure interactions because the hexahedra are not aligned with the surface and there is no recognition of scale. A promising alternative is the isogeometric approach pioneered by Hughes et al. [15
], and adapted for fluid-structure analysis of treed biomedical geometries by Zhang et al. [36
]. This method has the dual advantages of discretizing volumes with quality hexahedra and being fully integrated with the numerical solver. However, it is unsuited to domains that cannot be reduced to a one-dimensional network of curves. Thus, it is well suited to arterial trees, but unsuited to the entire cardiac geometry.
In contrast, our method uses a gradient-limited feature size (GLFS) to guide the adaptation of the surface mesh, the generation of layered tetrahedra, the generation of the prismatic mesh at the boundary layer and the generation of tetrahedral mesh in the interior. This achieves scale-invariance such that even the smallest resolved coronary vessel is guaranteed to have a preset number of layers across its diameter. Furthermore, guided by the GLFS, we are able to laminate the cardiac tissue, seamlessly transitioning from myocaridum, to cardiac valve to coronary lumen. In addition, this approach enables the automatic adaption of the boundary layer by scale. Our method is also automatic in the sense that the user must only specify a few meaningful parameters: number of layers, degree of anisotropy, growth ratio for the prisms, and minimum and maximum resolved feature sizes.