A practical framework for generating volumetric meshes of subject-specific soft tissue

Abstract

Studying human motion using musculoskeletal models is a common practice in the field of biomechanics. By using such models, recorded subject’s motions can be analyzed in successive steps from kinematics and dynamics to muscle control. However simulating muscle deformation and interaction is not possible, but other methods such as a finite element (FE) simulation are very well suited to simulate deformation and interaction of objects. In this paper we present a practical framework for the automatic generation of FE ready meshes based on subject-specific segmented MRI data. The proposed method resolves several types of data inconsistencies: noise, an incomplete data set and self-intersections. This paper shows the different steps of the method, such as solving overlaps in the segmented surfaces, generating the volume mesh and the connection to a musculoskeletal simulation.

Appendices

Appendix A

1.1 A.1 Push-based method

During this research, we tried several methods of approaching the overlap problem. We first developed a method where two objects would push each other away. The reason for using a push-based method is that the topology of the object remains intact, and the attachment and tendon information are transformed together with the surface. The surfaces that have smaller volumes have a higher priority, meaning they will be the last to be deformed, because their changes will most likely have the most impact on the shape of the object. The bones will have the highest priority since they have the best segmentation quality. First we apply an iterative smoothing algorithm on the pushed surface to represent the general trend that the surface has. We then use this information to determine the direction whereto each overlapping vertex should move when pushed away by the other surface. As we also want to take into account the geometry of the pushing surface, this direction is combined with the average vertex normal of the closest vertices in the pushing surface. This lookup is done using a kd-tree, for which the implementation from the CGAL library is used [26, 33]. We take the 10 closest vertices, calculate the vertex normals of those and average them. This gives us a second pushing direction vector which is combined to the smoothed direction with an equal weight. A vertex is pushed away in that direction until it does not lie inside the other surface anymore. Unfortunately, this method fails if a surface lies too far inside another one, because pushing directions cannot be accurately determined anymore. Furthermore, the convergence of this kind of technique will be difficult to ensure when the data set is made of many objects e.g. with all muscles, fat-tissue and skin layers.

1.2 A.2 Shrink-based method

We also looked at methods that can make two surfaces push each other away by both moving vertices to a shrunk copy of the surface. For this, we need an algorithm that can give for a surface S a shrunk version S′. Shrinking an object iteratively gives a sequence of shrunk surfaces S ¹,…,S ⁿ. Similarly, each vertex v _i ∈V, has a sequence of shrunk positions $v_{i}^{1}, \ldots, v_{i}^{n}$. By applying the shrinking algorithm iteratively on the overlapping surfaces, they will naturally move apart. We tried the following two shrinking algorithms.

Shrinking by smoothing

We used a shrinking algorithm based on standard Laplacian smoothing [22]. Each vertex is moved in the average direction of its neighbors.

$$ v_i^{t+1} = v_i^{t} + \lambda \frac{\sum_{j \in v_i^{av}}{v_j^{t} - v_i^t }}{|v_i^{av}|} $$

(9)

Here λ is the scaling factor of the displacement done at each shrinking step. This method produces smooth shrinking, and can be done in arbitrarily small steps. The problem of this method is that globally the surface will shrink, but locally the surface can expand. Specifically, concave parts of the surface will expand, and convex parts will shrink.

Shrinking to skeleton

The second shrinking algorithm we developed moves the vertices of the surface to a precalculated ‘skeleton,’ which represents the global minimal structure of an object. We determine this skeleton using the method of Au et al. [34]. The algorithm takes a surface S={V,F} and produces a curve-skeleton K=(U,E) with skeleton nodes U and edges E, where $U = \{ u_{1}, u_{2}, \ldots, u_{n_{U}} \}$ are the node positions. The algorithm also produces a mapping coll(v _i )=u _k from vertices to skeletal nodes, which gives for each vertex v _i the node u _k to which it was collapsed in the skeleton. The shrinking algorithm we developed moves the vertices towards an attraction point on the skeleton. The attraction points are dynamically determined by the mapping coll(v _i ), where the influence of the parent node of each vertex increases as the vertex moves closer to the skeleton. The smoothness of the shrinking could be improved by a different strategy of moving towards the skeleton. Unfortunately, the algorithm has a bigger defect: the algorithm produces for some muscles a skeleton that lies outside of the surface shape, which is not tolerable for our application, since the vertices can never go outside of the original shape.

Appendix B

Static MRI scans of the whole lower limb have been performed. The subject was lying in a resting pose. In close collaboration with radiologists, adequate protocols for the imaging of soft and bony tissues was defined. The following protocol was used: 1.5T Philips Medical Systems machine, axial 2D T1 Turbo Spin Echo, TR/TE=578/18 ms, FOV/FA=40 cm/90^∘, matrix/resolution=512×512/0.78×0.78 mm, thickness: from 2 mm (near joints) to 10 mm (long bones).

Not all algorithms in the pipeline can be used without configuration. We describe hereafter the most significant parameters used in the pipeline. The smoothing algorithm (Sect. 3.3) has three parameters: λ factor, μ factor and number N. The surfaces in our data set are smoothed with the values λ=0.33, μ=−0.331 and N=1000 iterations.

In the tendon generation step (Sect. 3.4), we add regularly spaced vertices to fill the gap between the end of the tendon and the start of the muscle. The space between these rows of vertices is set to 2 mm.

The degenerate triangle method (Sect. 3.5) needs two ϵ values, ϵ _e and ϵ _α , that respectively indicate the minimal edge length and minimal angle. In the experiment, we used ϵ _e =0.5 mm and ϵ _α =0.01 degrees.

The overlap solving algorithm (Sect. 3.6) creates a raw-offset of the surface that will be subtracted. The parameter λ indicating the distance of the offset is set to λ=1 mm.

The volume mesh generation step has several parameters determining the final density and accuracy of the output mesh. Table 1 lists the parameters used on our data set. For both muscle and bone objects, we choose the same values for angular bound and radius-edge bound. This is because these values do not influence the amount of elements created in the output mesh, but only influence the runtime of the mesh generation. Facet angular bound is set to 30 degrees because the CGAL algorithm will guarantee a solution for this bound. The radius-edge ratio is set to 1.25, because we want to avoid unbalanced tetrahedra. For the radius bounds, we choose for the muscle a lower value since it will deform during the simulation and therefore needs a higher resolution. The bones are not going to deform so the size of the tetrahedra inside of the bones is not important, hence the high value.

Table 1 CGAL mesh generation parameters