Comparison of Displacement-Based and Force-Based Mapped Meshing

The finite element (FE) method is a powerful tool for the study of biomechanics. One of the limiting factors in transitioning this tool into the clinic is the time required to generate high quality meshes for analysis. Previously, we developed a mapped meshing technique that utilized force control and a finite element solver to warp a template mesh onto subject specific surfaces. This paper describes a displacement based method that directly warps the template mesh onto subject specific surfaces using distance as the driving measure for the deformable registration. The resulting meshes were evaluated for mesh quality and compared to the force based method. An initial evaluation was performed using a mathematical phantom. The algorithm was then applied to generate meshes for the phalanx bones of the human hand. The algorithm successfully mapped the template bone to all of the bony surfaces, with the exception of the distal phalanx bone. In this one case, significant differences existed between the geometries of the template mesh and the distal phalanx. Further refinement of the algorithm may allow the algorithm to successfully generate meshes even in the presence of large geometric shape differences. In this paper we present an algorithm for the generation of finite element meshes using a mapped meshing approach. The approach computes the distance between the template mesh and the subject specific surface. The position of the surface nodes of the template mesh are incrementally projected towards the surface of interest. The number of iterations used in this mapping is specified by the user. The interior nodes are then mapped to the new representation using a thin-plate spline transformation. This algorithm is initially evaluated using a mathematical phantom data set consisting of a cube template mesh and spherical surface. Next, the algorithm is applied to map a template mesh to the


Introduction
The finite element (FE) method is a powerful tool, widely used in the field of biomedical engineering. One of the limiting factors in transitioning this tool into the clinic is the time and manual effort required to generate anatomic models. The focus of our work is the coupling of mesh generation to medical imaging data to create a pipeline for the rapid generation of finite element models. The ultimate goal is to develop tools that are able to create high quality hexahedral meshes on a subject-specific basis. Our previous work has included the development of a novel building block approach [1] and a deformable registration algorithm that utilized forces to drive the registration [2]. The building block approach is the core of IA-FEMesh, an interactive meshing tool for the generation of high quality hexahedral meshes [3]. Using this tool, we have been able to generate models of several long bones in the human body relatively quickly (on the order of minutes) as well as more complex structures such as the spine. This tool has reduced the time required for mesh generation by an order of magnitude over traditional meshing tools. Once a mesh of high quality is established, especially for complex geometries, it would be ideal to map it to surface representations of similar size and shape, as opposed to redefine the mesh.
Since Couteau et al. [4] proposed a mesh-matching algorithm for automatic grid-based finite element model generation, a number of researchers have published methods for warping a template mesh onto a structures of interest. Jaume et al. [5] labeled the brain surface using a deformable multi-resolution mesh. Gibson et al. [6] generated a finite element mesh from the surface of an adult head subsequently fit to the surface of a neonatal head. We have previously explored the use of a deformable registration based on the finite element method [2]. In this algorithm, the surface and the template were first aligned using an affine transformation to account for differences in the position, orientation, and scale between the datasets. A hierarchal approach was then used to deform the template mesh onto the subject surface by applying forces to the surface nodes.
Using a mapped meshing approach, a convergence study could be conducted on a bone type of interest (e.g., proximal phalanx bone of the hand) to establish the mesh density appropriate for this structure in only a single subject. Once the template mesh is created for this subject, it can then be mapped onto a subject specific surface. An appropriate mapping is one that can accurately align the mesh to the subject surface while maintaining high quality elements in the resulting mesh. In this work, we present a displacement based mapped meshing solution and compare the results to the force based solution that has been previously reported. The main impetus for this work was to decrease the time required for the mapped meshing algorithm. Our apriori hypothesis was that the displacement based solution would be faster, but would be more prone to the generation of distorted elements in the resulting mesh.

Displacement Based Mapped Meshing
The input to the mapped meshing algorithm is the template mesh and a subject-specific surface. The initial registration defines an affine transformation to bring the template mesh into close correspondence with the target surface; thereby accounting for differences in spatial orientation, size, and position.
Thereafter, the registration algorithm is used to define a nonlinear mapping to locally register the surface nodes of the template mesh onto the bony surface. The registration algorithm implemented is based on the distances between the surface nodes of the template mesh and corresponding points on the target surface, identified by intersecting the normals from the template mesh with the triangulated target surface. Each node is repositioned toward the bony surface along the point normal. Consequently, if the node is positioned outside the target surface the distance is negated. The nodal positions are then updated according to the following equation: Where x′ is the updated nodal position, x is the original nodal position, n is the normal direction, d is the distance between the nodal position and the subject surface, N is the total number of iterations, and i is the increment number (i = 1…N). If the current iteration is not the last iteration, the resulting deformed mesh is smoothed using Laplacian smoothing. This was done to maintain a smooth representation of the deforming mesh while eliminating large local changes in the mesh resulting from the projection techniques. This process is repeated for the number of iterations specified by the user. The incremental adjustment of the nodal positions was implemented to allow for greater variation between the template and surface without generating distorted representations of the surface faces.
Once the surface nodes have been mapped using the algorithm described above, the next step is to recompute the distribution of the interior nodes. The internal nodes are repositioned using a thin-plate spline transform [7]. The thin-plate spline (TPS) transform is created via the surface nodes of the original template mesh as the source landmark positions, while the final position of the mapped nodal point positions are used as the target landmarks. The resulting transform defines a smooth mapping for the interior nodal positions from their original position to their new position in the mapped mesh. The user has the ability to specify the density of the points used to define the TPS. The complete displacement controlled mapped meshing algorithm is summarized in Figure 1.

Evaluation of the Displacement Based Mapped Meshing Algorithm
Two data sets were used to evaluate the displacement-based mapped meshing algorithm. The first data set consisted of a cube representing the template mesh and a sphere for the subject surface. The cube had dimensions of 22.5 mm × 22.5 mm × 22.5 mm and consisted of 1728 elements (12×12×12). The sphere had a diameter of 20 mm and contained 19,600 triangles.
The second dataset was used to evaluate the applicability of the method to irregular anatomic structures. Two cadaveric specimens were acquired from the Anatomy Gifts Registry located in Hanover, Maryland. CT datasets were collected at The University of Iowa. The specimens were imaged in the axial plane on a Siemens Sensation 64 CT scanner (Matrix = 512×512, FOV = 172×172 mm, KVP = 120, Current = 94 mA, Exposure = 105 mA) with an in-plane resolution of 0.34 mm and a slice thickness of 0.40 mm. The regions of interests defining the phalanx bones of the index, middle, ring, and little finger were defined by a manual rater using the BRAINS image analysis suite [8]. The resulting regions of interest were converted into a surface representation. The proximal phalanx bone of the index finger from the first specimen was meshed using the building block approach in IA-FEMesh. A mesh convergence study was performed to determine the optimal mesh density for a static FE analysis. The mesh density was increased until minimal change in the internal stress values was observed with 6765 elements. This represented an average element edge length of 1.0mm. This mesh was then used as the template mesh for this study. In order to evaluate the ability of the algorithm to mesh structures of similar shape, the template mesh was mapped to the proximal phalanx bones of the index, middle, ring, and little fingers of the second cadaveric specimen. This experiment was undertaken to simulate the variation in bone size that exists across subjects. To evaluate algorithm's performance when significant variations exist in geometries between the template mesh and subject surface, the template mesh of the proximal phalanx bone was mapped to the middle and distal phalanx of the index finger. For all of the experiments, the registration algorithm was run for 5 iterations, 1000 iterations of Laplacian smoothing were used after each step, and all of the surface points were used in the TPS based interpolation.
We previously reported on the application of a finite element based mapped meshing algorithm to the same dataset [3]. In this previous work, forces were applied to the template mesh to drive it into correspondence with a subject surface. The two mapped meshing algorithms are compared in this paper using the metrics described below.

Evaluation Metrics
In order to valuate the resulting registration, the quality of the mesh was checked using an in-house program based on the Verdict library [9]. The mesh quality was evaluated in terms of element volume and distortion. The goal was not to introduce any zero-volume or distorted elments to the mesh, as a result of the registration algorithm. An element with an angle between its isoparametric lines either less than 45° or greater than 135° was considered distorted.

Software
The algorithm was written in C++ and utilized the Visualization Toolkit (www.vtk.org). The software is available as a command line tool within IA-FEMesh (http://www.ccad.uiowa.edu/mimx/IA-FEMesh). IA-FEMesh is an open source software toolkit for rapid anatomic finite element model development. The core of IA-FEMesh employs a multiblock meshing scheme aimed at hexahedral mesh development. An emphasis has been placed on making the tools interactive, in an effort to create a user friendly environment. This tool provides efficient, reliable and valid methods for model development, visualization and mesh quality evaluation. Figure 2 demonstrates the iterative nature of the displacement-based registration process as the cube is mapped onto the sphere. As illustrated, the initial configuration is accompanied by mesh definitions following the first, second, and final (fifth) iterations, respectively. The corner points undergo the greatest displacement during the first iteration due to the fact that they initially reside furthest away from the target surface. Figure 3 illustrates the element volumes for the resulting mapped mesh. None of the resulting elements had a negative volume that would preclude this mesh from a finite element analysis. The resulting element volumes ranged from 0.11 to 5.93mm 3 , as compared to the original element volume of 6.59mm 3 . The resulting mesh did have 584 distorted elements having an angle smaller than 45 degrees or larger than 135 degrees.

Results
Mapping the template mesh to the bony surfaces using the displacement-based technique was successful in all cases except the distal phalanx bone. Figure 4 illustrates the bony surface and resulting mesh for the proximal and middle phalanx bones of the index finger. There is a slight variation in the geometry between the subject surface and the mapped template for the middle phalanx bone. This is seen in the left proximal aspect of the middle phalanx as shown in Figure 4c. The original template mesh consisted of 340 distorted elements or approximately 5% of the total elements. More distorted element existed in the resulting mapped meshes as compared to the template mesh. The number of distorted elements increased as the variation between the template and target geometry increased. The percentage increase in the number of distorted elements is shown in Table 1. While the percentage increase is fairly large (70-159%), this still represents a small fraction of the total elements (8.5-13.0%). The resulting distal phalanx mesh resulted in 584 elements with zero, or negative volume.
The FE-based solution that we had published previously [2] showed similar trends as the displacement based mapped meshing algorithm. The number of distorted elements tends to increase as the geometry differences between the surface and template mesh increase. The finite element method resulted in fewer distorted elements. For FE based algorithm, we found a maximum increase of 5.6% for the number of distorted elements in mapping the proximal phalanx bones across the index, middle, ring, and little fingers. The number increased to 10.6% for the middle phalanx of the index finger. This algorithm also generated a mesh with zero or negative element volumes for the distal phalanx bone.

Discussion
The mapped meshing algorithm based on a deformable registration was successfully applied to the phalanx bones of the hand. In this work, we mapped the phalanx bones of the hand to evaluate the effectiveness of the algorithm on a variety of bone sizes. In evaluating the algorithm across different geometries, the algorithm successfully mapped to the middle phalanx bone of the index finger, but generated a mesh that had several hundred zero volume elements in the resulting mesh. We had previously evaluated a FE-based force control mapped mesh algorithm that tended to produce fewer distorted elements and generated only a single zero volume element in the distal phalanx bone. In general, the finite element solution is less sensitive to differences in geometry between the template mesh and the subject surface. Both algorithms would benefit from a higher mesh density within this region. With more nodal points in this region the algorithms would be able to account for more variability between the template mesh and the subject surface.
One advantage of this algorithm was that we were able to utilize a single ICP based algorithm for the all of the surfaces. The FE-based method required the use of an additional algorithm (i.e., Procrustes) to attain initial alignment of the template with the subject surface. Although this alignment was not required, nor used, during the deformation-based approach, future implementations therefore may benefit from an improved initial alignment. Additional consistency in alignment may improve the overall reliability of the method.
An additional advantage of this algorithm was that the resulting mesh could be generated in about half of the time required by the FE based method. Finally, the displacement based mapped meshing eliminated the need for hierarchal mesh refinement that was used in the FE based method to improve performance. The displacement control method presented here did suffer from an increase in the number of distorted elements. It is anticipated that increasing the number of iterations used for registration, along with an improved initial alignment, will help to minimize the mesh distortion.
In summary, the displacement-based solution presented here was able to successfully map a template mesh to a subject surface. In the presence of large geometric shape differences between the template mesh and individual surface, the resulting mesh may contain zero volume elements. It is anticipated that a better initial alignment and increasing the number of iterations will minimize the generation of poor quality elements.  The cube template mesh shown being registered to a sphere. (A) Before registration, (B) after the first iteration, (C) after the second iteration, and (D) after the final iteration.  Results showing the warped template mesh. The warped template mesh is shown in white and the subject surface in red for the index finger proximal phalanx (A) and middle phalanx (C). The volume mesh quality metric is shown for proximal phalanx (B) and middle phalanx (D).