Efficient tetrahedralization of multi-material images with quality, fidelity, and topological guarantees

Highlights

• The mesh offers a faithful topology to the materials.

• Elements with arbitrarily small angles do not appear in the mesh.

• The mesh offers a close representation (fidelity) of the underlying materials.

• The mesh contains a small number of elements that comply with the guarantees above.

• The mesh can be constructed within tight time constraints.

Abstract

An algorithm is proposed for generating three-dimensional unstructured tetrahedral meshes of multi-material images. The algorithm produces a mesh whose boundary is proved to be homeomorphic to the object surface. In addition, it provides a guaranteed dihedral angle bound of up to $19.47^{\circ }$ for the output tetrahedra. Moreover, it allows for user-specified guaranteed bounds on the distance between the boundaries of the mesh and the boundaries of the materials. It produces a small number of mesh elements that comply with these guarantees, and is compatible in performance with other software. The algorithm offers all of the following properties: bounded dihedral angles, customized two-sided Hausdorff distance, topological faithfulness, a small number of elements, and run-time performance. The experimental evaluation on synthetic and real medical data illustrates the efficiency and effectiveness of the method.

Introduction

The algorithm proposed in this article is an improved and extended version of the preliminary results published in the conference paper [1]. The algorithm constructs tetrahedral volume meshes from segmented multi-material images, and it is the first to generate a mesh which satisfies all the following criteria:

1. The mesh offers a faithful topology to the materials, i.e., there is a homeomorphism between the boundary of the mesh and the boundary of the materials. The term homeomorphism, also called a continuous transformation, is a topological notion of equivalence. Since space-tiling schemes can miss important details of the surface, correct topology usually needs extra efforts to be guaranteed. For implicit surfaces, tools such as critical points [2], [3], Lipschitz conditions [3], [4], [5] and interval arithmetic [6] can be used to capture the topological details. For images, a sufficient condition (single manifold condition, see Section 2.1) is given for the approximation to offer homeomorphism.

2. Elements with arbitrarily small angles which cause the stiffness matrix in finite element analysis to be ill-conditioned do not appear in the mesh. Specifically, all dihedral angles are guaranteed to be above a user-specified lower bound which can be set to any value up to $19.47°$. The algorithm exposes the parameter that allows for a trade-off between the minimum dihedral angle and the final number of elements. In contrast, the state-of-the-art guaranteed quality Delaunay method only allows for a bound on the circumradius-to-shortest-edge ratio of tetrahedra, which in three dimensions is not equivalent to a bound on the smallest dihedral angles.

3. The mesh offers a close representation (fidelity) of the underlying materials. In particular, the two-sided Hausdorff distance between the boundaries of the materials and the boundaries of the mesh respects the user specified fidelity bounds. The Hausdorff distance measures how far two subsets of a metric space are from each other (see the definition in Section 2.2.2). This bound can be zero, which means that the mesh is an exact match to individual pixels along the boundaries of the object. However, a strict matching to image boundaries will produce a large number of elements that will slow down the solver. In some cases, the analysis requires a loose fidelity bound but a faster solver, for example, for image registration [7]. The presented solution exposes parameters that allow for a trade-off between the fidelity and the final number of elements.

4. The mesh contains a small number of elements that comply with the three guarantees above. The reason for this criterion is that the cost of assembling and solving a sparse system of linear equations in the finite element or similar method directly depends on the number of tetrahedra [8], [9]. To achieve this goal, the proposed method is able to offer the ability to transition from small to large elements over a relatively short distance, not only in the volume but also on the surface. A specialized post-processing procedure was implemented to decimate the mesh such that the number of tetrahedra in the mesh is significantly reduced.

5. The mesh can be constructed within tight time constraints enforced by clinical simulation and applications that require interactive re-meshing such as surgical simulations and image-guided interventions. The efficient implementation described below is compatible in performance with other software such as a state-of-the-art Delaunay code.

The definition of the input domain plays an important role in the design of the algorithm and in the guarantees it can offer. Meshing techniques have been proposed to mesh the domain bounded by a Piecewise Linear Complex (PLC). A PLC can be composed of polygons of any shape. The challenge is that the quality of the input PLC affects the quality of the final mesh, because the mesh has to match exactly to the boundaries of the model. For images, the faces of each boundary voxel can be considered as the input PLC, since these faces meet at angles of 90° or 180°, thus alleviating this challenge. The most popular technique for generating tetrahedral meshes in this category of techniques is Delaunay refinement [10], [11], [12], [13], [14], [15]. However, the problem with Delaunay refinement is that it only satisfies a bound on the circumradius-to-shortest-edge ratio, which works well in two dimensions, but in three dimensions it does not imply a bound on dihedral angles. As a consequence, tetrahedra with arbitrarily small dihedral angles called slivers can survive. Sliver exudation [13], [16], [17] or other optimization post-processing techniques [10], [12], [18], [19], [20] can be used to eliminate degenerate elements, but those methods are unable to guarantee meaningful dihedral angle bounds(they exist in theory but too small to even be computed).

Meshing techniques also have been proposed to mesh the domains specified implicitly by a function $f:R^3\to R$. One guaranteed-quality technique in this category is based on the Delaunay refinement [2], [21], [22], [23], where a piecewise-linear approximation of a surface is generated from a finite set of sufficiently dense sample points [4], [5], [24]. Another guaranteed-quality technique in this category employs a space-tiling background grid to guide the creation of a mesh [25], [26], [27], [28], [29], the focus of this article.

Isosurface Stuffing [25] approximates general surfaces with guaranteed-quality using interval arithmetic to account for vertex movement. If the surface is a smooth manifold with bounded curvature, it provides the two-sided Hausdorff distance guarantee between the model and the mesh. However, the presented algorithm is able to achieve a minimum dihedral angle bound of $19.47°$, higher than the angle bound $10.7°$ the Isosurface Stuffing provided. Expanding on the ideas from Isosurface Stuffing, Walker [26] proposed a method for generating quasi-uniform tetrahedral meshes of solids whose boundary is a smooth surface using edge rearrangement. If the lattice spacing is smaller than the local feature size, then the dihedral angles are bounded above $11.4°$. The proposed algorithm is able to achieve a higher minimum dihedral angle bound even with a fully graded case using an octree decomposition. Foteinos et al. [23] presented a Delaunay meshing algorithm with several mathematical guarantees. Using a strategy called $\epsilon$-sampling [4], [5], the mesh boundary is proved to be ambient isotopic to the object surface. Later Feng et al. [30], [31] extended this theory to hybrid parallel Delaunay image-to-mesh conversion algorithm on both distributed-memory and distributed shared memory architectures. Although the circumradius-to-shortest-edge ratio of the tetrahedra in the output mesh is proved to be less than 1.93, they did not provide any dihedral angle bound. Even if the radius-edge ratio is very small, it cannot avoid the almost flat tetrahedra. Furthermore, the two-sided Hausdorff distance between the object surface and mesh boundary is bounded by $O\left({\delta }^2\right)$, where $\delta$ is a small constant 0.0168. However, if the constant implied in Big-O is very large, the two-sided Hausdorff distance bound could also be very large.

The rest of the paper is organized as follows. In Section 2 the proposed algorithm is described along with the necessary definitions and preliminaries. The implementation details and evaluation are presented in Section 3. Section 4 concludes the paper.