Automated segmentation of swine skulls for finite element model creation using high resolution µ-CT images

Automated segmentation of skulls for finite element model creation high Abstract: Finite element (FE) analysis has been widely used to investigate bone responses to mechanical loading. Research in long bones has been straight forward because modeling of these bones requires only two material properties. Such an FE model may provide an adequate approximation of the anatomy for many cases. However, a more detailed model of skull bones is needed to accurately capture its complex structure of multiple bone pieces and the various mineral densities distributed throughout these bone pieces. Unfortunately, FE model development incorporating both complex geometries and anatomically accurate material properties is both computationally and labor intensive. In this study, a method is proposed to automatically segment micro-computed tomography (µ-CT) scan images of bone pieces to build an FE model of a full swine hemi-skull. Using the Digital Imaging and Communications in Medicine (DICOM) files from scanned bones, the complete geometry of each bone piece is recreated through seven customized processing algorithms. After assembling the bone pieces to form the skull, experimentally derived Young's modulus values are correlated to grayscale values to produce a detailed FE model for accurate simulation. This detailed skull model can be used to predict strain/stress patterns in response to various loading regimes to facilitate research questions in fracture healing and growth, as well as bone tissue engineering and bone mineral density loss (e.g. osteoporosis). Abstract Finite element (FE) analysis has been widely used to investigate bone responses to mechanical loading. Research in long bones has been straight forward because modeling of these bones requires only two material properties. Such an FE model may provide an adequate approximation of the anatomy for many cases. However, a more detailed model of skull bones is needed to accurately capture its complex structure of multiple bone pieces and the various mineral densities distributed throughout these bone pieces. Unfortunately, FE model development incorporating both complex geometries and anatomically accurate material properties is both computationally and labor intensive. In this study, a method is proposed to automatically segment micro-computed tomography (µ-CT) scan images of bone pieces to build an FE model of a full swine hemi-skull. Using the Digital Imaging and Communications in Medicine (DICOM) files from scanned bones, the complete geometry of each bone piece is recreated through seven customized processing algorithms. After assembling the bone pieces to form the skull, experimentally derived Young’s modulus values are correlated to grayscale values to produce a detailed FE model for accurate simulation. This detailed skull model can be used to predict strain/stress patterns in response to various loading regimes to facilitate research questions in fracture healing and growth, as well as bone tissue engineering and bone mineral density loss ( e.g. osteoporosis).

Abstract: Finite element (FE) analysis has been widely used to investigate bone responses to mechanical loading. Research in long bones has been straight forward because modeling of these bones requires only two material properties. Such an FE model may provide an adequate approximation of the anatomy for many cases. However, a more detailed model of skull bones is needed to accurately capture its complex structure of multiple bone pieces and the various mineral densities distributed throughout these bone pieces. Unfortunately, FE model development incorporating both complex geometries and anatomically accurate material properties is both computationally and labor intensive. In this study, a method is proposed to automatically segment microcomputed tomography (µ-CT) scan images of bone pieces to build an FE model of a full swine hemi-skull. Using the Digital Imaging and Communications in Medicine (DICOM) files from scanned bones, the complete geometry of each bone piece is recreated through seven customized processing algorithms. After assembling the bone pieces to form the skull, experimentally derived Young's modulus values are correlated to grayscale values to produce a detailed FE model for accurate simulation. This detailed skull model can be used to predict strain/stress patterns in response to various loading regimes to facilitate research questions in fracture healing and growth, as well as bone tissue engineering and bone mineral density loss (e.g. osteoporosis).

Abstract
Finite element (FE) analysis has been widely used to investigate bone responses to mechanical loading. Research in long bones has been straight forward because modeling of these bones requires only two material properties. Such an FE model may provide an adequate approximation of the anatomy for many cases. However, a more detailed model of skull bones is needed to accurately capture its complex structure of multiple bone pieces and the various mineral densities distributed throughout these bone pieces. Unfortunately, FE model development incorporating both complex geometries and anatomically accurate material properties is both computationally and labor intensive. In this study, a method is proposed to automatically segment micro-computed tomography (µ-CT) scan images of bone pieces to build an FE model of a full swine hemi-skull. Using the Digital Imaging and Communications in Medicine (DICOM) files from scanned bones, the complete geometry of each bone piece is recreated through seven customized processing algorithms. After assembling the bone pieces to form the skull, experimentally derived Young's modulus values are correlated to grayscale values to produce a detailed FE model for accurate simulation. This detailed skull model can be used to predict strain/stress patterns in response to various loading regimes to facilitate research questions in fracture healing and growth, as well as bone tissue engineering and bone mineral density loss (e.g. osteoporosis).

Introduction
Mechanical stimulation alters the patterns of growth, remodeling, and fracture repair in living bones. (Hambli, 2014) Investigations aimed at elucidating the mechanisms of bone responses have utilized finite element (FE) analyses (Liu and Quek, 2013) to predict strain patterns resulted from mechanical stimulations (Crawford et al., 2003;Gautam et al., 2007;Geng et al., 2001;Metzger et al., 2005) (Hambli, 2014;Pottecher et al., 2016;Yousif and Aziz, 2012). Limitations have surfaced, however, as differences in the assumed model parameters (Kim et al., 2012;Rho et al., 1998;Rho et al., 1997) and conflicts between model and experimental results have been reported (Crawford et al., 2003). FE models will fail to faithfully represent in vivo strains if the models oversimplify the geometry and/or incorrectly assume mechanical properties (e.g., elasticity) of the bones. This is most challenging for complex bone structures such as skulls, where loading directions and magnitudes are multifactor, geometries are intricate, and internal distributions of material properties can vary greatly with anatomic location (Morgan et al., 2003).
Previous FE models in skulls have taken one of two strategies when balancing the competing factors of biological details, time for data collection, and computational efforts. Frequently, a representation of the bone is built from CT scans with 5-25 millimeter resolutions. This limits detail for materials with properties on the micrometer scale (e.g., trabeculae), forcing models to assume unrealistic anatomic characteristics (i.e., a skull built with a 1 mm thick shell of cortical bone; Field et al., 2010;Song et al., 2015). On the other hand, when including micrometer-level details (i.e., the distribution of cortical and trabecular bone), models accurately incorporate the different mechanical properties of these tissue types, but face other limitations. These limitations include: 1) time-consuming and tedious segmentation protocols needed to produce 3D reconstructions that give accurate distributions of bone types; 2) limited scan capacity and increased expense of µ-CT; and 3) a large computational burden to process FE models with possibly hundreds of millions of elements.
Another source of error, generally, stems from the widely-used assumption that skeletal structures will have only one or two material properties. Literature-derived or experimentally determined elasticities for "cortical" and "trabecular" bone are typically assigned, with no regard for the complex distributive nature of bone properties. In a few applications, apparent bone density has been approximated through the use of grayscale in CT scans (Zannoni et al., 1999), but this has not been done for large, complex skeletal elements, nor at resolutions that would demonstrate the distribution at the micrometer level (Kim et al., 2012), and lastly, not in a way that allows for a continual scale of material properties related to the amount of bone mineral present in a sample. This is important because 'bone type' is a categorical description of a continuous property and weakly mineralized cortical bone can have a less than or equal stiffness to a more densely packed and well-mineralized sample of trabecular bone (Morgan et al., 2003). It would be ideal if FE models could incorporate more of this continual spectrum rather than artificially categorizing all bones into just two types of materials.
In this paper, we propose an automated segmentation protocol and FE model creation technique that allow us to maintain micrometer scale resolution (as captured in µ-CT scans) for a large skeletal element: a juvenile porcine hemi-skull. We provide detailed instructions on how to assemble skeletal pieces in silico for use in further modeling, and outline the meshing protocols that provide a complete representation of 3D bone configuration. Algorithms were developed to automatically assign elastic modulus properties to elements in a distributive nature, based upon several experimentally derived Young's modulus measures from multiple anatomic locations. The presented techniques and algorithms will be applicable to skeletal research generally, and also for understanding the biological mechanisms underpinning repair rates in our own tissue engineering research (Caballero et al., 2015). The central hemi-skull was deconstructed into pieces to accommodate the bed limitations of the scanner. The frontal (F), vomer (R), turbinate (T), premaxilla (P), maxilla (X), ethmoid (E), lacrimal (L), and malar (M) bones were disarticulated along the sutures. The brain encasing portion of the skull was cut (red line) into dorsal (D) and ventral (V) pieces. The styloid process was not included for the completed model.

Acquisition of skeletal images
The right half of the skull of a cadaveric Yorkshire farm pig (aged 6 weeks) was skeletonized and separated into 10 pieces ( Figure 1) so that it could physically fit into the bed of a µ-CT scanner. Only the right side of the skull was used, because significant bilateral differences were not expected. High resolution images were acquired in air by the Imaging Research Core at the Cincinnati Children Medical Center (CCHMC) using an ImTek Micro CT with a resolution of 52.06 µm. DICOM files for each piece were saved, which will be used as inputs for our processing to create the FE model.

2.2
Algorithms for processing grayscale data Figure 2. Algorithms applied on a DICOM file and the procedure of the proposed image processing. A) Original DICOM image from a portion of the parietal bone, demonstrating the complex internal architecture of the mineralized tissue. B-H) Sequential images produced by our processing algorithms applied to (A), based on the procedures described in the text.
Each DICOM file was imported into python using an open source package (PYDICOM http://pydicom.readthedocs.io/en/stable/getting_started.html) as a two-dimensional (2D) grayscale image, which clearly illustrates the distribution of mineralized tissue in the slice ( Figure 2A). To convert the grayscale values in the DICOM files to an outline for production of a meshed bone piece, images were processed with sequential algorithms. Smoothing algorithms were designed to provide a balance between capturing the complexity of bone and producing a surface that would require reduced computational time. The images were first subjected to grayscale dilation, a trimmed version was next produced, followed by a Gaussian smoothing operation. Then pixel thresholds were determined, converted to binary values, and the resultant image was subjected to a final closing algorithm to ensure the outline was free of gaps.
As a first step to remove gaps in the external surface and reduce image noise from nonbone regions, a grayscale dilation was performed using (eq. 1), a clipping algorithm produced a trimmed version of the image through operation defined in (eq. 2), and a Gaussian smoothing algorithm (eq. 3) was used to simplify the space between trabeculae (Sternberg, 1986).
and represents grayscale dilatation and erosion, is a function of the grayscale values of the image with respect to different coordinates in 2D (x, y) space, and is the function of structuring element.
is the domain of the image, is the domain where is defined, and represents grayscale closing.
The structuring element needs to be adjusted to fill the spaces between trabeculae while not expanded to the point that the outline of the image is altered. In this case ( Figure 2B), the discrete structuring element was a matrix with all elements equal to 1. The initial result demonstrates areas with low grayscale values that may cause the outline to be inaccurate. To account for this, the Gaussian smoothing algorithm calculates a new grayscale number with a weighted average around its neighbor for each pixel by convolving the images. The discrete Gaussian smoothing operator (eq. 3) can be described as a matrix form: ⑷ where x is the distance from the center point of the Gaussian distribution in rows, y is the distance from the center point of the Gaussian distribution in columns, and is the standard deviation. Larger values of produce more blurred images. In theory, the size of the Gaussian smoothing operator matrix is infinite, and all the values are positive. To ensure this in our code, any number less than in the matrix was truncated. After applying the Gaussian smoothing operator, the image becomes more uniform and is acceptable for the next step (Error! Reference source not found.C).
The grayscale closing algorithm (eq 2) is applied ( Figure 2D) to trim the image to provide the positive relief needed to produce a surface outline. A subtracting algorithm is applied:

⑸
The result (Error! Reference source not found.E) now requires a quantified threshold to determine which pixels should be preserved as outline node information. Attention to detail is needed in determine the threshold value. If the threshold is set too high, the outline of low grayscale areas may be neglected leading to gaps in the outline. If the threshold is too low, excess pixels will be selected, complicating how the images are stacked together for the meshes in the subsequent steps. Another complication relates to how the grayscales in each of the image slices vary across different mineral densities. Due to this, it is necessary to define a unique threshold for each image. Therefore, an additional simplification converted the grayscale images into binary images (0 or 1) to generate an outline ( Figure 2F) using: ⑹ ⑺ where and are the user defined thresholds ( and ), is the maximum grayscale in each image, whereas is the minimum grayscale in the image. Finally, a thinner outline, the closing algorithm (eq. 2) was applied ( Figure 2G) and produced an image with a thinner node cloud needed for the completion of each bone piece model.

Modeling each bone piece
The 2D node clouds for each bone piece were stacked sequentially at a distance of along the z-axis to match the resolution of the µ-CT scanner. Then, all the node coordinates were imported to Geomagic (3D Systems, Inc.) and image noise and unnecessary support frames were manually removed prior to constructing the modeled bone piece. Detailed trabecular structures were connected where possible and minor canals for blood vessels and nerves of less than were smoothed. Files composed of the stacked grayscale elements were exported as Initial Graphics Exchange Specification (iges) files and exported into Hypermesh (Altair Engineering, Inc.). Meshes were built for each bone piece using tetrahedral elements ( Figure 3A) and the grayscale of each element was calculated by averaging all the grayscale values inside the tetrahedron within the 2D images ( Figure 3B).

Assembling the bone pieces
Rigid body rotation and transformation for each meshed bone piece were calculated, because they needed to be combined into one integrated skull, and all the meshed nodes of each piece needed to be transferred into the global coordinates shown in Error! Reference source not found.C. This procedure is described in the following paragraph. All the bone pieces were imported into Geomagic and transformed to their original places using the Object Mover function. After finishing the translation and rotation of each bone piece in stl file form, all the suture lines should be merged, and exported to an iges file for further meshing. In addition, all the transformation and rotation matrices were exported from Geomagic. Then, the integrated model was meshed with a total of 176,058 T4 (tetrahedrons with 4 nodes) elements with an element size of in Hypermesh, as shown in Figure 3D. For assigning the grayscale, the nodes and their grayscale data from each individual piece were also needed to form the entire skull, as a new node data file corresponding to the meshed model. For acquiring the new global coordinates of the nodes, a translation vector and a rotation matrix should be applied on the original coordinates, which can be expressed by the following equation: represents the node coordinates when they are meshed, and stands for the new global coordinates of each nodes. R is the rotation matrix for each bone piece, and is the translation vector.

Calculation of full grayscale model
For acquiring the grayscale of the integrated model, for each element, the grayscale of node data within each tetrahedron was averaged. The result is shown in Tecplot (Tecplot, Inc.). From Error! Reference source not found.-E, we can observe that teeth are at the highest position of the scale, which means they have the highest density. Comparatively, the dark blue parts, where openings in bone or cartilage, are present, have the lowest density.

Mechanical testing of bone (Young's Modulus)
The relationship between the bone's apparent density and grayscale average in the region of interest can be calculated using the following formula. (Morgan et al., 2003) ⑼ where is apparent bone density, is the average grayscale value from a CT (or µ-CT) image, and and are constant parameters depending on the exposure conditions. The density correlation to Young's Modulus is . ⑽ where and depend on the bone locations. (Zannoni et al., 1999) Combining these two equations, the relationship between Young's Modulus and the grayscale value can be expressed as a polynomial. Such a polynomial curve can be fitted using experimental data. In our case, a 2 nd power polynomial was used in the curve fitting. An indent experiment of measuring the Young's Modulus was conducted by both the left and right side of the teeth, hard palate and parietal bone in the front, middle and back. The average elastic modulus for tissue types were significantly different (p < 0.0001), with teeth equal to 47311.47 ± 7030 MPa; hard palate was 11580.94 ± 3910.30 MPa; parietal bone with 3 different locations were 13095.21 ± 4764.90 MPa (frontal), 11589.04 ± 805.75 MPa (middle) and 6935.10 ± 1946.38 MPa (back). A post-hoc comparison of tissue type indicated no difference between the frontal bone and the hard palate (p = 0.127), but that both of these tissues were significantly less elastic than the teeth (p < 0.0001). The general linear model demonstrated no significant differences among pigs (p = 0.546) or between the left and right sides (p = 0.567). For the curve fitting, the grayscales were sampled multiple-times in the corresponding location where they were tested. The Young's Modulus data and corresponding grayscales are summarized in the following table.

Gray scale modeling
In µ-CT scans, x-ray attenuation is largely dependent upon tissue density and is represented visually through a black-to-white scale by assigning pixels a gray color based upon the percentage of x-rays absorbed at that location. At a micrometer scale, a structural difference between trabecular and cortical bone is related to the distribution of mineralized and nonmineralized tissue: trabecular bone has dense trabeculae in a soft tissue matrix and cortical bone is more uniformly mineralized. Our protocol generates elements (200 -400 µm in size) based upon the average grayscale in pixels (20 -150 pixels, depending upon the resolution of the scan, Figure 2). After the elements are defined, a more visually interpretable heat map is generated from the mean grayscales where tissues like cartilage and fenestra (for neural-vascular bundles) expressed by the lowest grayscales in blue, while dense tissues, like teeth, are red. This was accomplished by generating a meshed model for each bone segment to provide a framework for defining grayscale. The original scans from which these elements are defined represent the 3D characteristics of each bone, and therefore provide a realistic model of the internal architecture representing the distribution of both trabecular and cortical bone, in addition to any endosteal spaces, with greater detail than previously available (Figure 3).
To overcome the bed size limitations of the high resolution µ-CT scanner, it was necessary to physically separate smaller pieces of the skull for initial data collection. Our work represents the first time a skull or areas of bone larger than 50mm (Fajardo and Müller, 2001) (Tuan andHutmacher, 2005) (de Crespigny et al., 2008) have been fully modeled to capture the distribution of cortical and trabecular bone. This was accomplished through use of a novel script to integrate bone pieces to overcome the limitation of sample size.

Calculation of Young's modulus
We have generated a mathematically determined correlation between measured mechanical properties of bone and the modeled gray scale (Figure 4). The relationship between measured elasticity and grayscale is described through a 2 nd order polynomial that shows correspondence with others' work. (Morgan et al., 2003) (Cong et al., 2011;Roberts and Garboczi, 2000) The polynomial is shown as: ⑾ is described as the Young's Modulus and is grayscale. The coefficient is 0.991.

Mechanical properties of full skull
The distribution of the elastic Young's modulus in the pig skull was determined using our automated procedure. It provides a more detailed and biologically valid representation of the skeletal anatomy, including internal architecture previously missing in published models ( Figure  5) (Buie et al., 2007). Cortical bones are not only shown as a 2 mm homogeneous cortical shell, but a detailed distribution of Young's modulus, which means that the Young's Modulus should be different with respect to different areas in cortical layer. (Bright and Rayfield, 2011) Cavities where nerves and cartilage should be present can also be described properly in this model by a low Young's Modulus. Teeth can be shown as the highest Young's Modulus in this model, which also have different mechanical properties on the enamel and dentin areas.

Discussion
Commercial software (Fischerauer et al., 2013;Jaecques et al., 2004) and some in-house software (Bouxsein et al., 2010;Hildebrand and Rüegsegger, 1997;Rüegsegger et al., 1996) (Niebur et al., 2000) have been used for creating FE models using µ-CT scans. These techniques are difficult to use for large-sized bones, such as swine skulls, because they cannot fit into the small bed of µ-CT scanners. In this proposed method, we cut the large skull into many smaller pieces of bones, each of which can be hosted in a µ-CT scanner. Each bone piece is then scanned and stored in a DOCOM file. An image processing technique was developed to assemble all of these pieces to construct an FE model with the desired number of elements. The proposed algorithm ensures that each of the elements will be assigned with a proper grayscale value (and then a Young's Modulus value) from the high resolution µ-CT scan images. The parameters and thresholds used in our algorithms were optimized for complex swine skulls. In addition, the threshold shown in eq-6 was set up and optimized to capture trabeculae on the surface in which the grayscale is lower than normal, and to clear background noise using the created model. To improve our algorithms, we also tested the following imaging processing techniques.
First, we tested another common method, known as the Sobel operator, for calculating the boundary lines of the bone pieces. The Sobel operator has usually been used in a variety of situations to determine edges of pictures (Kanopoulos et al., 1988). The normal Sobel operator can be expressed as follows, ⑿ ⒀ in which A represents the origin plot data, * is a convolution operator, is the gradient in the x direction, and is the gradient in the y direction. The magnitude value of the gradient in each pixel is then calculated by . ⒁ Figure 6. The image processing result using Sobel operator instead of image subtracting.
If the subtraction procedure is substituted by the Sobel operator to determine the outline, and the same procedure is followed to the end of image processing, the result is shown in Figure  6, which elucidates greater detail of the inner structure of the bone slice. Thus, the Sobel operator can also be used to improve calculations.
Secondly, when using the Gaussian filter, the outlines of each bone piece may be expanded by a certain number of pixels, depending on the order of the Gaussian filter used. In our case, the resolution of the scans (distance between two adjoining pixels) is , and outlines can be expanded by 10 pixels. The grayscale value computed in an element on the external surface of bone pieces will be lower than the real bone surface because of the artificial expansion. Users of our method need to take note of this error.
Third, we noted that the grayscales in the µ-CT files represent the resorption capacity of calcium, which correlates to the real densities of trabecular or cortical bone at their specific positions. Hounsfield Units (HU) are used to express these attenuation values, which define water with an attenuation value (HU) of 0 (Hsieh, 2009). In some areas, the grayscale of trabecular bone is higher than the nearby cortical bone, which means that the real density of that specific trabeculae is higher than the cortical bone. However, due to the low volume of the trabeculae, the apparent density and Young's Modulus is usually lower than cortical bone. In this case, the apparent density of each element is acquired from averaging all the pixels of real bone grayscales and the intervals. Some research groups use two different equations to relate HU and apparent density depending on cortical or trabecular bone (Cattaneo et al., 2001). HU is used to determine a threshold to differentiate cortical and trabecular bone. However, as the CT scan resolution is increased with µ-CT data, HU is not a criterion to distinguish cortical and trabecular bone. Note that the assumptions of cortical and trabecular types (Bright and Rayfield, 2011) and a uniform thickness of cortical shells can result in much stiffer models. In our study, assigning Young's Modulus only depends on the apparent density of each element. There is no need to consider whether the bone is cortical or not, and we base our predictions on the grayscale that may vary continuously. Therefore, our model is a much closer approximation of reality due to the high resolution µ-CT images.
Fourth, we note that all the elements are simplified as isotropic in this study. Generally speaking, the mechanical properties of bones are anisotropic.
Fifth, Young's Modulus obtained experimentally can vary significantly due to the complex nature of bone structure and many other unquantifiable factors. This causes large standard deviations in the input material data for FE models and may lead to inaccuracies when fitting the curve with polynomials. To solve this problem, the average Young's Modulus is used and correlated to grayscale.
Finally, we note that our algorithms are applicable not only to swine skull models, but also to the skulls of other similar species (Sigal et al., 2008). In addition, we may be able to use the FE model of one pig to estimate other pigs, by morphing the technique using least square optimization, to increase the efficiency of specimen-specific analysis. (Taddei et al., 2007) Manual landmarks have historically been needed for morphing and mapping skull distribution. Once we have determined the Young's Modulus distribution of bone, corresponding nodes can be selected as landmarks from both the meshed model and target model. Using morphing and warping algorithms to change the meshed model into the shape of the target model, we will be able to preserve the Young's Modulus information in each element.

Conclusions
We are now able to draw the following conclusions: 1) We have presented a novel automated method with a set of computed algorithms for image processing to build a detailed whole hemi-skull FE model with µ-CT level ( ) resolution, which was able to overcome the limitation of scanner bed size. 2) The Young's Modulus used in the elements are derived directly from the grayscale value provided in the DICOM file of the µ-CT image. In theory, we can even assign each individual element a distinct Young's Modulus. Therefore, we no longer need to artificially classify bone as just two types, trabecular bone and a cortical shell.
3) The process is automated, easy to use, and can create FE models of complicated fullsized skulls. 4) The created FE model becomes more realistic, both in terms of geometry representation and material property assignment, because it makes the fullest use of high resolution µ-CT image data in an automated manner.
In future experiments, we plan to study in greater detail the correlations of grayscale and Young's Modulus in swine skulls. Using the FE model and the Grayscale-Young's Modulus correlations, we will also conduct FE analysis by adding muscle loads. We will also further develop bone remodeling algorithms. In addition, once our pig skull studies are concluded, we will perform a similar set of studies on human skulls.