- Email: [email protected]
Computational biomechanical modelling of the lumbar spine using marching-cubes surface smoothened finite element voxel meshing
Computer Methods and Programs in Biomedicine (2005) 80, 25—35
Computational biomechanical modelling of the lumbar spine using marching-cubes surface smoothened finite element voxel meshing Z.L. Wang a, J.C.M. Teo a, C.K. Chui a, S.H. Ong b,c, C.H. Yan b, S.C. Wang d, H.K. Wong e, S.H. Teoh a,∗ a
Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, Singapore 117576, Singapore b Department of Electrical & Computer Engineering, National University of Singapore, Singapore c Division of Bioengineering, National University of Singapore, Singapore d Department of Diagnostic Radiology, National University Hospital, Singapore e Department of Orthopaedic Surgery, National University Hospital, Singapore Received 14 January 2005 ; received in revised form 7 June 2005; accepted 7 June 2005 KEYWORDS Modelling; Finite element method (FEM); Voxel-based mesh; Marching-cubes; Volumetric mesh
Summary There is a need for the development of finite element (FE) models based on medical datasets, such as magnetic resonance imaging and computerized tomography in computation biomechanics. Direct conversion of graphic voxels to FE elements is a commonly used method for the generation of FE models. However, conventional voxel-based methods tend to produce models with jagged surfaces. This is a consequence of the inherent characteristics of voxel elements; such a model is unable to capture the geometries of anatomical structures satisfactorily. We have developed a robust technique for the automatic generation of voxel-based patientspecific FE models. Our approach features a novel tetrahedronization scheme that incorporates marching-cubes surface smoothing together with a smooth-distortion factor (SDF). The models conform to the actual geometries of anatomical structures of a lumbar spine segment (L3). The resultant finite element analysis (FEA) at the surfaces is more accurate compared to the use of conventional voxel-based generated FE models. In general, models produced by our method were superior compared to that obtained using the commercial software ScanFE. © 2005 Elsevier Ireland Ltd. All rights reserved.
1. Introduction ∗ Corresponding author. Tel.: +65 6874 6345; fax: +65 6777 3537. E-mail address: [email protected] (S.H. Teoh).
There is a current trend towards the use of computational biomechanics in clinical studies. Prior to biomechanical computation, there is a need for the development of finite element (FE) models based
0169-2607/$ — see front matter © 2005 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2005.06.006
26 on medical datasets, such as magnetic resonance imaging (MRI) and computerized tomography (CT). After the domain boundaries have been defined for automatic finite element meshing of anatomical structures, algorithms are available to fill the domain with tetrahedral elements or hexahedral elements. Tetrahedral elements are used most frequently as this element type is geometrically versatile and can conform to the complex boundaries of anatomical structures. Linear tetrahedral elements usually have the problem of being overly stiff due to their ‘‘odd’’ shape relative to the hexahedral element in theory [1], but it is unclear from published literature that hexahedral elements are better than tetrahedral elements in practice. Tetrahedral elements have the advantage of being able to better model the complex shape of the object. The voxel-based FE meshing method directly converts computer images from 3D voxels to eightnode hexahedral finite elements and fills the domain with hexahedral elements. Despite the intuitive value of such a direct voxel-based method of producing FE meshing models, this method is unable to capture surface boundaries satisfactorily. Jagged or digitized surfaces are unavoidable in voxel-based FE models, which in turn leads to inaccurate finite element analysis (FEA) results. They, however, are still employed for computational biomechanical investigations [2—4]. One reason is that voxel-based modelling also has the ability to automatically assign orthotropic bone mechanical properties based on greyscale intensities [5], thus allowing heterogeneous properties throughout the FE model of the anatomical structure. Voxel-based FE models have been used in biomechanical studies of human distal radius [6], femur [7] and vertebral body [5,8] with some success. Current software that provides computational biomechanical modelling based on medical datasets includes ScanFE and Mimics. ScanFE (Simpleware Ltd., Exeter, UK) has the capability to manipulate medical datasets and obtain FE models automatically. It also produces a voxel-based mesh of the interior together with a tetrahedral surface. Mimics (Materialise, Leuven, Belgium) can produce a surface mesh, which then has to be exported into other FE pre-processors for the generation of a volumetric all tetrahedral FE mesh. Because of the surface irregularities inherent in the models produced by conventional voxel-based methods, Grosland et al. [9,10] proposed a technique to smoothen the irregularities for the FE contact stress analysis of articular joints. Their method generates smooth surfaces at plausible con-
Z.L. Wang et al. tact regions with local refinement. However, due to the use of multiple point constraints (MPCs), the object is modelled in two resolutions and the noncontact regions with jagged surfaces still tend to produce inaccurate FEA results. Camacho et al. [11] also proposed a smoothening algorithm to eliminate the sharp edges of voxel models at surfaces and at the interfaces between materials without the introduction of new elements. This method was extended and implemented in three dimensions (3D) for the FE study of a vertebral body [4]. However, the reported performance of the finite element analysis is not significantly better than conventional voxel-based methods. We have developed a new approach that overcomes the above mentioned deficiencies. The marching-cubes method is employed to extract the smooth surfaces of anatomical structures and this followed by the use of a novel tetrahedronization scheme to generate a voxel-based mesh. Hence, the outer layer in our mesh model is composed of tetrahedral elements while the internal elements are all hexahedral. The two layers share the same set of nodes and are connected seamlessly, thus maintaining FE compliancy. Our proposed method is able to capture pathologies from images and also allows the assignment of heterogeneous mechanical properties to different structures. In addition, it permits automatic volumetric finite element modelling from medical images with satisfactory surface conformability and also requires less CPU time and computational resources compared to current techniques.
1.1. Marching-cubes algorithm The marching-cubes algorithm [12] is conventionally used in the field of computer graphics for surface extraction from a stack of images. The input image dataset can range from medical images to geological scans. Here we regard voxels as points and define a cube by the neighbouring voxels at the eight vertices. If some voxels of a cube have values (greyscale) less than a user-specified value, the isovalue, while others have values greater than this value, the cube must contribute some component of a surface, correspondingly called the isosurface. From the eight cubic vertex voxels, those that are found to be greater than or equal to the isovalue are considered to be inside the surface, while those that are less than the isovalue are considered to be outside the surface. By determining which edges of the cube are intersected by the isosurface, triangular patches/facets are created that divide the cube into two portions, one within the surface and one outside. These triangular patches are stitched
Computational biomechanical modelling of the lumbar spine together to obtain a surface representation of the object.
2. Methods and materials 2.1. Solution strategy The marching-cubes algorithm can provide a reasonably smooth surface. However, for volumetric filling, a new method is required to generate a voxel mesh that can incorporate with the marching-cubes surface. The input to our modelling method includes a 3D medical dataset and a user-specified value (isovalue) that is used for the generation of the marching-cubes surface. We illustrate the underlying idea of the algorithm in 2D with the example in Fig. 1. The 2D grid represents an image, where black dots are pixels with values greater than or equal to the isovalue, while the small circles are pixels with values less than the isovalue. Based on the marching-cubes algorithm, the 2D surfaces are inserted as solid lines. We define a cell composed of four neighbouring pixels (square) in 2D and eight voxels (cube) in 3D, respectively. The cells are divided into the following three groups in both 2D and 3D: Group 1: Cells completely inside the marchingcubes surface. Group 2: Cells completely outside the marchingcubes surface. Group 3: Cells cut through by the marching-cubes surface. The first group of cells (e.g., the shadowed cells in Fig. 1) are naturally FEM compliant elements in which the vertex voxels/pixels are used directly
Fig. 1 A 2D example of our meshing technique.
27
as nodes. The second group of cells will not be included in the model at all. The third group of cells are the interfaces between internal elements and the marching-cubes surface. Unfortunately, the remaining portions of the group 3 cells cut by the marching-cubes surface are not all FEM compliant. Most of them cannot be used as FE elements directly and therefore need to be further divided into FEM compliant primitives. As seen in the example of Fig. 1, we split the non-compliant group 3 cells into triangles using dashed oblique lines. In 3D, a similar process tetrahedronizes the non-compliant group 3 cells or cubes. Together with the new nodes introduced by the intersection points between marching-cubes surface and the grid edges, a smoothened voxel-based mesh is produced. Since a cube has eight vertices and each vertex has two states, inside and outside, there are 28 = 256 ways a surface can intersect the cube; thus, there are 256 possible patterns that the cube may possess. By enumerating and tetrahedronizing all the 256 cases, the smoothened voxel-based mesh can be achieved with the incorporation of the marching-cubes surface. However, this task is tedious and error-prone. Therefore, by taking rotational symmetry into account, we first consider only the 23 patterns (Fig. 2), which include the initial fifteen patterns proposed by Lorensen and Cline [12] and the eight improved complementary cases proposed by Shoeb [13]. We tetrahedronized the 23-pattern cubes that are cropped by the triangular patches of the marching cubes. In Fig. 2, voxels are located at the cube vertices, and the dashed lines indicate the boundaries between tetrahedrons. Cube vertices with black dots represent voxels inside the marchingcubes surface, while those without black dots are outside the marching-cube surface. Based on the state of the vertex, an eight-bit index containing one bit for each vertex/voxel is created to uniquely identify each of the 256 cases; this index serves as an entry pointer into a lookup table that contains all the tetrahedronization solutions. A tetrahedron generation module was developed to derive the solutions of all the 256 cases from the tetrahedronization solution of the 23 patterns. The algorithm is based on a pattern-rotation method to enumerate the symmetric cases of basic patterns. It can be seen in Fig. 2 that some patterns are naturally tetrahedral (pattern 2) or do not need to be tetrahedronized (patterns 1 and 23), some patterns (4, 5, 7, 8, 10, 11, 12 and 18) are combinations of other simple patterns, and some (3, 6, 9, 13, 14, 15, 16, 17, 19, 20, 21 and 22) have been manually tetrahedronized. With this information, the
28
Z.L. Wang et al.
Fig. 2 Tetrahedronization of the 23 marching-cubes patterns.
module can automatically identify all the symmetric case indices for any given pattern and automatically derive the tetrahedronization solution from either a tetrahedronized pattern or a combination pattern. Apart from the improved performance in implementation, we reduce the likelihood of introducing errors. This meshing technique is sufficiently flexible to be adapted to any variations or improvements due to its full automation. Table 1 shows the number of symmetric cases and the number of tetrahedrons contained for each pattern.
Numbers ranging from 0 to 19 are used to name the nodes in the cubes as shown in Fig. 3, where nodes 0—7 are voxels at cube vertices and others are the new nodes introduced by the intersection of the marching-cubes surface and the cube edges. It can be seen from the above discussion that we have derived a 256-entry lookup table that gives the tetrahedronization solutions for all 256 cases accordingly. In the table, the case index is used as the entry pointer. Part of the table is shown in Table 2 as an illustration. In each entry, the num-
Computational biomechanical modelling of the lumbar spine
Table 1 Symmetric case number and tetrahedron number of 23 patterns Pattern
Symmetric case no.
Tetrahedron no.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
1 8 12 12 4 24 8 24 6 2 24 6 8 12 12 24 24 8 12 12 4 8 1
0 1 3 2 2 6 3 4 6 4 7 6 7 7 7 8 9 8 10 9 10 10 0
29
Table 2
Tetrahedronization lookup table
Index
Tetrahedron number
0 1 2 3
0 1 1 3
.. .
.. .
254
10
255
0
ber of the tetrahedron included is given first, followed by a series of node numbers. Node numbers are divided into group numbers as described above, with each group containing four nodes that form a tetrahedron.
2.2. Smoothening enhancement In the standard marching-cubes algorithm, a linear interpolation algorithm is normally used to calculate the coordinates of the newly inserted nodes on the edges of the cubes. It is possible that such interpolation could produce nodes very close to the cube vertices, and consequently lead to severely
NIL 0, 16, 8, 11 1, 17, 9, 8 0, 11, 16, 9, 1, 17, 9, 16, 1, 0, 16, 9 .. . 6, 2, 5, 4, 5, 8, 7, 6, 7, 5, NIL
3, 11, 6, 16, 16, 7, 3, 11, 6, 16,
3, 3, 6, 7, 8,
1, 5, 8, 11, 8, 7, 16, 11, 7, 3, 5, 8, 11, 16, 7
distorted tetrahedral elements in FEM. However, an intermediate interpolation of the new nodes is able to reduce the distortion significantly at the cost of smoothness. Therefore, a smooth-distortion factor (SDF), denoted by , was developed and implemented in our method as a parameter for distortion controllable meshing. The new node, PC , is inserted according to the equations PC = (PL − C) + C V1 + V2 C= 2
Fig. 3 Node naming.
Tetrahedrons
∈ [0, 1]
(1)
where V1 and V2 are the coordinates of two neighbouring vertices/voxels, C the mid-point of the line joining them and PL is the linear-interpolated point. The degree of smoothness and element distortion are controlled by varying in the range [0, 1]. The extreme values of PC are C =0 PC = (2) PL = 1 Thus, when = 0, PC is the intermediateinterpolated point and the FE mesh will have the lowest distortion; when = 1, the node is the linearinterpolated point and the FE mesh will have the greatest smoothness.
2.3. Segmentation The model development time and accuracy will be reduced when the input images consist solely of the bone of interest. Our region of interest is the vertebral body and the surrounding soft tissue. The goal of segmentation is to separate this region of interest from the surrounding. In CT segmentation, it is not possible to segment out the vertebral body by
30
Z.L. Wang et al.
using a single threshold value due to partial volume effect and noise. A low threshold value is not sufficient to separate bones from the tissue and too high a threshold will misclassify bone regions that have low CT numbers as background. We developed and implemented a 3D adaptive threshold algorithm to provide an accurate segmentation under these conditions. After segmentation, the results can be used as input for marching-cube surface extraction which is then used to build the finite element model of the spine. Details of the adaptive threshold algorithm can be found in Ref. [14].
3. Results A CT scan of a human lumbar L3 vertebra was employed to test our meshing technique. The CT dataset comprised 53 1-mm slices with an imaging matrix of 256 pixel × 256 pixel and an in-plane resolution of 0.664 mm × 0.664 mm. This dataset was thus non-isotropic. The experiment was conducted on a PC with a Pentium IV 2.4GHz processor and 512 MB system memory. Our method was implemented with Visual C++. When an SDF () of 1 was used for full smoothening, all 292,081 mesh elements were generated in less than 6 seconds including the file output time. Table 3 and Fig. 4 compare the error and warning elements associated with each SDF value and selected SDF value, respectively. Note that the error elements are indicated in pink and warning elements in yellow in Fig. 4. In this context, the warning element refers to only a distorted element. There are no error elements and the warning elements are tetrahedral elements found only on the surface of the model. The number of warning elements is small relative to the size of the model— –accounting for less than 10% of the total number Table 3 Number of error and warning elements with different SDF values SDF
Number of error elements
Number of warning elements
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0 0 0 0 0 0 0 0 0 0 0
850 661 639 764 1052 2782 4745 8832 15310 22537 28840
(0.29%) (0.23%) (0.22%) (0.26%) (0.36%) (0.95%) (1.62%) (3.02%) (5.24%) (7.72%) (9.875%)
Fig. 4 Comparison of distorted elements with selected SDF values.
of elements. It is also clear that the number of possibly distorted elements is highest with SDF = 1.0. This number is smallest when SDF = 0.2. Fig. 5a shows the FE model generated with our marching-cube surface smoothened voxel meshing technique; for comparison, a conventional voxel mesh is shown in Fig. 5b. A 3D volume rendered image of the vertebra from the same CT scan dataset is shown in Fig. 5c. The volume rendering technique is capable of displaying both surface and internal details of the actual anatomical structures directly from the 3D images. Close similarity can be observed between our model and the volumetric model. We have cropped the smoothened model to reveal the internal mesh structures. A magnified view of part of the vertebral body is shown in Fig. 5d so that the tetrahedral and hexahedral elements and their distributions can be clearly identified. To further validate the performance of our method in representing the geometry of the structure, we
Computational biomechanical modelling of the lumbar spine
31
Fig. 5 Biomechanical modelling of a human lumbar L3 vertebra (a) FE model generated with marching-cube surface smoothened voxel meshing technique ( = 1), (b) FE model generated with the conventional voxel meshing method, (c) 3D volume rendered image generated directly from the CT scan dataset and (d) magnified view of the vertebra body of the cropped smoothened FE model.
compared the corresponding views of our model and CT images reconstructed using volume rendering in Fig. 6. The observed close similarity is an indication of the superior performance of our method compared to the conventional method in the genTable 4
eration of models conforming to actual geometries of anatomical structures. Subsequently, we compared our mesh to those generated by the commercial software ScanFE whose mesh output is also ABAQUS (Hibbit, Karlsson
Quantitative comparison between ScanFE and the new method
Mesh resolution ScanFE 64 × 64 128 × 128 256 × 256 New method 64 × 64 128 × 128 256 × 256
No. of elements
23688 106829 309964 19477 (17.8%) 98000 (8.3%) 300685 (3.0%)
No. of nodes
8567 45982 154151 6590 (23.1%) 40027 (13.0%) 142847 (7.3%)
CPU time (s)
163 6285 12236 113 (30.7%) 5861 (6.8%) 7533 (38.4%)
Percentage reduction as compared to ScanFE models are illustrated in parentheses.
Memory required (MB)
Disk space required (MB)
66.66 347.29 1170.00
75.83 721.03 3590.00
58.60 (12.1%) 313.51 (9.7%) 1040.00 (11.1%)
55.76 (26.5%) 610.47 (15.3%) 3360.00 (6.4%)
32
Z.L. Wang et al.
Fig. 6 Comparison of corresponding views of marching-cube surface smoothened FE model and 3D reconstructed CT images.
and Sorensen, Inc., USA) compliant. Using the same dataset as mentioned above, we sub-sampled the in-plane resolution to 128 × 128 and 64 × 64 to better evaluate the meshes. All three datasets were used as inputs to our technique and ScanFE. The FE meshes generated were subjected to the same analysis which was chosen to simulate experimental testing to determine vertebral compression fractures. With the inferior end of the vertebral body fixed, the superior end was vertically displaced by 1 mm downwards (Fig. 7). Our mesh evaluation is based on the comparison of the meshes generated from our method and ScanFE. The following characteristics are considered: visual inspection, continuity of stress distribution (Fig. 8), number of nodes/elements, CPU
Fig. 7 Loading and boundary conditions for the finite element analysis chosen for mesh evaluation.
Computational biomechanical modelling of the lumbar spine
33
Fig. 8 Visual Inspection (left column) and Von Mises stress distribution (right column) of meshes generated using ScanFE and our new method.
time, physical memory requirements and disk space requirements (Table 4). In general, all the models exhibit a similar stress pattern which extends from the anterior-superior region diagonally down towards the posteriorinferior region of the vertebral body (Fig. 8). From closer visual inspection, we can see that the meshes from ScanFE consist of both hexahedral and tetrahedral elements at the surface, which we expect to lead to undesirable results. This is because, mathematically, tetrahedral elements are stiffer than hexahedral elements. This difference in stiffness
across adjacent finite elements will lead to undesirable stress discontinuity. This is demonstrated clearly in the stress distribution of models from ScanFE, and is most obvious in the 64 × 64 dataset (Fig. 8a, left). The discontinuity effects are less severe at increasing dataset resolutions (Fig. 8c and e). Meshes from our technique do not exhibit these discontinuities even in the 64 × 64 dataset (Fig. 8a, right). In terms of the number of elements, number of nodes, CPU time required, memory and disk space required, the models generated using our technique
34 performed better (Table 4). For 64 × 64, 128 × 128 and the 256 × 256 resolutions, our method resulted in a 17.8%, 8.3% and 3.0% reduction, respectively, in the number of elements. This contributed to the reduction in CPU time required to run the model analyses regardless of the dataset resolution. The computation time was reduced by 30.7% for the 64 × 64 dataset, 6.8% for the 128 × 128 dataset, and by 38.4% for the 256 × 256 dataset. On average, about 11% less memory and 16% less disk space were required by our method.
4. Discussion and conclusion We have proposed a new method for the automatic generation of patient-specific volumetric voxelbased finite element model from medical images. Such a model appears to conform well to actual geometries of anatomical structures. The method comprises of a novel tetrahedronization scheme incorporating with marching-cubes surface. The method is further enhanced by the introduction of a smooth-distortion factor for distortion controllable meshing. From our testing, this new meshing technique presents superior performance over existing methods in terms of both geometry conformability and FEA performance. The innovative tetrahedronization scheme that incorporates a marching-cube surface is different from existing voxel meshing techniques. We take advantage of the robustness of the marching-cube method and its resultant smooth surface model to guide the meshing process. The internal hexahedral elements can be tetrahedronized if an all tetrahedral element FE model is desired. Similarly, the tetrahedral elements on the outer surface can be combined to form hexahedral elements. The FEA experiments described in the previous section confirmed that the combination of tetrahedral and hexahedral elements gives good results. Nevertheless, it is important to compare the performance using a model comprising all-hexahedral, all-tetrahedral and mixed element types. This is not within the scope of this paper. Distortion controllable meshing using a smoothdistortion factor is another feature of this method. In the previous section, we have shown that when SDF = 1, the number of possible distorted elements is almost 45 times of that when SDF = 0.2. Our study also reveals that the optimal SDF is not necessary the smallest possible number. We demonstrated the ability of our meshing method to conform to the actual geometry of the anatomical structure by comparing the corresponding rendered views of FE model and 3D recon-
Z.L. Wang et al. structed CT images. They are clearly very similar. A more quantitative approach to assess the similarity would be desirable; we are currently investigating the possibilities of using regression techniques to determine the differences quantitatively. We use extracted surfaces from the conventional meshing cube method as the input to our voxelbased mesh generation algorithm. However, our method is certainly not restricted to marching-cube surfaces. There are newer methods that will enable interactive geometry triangular surface meshing [15,16]. These methods in the computer graphics community have produced smoother and more efficient triangular meshes. They are likely to improve the accuracy of our FE meshing method if adopted. There are limitations with our methods. Being a voxel-based method, the resultant number of finite elements is significant. Hence, the method is best for computational application, such as simulation of vertebroplasty where a bone needle is inserted into a vertebra on lumbar spine. However, this FE meshing method is applicable to the whole spine. We plan to further evaluate our method by comparing the results obtained from finite element analysis to experimental studies. More accurate analysis will also be done with our models with the assignment of the patient-specific heterogeneous mechanical properties.
Acknowledgement The project is funded by Virtual Spine Workstation grant (R-265-000-184-305) from the Agency for Science, Technology and Research (A*STAR), Singapore.
References [1] ABAQUS Version 6.4 Analysis User’s Manual Vol IV: Elements (2004). [2] D.P. Fyhrie, M.S. Hamid, R.F. Kuo, S.M. Lang, Direct threedimensional finite element analysis of human vertebral cancellous bone, Trans. Ann. Meeting Orthop. Res. Soc. 17 (1992) 551. [3] B. Rietbergen, H. Weinans, R. Huiskes, A. Odgaard, A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models, J. Biomech. 28 (1995) 69—81. [4] J.C.M. Teo, Z.L. Wang, S.H. Teoh, Effects of geometric smoothening to voxel-based finite element models, in: Proceedings of the Scientific Meeting of Biomedical Enginnering Society of Singapore, 2004. [5] R.P. Crawford, C.E. Cann, T.M. Keaveny, Finite element models predict in vitro vertebral body compressive strength better than quantitative computed tomography, Bone 33 (4) (2003) 744—750.
Computational biomechanical modelling of the lumbar spine [6] E. Pistoia, B. van Rietbergen, E.M. Lochmuller, A.C. Lill, F. Eckstein, P. Rugsegger, Estimation of distal radius failure load with micro-finite element analysis models based on three-dimensional peripheral quantitative computed tomography images, Bone 30 (6) (2002) 842—848. [7] D.D. Cody, G.J. Gross, F.J. Hou, H.J. Spencer, S.A. Goldstein, D.P. Fyhrie, Femoral strength is better predicted by finite element models than QCT and DXA, J. Biomech. 32 (1999) 1013—1020. [8] J.H. Keyak, J.M. Meagher, H.B. Skinner, C.D. Mote, Automated three-dimensional finite element modelling of bone: a new method, J. Biomed. Eng. 12 (1990) 389—797. [9] N.M. Grosland, D.R. Pedersen, T.D. Brown, An anatomical voxel-based FE contact analysis formulation, in: Proceedings of the Annual Meeting of the American Society of Biomechanics, 2001. [10] N.M. Grosland, T.D. Brown, A voxel-based formulation for contact finite element analysis, Comput. Methods Biomech. Biomed. Eng. 5 (2002) 21—32.
35
[11] D.L. Camacho, R.H. Hopper, G.M. Lin, B.S. Myers, An improved method for finite element mesh generation of geometrically complex structures with application to the skullbase, J. Biomech. 30 (1997) 1067— 1070. [12] W.E. Lorensen, H.E. Cline, Marching cubes: a high resolution 3D surface construction algorithm, Comput. Graph. 21 (4) (1987) 163—169. [13] Shoeb, Improved Marching Cubes, http://enuxsa.eas.asu. edu/∼shoeb/graphics/improved.html (1998). [14] C.H. Yan, S.H. Ong, Y. Ge, S.H. Teoh, C.K. Chui, Accurate and precise 3D surface based registration, in: Proceedings of the IASTED Signal and Image Processing’04, 2004. [15] P. Alliez, M. Meyer, M. Desbrun, Interactive geometry remeshing, in: Proceedings of the ACM SIGGRAPH, 2002, pp. 347—354. [16] I. Guskov, K. Vidimce, W. Sweldens, P. Schroder, Normal meshes, in: Proceedings of the ACM SIGGRAPH, 2000, pp. 95—102.
Computational biomechanical modelling of the lumbar spine using marching-cubes surface smoothened finite element voxel meshing Z.L. Wang a, J.C.M. Teo a, C.K. Chui a, S.H. Ong b,c, C.H. Yan b, S.C. Wang d, H.K. Wong e, S.H. Teoh a,∗ a
Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, Singapore 117576, Singapore b Department of Electrical & Computer Engineering, National University of Singapore, Singapore c Division of Bioengineering, National University of Singapore, Singapore d Department of Diagnostic Radiology, National University Hospital, Singapore e Department of Orthopaedic Surgery, National University Hospital, Singapore Received 14 January 2005 ; received in revised form 7 June 2005; accepted 7 June 2005 KEYWORDS Modelling; Finite element method (FEM); Voxel-based mesh; Marching-cubes; Volumetric mesh
Summary There is a need for the development of finite element (FE) models based on medical datasets, such as magnetic resonance imaging and computerized tomography in computation biomechanics. Direct conversion of graphic voxels to FE elements is a commonly used method for the generation of FE models. However, conventional voxel-based methods tend to produce models with jagged surfaces. This is a consequence of the inherent characteristics of voxel elements; such a model is unable to capture the geometries of anatomical structures satisfactorily. We have developed a robust technique for the automatic generation of voxel-based patientspecific FE models. Our approach features a novel tetrahedronization scheme that incorporates marching-cubes surface smoothing together with a smooth-distortion factor (SDF). The models conform to the actual geometries of anatomical structures of a lumbar spine segment (L3). The resultant finite element analysis (FEA) at the surfaces is more accurate compared to the use of conventional voxel-based generated FE models. In general, models produced by our method were superior compared to that obtained using the commercial software ScanFE. © 2005 Elsevier Ireland Ltd. All rights reserved.
1. Introduction ∗ Corresponding author. Tel.: +65 6874 6345; fax: +65 6777 3537. E-mail address: [email protected] (S.H. Teoh).
There is a current trend towards the use of computational biomechanics in clinical studies. Prior to biomechanical computation, there is a need for the development of finite element (FE) models based
0169-2607/$ — see front matter © 2005 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2005.06.006
26 on medical datasets, such as magnetic resonance imaging (MRI) and computerized tomography (CT). After the domain boundaries have been defined for automatic finite element meshing of anatomical structures, algorithms are available to fill the domain with tetrahedral elements or hexahedral elements. Tetrahedral elements are used most frequently as this element type is geometrically versatile and can conform to the complex boundaries of anatomical structures. Linear tetrahedral elements usually have the problem of being overly stiff due to their ‘‘odd’’ shape relative to the hexahedral element in theory [1], but it is unclear from published literature that hexahedral elements are better than tetrahedral elements in practice. Tetrahedral elements have the advantage of being able to better model the complex shape of the object. The voxel-based FE meshing method directly converts computer images from 3D voxels to eightnode hexahedral finite elements and fills the domain with hexahedral elements. Despite the intuitive value of such a direct voxel-based method of producing FE meshing models, this method is unable to capture surface boundaries satisfactorily. Jagged or digitized surfaces are unavoidable in voxel-based FE models, which in turn leads to inaccurate finite element analysis (FEA) results. They, however, are still employed for computational biomechanical investigations [2—4]. One reason is that voxel-based modelling also has the ability to automatically assign orthotropic bone mechanical properties based on greyscale intensities [5], thus allowing heterogeneous properties throughout the FE model of the anatomical structure. Voxel-based FE models have been used in biomechanical studies of human distal radius [6], femur [7] and vertebral body [5,8] with some success. Current software that provides computational biomechanical modelling based on medical datasets includes ScanFE and Mimics. ScanFE (Simpleware Ltd., Exeter, UK) has the capability to manipulate medical datasets and obtain FE models automatically. It also produces a voxel-based mesh of the interior together with a tetrahedral surface. Mimics (Materialise, Leuven, Belgium) can produce a surface mesh, which then has to be exported into other FE pre-processors for the generation of a volumetric all tetrahedral FE mesh. Because of the surface irregularities inherent in the models produced by conventional voxel-based methods, Grosland et al. [9,10] proposed a technique to smoothen the irregularities for the FE contact stress analysis of articular joints. Their method generates smooth surfaces at plausible con-
Z.L. Wang et al. tact regions with local refinement. However, due to the use of multiple point constraints (MPCs), the object is modelled in two resolutions and the noncontact regions with jagged surfaces still tend to produce inaccurate FEA results. Camacho et al. [11] also proposed a smoothening algorithm to eliminate the sharp edges of voxel models at surfaces and at the interfaces between materials without the introduction of new elements. This method was extended and implemented in three dimensions (3D) for the FE study of a vertebral body [4]. However, the reported performance of the finite element analysis is not significantly better than conventional voxel-based methods. We have developed a new approach that overcomes the above mentioned deficiencies. The marching-cubes method is employed to extract the smooth surfaces of anatomical structures and this followed by the use of a novel tetrahedronization scheme to generate a voxel-based mesh. Hence, the outer layer in our mesh model is composed of tetrahedral elements while the internal elements are all hexahedral. The two layers share the same set of nodes and are connected seamlessly, thus maintaining FE compliancy. Our proposed method is able to capture pathologies from images and also allows the assignment of heterogeneous mechanical properties to different structures. In addition, it permits automatic volumetric finite element modelling from medical images with satisfactory surface conformability and also requires less CPU time and computational resources compared to current techniques.
1.1. Marching-cubes algorithm The marching-cubes algorithm [12] is conventionally used in the field of computer graphics for surface extraction from a stack of images. The input image dataset can range from medical images to geological scans. Here we regard voxels as points and define a cube by the neighbouring voxels at the eight vertices. If some voxels of a cube have values (greyscale) less than a user-specified value, the isovalue, while others have values greater than this value, the cube must contribute some component of a surface, correspondingly called the isosurface. From the eight cubic vertex voxels, those that are found to be greater than or equal to the isovalue are considered to be inside the surface, while those that are less than the isovalue are considered to be outside the surface. By determining which edges of the cube are intersected by the isosurface, triangular patches/facets are created that divide the cube into two portions, one within the surface and one outside. These triangular patches are stitched
Computational biomechanical modelling of the lumbar spine together to obtain a surface representation of the object.
2. Methods and materials 2.1. Solution strategy The marching-cubes algorithm can provide a reasonably smooth surface. However, for volumetric filling, a new method is required to generate a voxel mesh that can incorporate with the marching-cubes surface. The input to our modelling method includes a 3D medical dataset and a user-specified value (isovalue) that is used for the generation of the marching-cubes surface. We illustrate the underlying idea of the algorithm in 2D with the example in Fig. 1. The 2D grid represents an image, where black dots are pixels with values greater than or equal to the isovalue, while the small circles are pixels with values less than the isovalue. Based on the marching-cubes algorithm, the 2D surfaces are inserted as solid lines. We define a cell composed of four neighbouring pixels (square) in 2D and eight voxels (cube) in 3D, respectively. The cells are divided into the following three groups in both 2D and 3D: Group 1: Cells completely inside the marchingcubes surface. Group 2: Cells completely outside the marchingcubes surface. Group 3: Cells cut through by the marching-cubes surface. The first group of cells (e.g., the shadowed cells in Fig. 1) are naturally FEM compliant elements in which the vertex voxels/pixels are used directly
Fig. 1 A 2D example of our meshing technique.
27
as nodes. The second group of cells will not be included in the model at all. The third group of cells are the interfaces between internal elements and the marching-cubes surface. Unfortunately, the remaining portions of the group 3 cells cut by the marching-cubes surface are not all FEM compliant. Most of them cannot be used as FE elements directly and therefore need to be further divided into FEM compliant primitives. As seen in the example of Fig. 1, we split the non-compliant group 3 cells into triangles using dashed oblique lines. In 3D, a similar process tetrahedronizes the non-compliant group 3 cells or cubes. Together with the new nodes introduced by the intersection points between marching-cubes surface and the grid edges, a smoothened voxel-based mesh is produced. Since a cube has eight vertices and each vertex has two states, inside and outside, there are 28 = 256 ways a surface can intersect the cube; thus, there are 256 possible patterns that the cube may possess. By enumerating and tetrahedronizing all the 256 cases, the smoothened voxel-based mesh can be achieved with the incorporation of the marching-cubes surface. However, this task is tedious and error-prone. Therefore, by taking rotational symmetry into account, we first consider only the 23 patterns (Fig. 2), which include the initial fifteen patterns proposed by Lorensen and Cline [12] and the eight improved complementary cases proposed by Shoeb [13]. We tetrahedronized the 23-pattern cubes that are cropped by the triangular patches of the marching cubes. In Fig. 2, voxels are located at the cube vertices, and the dashed lines indicate the boundaries between tetrahedrons. Cube vertices with black dots represent voxels inside the marchingcubes surface, while those without black dots are outside the marching-cube surface. Based on the state of the vertex, an eight-bit index containing one bit for each vertex/voxel is created to uniquely identify each of the 256 cases; this index serves as an entry pointer into a lookup table that contains all the tetrahedronization solutions. A tetrahedron generation module was developed to derive the solutions of all the 256 cases from the tetrahedronization solution of the 23 patterns. The algorithm is based on a pattern-rotation method to enumerate the symmetric cases of basic patterns. It can be seen in Fig. 2 that some patterns are naturally tetrahedral (pattern 2) or do not need to be tetrahedronized (patterns 1 and 23), some patterns (4, 5, 7, 8, 10, 11, 12 and 18) are combinations of other simple patterns, and some (3, 6, 9, 13, 14, 15, 16, 17, 19, 20, 21 and 22) have been manually tetrahedronized. With this information, the
28
Z.L. Wang et al.
Fig. 2 Tetrahedronization of the 23 marching-cubes patterns.
module can automatically identify all the symmetric case indices for any given pattern and automatically derive the tetrahedronization solution from either a tetrahedronized pattern or a combination pattern. Apart from the improved performance in implementation, we reduce the likelihood of introducing errors. This meshing technique is sufficiently flexible to be adapted to any variations or improvements due to its full automation. Table 1 shows the number of symmetric cases and the number of tetrahedrons contained for each pattern.
Numbers ranging from 0 to 19 are used to name the nodes in the cubes as shown in Fig. 3, where nodes 0—7 are voxels at cube vertices and others are the new nodes introduced by the intersection of the marching-cubes surface and the cube edges. It can be seen from the above discussion that we have derived a 256-entry lookup table that gives the tetrahedronization solutions for all 256 cases accordingly. In the table, the case index is used as the entry pointer. Part of the table is shown in Table 2 as an illustration. In each entry, the num-
Computational biomechanical modelling of the lumbar spine
Table 1 Symmetric case number and tetrahedron number of 23 patterns Pattern
Symmetric case no.
Tetrahedron no.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
1 8 12 12 4 24 8 24 6 2 24 6 8 12 12 24 24 8 12 12 4 8 1
0 1 3 2 2 6 3 4 6 4 7 6 7 7 7 8 9 8 10 9 10 10 0
29
Table 2
Tetrahedronization lookup table
Index
Tetrahedron number
0 1 2 3
0 1 1 3
.. .
.. .
254
10
255
0
ber of the tetrahedron included is given first, followed by a series of node numbers. Node numbers are divided into group numbers as described above, with each group containing four nodes that form a tetrahedron.
2.2. Smoothening enhancement In the standard marching-cubes algorithm, a linear interpolation algorithm is normally used to calculate the coordinates of the newly inserted nodes on the edges of the cubes. It is possible that such interpolation could produce nodes very close to the cube vertices, and consequently lead to severely
NIL 0, 16, 8, 11 1, 17, 9, 8 0, 11, 16, 9, 1, 17, 9, 16, 1, 0, 16, 9 .. . 6, 2, 5, 4, 5, 8, 7, 6, 7, 5, NIL
3, 11, 6, 16, 16, 7, 3, 11, 6, 16,
3, 3, 6, 7, 8,
1, 5, 8, 11, 8, 7, 16, 11, 7, 3, 5, 8, 11, 16, 7
distorted tetrahedral elements in FEM. However, an intermediate interpolation of the new nodes is able to reduce the distortion significantly at the cost of smoothness. Therefore, a smooth-distortion factor (SDF), denoted by , was developed and implemented in our method as a parameter for distortion controllable meshing. The new node, PC , is inserted according to the equations PC = (PL − C) + C V1 + V2 C= 2
Fig. 3 Node naming.
Tetrahedrons
∈ [0, 1]
(1)
where V1 and V2 are the coordinates of two neighbouring vertices/voxels, C the mid-point of the line joining them and PL is the linear-interpolated point. The degree of smoothness and element distortion are controlled by varying in the range [0, 1]. The extreme values of PC are C =0 PC = (2) PL = 1 Thus, when = 0, PC is the intermediateinterpolated point and the FE mesh will have the lowest distortion; when = 1, the node is the linearinterpolated point and the FE mesh will have the greatest smoothness.
2.3. Segmentation The model development time and accuracy will be reduced when the input images consist solely of the bone of interest. Our region of interest is the vertebral body and the surrounding soft tissue. The goal of segmentation is to separate this region of interest from the surrounding. In CT segmentation, it is not possible to segment out the vertebral body by
30
Z.L. Wang et al.
using a single threshold value due to partial volume effect and noise. A low threshold value is not sufficient to separate bones from the tissue and too high a threshold will misclassify bone regions that have low CT numbers as background. We developed and implemented a 3D adaptive threshold algorithm to provide an accurate segmentation under these conditions. After segmentation, the results can be used as input for marching-cube surface extraction which is then used to build the finite element model of the spine. Details of the adaptive threshold algorithm can be found in Ref. [14].
3. Results A CT scan of a human lumbar L3 vertebra was employed to test our meshing technique. The CT dataset comprised 53 1-mm slices with an imaging matrix of 256 pixel × 256 pixel and an in-plane resolution of 0.664 mm × 0.664 mm. This dataset was thus non-isotropic. The experiment was conducted on a PC with a Pentium IV 2.4GHz processor and 512 MB system memory. Our method was implemented with Visual C++. When an SDF () of 1 was used for full smoothening, all 292,081 mesh elements were generated in less than 6 seconds including the file output time. Table 3 and Fig. 4 compare the error and warning elements associated with each SDF value and selected SDF value, respectively. Note that the error elements are indicated in pink and warning elements in yellow in Fig. 4. In this context, the warning element refers to only a distorted element. There are no error elements and the warning elements are tetrahedral elements found only on the surface of the model. The number of warning elements is small relative to the size of the model— –accounting for less than 10% of the total number Table 3 Number of error and warning elements with different SDF values SDF
Number of error elements
Number of warning elements
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0 0 0 0 0 0 0 0 0 0 0
850 661 639 764 1052 2782 4745 8832 15310 22537 28840
(0.29%) (0.23%) (0.22%) (0.26%) (0.36%) (0.95%) (1.62%) (3.02%) (5.24%) (7.72%) (9.875%)
Fig. 4 Comparison of distorted elements with selected SDF values.
of elements. It is also clear that the number of possibly distorted elements is highest with SDF = 1.0. This number is smallest when SDF = 0.2. Fig. 5a shows the FE model generated with our marching-cube surface smoothened voxel meshing technique; for comparison, a conventional voxel mesh is shown in Fig. 5b. A 3D volume rendered image of the vertebra from the same CT scan dataset is shown in Fig. 5c. The volume rendering technique is capable of displaying both surface and internal details of the actual anatomical structures directly from the 3D images. Close similarity can be observed between our model and the volumetric model. We have cropped the smoothened model to reveal the internal mesh structures. A magnified view of part of the vertebral body is shown in Fig. 5d so that the tetrahedral and hexahedral elements and their distributions can be clearly identified. To further validate the performance of our method in representing the geometry of the structure, we
Computational biomechanical modelling of the lumbar spine
31
Fig. 5 Biomechanical modelling of a human lumbar L3 vertebra (a) FE model generated with marching-cube surface smoothened voxel meshing technique ( = 1), (b) FE model generated with the conventional voxel meshing method, (c) 3D volume rendered image generated directly from the CT scan dataset and (d) magnified view of the vertebra body of the cropped smoothened FE model.
compared the corresponding views of our model and CT images reconstructed using volume rendering in Fig. 6. The observed close similarity is an indication of the superior performance of our method compared to the conventional method in the genTable 4
eration of models conforming to actual geometries of anatomical structures. Subsequently, we compared our mesh to those generated by the commercial software ScanFE whose mesh output is also ABAQUS (Hibbit, Karlsson
Quantitative comparison between ScanFE and the new method
Mesh resolution ScanFE 64 × 64 128 × 128 256 × 256 New method 64 × 64 128 × 128 256 × 256
No. of elements
23688 106829 309964 19477 (17.8%) 98000 (8.3%) 300685 (3.0%)
No. of nodes
8567 45982 154151 6590 (23.1%) 40027 (13.0%) 142847 (7.3%)
CPU time (s)
163 6285 12236 113 (30.7%) 5861 (6.8%) 7533 (38.4%)
Percentage reduction as compared to ScanFE models are illustrated in parentheses.
Memory required (MB)
Disk space required (MB)
66.66 347.29 1170.00
75.83 721.03 3590.00
58.60 (12.1%) 313.51 (9.7%) 1040.00 (11.1%)
55.76 (26.5%) 610.47 (15.3%) 3360.00 (6.4%)
32
Z.L. Wang et al.
Fig. 6 Comparison of corresponding views of marching-cube surface smoothened FE model and 3D reconstructed CT images.
and Sorensen, Inc., USA) compliant. Using the same dataset as mentioned above, we sub-sampled the in-plane resolution to 128 × 128 and 64 × 64 to better evaluate the meshes. All three datasets were used as inputs to our technique and ScanFE. The FE meshes generated were subjected to the same analysis which was chosen to simulate experimental testing to determine vertebral compression fractures. With the inferior end of the vertebral body fixed, the superior end was vertically displaced by 1 mm downwards (Fig. 7). Our mesh evaluation is based on the comparison of the meshes generated from our method and ScanFE. The following characteristics are considered: visual inspection, continuity of stress distribution (Fig. 8), number of nodes/elements, CPU
Fig. 7 Loading and boundary conditions for the finite element analysis chosen for mesh evaluation.
Computational biomechanical modelling of the lumbar spine
33
Fig. 8 Visual Inspection (left column) and Von Mises stress distribution (right column) of meshes generated using ScanFE and our new method.
time, physical memory requirements and disk space requirements (Table 4). In general, all the models exhibit a similar stress pattern which extends from the anterior-superior region diagonally down towards the posteriorinferior region of the vertebral body (Fig. 8). From closer visual inspection, we can see that the meshes from ScanFE consist of both hexahedral and tetrahedral elements at the surface, which we expect to lead to undesirable results. This is because, mathematically, tetrahedral elements are stiffer than hexahedral elements. This difference in stiffness
across adjacent finite elements will lead to undesirable stress discontinuity. This is demonstrated clearly in the stress distribution of models from ScanFE, and is most obvious in the 64 × 64 dataset (Fig. 8a, left). The discontinuity effects are less severe at increasing dataset resolutions (Fig. 8c and e). Meshes from our technique do not exhibit these discontinuities even in the 64 × 64 dataset (Fig. 8a, right). In terms of the number of elements, number of nodes, CPU time required, memory and disk space required, the models generated using our technique
34 performed better (Table 4). For 64 × 64, 128 × 128 and the 256 × 256 resolutions, our method resulted in a 17.8%, 8.3% and 3.0% reduction, respectively, in the number of elements. This contributed to the reduction in CPU time required to run the model analyses regardless of the dataset resolution. The computation time was reduced by 30.7% for the 64 × 64 dataset, 6.8% for the 128 × 128 dataset, and by 38.4% for the 256 × 256 dataset. On average, about 11% less memory and 16% less disk space were required by our method.
4. Discussion and conclusion We have proposed a new method for the automatic generation of patient-specific volumetric voxelbased finite element model from medical images. Such a model appears to conform well to actual geometries of anatomical structures. The method comprises of a novel tetrahedronization scheme incorporating with marching-cubes surface. The method is further enhanced by the introduction of a smooth-distortion factor for distortion controllable meshing. From our testing, this new meshing technique presents superior performance over existing methods in terms of both geometry conformability and FEA performance. The innovative tetrahedronization scheme that incorporates a marching-cube surface is different from existing voxel meshing techniques. We take advantage of the robustness of the marching-cube method and its resultant smooth surface model to guide the meshing process. The internal hexahedral elements can be tetrahedronized if an all tetrahedral element FE model is desired. Similarly, the tetrahedral elements on the outer surface can be combined to form hexahedral elements. The FEA experiments described in the previous section confirmed that the combination of tetrahedral and hexahedral elements gives good results. Nevertheless, it is important to compare the performance using a model comprising all-hexahedral, all-tetrahedral and mixed element types. This is not within the scope of this paper. Distortion controllable meshing using a smoothdistortion factor is another feature of this method. In the previous section, we have shown that when SDF = 1, the number of possible distorted elements is almost 45 times of that when SDF = 0.2. Our study also reveals that the optimal SDF is not necessary the smallest possible number. We demonstrated the ability of our meshing method to conform to the actual geometry of the anatomical structure by comparing the corresponding rendered views of FE model and 3D recon-
Z.L. Wang et al. structed CT images. They are clearly very similar. A more quantitative approach to assess the similarity would be desirable; we are currently investigating the possibilities of using regression techniques to determine the differences quantitatively. We use extracted surfaces from the conventional meshing cube method as the input to our voxelbased mesh generation algorithm. However, our method is certainly not restricted to marching-cube surfaces. There are newer methods that will enable interactive geometry triangular surface meshing [15,16]. These methods in the computer graphics community have produced smoother and more efficient triangular meshes. They are likely to improve the accuracy of our FE meshing method if adopted. There are limitations with our methods. Being a voxel-based method, the resultant number of finite elements is significant. Hence, the method is best for computational application, such as simulation of vertebroplasty where a bone needle is inserted into a vertebra on lumbar spine. However, this FE meshing method is applicable to the whole spine. We plan to further evaluate our method by comparing the results obtained from finite element analysis to experimental studies. More accurate analysis will also be done with our models with the assignment of the patient-specific heterogeneous mechanical properties.
Acknowledgement The project is funded by Virtual Spine Workstation grant (R-265-000-184-305) from the Agency for Science, Technology and Research (A*STAR), Singapore.
References [1] ABAQUS Version 6.4 Analysis User’s Manual Vol IV: Elements (2004). [2] D.P. Fyhrie, M.S. Hamid, R.F. Kuo, S.M. Lang, Direct threedimensional finite element analysis of human vertebral cancellous bone, Trans. Ann. Meeting Orthop. Res. Soc. 17 (1992) 551. [3] B. Rietbergen, H. Weinans, R. Huiskes, A. Odgaard, A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models, J. Biomech. 28 (1995) 69—81. [4] J.C.M. Teo, Z.L. Wang, S.H. Teoh, Effects of geometric smoothening to voxel-based finite element models, in: Proceedings of the Scientific Meeting of Biomedical Enginnering Society of Singapore, 2004. [5] R.P. Crawford, C.E. Cann, T.M. Keaveny, Finite element models predict in vitro vertebral body compressive strength better than quantitative computed tomography, Bone 33 (4) (2003) 744—750.
Computational biomechanical modelling of the lumbar spine [6] E. Pistoia, B. van Rietbergen, E.M. Lochmuller, A.C. Lill, F. Eckstein, P. Rugsegger, Estimation of distal radius failure load with micro-finite element analysis models based on three-dimensional peripheral quantitative computed tomography images, Bone 30 (6) (2002) 842—848. [7] D.D. Cody, G.J. Gross, F.J. Hou, H.J. Spencer, S.A. Goldstein, D.P. Fyhrie, Femoral strength is better predicted by finite element models than QCT and DXA, J. Biomech. 32 (1999) 1013—1020. [8] J.H. Keyak, J.M. Meagher, H.B. Skinner, C.D. Mote, Automated three-dimensional finite element modelling of bone: a new method, J. Biomed. Eng. 12 (1990) 389—797. [9] N.M. Grosland, D.R. Pedersen, T.D. Brown, An anatomical voxel-based FE contact analysis formulation, in: Proceedings of the Annual Meeting of the American Society of Biomechanics, 2001. [10] N.M. Grosland, T.D. Brown, A voxel-based formulation for contact finite element analysis, Comput. Methods Biomech. Biomed. Eng. 5 (2002) 21—32.
35
[11] D.L. Camacho, R.H. Hopper, G.M. Lin, B.S. Myers, An improved method for finite element mesh generation of geometrically complex structures with application to the skullbase, J. Biomech. 30 (1997) 1067— 1070. [12] W.E. Lorensen, H.E. Cline, Marching cubes: a high resolution 3D surface construction algorithm, Comput. Graph. 21 (4) (1987) 163—169. [13] Shoeb, Improved Marching Cubes, http://enuxsa.eas.asu. edu/∼shoeb/graphics/improved.html (1998). [14] C.H. Yan, S.H. Ong, Y. Ge, S.H. Teoh, C.K. Chui, Accurate and precise 3D surface based registration, in: Proceedings of the IASTED Signal and Image Processing’04, 2004. [15] P. Alliez, M. Meyer, M. Desbrun, Interactive geometry remeshing, in: Proceedings of the ACM SIGGRAPH, 2002, pp. 347—354. [16] I. Guskov, K. Vidimce, W. Sweldens, P. Schroder, Normal meshes, in: Proceedings of the ACM SIGGRAPH, 2000, pp. 95—102.