2.3 Voxel Conforming Mesh Generation

In this step, a tetrahedral mesh that conforms exactly to the voxelized regions of interest is created. Every 1-by-1 square on the boundary is composed of two tetrahedral faces. On each side of the boundary, two tetrahedra, taken together, form a pyramid, whose base is the 1-by-1 square of the boundary and whose apex is exactly at the center of a unit cube formed from a single voxel. In this way, no tetrahedra are permitted to have all 4 nodes on the boundary.

For the interior tetrahedra, an adaptive body-centered cubic (BCC) mesh based on

Molina et al. (2003) is used. It begins by creating large tetrahedra and refines each of these into 8 identical smaller tetrahedra. This process is continued in an iterative fashion until the degree of fit is satisfied. The signed distance function (

Osher and Sethian 1988) is used to determine how many refinements are made. Tetrahedra near the boundary are refined many times, making them smaller, while the interior tetrahedra are refined less and so will tend to be larger. T-junctions can be formed when a node of one tetrahedron lies on the edge of another tetrahedra and results in an unusable mesh. These are identified and systematically removed allowing interior tetrahedral elements to be adaptively sized and have good quality.

Most meshing methods allow for adaptively sized elements on region boundaries. This is an extremely useful property for meshing an object of arbitrary detail that is flat or nearly in some areas. For example, a mesh of an airplane might only need large tetrahedra for the wings but smaller tetrahedra would be needed for more detailed areas such as the engine nacelle. Medical images, such as those of the brain, usually do not have large flat areas where large tetrahedra can be effectively used. Two other image-to-mesh methods referenced in the introduction,

Zhang et al. (2005) and

Fang and Boas (2009), employ tetrahedral generation methods with the capacity for adaptively sized boundaries, but the boundaries of segmented brain regions end up being universally small on the boundaries. In contrast, the methods that mesh arbitrarily detailed boundaries of man-made objects make effective effect use of adaptively sized boundary tetrahedra. Additionally, many image processing applications like registration are concerned with the pooling of data to assess between-group differences. A smooth boundary in one mesh should be able to be registered to a less smooth boundary in another mesh or image. Keeping the boundaries at a uniformly high level of detail allows for more types of deformation at the level of detail allowed by the image resolution.

2.4 Boundary Smoothing Method

The effects of voxelization are removed while maintaining the overall geometry of the regions. For the notation of the equations below, we utilize symbols that are standard in continuum mechanics (such as

*x*,

*X, F*) and some that are standard in level set applications (such as δ, Ω) when possible. The initial voxel conforming state is given by

*X*, while the deformed state is denoted by

*x*, and the displacement vector field, ν, is given by:

The smooth mesh is obtained by finding a

*v* that minimizes an energy

*G(ν). G(ν)* has three additive terms, an elastic term,

*E(ν)*, a smoothness term,

*S(ν)*, and a fidelity term

*B(ν)*. The formula for the energy is:

The elastic energy prevents dihedral angles from becoming too small and tetrahedral elements from collapsing. In places where the voxelized region is smooth and well shaped, the elastic term has no effect on the location of the boundary nodes. It is undesirable to have the boundary moved unnecessarily to improve element quality. In poorly shaped regions (most likely from errors in the automatic segmentation), the smoothing should be done as much as possible, but the elasticity should still maintain element quality and affect the location of the boundary nodes. These features allow to the proposed meshing method to succeed in producing a mesh with uniformly good element quality and smoothness and/or accuracy consistent with the accuracy of the initial segmented input image. We employ techniques from nonlinear elasticity (

Hughes 2000), since these allow larger and smoother deformations.

The displacement vector field is related to measures of deformation

*F* and

*C*. Here,

*I* refers to the identity matrix.

*F* is often called the deformation gradient and

*C* is referred to as the right Cauchy-Green Deformation Tensor (

Gonzalez and Stuart 2008).

As done with many existing elasticity models such as the Mooney-Rivlin elasticity (2000), an elasticity penalty is computed using the three invariants of

*C*. If λ

_{1}, λ

_{2}, λ

_{3} are the eigenvalues of

*F*:

The first invariant is most affected by large amounts of stretching, the second by large amounts of shearing, and the third by significant changes in volume. The elasticity will provide no resistance to small deformations, but will resist large stretch, shear, and volume changes. While this is not a physically realistic model for solids, it is effective for spreading deformations evenly throughout the mesh. This allows the elasticity to have as little effect as possible on the boundary location while enforcing uniformly high element quality. We first introduce the following functions (as part of the nonlinear elastic component),

and the elastic penalty term is given by:

where

*k*_{1},

*k*_{2},

*k*_{3},

*e*_{1},

*e*_{2},

*e*_{3} are parameters of the meshing model. The function

*W*_{3} will resist a change in the volume. The number

*e*_{3} when greater than 0 will allow the volume to decrease before any penalty is applied. The functions

*W*_{1} and

*W*_{2} will apply a penalty if too much stretching or shearing, respectively is applied relative to the volume change. They are initially calibrated to be 0, regardless of

*e*_{1} and

*e*_{2}, if the deformation is equal in all directions,

*e.g.* if the three eigenvalues of

*F* are all equal. The values of

*e*_{1} and

*e*_{2} when greater than 0 allow modest stretching and shearing to occur before any penalty is applied. All parameter values used are provided in a detailed listing in

Appendix A. This elasticity,

*E*, is a function of all three invariants of the deformation and should resist all possible types of changes in the shape of a tetrahedron. As discussed in the results section, in practice no tetrahedra were noted to approach becoming inverted and a variety of measures show the elements to be of good quality.

Another possibility for maintaining element quality is to directly optimize common quality measures, such as a tetrahedron's dihedral angles (

Freitag and Ollivier-Gooch 1997). However, in these cases, it may be difficult to adapt a general quality measure to limit its effect on the boundary (if the boundary is not fixed). The proposed elasticity is also well suited for optimization where as differentiating a general quality measure may be more susceptible to local minima.

The second term of

*G(ν)* is a smoothing term,

*S(ν)*, involving level sets,

_{0,m}. Each region 1, ‥,

*m*, …

*N*, has its own associated level set. No topological changes are permitted at this stage and the level sets themselves are not evolved. Instead, a map ν is morphed, as performed in

Le Guyader and Vese (2009). We define the smoothing term using the harmonic energy, as follows,

where we recall the Heaviside function and its distributional derivative delta function,

Recalling that

_{0}(

*x*) is the distance function to the union of all boundaries of segmented regions, we use in the above smoothing term another function, defined for each

*m* between 1 and

*N*, as the distance function inside region

*m*, and zero outside:

An approximation to the delta function is typically used when functionals of this type are used in finite difference applications. In the proposed method, only the boundary is moved to minimize the above functional and the interior nodes are allowed to position themselves in a way that maximizes element quality.

The third term of

*G(v)* is a data fidelity term,

*B(ν)*, similar to ones used in

Chan-Vese (2001) and

Le Guyader and Vese (2009). Inside the integral for each region, there are three factors multiplied together. The product of the first two factors is either just, at each

*x*, 0, if the mesh label and image label match, and 1 if the mesh label and image label do not match. The third multiplicative factor increases the penalty as a node moves farther away from its corresponding voxelized region. The proposed data fidelity term is thus given by

where we defined for each

*m* between 1 and

*N*, the label function of each region,

and

As with the smoothing term, only the boundary nodes are moved to minimize the above functional. However, to avoid overly strict enforcement of the voxelization, the value of the Heaviside function on the boundary nodes is taken to be a weighted average over nodes in the surrounding volume using the connectivity of the mesh. If the constant *c* is 0.5, no penalty is applied until a boundary node is half a voxel away from its corresponding region.

The minimum of the above energy is found by solving the Euler-Lagrange equation via gradient descent:

At each time step, the

*L*^{2} gradient and Laplacian of the level set function are computed independently using linear basis functions and a straightforward application of the finite element method (

Bank 1996;

Mezger, Thamszewski *et al.* 2009). Computationally, this can be done quickly and does not require solving a system of equations at each time step. In contrast, a strict use of the finite element method for non-linear differential equations frequently requires solving a nonlinear system of equations, which can be computationally expensive.

To further improve the speed of the numerical implementation, the Sobolev Gradient Method (1997; 2005) is employed by solving

The Sobolev Gradient implementation requires a linear system to be solved at each time step where *K>0* is a constant step-size parameter.