Bio-Medical Materials and Engineering 24 (2014) 1469–1484 DOI 10.3233/BME-130951 IOS Press
1469
Voxel-based approach to generate entire human metacarpal bone with microscopic architecture for finite element analysis C.Y. Tang a,∗ , C.P. Tsui a , Y.M. Tang a , L. Wei a , C.T. Wong a , K.W. Lam b , W.Y. Ip b , W.W.J. Lu b and M.Y.C. Pang c a
Department of Industrial and Systems Engineering, The Hong Kong Polytechnic University, Hong Kong, China b Department of Orthopaedics and Traumatology, Li Ka Shing Faculty of Medicine, The University of Hong Kong, Hong Kong, China c Department of Rehabilitation Sciences, The Hong Kong Polytechnic University, Hong Kong, China Received 5 June 2012 Accepted 31 July 2013 Abstract. With the development of micro-computed tomography (micro-CT) technology, it is possible to construct threedimensional (3D) models of human bone without destruction of samples and predict mechanical behavior of bone using finite element analysis (FEA). However, due to large number of elements required for constructing the FE models of entire bone, this demands a substantial computational effort and the analysis usually needs a high level of computer. In this article, a voxel-based approach for generation of FE models of entire bone with microscopic architecture from micro-CT image data is proposed. To enable the FE analyses of entire bone to be run even on a general personal computer, grayscale intensity thresholds were adopted to reduce the amount of elements. Human metacarpal bone (MCP) bone was used as an example for demonstrating the applicability of the proposed method. The micro-CT images of the MCP bone were combined and converted into 3D array of pixels. Dual grayscale intensity threshold parameters were used to distinguish the pixels of bone tissues from those of surrounding soft tissues and improve predictive accuracy for the FE analyses with different sizes of elements. The method of selecting an appropriate value of the second grayscale intensity threshold was also suggested to minimize the area error for the reconstructed cross-sections of a FE structure. Experimental results showed that the entire FE MCP bone with microscopic architecture could be modeled and analyzed on a personal computer with reasonable accuracy. Keywords: Voxel-based, finite element model, micro-computed tomography (micro-CT), microscopic architecture, human metacarpal bone, grayscale intensity threshold
1. Introduction The Finite Element method (FEM) has been shown to be a powerful tool for determining the solution of the physical mechanical behavior of the human skeleton in the field of orthopedic biomechanics [1,2]. However, bone is an osseous tissue with very complex organization consisting of an intricate cortical * Address for correspondence: C.Y. Tang, Department of Industrial and Systems Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, China. Tel.: +852 27666608; Fax: +852 23625267; E-mail: [email protected].
0959-2989/14/$27.50 © 2014 – IOS Press and the authors. All rights reserved
1470
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
layer and a randomized trabecular pattern of internal microstructures. The difficulties of capturing fine detailed geometrical data and generation of three-dimensional (3D) finite element (FE) meshes complicate the implementation of finite element analysis [3]. With the recent advancement in computational power and medical imaging techniques, FEM is now able to nondestructively model stress, strain and deformation in complex 3D bony structures [4]. 3D FE models are usually generated from computed tomography (CT) scan datasets through various approaches [5–7], because the volume and surface information of either a whole bone or bone segments can be acquired from the CT scan. Based on anatomically realistic bone geometry obtained from a CT scan dataset, Bandak [5] developed a 3D FE model of the human lower extremity for studying mechanisms of impact injury to the lower leg. The method requires the CT scan dataset to be pre-processed to extract the geometry of the bone and through creation of inner and outer surfaces of the bone by using spline interpolation between the virtual data points from the CT image. A grayscale intensity threshold was used to isolate the bone region from soft tissues and the inner medullary surface [6,8]. Another commonly used method for creation of the CT-based 3D bone model is to firstly segment bone and inner medullary surfaces into triangles. The segmented results was saved in the Stereo-Lithography Interface Format (STL) for screen visualization and rapid prototyping [6,9]. The surface meshes in STL format after mesh optimization could be transformed into volume meshes using software such as MIMICS and 3D Slicer, the output of which could be subsequently used for finite element analysis using a software package such as ABAQUS. However in these methods, most microscopic architecture of bone was removed during the transformation of the surface meshes into the volume ones. Neglecting this important structural characteristic of bone may lead to difficulty in accurate prediction of mechanical properties of a bone structure in subsequent finite element analysis. Although some of the microscopic bone features may be reconstructed manually or semi-automatically using image processing software, substantial time and substantial effort are required for remeshing. The voxel-based meshing method was developed by Keyak et al. [10–12] to directly generate “voxel mesh” from a dataset of stacked sliced CT image data for eliminating the geometry extraction step. With a high spatial resolution of microfocus computed tomography (micro-CT) scanner, detailed characterization of microscopic bone architecture in three dimensions can be performed [13–16]. Keyak’s method has been adopted to build microstructural finite element (micro-FE) models of microscopic structures, such as small regions of trabecular bone using micro-sized elements usually around 50–60 µm per side [17–19]. In the work of van Rietbergen et al. [18], a 3D reconstruction of a small 5 mm trabecular bone cube was firstly obtained by digitizing micro-CT image cross-sections, and then converted into a microFE model with 296,679 equally sized linear brick elements from the voxels size of ∼50 µm per side. Instead of using the micro-CT images for constructing the micro-FE model, high resolution (20 µm) images of a small trabecular bone sample which was embedded in methyl-methacrylate and then milled off a layer approximately 20 µm thick could be captured with a digital camera [20,21]. Although the microscopic architecture of the bone could be analyzed by using micro-FE models, a huge amount of elements is required even for a tiny part of bone, limiting the application of this method for the FE analysis of small pieces of bone rather than a whole bone segment with different types of bone such as trabecular and cortical ones. The micro-FE model was also applied to investigate the post-yield behavior of trabecular bone tissue [19]. In order to reduce computational time, the region averaging technique [19] was used to scale down the image resolution from 20 µm to 60 µm in their studies for subsequent finite element mesh generation. However, the effects of a varying number of image pixels used in forming an element on predictive accuracy of their FE analyses were not reported.
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1471
Fig. 1. Human metacarpal bone obtained using the micro-CT images.
The purpose of the present work is to develop a method for automated generation of 3D FE models of entire bone with microscopic architecture from micro-CT scan images. A metacarpal (MCP) bone was adopted to demonstrate the proposed method. The novelty of the proposed method is that FE analyses for the entire bone with microscopic architecture could be made possible even on a general personal computer by using grayscale intensity thresholds to reduce the amount of elements. Dual grayscale intensity threshold parameters were incorporated into the meshing method to distinguish the bone tissues from the surrounding soft tissues and improve predictive accuracy for the FE analyses with different sizes of elements. The effects of varying element size and the second grayscale threshold parameter on the geometry of the reconstructed bone as well as a FE analysis performed on a personal computer and validation of the methodology were investigated. 2. Materials and method In order to automatically construct a FE model of MCP bone with internal microscopic architecture from the micro-CT sliced images of the bone, a computer program was developed using MATLAB software. The FE model of the MCP bone with internal microscopic architecture was created in the form of ABAQUS input file. 2.1. Image acquisition and pre-processing The micro-CT images that used to construct the MCP bone shown in Fig. 1 were obtained by using a SkyScan 1076 machine (from SKYSCAN, Kontich, Belgium) at 100 kV source voltage and 100 µA current by line scanning at 80 ms intervals. During the scanning, a 1 mm thick aluminum filter was used to reduce the artifacts from the rotary scanning at the rotational step of 0.7 degree. The scanning images were set at 16-bit grayscale with 616 pixels squared with each pixel size of 34.66 µm per side. The spatial resolution of this device could be as low as 9 µm pixel size so that the inter-slices spacing of these set of images were set at 34.66 µm with no overlapping on adjacent images for the ease of handling at the phase of 3D reconstruction. A total of 1780 sections covering the complete bone length along the axis of the bone shaft from which images at cross-section no. 300, 890 and 1480 were acquired as shown in Fig. 2. The micro-CT sliced images were firstly segmented by the grayscale intensity method. The segmentation process makes use of the first grayscale intensity threshold T1 to remove the pixels corresponding to the non-bony tissues in the micro-CT images [22]. Those pixels with the grayscale intensity below a threshold value T1 were considered as the region of the non-bony tissues like soft tissues, and then removed from the corresponding micro-CT sliced images of the MCP bone by changing the grayscale intensity of these pixels to zero.
1472
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 2. Samples of micro-CT sliced images at frame numbers 300, 890 and 1480.
To enhance the image processing in the following steps, the size of all obtained micro-CT images was reduced from 616 pixel × 616 pixel to 500 pixel × 500 pixel by cropping the black borders, while the pixel size remained unchanged. The pre-processed micro-CT scan dataset were converted into a threedimensional array of pixels describing their identification number, position and grayscale intensity. 2.2. FE meshing Each 3D voxel was treated as a cube containing one pixel or an array of pixels, and converted into an equally sized three-dimensional linear eight-node finite element, no matter how many pixels are contained in the voxel. As the pixel size is 34.66 µm for the pre-processed micro-CT images of the MCP bone, a voxel containing only one pixel and the corresponding element would have a length of 34.66 µm per side. If a voxel could contain up to 27 pixels (i.e., a 3D array of 3 × 3 × 3 pixels), the element would become 103.98 µm per side, which is three times bigger in size than the element built from the voxel of 34.66 µm per side. It should be noted that, for the region with both bony and non-bony tissues, the pixels corresponding to these non-bony materials would be deleted using the T1 value during the image pre-processing step. Therefore, the number of pixels corresponding to the bony tissue in a voxel could be lower than 27 or even only one for the element of 103.98 µm per side. 2.3. Element reduction For each element, there is a corresponding voxel containing from one pixel up to a maximum number of pixels depending on the size of the element used. Therefore, an average grayscale intensity of the element was calculated by taking the average of the grayscale intensity at all pixels for the bony region within the voxel. For the element with low average grayscale intensity, the element could correspond to the voxel containing either only a few pixels or the pixels with very low grayscale intensity. If the average grayscale intensity of a certain element is lower than the second grayscale intensity threshold, T2 , these elements would have less effect on the structural response of the bone and hence be deleted from the FE mesh of the bone. Figure 3 demonstrates how to remove the less important elements using the value T2 . Suppose a simple cylinder with a diameter, D of 160 mm and the image of its cross-section consisting of 192 × 192 pixels as shown in Fig. 3(a) was used. Figure 3(b)–(d) showed a comparison between the original boundary of the cylinder and the boundaries of the FE models of the reconstructed cylinder cross-section using different sizes of elements with side lengths of 3, 6 and 9 mm. It could be apparently found that the size of the cylinder cross-section reconstructed by using the larger sized elements as shown in Fig. 3(d) was much bigger than that with the smaller sized elements as shown in Fig. 3(b). From Fig. 4, the effects of varying values of T2 from 0 to 90 on the geometry of the cylinder cross-sections reconstructed by
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1473
Fig. 3. Effect of element sizes on the geometry of the reconstructed cylinder cross-sections with a fixed value of T2 , (a) an image of cylinder cross-section, (b) element side length = 3 mm, T2 = 0, (c) element side length = 6 mm, T2 = 0 and (d) element side length = 9 mm, T2 = 0.
Fig. 4. Effect of varying T2 on the geometry of the reconstructed cylinder cross-sections, (a) element size = 9 mm, T2 = 0, (b) element size = 9 mm, T2 = 30, (c) element size = 9mm, T2 = 60 and (d) element size = 9 mm, T2 = 90.
using the elements of the size of 9 mm. It could be observed that the size of the reconstructed cylinder cross-section without the use of T2 as shown in Fig. 4(a) was much bigger and has a large deviation in shape from the original one as compared to that with considering T2 value of 90 as shown in Fig. 4(d). We aimed to determine the value of T2 such that the area error is closed to zero. An area error, ea , which was defined as the relative difference between the reconstructed cross-sectional area from the FE model, Ab and the original cross-sectional area of the cylinder, Af , was calculated by ea = c − 1 =
Ab − Af , Af
(1)
where c is the relative area coefficient (c = Ab /Af ). The area, Ab , is calculated by determining the product of number of elements, n and area of one element, ab . From Eq. (1), the area error is minimized when the value of c is close to one. Therefore, the value of c close to one for minimizing the area error were used to serve as the criterion for determining the value of T2 in reconstructing FE model of human bone in this work. Finally, the disconnected elements in the FE model of the reconstructed structure were deleted by recording the neighbouring information of all elements in the FE model. Any elements without any neighbouring elements and any isolated group of elements were removed from the FE mesh as shown in Fig. 5. 2.4. Validation The accuracy of the constructed FE models using our proposed methodology was validated by conducting a tensile test to compare the results obtained from the real experiment and FE simulation.
1474
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 5. Deletion of the isolated elements in the FE models, (a) sliced images, (b) FE model without deletion of the isolated elements, (c) FE model after deletion of the isolated elements. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
Fig. 6. Loading and boundary conditions on the human metacarpal bone for FE analysis. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
2.4.1. Setup of FE simulation The deformation behavior of MCP bone under an applied load was simulated using ABAQUS software. Elements with different sizes of the corresponding threshold value T2 were used in constructing the FE meshes. Figure 6 shows the positions of the loading and boundary conditions. Since the present analysis is mainly focused on the central portion of the bone model, regions A and B were not considered in the simulation to reduce the computation time. The left end of the bone model after removal of the region A was fixed along z-direction while a concentrated force is applied on the right end after removal of the region B along z-direction. 2.4.2. Real experimental setup In the real experiment, tensile tests of MCP bone specimen were conducted on a material testing system (Instron 5565) at the rate of 100 N/min up to the required load at the room temperature. Figure 7 shows the setup of the experiment. Both ends of the MCP bone specimen with an average diameter of 9.8 mm and a full length of 86 mm were completely wrapped with an epoxy-based adhesive and separately fixed in aluminum alloy specimen holders as shown in Fig. 8. Upper and lower specimen support fixtures with four adjustable bolts as shown in Fig. 8 were specially designed to ensure concentricity of the applied load and parallelism of the interlocking surfaces between the specimen holder and the fixture. Two additional bolts were incorporated into the fixtures to allow for the specimen to be mounted
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1475
Fig. 7. Setup for the tensile testing of the MCP bone specimen. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
Fig. 8. MCP bone specimen mounted in custom-made support fixtures. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
onto the grippers of the tensile tester. The displacements of the specimen were measured by a calibrated extensometer with a gauge length of 25 mm.
1476
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
3. Results In order to compare the results of varying the element size on the geometry of the reconstructed bone model and stress/strain results of FE analysis, five different sizes of the 3D brick element, Se (A–E) with dimensions and its relationship with the number of elements in a voxel as listed in Table 1 were generated. The 3D micro-CT image and an example of the generated 3D FE mesh of the human MCP bone with element size “D” were shown in Fig. 9(a) and (b) respectively. 3.1. Effect of element size It has been well known that the computational time for a FE analysis is highly depends on the number of elements used in the FE model. If the element size is too small, the number of elements required for constructing the FE model of bone increases sharply and may exceed computational capability. For the case of using too large element, the internal microscopic architecture of bone may not be correctly reconstructed in the FE model. By using element sizes of “A” to “E” in building FE models of the reconstructed human MCP bone, a different resolution of micro-architecture in the bone at the cross-sections corresponding to microCT sliced image numbered as 300, 890 and 1480 were reconstructed in Fig. 10. It can be apparently observed from Fig. 10 that the area of the bony region in the FE mesh of the human MCP bone created using different element sizes was different for the same cross-section of the bone. In order to study their differences quantitatively, the relative area coefficient, c for different element sizes for three different cross-sections of the bone, Ncs derived from the micro-CT scan was calculated as listed in Table 2. As Table 1 Relationship between element dimensions and number of elements in an image voxel Element size (Se ) A B C D E
Dimensions of the element 34.66 µm × 34.66 µm × 34.66 µm 103.98 µm × 103.98 µm × 103.98 µm 138.64 µm × 138.64 µm × 138.64 µm 173.30 µm × 173.30 µm × 173.30 µm 242.62 µm × 242.62 µm × 242.62 µm
Number of elements in an image voxel 1 1–27 1–64 1–125 1–343
Fig. 9. Human metacarpal bone, (a) 3D micro-CT image and (b) 3D FE model. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1477
Table 2 Relationship between element sizes and relative area coefficient c at three different cross-sections of the human MCP bone Cross-section number (Ncs ) 300
890
1480
Element size (Se ) A B C D E A B C D E A B C D E
Number of elements used (n) 26,029 5288 3455 2379 1468 20,057 2585 1494 978 542 24,040 4337 2760 1905 1141
Area of an element (ab (mm2 )) 1af 9af 16af 25af 49af 1af 9af 16af 25af 49af 1af 9af 16af 25af 49af
Area of the bony region (Ab (mm2 )) Af = 26,029af 47,592af 55,280af 59,475af 71,932af Af = 20,057af 23,265af 23,904af 24,450af 26,558af Af = 24,040af 39,033af 44,160af 47,625af 55,909af
Relative area coefficient (c) 1.00 1.83 2.12 2.28 2.76 1.00 1.16 1.19 1.22 1.32 1.00 1.62 1.84 1.98 2.32
the element of size “A” was built from a voxel containing only one pixel, the area of the bony region in this element is regarded as the reference area, Af in which the area of one element was treated as a reference value, af (i.e., af = 0.034662 mm2 ). In this study, the base material properties of the trabecular and cortical bone tissues were assumed to be the same [20]. With the increase of the element size, the FE models with elements of sizes “D” and “E” as shown in Fig. 10(d) and (e) became unable to accurately reflect the microscopic architecture of the MCP bone because the relative area coefficient, c for these FE meshes were much higher than 1 as shown in Table 2. Thus, the stress and strain results of the corresponding FE analyses would be affected. In order to reduce the FE analysis error, a second grayscale intensity threshold T2 was introduced to eliminate the elements which had reduced the area coefficient c close to one. 3.2. Effect of varying the value of T2 The effect of varying second grayscale intensity threshold values T2 on the relative area coefficient c at three different cross-sections of the human MCP bone are determined as listed in Table 3. From this table, for the elements sizes “B” and “C”, the most appropriate value of T2 is 65 in order to reduce the value of c approaching to one and hence minimize the area errors of the reconstructed cross-sections of the bone, while the most appropriate value of T2 is 60 for the element size “D”. By using the relative area coefficient c with the value of 1, the total elemental volume of FE meshes constructed from the elements with sizes “C”–“E” could be kept to vary only within 5%. In order to compare the microscopic geometry of the human MCP bone generated by commercially available 3D reconstruction tool “Mimics” and the proposed method based on the obtained micro-CT data, three different regions containing trabecular bone, cortical bone,and both of them are selected and labeled as “P”, “Q” and “R”, respectively in Fig. 11(a). By using different element sizes and corre-
1478
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 10. FE meshes of the reconstructed human MCP bone with varying element sizes (A–E) showing micro-architecture of the bone at three different cross-sections corresponding to micro-CT sliced image numbered as 300 (left), 890 (middle) and 1480 (right), (a) 34.66 µm×34.66 µm×34.66 µm, (b) 103.98 µm×103.98 µm×103.98 µm, (c) 138.64 µm×138.64 µm×138.64 µm, (d) 170.30 µm×170.30 µm×173.30 µm, and (e) 242.62 µm×242.62 µm×242.62 µm. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
sponding values of T2 , various FE meshes for three different regions of the bone were created as shown in Fig. 11(b–d). The FE mesh created using the element size “A” shown in Fig. 11(b) was used to serve as a reference for comparison with other meshes. It can be observed from Fig. 11 that the resolution of the reconstructed bone regions was reduced with increasing element size. Therefore, too large element would not be recommended to be used for local microscopic stress analysis of the bone. Although it was the case, it can be observed that the overall geometry of some bone regions reconstructed using the larger sized elements of “B” and “C” (Fig. 11(c) and (d)) is similar from those with the element size “A”, because the shape errors of these bone regions have been minimized by using the T2 value.
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1479
Table 3 Effect of varying second grayscale intensity threshold values T2 on the relative area coefficient c at three different cross-sections of the human MCP bone Cross-section number (Ncs ) 300
890
1480
Element size (Se ) A B B B B B C D E A B B B B B C D E A B B B B B C D E
Second grayscale intensity threshold parameter (T2 ) 0 0 20 40 60 65 65 60 60 0 0 20 40 60 65 65 60 60 0 0 20 40 60 65 65 60 60
Number of elements used (n) 26,029 5288 4121 3391 2917 2836 1592 1019 520 20,057 2585 2433 2332 2265 2250 1272 839 448 24,040 4337 3591 3114 2737 2691 1533 1021 540
Area of an element (ab (mm2 )) 1af 9af 9af 9af 9af 9af 16af 25af 49af 1af 9af 9af 9af 9af 9af 16af 25af 49af 1af 9af 9af 9af 9af 9af 16af 25af 49af
Area of the bony region (Ab (mm2 )) Af = 26,029af 47,592af 37,089af 30,519af 26,253af 25,524af 25,472af 25,475af 25,480af Af = 20,057af 23,265af 21,897af 20,988af 20,385af 20,250af 20,352af 20,975af 21,952af Af = 24,040af 39,033af 32,319af 28,026af 24,633af 24,219af 24,528af 25,525af 26,460af
Relative area coefficient (c) 1.00 1.83 1.42 1.17 1.01 0.98 0.98 0.98 0.98 1.00 1.16 1.09 1.05 1.01 1.01 1.01 1.05 1.09 1.00 1.62 1.34 1.17 1.02 1.01 1.02 1.06 1.10
3.3. FE analysis In order to illustrate the applicability of the proposed method, a series of FE analyses of the MCP bone were conducted using ABAQUS software with a general personal computer (Mainboard: MSI P55GD55, CPU: Intel Core i5-750 2.67 GHz, RAM: 8 GB DDR3 1333). The effects of varying element sizes on the stress distribution and force–displacement relation of the bone under the same boundary and load conditions were investigated. Elements of three different sizes “C”, “D” and “E” with the corresponding threshold value T2 were used in constructing the FE meshes. For the FE model composed of the elements of size “C”, a total of 496,149 elements and 829,443 nodes are used. The corresponding results were plotted in Figs 12 and 13. It can be observed from Fig. 12 that the pattern of Von Mises stress distribution for different FE meshes is similar, although the maximum stress was higher in the FE model with the elements size “C”. From Fig. 13, the force–displacement data predicted by using the element size “C” lied between those with larger-size elements “D” and “E”, but closer to that with the element of size D. This implies that larger-sized elements could be used for the macroscopic stress analysis of the bone, but
1480
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 11. 3D reconstruction of three particular regions “P”, “Q” and “R” in human MCP bone based on micro-CT data of the bone, (a) 3D images in Mimics, (b)–(d) 3D FE meshes generated with variation of the values of Se , T2 and n: (b) Se = A, T2 = 0, (c) Se = B, T2 = 65, (d) Se = C, T2 = 65. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1481
Fig. 12. Von Mises stress contour plots of FE models of the MCP bone with three different sizes of elements and corresponding values of T2 , (a) Se = C, T2 = 65, n = 496,149, (b) Se = D, T2 = 60, n = 265,452 and (c) Se = E, T2 = 60, n = 96,917. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
too large element like the size “E” is not recommended due to the large deviation in the displacement as compared with the smaller sized elements. 3.4. Model validation A tensile test that compares the force–displacement relationship of simulated and real experiments was conducted to validate the FE MCP bone model reconstructed by using our proposed methodology in reducing the number of elements. The trabecular and cortical bone tissues are assumed to have the same material properties. The Young’s modulus (E) and Poisson’s ratio (ν) of the bone tissues are assumed to be 18.7 GPa and 0.3, respectively [20]. Based on the results from Fig. 13, the element size “D” was used and the applied load was 40 N. As the cross-sectional area of the bone specimen is non-uniform, two additional measurement of the force–displacement relation of the specimen is made by rotating the specimen along the loading axis by 120◦ and 240◦ . The force–displacement relationship was illustrated in Fig. 14. With the increase in the displacement, the load generally showed an increasing trend for the three different testing positions. It can be observed that the trends of the predicted results were in good agreement with the experimental data, demonstrating the general applicability of the proposed method to the model validation.
1482
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 13. Displacement–force relation of the MCP bone with three different sizes of elements and corresponding values of T2 . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
Fig. 14. Load–displacement plots for tensile testing of the MCP bone specimen.
4. Conclusion A methodology for automated generation of three-dimensional FE models of entire bone with microscopic architecture from micro-CT image data is successfully developed. The applicability of the
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1483
proposed method is demonstrated using human MCP bone. The dual grayscale intensity threshold parameters are adopted to distinguish the bone tissues from the surrounding soft tissues and improve predictive accuracy for the FE analyses for different sizes of elements. With the proposed method, FE analyses of an entire bone can be run on a personal computer by using a smaller number of larger-size elements with maintaining certain predictive accuracy. Experiments are performed to simulate different sizes of FE element models. In addition, we attempted to validate the methodology by comparing the results from the FE analysis and the tensile test data. The simulation results show good agreement with the experimental results. It is believed that the proposed method can also be extended to model the bone with variation of mechanical properties internally and functional graded structures based on the average grayscale intensity within elements. Acknowledgement The work described in this paper is supported by a grant from the Research Committee of the Hong Kong Polytechnic University (Project No. G-YG18). References [1] R. Huiskes and E.Y.S. Chao, A survey of finite element analysis in orthopedic biomechanics: the first decade, J. Biomech. 16 (1983), 385–409. [2] T. Nomura, M.P. Powers, J.L. Katz and C. Saito, Finite element analysis of a transmandibular implant, J. Biomed. Mater. Res. B – Appl. Biomater. 80 (2007), 370–376. [3] M. Viceconti, L. Bellingeri, L. Cristofolini and A. Toni, A comparative study on different methods of automatic mesh generation of human femurs, Med. Eng. Phys. 20 (1998), 1–10. [4] O. Panagiotopoulou, N. Curtis, P.O. Higgins and S.N. Cobb, Modelling subcortical bone in finite element analyses: A validation and sensitivity study in the macaque mandible, J. Biomech. 43 (2009), 1603–1611. [5] F.A. Bandak, R.E. Tannous and T. Toridis, On the development of an osseo-ligamentous finite element model of the human ankle joint, Intl. J. Solids Struct. 38 (2001), 1681–1697. [6] P. Messmer, F. Matthews, A.L. Jacob, R. Kikinis, P. Regazzoni and H. Noser, A CT database for research, development and education: concept and potential, J. Digital Imaging 20 (2007), 17–22. [7] J.A. MacNeil and S.K. Boyd, Bone strength at the distal radius can be estimated from high-resolution peripheral quantitative computed tomography and the finite element method, Bone 42 (2008), 1203–1213. [8] D.S. Barker, D.J. Netherway, J. Krishnan and T.C. Hearn, Validation of a finite element model of the human metacarpal, Med. Eng. Phys. 27 (2005), 103–113. [9] A. Kato and N. Ohno, Construction of three-dimensional tooth model by micro-computed tomography and application for data sharing, Clin. Oral Invest. 13 (2009), 43–46. [10] J.H. Keyak, J.M. Meagher, H.B. Skinner and C.D. Mote, Automated three-dimensional finite element modelling of bone: a new method, J. Biomed. Eng. 12 (1990), 389–397. [11] J.H. Keyak and H.B. Skinner, Three dimensional finite element modeling of bone: effect of element size, J. Biomed. Eng. 14 (1992), 483–489. [12] J.H. Keyak, M.G. Fourkas, J.M. Meagher and H.B. Skinner, Validation of an automated method of three dimensional finite-element modeling of bone, J. Biomed. Eng. 15 (1993), 505–509. [13] G. Kochi, S. Sato, T. Fukuyama, C. Morita, K. Honda, Y. Arai and K. Ito, Analysis on the guided bone augmentation in the rat calvarium using a microfocus computerized tomography analysis, Oral Surg. Oral Med. Oral Pathol. Oral Radiol. Endod. 107 (2009), e42–e48. [14] C. Efeoglu, S.E. Fisher, S. Erturk, F. Oztop, S. Gunbay and A. Sipahi, Quantitative morphometric evaluation of critical size experimental bone defects by microcomputed tomography, Brit. J. Oral. Maxillofacial Surg. 45 (2007), 203–207. [15] G. Bevill and T.M. Keaveny, Trabecular bone strength predictions using finite element analysis of micro-scale images at limited spatial resolution, Bone 44 (2009), 579–584. [16] J. Peng, C.Y. Wen, A.Y. Wang, Y. Wang, W.J. Xu, B. Zhao, L. Zhang, S.B. Lu, L. Qin, Q.Y. Guo, L.M. Dong and J.M. Tian, Micro-CT-based bone ceramic scaffolding and its performance after seeding with mesenchymal stem cells for repair of load-bearing bone defect in canine femoral head, J. Biomed. Mater. Res. B – Appl. Biomater. 96 (2011), 316–325.
1484
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
[17] S.J. Hollister, J.M. Brennan and N. Kikuchi, A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress, J. Biomech. 27 (1994), 433–444. [18] B.V. Rietbergen, H. Weinans, R. Huiskes and A. Odgaard, A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models, J. Biomech. 28 (1995), 69–81. [19] E. Verhulp, B.V. Rietbergen, R. Muller and R. Huiskes, Indirect determination of trabecular bone effective tissue failure properties using micro-finite element simulations, J. Biomech. 41 (2008), 1479–1485. [20] G.L. Niebur, J.C. Yuen, A.C. Hsia and T.M. Keaveny, Convergence behavior of high-resolution finite element models of trabecular bone, J. Biomech. Eng. 121 (1999), 629–635. [21] G.L. Niebur, M.J. Feldstein, J.C. Yuen, T.J. Chen and T.M. Keaveny, High-resolution finite element models with tissue strength asymmetry accurately predict failure of trabeclar bone, J. Biomech. 33 (2000), 1575–1583. [22] S.V.N. Jaecques, L. Muraru, C.V. Lierde, E.D. Smet, H.V. Oosterywck, M. Wevers, I. Naert and J.V. Sloten, In vivo micro-CT-based FE models of guinea pigs with titanium implants: an STL-based approach, Inter. Cong. Ser. 1268 (2004), 579–583.
Copyright of Bio-Medical Materials & Engineering is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
1469
Voxel-based approach to generate entire human metacarpal bone with microscopic architecture for finite element analysis C.Y. Tang a,∗ , C.P. Tsui a , Y.M. Tang a , L. Wei a , C.T. Wong a , K.W. Lam b , W.Y. Ip b , W.W.J. Lu b and M.Y.C. Pang c a
Department of Industrial and Systems Engineering, The Hong Kong Polytechnic University, Hong Kong, China b Department of Orthopaedics and Traumatology, Li Ka Shing Faculty of Medicine, The University of Hong Kong, Hong Kong, China c Department of Rehabilitation Sciences, The Hong Kong Polytechnic University, Hong Kong, China Received 5 June 2012 Accepted 31 July 2013 Abstract. With the development of micro-computed tomography (micro-CT) technology, it is possible to construct threedimensional (3D) models of human bone without destruction of samples and predict mechanical behavior of bone using finite element analysis (FEA). However, due to large number of elements required for constructing the FE models of entire bone, this demands a substantial computational effort and the analysis usually needs a high level of computer. In this article, a voxel-based approach for generation of FE models of entire bone with microscopic architecture from micro-CT image data is proposed. To enable the FE analyses of entire bone to be run even on a general personal computer, grayscale intensity thresholds were adopted to reduce the amount of elements. Human metacarpal bone (MCP) bone was used as an example for demonstrating the applicability of the proposed method. The micro-CT images of the MCP bone were combined and converted into 3D array of pixels. Dual grayscale intensity threshold parameters were used to distinguish the pixels of bone tissues from those of surrounding soft tissues and improve predictive accuracy for the FE analyses with different sizes of elements. The method of selecting an appropriate value of the second grayscale intensity threshold was also suggested to minimize the area error for the reconstructed cross-sections of a FE structure. Experimental results showed that the entire FE MCP bone with microscopic architecture could be modeled and analyzed on a personal computer with reasonable accuracy. Keywords: Voxel-based, finite element model, micro-computed tomography (micro-CT), microscopic architecture, human metacarpal bone, grayscale intensity threshold
1. Introduction The Finite Element method (FEM) has been shown to be a powerful tool for determining the solution of the physical mechanical behavior of the human skeleton in the field of orthopedic biomechanics [1,2]. However, bone is an osseous tissue with very complex organization consisting of an intricate cortical * Address for correspondence: C.Y. Tang, Department of Industrial and Systems Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, China. Tel.: +852 27666608; Fax: +852 23625267; E-mail: [email protected].
0959-2989/14/$27.50 © 2014 – IOS Press and the authors. All rights reserved
1470
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
layer and a randomized trabecular pattern of internal microstructures. The difficulties of capturing fine detailed geometrical data and generation of three-dimensional (3D) finite element (FE) meshes complicate the implementation of finite element analysis [3]. With the recent advancement in computational power and medical imaging techniques, FEM is now able to nondestructively model stress, strain and deformation in complex 3D bony structures [4]. 3D FE models are usually generated from computed tomography (CT) scan datasets through various approaches [5–7], because the volume and surface information of either a whole bone or bone segments can be acquired from the CT scan. Based on anatomically realistic bone geometry obtained from a CT scan dataset, Bandak [5] developed a 3D FE model of the human lower extremity for studying mechanisms of impact injury to the lower leg. The method requires the CT scan dataset to be pre-processed to extract the geometry of the bone and through creation of inner and outer surfaces of the bone by using spline interpolation between the virtual data points from the CT image. A grayscale intensity threshold was used to isolate the bone region from soft tissues and the inner medullary surface [6,8]. Another commonly used method for creation of the CT-based 3D bone model is to firstly segment bone and inner medullary surfaces into triangles. The segmented results was saved in the Stereo-Lithography Interface Format (STL) for screen visualization and rapid prototyping [6,9]. The surface meshes in STL format after mesh optimization could be transformed into volume meshes using software such as MIMICS and 3D Slicer, the output of which could be subsequently used for finite element analysis using a software package such as ABAQUS. However in these methods, most microscopic architecture of bone was removed during the transformation of the surface meshes into the volume ones. Neglecting this important structural characteristic of bone may lead to difficulty in accurate prediction of mechanical properties of a bone structure in subsequent finite element analysis. Although some of the microscopic bone features may be reconstructed manually or semi-automatically using image processing software, substantial time and substantial effort are required for remeshing. The voxel-based meshing method was developed by Keyak et al. [10–12] to directly generate “voxel mesh” from a dataset of stacked sliced CT image data for eliminating the geometry extraction step. With a high spatial resolution of microfocus computed tomography (micro-CT) scanner, detailed characterization of microscopic bone architecture in three dimensions can be performed [13–16]. Keyak’s method has been adopted to build microstructural finite element (micro-FE) models of microscopic structures, such as small regions of trabecular bone using micro-sized elements usually around 50–60 µm per side [17–19]. In the work of van Rietbergen et al. [18], a 3D reconstruction of a small 5 mm trabecular bone cube was firstly obtained by digitizing micro-CT image cross-sections, and then converted into a microFE model with 296,679 equally sized linear brick elements from the voxels size of ∼50 µm per side. Instead of using the micro-CT images for constructing the micro-FE model, high resolution (20 µm) images of a small trabecular bone sample which was embedded in methyl-methacrylate and then milled off a layer approximately 20 µm thick could be captured with a digital camera [20,21]. Although the microscopic architecture of the bone could be analyzed by using micro-FE models, a huge amount of elements is required even for a tiny part of bone, limiting the application of this method for the FE analysis of small pieces of bone rather than a whole bone segment with different types of bone such as trabecular and cortical ones. The micro-FE model was also applied to investigate the post-yield behavior of trabecular bone tissue [19]. In order to reduce computational time, the region averaging technique [19] was used to scale down the image resolution from 20 µm to 60 µm in their studies for subsequent finite element mesh generation. However, the effects of a varying number of image pixels used in forming an element on predictive accuracy of their FE analyses were not reported.
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1471
Fig. 1. Human metacarpal bone obtained using the micro-CT images.
The purpose of the present work is to develop a method for automated generation of 3D FE models of entire bone with microscopic architecture from micro-CT scan images. A metacarpal (MCP) bone was adopted to demonstrate the proposed method. The novelty of the proposed method is that FE analyses for the entire bone with microscopic architecture could be made possible even on a general personal computer by using grayscale intensity thresholds to reduce the amount of elements. Dual grayscale intensity threshold parameters were incorporated into the meshing method to distinguish the bone tissues from the surrounding soft tissues and improve predictive accuracy for the FE analyses with different sizes of elements. The effects of varying element size and the second grayscale threshold parameter on the geometry of the reconstructed bone as well as a FE analysis performed on a personal computer and validation of the methodology were investigated. 2. Materials and method In order to automatically construct a FE model of MCP bone with internal microscopic architecture from the micro-CT sliced images of the bone, a computer program was developed using MATLAB software. The FE model of the MCP bone with internal microscopic architecture was created in the form of ABAQUS input file. 2.1. Image acquisition and pre-processing The micro-CT images that used to construct the MCP bone shown in Fig. 1 were obtained by using a SkyScan 1076 machine (from SKYSCAN, Kontich, Belgium) at 100 kV source voltage and 100 µA current by line scanning at 80 ms intervals. During the scanning, a 1 mm thick aluminum filter was used to reduce the artifacts from the rotary scanning at the rotational step of 0.7 degree. The scanning images were set at 16-bit grayscale with 616 pixels squared with each pixel size of 34.66 µm per side. The spatial resolution of this device could be as low as 9 µm pixel size so that the inter-slices spacing of these set of images were set at 34.66 µm with no overlapping on adjacent images for the ease of handling at the phase of 3D reconstruction. A total of 1780 sections covering the complete bone length along the axis of the bone shaft from which images at cross-section no. 300, 890 and 1480 were acquired as shown in Fig. 2. The micro-CT sliced images were firstly segmented by the grayscale intensity method. The segmentation process makes use of the first grayscale intensity threshold T1 to remove the pixels corresponding to the non-bony tissues in the micro-CT images [22]. Those pixels with the grayscale intensity below a threshold value T1 were considered as the region of the non-bony tissues like soft tissues, and then removed from the corresponding micro-CT sliced images of the MCP bone by changing the grayscale intensity of these pixels to zero.
1472
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 2. Samples of micro-CT sliced images at frame numbers 300, 890 and 1480.
To enhance the image processing in the following steps, the size of all obtained micro-CT images was reduced from 616 pixel × 616 pixel to 500 pixel × 500 pixel by cropping the black borders, while the pixel size remained unchanged. The pre-processed micro-CT scan dataset were converted into a threedimensional array of pixels describing their identification number, position and grayscale intensity. 2.2. FE meshing Each 3D voxel was treated as a cube containing one pixel or an array of pixels, and converted into an equally sized three-dimensional linear eight-node finite element, no matter how many pixels are contained in the voxel. As the pixel size is 34.66 µm for the pre-processed micro-CT images of the MCP bone, a voxel containing only one pixel and the corresponding element would have a length of 34.66 µm per side. If a voxel could contain up to 27 pixels (i.e., a 3D array of 3 × 3 × 3 pixels), the element would become 103.98 µm per side, which is three times bigger in size than the element built from the voxel of 34.66 µm per side. It should be noted that, for the region with both bony and non-bony tissues, the pixels corresponding to these non-bony materials would be deleted using the T1 value during the image pre-processing step. Therefore, the number of pixels corresponding to the bony tissue in a voxel could be lower than 27 or even only one for the element of 103.98 µm per side. 2.3. Element reduction For each element, there is a corresponding voxel containing from one pixel up to a maximum number of pixels depending on the size of the element used. Therefore, an average grayscale intensity of the element was calculated by taking the average of the grayscale intensity at all pixels for the bony region within the voxel. For the element with low average grayscale intensity, the element could correspond to the voxel containing either only a few pixels or the pixels with very low grayscale intensity. If the average grayscale intensity of a certain element is lower than the second grayscale intensity threshold, T2 , these elements would have less effect on the structural response of the bone and hence be deleted from the FE mesh of the bone. Figure 3 demonstrates how to remove the less important elements using the value T2 . Suppose a simple cylinder with a diameter, D of 160 mm and the image of its cross-section consisting of 192 × 192 pixels as shown in Fig. 3(a) was used. Figure 3(b)–(d) showed a comparison between the original boundary of the cylinder and the boundaries of the FE models of the reconstructed cylinder cross-section using different sizes of elements with side lengths of 3, 6 and 9 mm. It could be apparently found that the size of the cylinder cross-section reconstructed by using the larger sized elements as shown in Fig. 3(d) was much bigger than that with the smaller sized elements as shown in Fig. 3(b). From Fig. 4, the effects of varying values of T2 from 0 to 90 on the geometry of the cylinder cross-sections reconstructed by
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1473
Fig. 3. Effect of element sizes on the geometry of the reconstructed cylinder cross-sections with a fixed value of T2 , (a) an image of cylinder cross-section, (b) element side length = 3 mm, T2 = 0, (c) element side length = 6 mm, T2 = 0 and (d) element side length = 9 mm, T2 = 0.
Fig. 4. Effect of varying T2 on the geometry of the reconstructed cylinder cross-sections, (a) element size = 9 mm, T2 = 0, (b) element size = 9 mm, T2 = 30, (c) element size = 9mm, T2 = 60 and (d) element size = 9 mm, T2 = 90.
using the elements of the size of 9 mm. It could be observed that the size of the reconstructed cylinder cross-section without the use of T2 as shown in Fig. 4(a) was much bigger and has a large deviation in shape from the original one as compared to that with considering T2 value of 90 as shown in Fig. 4(d). We aimed to determine the value of T2 such that the area error is closed to zero. An area error, ea , which was defined as the relative difference between the reconstructed cross-sectional area from the FE model, Ab and the original cross-sectional area of the cylinder, Af , was calculated by ea = c − 1 =
Ab − Af , Af
(1)
where c is the relative area coefficient (c = Ab /Af ). The area, Ab , is calculated by determining the product of number of elements, n and area of one element, ab . From Eq. (1), the area error is minimized when the value of c is close to one. Therefore, the value of c close to one for minimizing the area error were used to serve as the criterion for determining the value of T2 in reconstructing FE model of human bone in this work. Finally, the disconnected elements in the FE model of the reconstructed structure were deleted by recording the neighbouring information of all elements in the FE model. Any elements without any neighbouring elements and any isolated group of elements were removed from the FE mesh as shown in Fig. 5. 2.4. Validation The accuracy of the constructed FE models using our proposed methodology was validated by conducting a tensile test to compare the results obtained from the real experiment and FE simulation.
1474
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 5. Deletion of the isolated elements in the FE models, (a) sliced images, (b) FE model without deletion of the isolated elements, (c) FE model after deletion of the isolated elements. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
Fig. 6. Loading and boundary conditions on the human metacarpal bone for FE analysis. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
2.4.1. Setup of FE simulation The deformation behavior of MCP bone under an applied load was simulated using ABAQUS software. Elements with different sizes of the corresponding threshold value T2 were used in constructing the FE meshes. Figure 6 shows the positions of the loading and boundary conditions. Since the present analysis is mainly focused on the central portion of the bone model, regions A and B were not considered in the simulation to reduce the computation time. The left end of the bone model after removal of the region A was fixed along z-direction while a concentrated force is applied on the right end after removal of the region B along z-direction. 2.4.2. Real experimental setup In the real experiment, tensile tests of MCP bone specimen were conducted on a material testing system (Instron 5565) at the rate of 100 N/min up to the required load at the room temperature. Figure 7 shows the setup of the experiment. Both ends of the MCP bone specimen with an average diameter of 9.8 mm and a full length of 86 mm were completely wrapped with an epoxy-based adhesive and separately fixed in aluminum alloy specimen holders as shown in Fig. 8. Upper and lower specimen support fixtures with four adjustable bolts as shown in Fig. 8 were specially designed to ensure concentricity of the applied load and parallelism of the interlocking surfaces between the specimen holder and the fixture. Two additional bolts were incorporated into the fixtures to allow for the specimen to be mounted
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1475
Fig. 7. Setup for the tensile testing of the MCP bone specimen. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
Fig. 8. MCP bone specimen mounted in custom-made support fixtures. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
onto the grippers of the tensile tester. The displacements of the specimen were measured by a calibrated extensometer with a gauge length of 25 mm.
1476
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
3. Results In order to compare the results of varying the element size on the geometry of the reconstructed bone model and stress/strain results of FE analysis, five different sizes of the 3D brick element, Se (A–E) with dimensions and its relationship with the number of elements in a voxel as listed in Table 1 were generated. The 3D micro-CT image and an example of the generated 3D FE mesh of the human MCP bone with element size “D” were shown in Fig. 9(a) and (b) respectively. 3.1. Effect of element size It has been well known that the computational time for a FE analysis is highly depends on the number of elements used in the FE model. If the element size is too small, the number of elements required for constructing the FE model of bone increases sharply and may exceed computational capability. For the case of using too large element, the internal microscopic architecture of bone may not be correctly reconstructed in the FE model. By using element sizes of “A” to “E” in building FE models of the reconstructed human MCP bone, a different resolution of micro-architecture in the bone at the cross-sections corresponding to microCT sliced image numbered as 300, 890 and 1480 were reconstructed in Fig. 10. It can be apparently observed from Fig. 10 that the area of the bony region in the FE mesh of the human MCP bone created using different element sizes was different for the same cross-section of the bone. In order to study their differences quantitatively, the relative area coefficient, c for different element sizes for three different cross-sections of the bone, Ncs derived from the micro-CT scan was calculated as listed in Table 2. As Table 1 Relationship between element dimensions and number of elements in an image voxel Element size (Se ) A B C D E
Dimensions of the element 34.66 µm × 34.66 µm × 34.66 µm 103.98 µm × 103.98 µm × 103.98 µm 138.64 µm × 138.64 µm × 138.64 µm 173.30 µm × 173.30 µm × 173.30 µm 242.62 µm × 242.62 µm × 242.62 µm
Number of elements in an image voxel 1 1–27 1–64 1–125 1–343
Fig. 9. Human metacarpal bone, (a) 3D micro-CT image and (b) 3D FE model. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1477
Table 2 Relationship between element sizes and relative area coefficient c at three different cross-sections of the human MCP bone Cross-section number (Ncs ) 300
890
1480
Element size (Se ) A B C D E A B C D E A B C D E
Number of elements used (n) 26,029 5288 3455 2379 1468 20,057 2585 1494 978 542 24,040 4337 2760 1905 1141
Area of an element (ab (mm2 )) 1af 9af 16af 25af 49af 1af 9af 16af 25af 49af 1af 9af 16af 25af 49af
Area of the bony region (Ab (mm2 )) Af = 26,029af 47,592af 55,280af 59,475af 71,932af Af = 20,057af 23,265af 23,904af 24,450af 26,558af Af = 24,040af 39,033af 44,160af 47,625af 55,909af
Relative area coefficient (c) 1.00 1.83 2.12 2.28 2.76 1.00 1.16 1.19 1.22 1.32 1.00 1.62 1.84 1.98 2.32
the element of size “A” was built from a voxel containing only one pixel, the area of the bony region in this element is regarded as the reference area, Af in which the area of one element was treated as a reference value, af (i.e., af = 0.034662 mm2 ). In this study, the base material properties of the trabecular and cortical bone tissues were assumed to be the same [20]. With the increase of the element size, the FE models with elements of sizes “D” and “E” as shown in Fig. 10(d) and (e) became unable to accurately reflect the microscopic architecture of the MCP bone because the relative area coefficient, c for these FE meshes were much higher than 1 as shown in Table 2. Thus, the stress and strain results of the corresponding FE analyses would be affected. In order to reduce the FE analysis error, a second grayscale intensity threshold T2 was introduced to eliminate the elements which had reduced the area coefficient c close to one. 3.2. Effect of varying the value of T2 The effect of varying second grayscale intensity threshold values T2 on the relative area coefficient c at three different cross-sections of the human MCP bone are determined as listed in Table 3. From this table, for the elements sizes “B” and “C”, the most appropriate value of T2 is 65 in order to reduce the value of c approaching to one and hence minimize the area errors of the reconstructed cross-sections of the bone, while the most appropriate value of T2 is 60 for the element size “D”. By using the relative area coefficient c with the value of 1, the total elemental volume of FE meshes constructed from the elements with sizes “C”–“E” could be kept to vary only within 5%. In order to compare the microscopic geometry of the human MCP bone generated by commercially available 3D reconstruction tool “Mimics” and the proposed method based on the obtained micro-CT data, three different regions containing trabecular bone, cortical bone,and both of them are selected and labeled as “P”, “Q” and “R”, respectively in Fig. 11(a). By using different element sizes and corre-
1478
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 10. FE meshes of the reconstructed human MCP bone with varying element sizes (A–E) showing micro-architecture of the bone at three different cross-sections corresponding to micro-CT sliced image numbered as 300 (left), 890 (middle) and 1480 (right), (a) 34.66 µm×34.66 µm×34.66 µm, (b) 103.98 µm×103.98 µm×103.98 µm, (c) 138.64 µm×138.64 µm×138.64 µm, (d) 170.30 µm×170.30 µm×173.30 µm, and (e) 242.62 µm×242.62 µm×242.62 µm. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
sponding values of T2 , various FE meshes for three different regions of the bone were created as shown in Fig. 11(b–d). The FE mesh created using the element size “A” shown in Fig. 11(b) was used to serve as a reference for comparison with other meshes. It can be observed from Fig. 11 that the resolution of the reconstructed bone regions was reduced with increasing element size. Therefore, too large element would not be recommended to be used for local microscopic stress analysis of the bone. Although it was the case, it can be observed that the overall geometry of some bone regions reconstructed using the larger sized elements of “B” and “C” (Fig. 11(c) and (d)) is similar from those with the element size “A”, because the shape errors of these bone regions have been minimized by using the T2 value.
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1479
Table 3 Effect of varying second grayscale intensity threshold values T2 on the relative area coefficient c at three different cross-sections of the human MCP bone Cross-section number (Ncs ) 300
890
1480
Element size (Se ) A B B B B B C D E A B B B B B C D E A B B B B B C D E
Second grayscale intensity threshold parameter (T2 ) 0 0 20 40 60 65 65 60 60 0 0 20 40 60 65 65 60 60 0 0 20 40 60 65 65 60 60
Number of elements used (n) 26,029 5288 4121 3391 2917 2836 1592 1019 520 20,057 2585 2433 2332 2265 2250 1272 839 448 24,040 4337 3591 3114 2737 2691 1533 1021 540
Area of an element (ab (mm2 )) 1af 9af 9af 9af 9af 9af 16af 25af 49af 1af 9af 9af 9af 9af 9af 16af 25af 49af 1af 9af 9af 9af 9af 9af 16af 25af 49af
Area of the bony region (Ab (mm2 )) Af = 26,029af 47,592af 37,089af 30,519af 26,253af 25,524af 25,472af 25,475af 25,480af Af = 20,057af 23,265af 21,897af 20,988af 20,385af 20,250af 20,352af 20,975af 21,952af Af = 24,040af 39,033af 32,319af 28,026af 24,633af 24,219af 24,528af 25,525af 26,460af
Relative area coefficient (c) 1.00 1.83 1.42 1.17 1.01 0.98 0.98 0.98 0.98 1.00 1.16 1.09 1.05 1.01 1.01 1.01 1.05 1.09 1.00 1.62 1.34 1.17 1.02 1.01 1.02 1.06 1.10
3.3. FE analysis In order to illustrate the applicability of the proposed method, a series of FE analyses of the MCP bone were conducted using ABAQUS software with a general personal computer (Mainboard: MSI P55GD55, CPU: Intel Core i5-750 2.67 GHz, RAM: 8 GB DDR3 1333). The effects of varying element sizes on the stress distribution and force–displacement relation of the bone under the same boundary and load conditions were investigated. Elements of three different sizes “C”, “D” and “E” with the corresponding threshold value T2 were used in constructing the FE meshes. For the FE model composed of the elements of size “C”, a total of 496,149 elements and 829,443 nodes are used. The corresponding results were plotted in Figs 12 and 13. It can be observed from Fig. 12 that the pattern of Von Mises stress distribution for different FE meshes is similar, although the maximum stress was higher in the FE model with the elements size “C”. From Fig. 13, the force–displacement data predicted by using the element size “C” lied between those with larger-size elements “D” and “E”, but closer to that with the element of size D. This implies that larger-sized elements could be used for the macroscopic stress analysis of the bone, but
1480
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 11. 3D reconstruction of three particular regions “P”, “Q” and “R” in human MCP bone based on micro-CT data of the bone, (a) 3D images in Mimics, (b)–(d) 3D FE meshes generated with variation of the values of Se , T2 and n: (b) Se = A, T2 = 0, (c) Se = B, T2 = 65, (d) Se = C, T2 = 65. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1481
Fig. 12. Von Mises stress contour plots of FE models of the MCP bone with three different sizes of elements and corresponding values of T2 , (a) Se = C, T2 = 65, n = 496,149, (b) Se = D, T2 = 60, n = 265,452 and (c) Se = E, T2 = 60, n = 96,917. (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
too large element like the size “E” is not recommended due to the large deviation in the displacement as compared with the smaller sized elements. 3.4. Model validation A tensile test that compares the force–displacement relationship of simulated and real experiments was conducted to validate the FE MCP bone model reconstructed by using our proposed methodology in reducing the number of elements. The trabecular and cortical bone tissues are assumed to have the same material properties. The Young’s modulus (E) and Poisson’s ratio (ν) of the bone tissues are assumed to be 18.7 GPa and 0.3, respectively [20]. Based on the results from Fig. 13, the element size “D” was used and the applied load was 40 N. As the cross-sectional area of the bone specimen is non-uniform, two additional measurement of the force–displacement relation of the specimen is made by rotating the specimen along the loading axis by 120◦ and 240◦ . The force–displacement relationship was illustrated in Fig. 14. With the increase in the displacement, the load generally showed an increasing trend for the three different testing positions. It can be observed that the trends of the predicted results were in good agreement with the experimental data, demonstrating the general applicability of the proposed method to the model validation.
1482
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
Fig. 13. Displacement–force relation of the MCP bone with three different sizes of elements and corresponding values of T2 . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/BME-130951.)
Fig. 14. Load–displacement plots for tensile testing of the MCP bone specimen.
4. Conclusion A methodology for automated generation of three-dimensional FE models of entire bone with microscopic architecture from micro-CT image data is successfully developed. The applicability of the
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
1483
proposed method is demonstrated using human MCP bone. The dual grayscale intensity threshold parameters are adopted to distinguish the bone tissues from the surrounding soft tissues and improve predictive accuracy for the FE analyses for different sizes of elements. With the proposed method, FE analyses of an entire bone can be run on a personal computer by using a smaller number of larger-size elements with maintaining certain predictive accuracy. Experiments are performed to simulate different sizes of FE element models. In addition, we attempted to validate the methodology by comparing the results from the FE analysis and the tensile test data. The simulation results show good agreement with the experimental results. It is believed that the proposed method can also be extended to model the bone with variation of mechanical properties internally and functional graded structures based on the average grayscale intensity within elements. Acknowledgement The work described in this paper is supported by a grant from the Research Committee of the Hong Kong Polytechnic University (Project No. G-YG18). References [1] R. Huiskes and E.Y.S. Chao, A survey of finite element analysis in orthopedic biomechanics: the first decade, J. Biomech. 16 (1983), 385–409. [2] T. Nomura, M.P. Powers, J.L. Katz and C. Saito, Finite element analysis of a transmandibular implant, J. Biomed. Mater. Res. B – Appl. Biomater. 80 (2007), 370–376. [3] M. Viceconti, L. Bellingeri, L. Cristofolini and A. Toni, A comparative study on different methods of automatic mesh generation of human femurs, Med. Eng. Phys. 20 (1998), 1–10. [4] O. Panagiotopoulou, N. Curtis, P.O. Higgins and S.N. Cobb, Modelling subcortical bone in finite element analyses: A validation and sensitivity study in the macaque mandible, J. Biomech. 43 (2009), 1603–1611. [5] F.A. Bandak, R.E. Tannous and T. Toridis, On the development of an osseo-ligamentous finite element model of the human ankle joint, Intl. J. Solids Struct. 38 (2001), 1681–1697. [6] P. Messmer, F. Matthews, A.L. Jacob, R. Kikinis, P. Regazzoni and H. Noser, A CT database for research, development and education: concept and potential, J. Digital Imaging 20 (2007), 17–22. [7] J.A. MacNeil and S.K. Boyd, Bone strength at the distal radius can be estimated from high-resolution peripheral quantitative computed tomography and the finite element method, Bone 42 (2008), 1203–1213. [8] D.S. Barker, D.J. Netherway, J. Krishnan and T.C. Hearn, Validation of a finite element model of the human metacarpal, Med. Eng. Phys. 27 (2005), 103–113. [9] A. Kato and N. Ohno, Construction of three-dimensional tooth model by micro-computed tomography and application for data sharing, Clin. Oral Invest. 13 (2009), 43–46. [10] J.H. Keyak, J.M. Meagher, H.B. Skinner and C.D. Mote, Automated three-dimensional finite element modelling of bone: a new method, J. Biomed. Eng. 12 (1990), 389–397. [11] J.H. Keyak and H.B. Skinner, Three dimensional finite element modeling of bone: effect of element size, J. Biomed. Eng. 14 (1992), 483–489. [12] J.H. Keyak, M.G. Fourkas, J.M. Meagher and H.B. Skinner, Validation of an automated method of three dimensional finite-element modeling of bone, J. Biomed. Eng. 15 (1993), 505–509. [13] G. Kochi, S. Sato, T. Fukuyama, C. Morita, K. Honda, Y. Arai and K. Ito, Analysis on the guided bone augmentation in the rat calvarium using a microfocus computerized tomography analysis, Oral Surg. Oral Med. Oral Pathol. Oral Radiol. Endod. 107 (2009), e42–e48. [14] C. Efeoglu, S.E. Fisher, S. Erturk, F. Oztop, S. Gunbay and A. Sipahi, Quantitative morphometric evaluation of critical size experimental bone defects by microcomputed tomography, Brit. J. Oral. Maxillofacial Surg. 45 (2007), 203–207. [15] G. Bevill and T.M. Keaveny, Trabecular bone strength predictions using finite element analysis of micro-scale images at limited spatial resolution, Bone 44 (2009), 579–584. [16] J. Peng, C.Y. Wen, A.Y. Wang, Y. Wang, W.J. Xu, B. Zhao, L. Zhang, S.B. Lu, L. Qin, Q.Y. Guo, L.M. Dong and J.M. Tian, Micro-CT-based bone ceramic scaffolding and its performance after seeding with mesenchymal stem cells for repair of load-bearing bone defect in canine femoral head, J. Biomed. Mater. Res. B – Appl. Biomater. 96 (2011), 316–325.
1484
C.Y. Tang et al. / Voxel-based approach to generate entire human metacarpal bone
[17] S.J. Hollister, J.M. Brennan and N. Kikuchi, A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress, J. Biomech. 27 (1994), 433–444. [18] B.V. Rietbergen, H. Weinans, R. Huiskes and A. Odgaard, A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models, J. Biomech. 28 (1995), 69–81. [19] E. Verhulp, B.V. Rietbergen, R. Muller and R. Huiskes, Indirect determination of trabecular bone effective tissue failure properties using micro-finite element simulations, J. Biomech. 41 (2008), 1479–1485. [20] G.L. Niebur, J.C. Yuen, A.C. Hsia and T.M. Keaveny, Convergence behavior of high-resolution finite element models of trabecular bone, J. Biomech. Eng. 121 (1999), 629–635. [21] G.L. Niebur, M.J. Feldstein, J.C. Yuen, T.J. Chen and T.M. Keaveny, High-resolution finite element models with tissue strength asymmetry accurately predict failure of trabeclar bone, J. Biomech. 33 (2000), 1575–1583. [22] S.V.N. Jaecques, L. Muraru, C.V. Lierde, E.D. Smet, H.V. Oosterywck, M. Wevers, I. Naert and J.V. Sloten, In vivo micro-CT-based FE models of guinea pigs with titanium implants: an STL-based approach, Inter. Cong. Ser. 1268 (2004), 579–583.
Copyright of Bio-Medical Materials & Engineering is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
