Finite Element Analysis of Infant
Skull Trauma using CT Images
Finite Element Analysis of Infant
Skull Trauma using CT Images
A R N A Ó S K A R S D Ó T T I R
Master’s Thesis in Scientific Computing (30 ECTS credits)Master Programme in Computer simulation for Science and Engineering 120 credits
Royal Institute of Technology year
2012Supervisor at KTH was Svein Kleiven
Examiner was Michael Hanke TRITA-MAT-E 2012:07 ISRN-KTH/MAT/E--12/07--SE
Abstract
Some cases of infant skull fracture fall under the category of forensic study where it is not obvious whether the head trauma happened due to an accident or abuse. To be able to determine the cause of the head trauma with sufficient accuracy, biomechanical analysis using finite element modeling of the infant cranium has been established. By simulating the trauma, one may be able to obtain the fracture propagation of the skull and from it determine if the scenario narrative is plausible. Geometry of skull, sutures, scalp and brain of a 2 month old infant head was obtained using CT images and meshed using voxel hexahedral meshing. Simulation of an impact to the head from a fall of 0.82 m height, to a rigid floor, was carried out in the non-linear finite element program LS-Dyna. Two scenarios were simulated: an impact to the occipital-parietal bones and an impact to the right occipital-parietal bone. The fracture propagation was obtained using the Chang-Chang Composite Failure Model as a constitutive model for the skull bones. The amount of material parameters gathered in the present study to predict fracture of the infant skull has not been obtained before, to the best knowledge of author.
Validation of the models’ ability to show relatively correct fracture propagation was carried out by comparing the obtained fracture pattern from the parietal-occipital impact against published fracture patterns of infant PMHS skulls from a free fall onto a hard surface. The fracture pattern was found to be in good compliance with the published data. The fracture pattern in the parietal bone from the impact was compared against a fracture pattern from a previously constructed model at STH. The patterns of the models show some similarities but improvements to the model and further validations need to be carried out.
Sammanfattning
Några skallskador hos spädbarn ger grund till kriminaltekniska studier där det inte är självklart om skallskadan skett på grund av en olycka eller misshandel. För att kunna fastställa orsaken till skallskadan med tillräcklig noggrannhet har biomekaniska analyser med finita element modeller av barns huvud genomförts. Genom att simulera traumat kan man kunna få sprickpropagering i skallbenet och från den avgöra om scenariot är rimligt. Geometrin för skallen, suturer, hårbotten och hjärnan hos ett 2 månader gammalt spädbarns huvud erhölls genom CT-bilder och Voxel hexahedermeshning. Simulering av påverkan på huvudet från ett fall på 0,82 m höjd mot ett hårt golv simulerades i det icke-linjära finita element programmet LS-Dyna. Två scenarier simulerades: ett islag mot nack-hjässbenet och ett mot det högra hjässbenet. Sprickpropagering simulerades med en Chang-Chang Composite konstitutiv frakturmodell för skallbenet. Den omfattande mängd materialparametrar som sammanfattades i denna studie för att prediktera skallbensfrakturer hos spädbarnets har, enligt författarens kännedom, inte erhållits tidigare. Validering av modellernas förmåga att visa relativt korrekt sprickpropagering genomfördes genom att jämföra det erhållna frakturmönstret från simuleringarna med publicerade frakturmönster från spädbarn för fritt fall mot en hård yta mot nack-hjässbenet. Frakturmönstret befanns vara i god överensstämmelse med publicerade data. Brottmönstret i hjässbenet jämfördes med frakturmönstret från en tidigare konstruerad modell på KTH. Brottmönstren från modellerna visar vissa likheter men förbättringar av modellen och ytterligare valideringar måste genomföras.
1.2 Background ... 2
1.3 Aim of the project ... 4
1.4 Anatomy of an infant cranium ... 4
1.4.1 The skull bones ... 5
1.4.2 Sutures and fontanelles ... 7
1.4.3 Dura mater ... 7
Biomechanical & material properties of infant skull and suture ... 9
2.1 Bone and suture elastic behavior ... 9
2.1.1 Isotropic material properties ... 10
2.1.2 Transversely isotropic material properties ... 12
2.2 Parameter selection... 14
2.2.1 Parameters obtained from literature ... 15
2.2.2 Parameters estimations ... 18
Creating Finite Element mesh from CT data ... 20
3.1 Introduction to computed tomography imaging ... 20
3.2 Image segmentation ... 21 3.2.1 3D Slicer ... 22 3.2.2 DeVIDE ... 25 3.2.3 Matlab ... 27 3.3 Hexahedral meshing ... 28 3.3.1 Dicer ... 28 3.3.2 Mapping ... 29 3.3.3 Mesh quality ... 32
Finite element Modeling ... 35
4.1 Finite element simulation in LS-Dyna ... 35
4.1.1 Governing equation ... 35
4.1.2 Discretization: Hexahedral element and its shape function ... 36
4.1.3 Time integration ... 38 4.1.4 Matrix calculation... 39 4.1.5 Volume integration ... 41 4.1.6 Hourglass control ... 42 4.2 Simulation in LS-Dyna ... 42 4.2.1 Model preparation ... 43
6.3 Models mesh ... 62
6.4 Model simulation ... 63
Conclusion ... 65
Appendix ... 66
Constitutive model for the brain ... 66
Chapter 1
Introduction
Studies on cranial impact have been done for different scenarios such as head molding during birth, bike and car accidents and also in forensic cases. Head injuries caused by falls are in majority when it comes to infants and juveniles [1] showing the importance of understanding the biomechanics of their head to predict the occurring of the injury [2]. Forensic studies on infant skull fracture may be done to determine if a head trauma was caused by an accident or abuse. To be able to determine the cause of the head trauma with sufficient accuracy, biomechanical analysis using finite element (FE) modeling of the infant cranium can be established. A model of the cranium can be constructed from the patient specific computed tomography (CT) images and used to simulate the head trauma. By simulating the trauma, one may be able to obtain the fracture propagation of the skull and consequently determine if the scenario narrative is plausible. Biomechanical studies of this kind can also help to determine the injury mechanism of infant skull bones from different impact scenarios.
1.1
Thesis overview
The thesis covers the whole process of simulating an infant skull trauma from a fall. Patient specific CT image data of a 2 month old infant cranium is provided. The image data is processed to obtain the complex geometry of scalp, skull, sutures and brain and further processed to construct a hexahedral mesh in the (un-commercial) program Dicer, resulting in a model suitable for FE simulation. The patient specific model is simulated in the non-linear Finite Element program LS-Dyna to obtain an analysis of a mechanical impact to the cranium from a fall of a certain height. Results from an impact to the parietal bone and parietal and occipital bones are analyzed. Validation of the models ability to show relatively correct fracture propagation is carried out by comparing the obtained fracture pattern from the parietal-occipital
code of the meshing program, Dicer, was used with the permission of its author. The commercial program LS-Dyna was accessed at STH.
The thesis material is organized as follows. Chapter 1.2 covers previous studies found in literature relevant for the background of this study. Chapter 1.3 goes through the aim of the project and chapter 1.4 covers the anatomy of the human infant skull, sutures and dura mater. Chapter 2 covers the biomechanical and material properties of infant skull and suture. Chapter 2.1 goes through the elastic behavior of skull and sutures and 2.2 explains how parameters for the skull bones were obtained and estimated. Chapter 3.1 goes through the image processing and chapter 3.2 the hexahedral meshing. Chapter 4 covers the simulation process in LS-Dyna. Chapter 5 shows results from impact simulations, chapter 6 covers discussions and chapter 7 lists the conclusions.
1.2
Background
Studies on the biomechanics of fetal heads go back at least to the year 1888 when Murray et al. studied the effects of compression on fetal skull during labor [4]. Until the year 1980 few other studies were done on the topic of the biomechanics of head molding during labor, giving more qualitative then quantitative results [5]. The mechanical complexity of the fetal heads structure was well understood at this time and explained by Holland at al. in 1922 as follows: “The fetal head consists of a non-rigid-shell, of a peculiar shape, composed of a loosely-joint plates of liable bone jointed on a rigid base and strengthened by the attachment of dura mater and its system of septa: the contents of this shell are partly solid and partly fluid” [6] [5].
McPherson and Kriewall, 1980, were the first ones to do mechanical studies on tissue properties of the fetal cranium. They investigated the fiber orientation of specimens of varying locations from 6 cranial cadavers ranging in gestational age. They found that the fiber pattern of each bone of the skull were radially oriented from the center of ossification [5] versus homogeneous grain structure in adults [7]. The elastic modulus was also determined by them using three point bending test for specimens with fiber direction running parallel and perpendicular to the long axis. Elastic modulus was found to differ significantly between parallel and perpendicular fiber direction and also to increase with increased gestational age [5]. Further research by Kriewall at al. confirmed that the elastic modulus for fibers running parallel to the long axis had higher elastic modulus than fibers running perpendicular to it. They also studied the stiffness of the parietal bone by applying load to the inner surface of the eminence with the bone placed concave upon a tripod. High correlation was found between the bone stiffness and gestational age and bone stiffness and fetal weight [8].
Mazuchowski and Thibault at al. studied the mechanical properties of human skull and porcine sutures in 1997 and constructed a FE biofidelic model to simulate the deformation and flexibility of skull and sutures respectively and its biomechanical function during head impact. The cranial vault was simplified to a hemi-elliptical shell model with sutures of 5 mm constant thickness and anterior fontanels as a 15 mm square plate. The cranial cavity was assumed to be filled with brain. They concluded that geometry such as varying bone thickness and mechanical properties such as elastic modulus of oriented fibers must be taken into account when developing models of infant cranium. They were not successful in their testing of sutures but
tolerance to impact due to their flexibility [2] . A continuation of this research was done in 2000 using sutures from infant porcines by Margulies at al. [9]. They constructed two FE models with the size of a one month old infant head in the form of ellipsoid using shell elements to model the skull, sutures, fontanel and foramen magnum and brick elements for the brain. Cranial bone and sutures were assumed to have homogeneous properties. The models, one with parameters for adults and fused fontanels and sutures and another with infant porcine and human properties for sutures and cranium respectively were simulated for the contact problem in LS-Dyna3D. Their models were fixed at the base of the cranial vault and impact load applied to the parietal bone. They concluded that the structural change from the one layered compact bones in infants to the three layered cortical bones in adults were an important factor when determining the mechanical response of an impact to the head. Elastic modulus and ultimate stress were found to increase with age both for the human cranial bone and porcine sutures [9] [10] [11].
In 2006 Coats and Margulies addressed the hypothesis if some mechanical properties of the infant skull were bending and tensile rate dependent. Their studies on elastic modulus and ultimate stress of infant cranial bone and infant sutures were obtained from specimens of 23 cadavers, taken with fiber direction perpendicular to the long axis, of age less the one year old. The bending rate of 2.4 m/sec or higher in three point bending test were chosen to resemble a skull response to fall from low heights [12]. Earlier studies on elastic modulus of infant cranial bones were done at a much lower rate, 8.3⋅10−6m/s, to simulate the skull response to uterine contraction during labor [5] [9]. They came to the conclusion that strain rate of bending and tension test did not have significant effect on the elastic modulus and ultimate stress of bone and sutures for the sample age [12]. This was in contradiction to earlier studies done on human adult cortical bone [13].They hypothesized that future studies on specimens, taken with fiber direction parallel to the long axis, would show rate dependence [12]. They found that ultimate stress and elastic modulus were larger for parietal bones than occipital bone and both increased with age. No relation of elastic modulus and ultimate stress was found for age and strain rate for the human sutures. Ultimate strain for cranial bone was found to be significantly lower than for sutures and 23 times stiffer. They also reported that human sutures were not comparable to porcine sutures as was hypothesized in their earlier studies [9]. They concluded that the skull should not be considered homogeneous, as a result of their studies on biomechanical properties, and emphasized the need to use relevant age specific material properties for the skull and sutures to obtain accurate results when simulating skull traumas using computational models [12]. Coats and Margulies constructed a FE model from CT images of a 1.5 month old infant in 1997 and simulated an impact to the occipital bone from a low fall height. As their FE model was to analyze skull fractures of cranium, brain was assumed isotropic and homogeneous. They concluded that too stiff brain (large shear modulus) could increase the stress of the cranial bone. Suture geometry was also of importance and sutures width bigger than 10 mm showed decrease in maximum stress resulting in decreased fracture, emphasizing once again the importance of accuracy of geometry in this kind of fracture analysis. Their model with obtained material
1.3
Aim of the project
The aim of this project is to generate a detailed FE model of an infant skull and sutures and study the skull fracture propagation from an impact from a fall of a certain height. The model is constructed using available CT images of an infant head of 2 month old age. It is of interest to evaluate the skull fracture propagation from models of skull with sutures using fine and coarse element mesh. The thesis aims to obtain the following:
• geometry of the infant skull, brain and scalp using CT images
• detailed sutures from CT images of unknown Hounsfield unit (HU) values • fine and coarse hexahedral mesh of skull, sutures, brain and scalp
• relevant material properties from literature and biomechanical relations for skull bones and sutures
• FE model of the head and simulation of its impact to a hard surface
• validation of the model ability to show relatively correct fracture propagation
If the model shows promising results of being a good estimation of fracture propagation one would like to construct models from CT images of infant heads from a range of age to have as comparison to future cases of skull traumas. The thesis aims to answer the following questions:
• Is it possible to obtain the sutures geometry from CT images? • Does fracture propagation change with mesh size?
• Does fracture propagation resemble results of earlier constructed model and fractures from Weber’s free fall study?
• Does the method used in constructing the age specific model and simulating the skull fracture demonstrate to be a good way to determine infant skull fracture?
1.4
Anatomy of an infant cranium
The geometry of the cranium varies between individuals in factors like age, gender, genetics and even environmental factors like sleeping position during the first years may influence its shape. The infant cranium is made up of the skull, sutures and fontanelles. The skull is formed from the following seven bones; two frontal bones, two parietal bones, two temporal bones and one occipital bone (Figure 1) in addition to facial bones and lower jaw. The bones of the skull are connected by four main sutures; metopic, coronal, lamboid and sagittal suture (Figure 4). The sutures intersect at two fontanelles which are spaces in between the bones of the infant skull. The posterior fontanelle is the junction of the parietal bones and occipital bone located at the back of the head and the anterior fontanelle is the junction of the frontal bones and parietal bones located at the top of the head [15].
Figure 1: Cranium of 8 month old infant showing frontal, parietal, occipital and temporal bones. Geometry is obtained from projects data and processed in the program Mimics from Materialize [16]
1.4.1 The skull bones
The skull can be divided into two parts: the neurocranium and the viscerocranium. The former one covers the brain and the latter one consists of lower face and jaw [15]. Here after the neurocranium will be referred to as skull. The structure of the infant and adult skull layer differs in composition. The infant skull has single layered compact bones (Figure 2) with varying porosity and fiber pattern radially oriented from the center of ossification [5] [8]. The adult skull layer consists of two flat compact bones and a spongy cancellous bone in the middle called the diploe (Figure 3). Bone thickness is not uniform throughout the skull [17] as can be seen in the plot of Figure 3 (right).
Figure 2: The left image shows a CT image of single layer compact frontal bone of 2 month old infant in sagittal view. The right image shows a Hounsfield unit (HU) distribution of the profile line (green arrow) across the skull bone in the left image. The smooth curve of the plot indicates that formation of spongy bone has not started.
The two frontal bones are located at the superior part of the face just above the eye sockets. Each of the two bones ossifies intramembranously from the end of eight weeks of gestation and is separated by metopic sutures. The two bones are usually completely fused together during an infant first year of life. The frontal bone is one of the first bones to obtain the morphology of an adult bone (Figure 3). Each frontal bone ossifies from a center which is located at the eminence of the bone and can be visualized as a rounded elevation in the middle of it [18] [19]. The eminence of a frontal bone can be seen in Figure 3 approximately where the lower profile line crosses the bone.
The two parietal bones are separated by the sagittal suture and lie in between the lamboid and coronal sutures. They ossify intramembranously from the bone’s center which is located at the parietal eminence. In infants and children the eminence is easily noticeable as an elevation approximately at the bone center. The parietal bone is very thin in infant’s first years and is more easily breakable than the surrounding skull bones [18] [19].
The occipital bone is located in between the parietal bones and lambdoid sutures. The ossification of the occipital bone is more complex than for the previously mentioned bones of the skull. The occipital bone consists initially of 5 separate parts at early fetal stage with only the center area, called the interparietal part, ossifying intramembranously. The other parts develop from preformed cartilage. After five gestational months the interparietal part has fused with its supraoccipital part but fusion of all occipital bones is usually completed at the end of infants first year [19].
Bones of infants are more compliant than adult bones as they have not been fully ossified. Most skull bones are formed from intramembranous ossification. The bones start as membranous layer of connective tissue which takes part in the formation of cells called osteoblasts. The osteoblasts transfer calcium to the bone layer from the surrounding blood supply to build up cartilage. The cartilage matrix is composed of calcium and protein. The proteins are partly built up of collagen which is flexible but strong material. As the cartilage matrix continues to receive calcium the osteoblasts change into osteocytes of the compact bone layer. The original connective tissue, now surrounding the external of the bone layer is called periosteum [20]
Figure 3: Beginning of formation of diploe layer in 8 month old frontal bone. Left image: Two profile lines cross the CT image of the frontal bone (green) in sagittal view. Right image: HU distribution of the profile lines. The dents in the plots indicate formation of spongy bone in-between two compact bones. The thickness of the bone can also be estimated from the plot of each profile line, giving 2.27 mm and 3.17 mm respectively. The thickness evaluation depends on the choice of the lower HU threshold.
1.4.2 Sutures and fontanelles
Figure 4: Cranium of 2 month old infant showing coronal, metopic, lambdoid and sagittal sutures. (The geometry is obtained from the project data and processed in the program Mimics [16]).
Ossification of the skull is incomplete at birth with sutures and fontanelles separating the individual bones (Figure 4). Sutures are membranous tissue built up of collagen fibers forming flexible joints which connect adjacent cranial bones [21] [22]. The morphology of sutures is easily distinguishable from the infant skull. Both interdigitated sutures and straight edged sutures exist [21] , for example between adjacent parietal bones and parietal and occipital bones respectively (Figure 4). The sutures and fontanelles are essential for molding of the head during normal birth [5] as well as following development of the growing and expanding brain [23]. The cranium can, because of the existence of sutures and fontanelles, be impacted from quasi-static or dynamic load without changing its shape significantly [2]. Mazuchowski et al. [2] did a study on the closing of sutures and fontanelles from MRI’s of infants of age from 30 week of gestation to 3 years old. He concluded that posterior fontanelles close between the age of 1 - 3 month old after birth and anterior fontanelles close somewhere in between the age of 12 to 24 month old [15] (Figure 5). Most sutures have grown together and fontanelles have closed by the age of three [2] [15] [24].
Figure 5: Infant age in months at which the fontanels start closing. Image is adapted from [2].
around the bone margin, especially in children. The periosteum surrounds the sutures and the bone edges [22] (Figure 6). Hereafter the endosteal dura mater will be referred to as dura mater.
Chapter 2
Biomechanical & material properties
of infant skull and suture
Infant skull response from an impact differs from the one of an adult skull. Material, geometrical properties and the structural composition of the skull affect the biomechanical behavior [9]. Few studies have been done to determine material properties of the infant skull and sutures and those which have been documented in literature are not always comparable due to difference in mechanical testing methods of the infants bone materials. The overall outcomes from these literatures of infant’s skull biomechanics are listed below:
• Fiber pattern of each bone of the infant skull is radially oriented from the center of ossification.
• Elastic modulus of infant skull bone is higher in the radial fiber direction from the center of ossification than perpendicular to the fiber direction.
• Elastic modulus and ultimate stress of bones increases with gestational age
• The one layer bone structure of the infant skull has different mechanical responses from an impact to a head than a three layered adult skull bone
• Ultimate stress and elastic modulus is larger for the parietal bones than the occipital bone
• Age specific material data should be used in age specific models
• Age and strain rate does not affect the elastic modulus and ultimate stress of infant sutures.
of author. For this reason the infant skull bone will be simplified to be homogeneous and have only elastic properties.
Bones in general are considered to be anisotropic with many of them having a mineral orientation in the direction of maximum load. Each bone has therefore directional properties according to its unique function. The skull bone does not follow the same principles as described above. Its properties are not based on the load direction but much more likely on stiffness and damage tolerance [25]. The infant skull bone has anisotropic material properties with fiber patterns radially oriented from the center of ossification. This means that the bones stiffness varies in response to force applied from different directions as opposed to isotropic material where the stiffness is the same for all directions [25].
The adult sutures are considered anisotropic with complex collagen fiber arrangement [21]. Studies of infant suture properties are very limited and as in the FE simulations done in those studies [9] [14], infant sutures will be assumed to be linearly elastic, homogeneous and incompressible with same density as the dura mater.
2.1.1 Isotropic material properties
Due to limited data in literature, sutures are assumed to be linear elastic. Hooks law describes the linear behavior of a material under load as:
ε
σ =E (1)
with σas stress, E as elastic modulus and εas strain.
For material of a linear elastic solid the stress strain relation is often expressed using the stiffness tensor (elastic modulus tensor), C, or the elastic compliance tensor, S, [26] such that:
kl ijkl ij C ε σ = i,j,k,l=1,2,3 (2) kl ijkl ij S σ ε = i,j,k,l=1,2,3 (3)
For isotropic materials the stiffness tensor has two components, C11 and C12, due to symmetric properties in every direction:
= 0 0 0 11 12 12 12 11 12 12 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C (4)
(
)
2 12 11 0 C C C = − (5)The components of the stiffness matrix C and the compliance matrix S are obtained in terms of the engineering constants in equations (6) and (7) respectively:
( )
(
il jk ik jl)
( )(
)
ij kl ijkl E E C δ δ ν ν δ δ δ δ ν 1 1 2 1 2 − + + + − = , (6)(
il jk ik jl)
ij kl ijkl E E S = +ν δ δ +δ δ +ν δ δ 2 1 (7)Note that C11≡C1111 and
C
12≡
C
1122=
C
2211. The same applies for the elastic compliance tensor S [26].Expansion of equation (2) to matrix form is:
(
)(
)
(
)
(
)
(
)
− − − − − − − + = 12 13 23 33 22 11 12 13 23 33 22 11 2 2 2 2 2 1 0 0 0 0 0 0 2 2 1 0 0 0 0 0 0 2 2 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 2 1 1 ε ε ε ε ε ε ν ν ν ν ν ν ν ν ν ν ν ν ν ν σ σ σ σ σ σ E (8)with E as the elastic modulus, ν the Poisson’s ratio, σ11=σ22=σ33 the stresses in x,y,z direction respectively, σ23 =σ13 =σ12 the shear stresses, ε11 =ε22 =ε33the strains in x,y,z direction respectively and ε23 =ε13 =ε12the shear strains.
The above matrix can be expressed in a more convenient form using index notation:
− + + =
∑
= 3 11 2 1 k ij kk ij ij E ε δ ν ν ε ν σ ∀i,j∈{1,2,3}, = = else j i if ij 0 1 δ (9)with the symmetry conditions:
ji ij ji ij σ ε ε
2.1.2 Transversely isotropic material properties
The material structure of a tree is often used as an explanation of an anisotropic structure that has symmetrical properties. The tree structure can be simplified to having orthotropic material properties with orthogonal planes of symmetry. The stiffness matrix has then nine components as it has three perpendicular symmetry planes with the principal axis directions going axial, radial, and tangential (circumferential) (Figure 7). Bone can be modeled as being orthotropic but is often simplified even further by assuming that it has a symmetry plane perpendicular to the fiber direction making it transversely isotropic. This reduces the parameters needed from nine, for orthotropic material, to five. For example, in long bones the fibers running perpendicular to the long axis can be assumed isotropic in its planes [25].
For infant skull bone it has been established that the elastic modulus parallel and perpendicular to the fiber direction is not the same. It is also known that the fiber direction is radially directed from the center of ossification. This means that the plane of the axial axis, representing the bone thickness, and the tangential axis is the plane of isotropy.
For transversely isotropic material of a linear elastic solid and isotropy in the 2nd-3rd plane, the stress strain relation can be represented using the stiffness matrix C:
= 55 55 44 22 23 12 23 22 12 12 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C (11) 13 12 C C = , C55=C66and
(
)
2 23 22 44 C C C = − (12)The compliance matrix, S, is a more convenient representation of the engineering constants E (elastic modulus), ν(Poisson’s ratio) and G(shear modulus) [26]:
Figure 7: Infant skull bone representation. Principal axis directions showing radial direction parallel to the fiber direction in infant skull bones, tangential direction perpendicular to the fiber direction and axial direction in the bone thickness.
− − − − − − = 12 12 23 3 2 23 1 31 3 23 2 1 21 3 13 2 12 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 G G G E E E E E E E E E S ν ν ν ν ν ν (13)
The independent material constants that need to be obtained are E1,E2,G12,ν12and ν23. The plane symmetry givesE2=E3,ν23=ν32 and ν21=ν31. The shear modulusG23is calculated in
equation (14):
(
23)
2 23 1 2 +ν = E G (14)The Poisson’s ratios are not symmetric but satisfy equation (15):
2 21 1 12 E E ν ν = (15)
The components of the stiffness matrix (equation (11)) can be derived from the engineering constants of the compliance matrix:
(
−)
Υ = 2 23 1 11 E 1 ν C (16)(
−)
Υ = = 33 2 12 21 22 C E 1 ν ν C (17)(
−)
Υ=(
−)
Υ = = 13 1 21 23 21 2 12 23 12 12 C E ν ν ν E ν ν ν C (18)(
−)
Υ = 2 23 12 21 23 E ν ν ν C (19)(
)
2 23 44 1 2 E v C = + (20) 12 55 1 G C = (21)(
12 21)
232 2(
12 21 23)
2 1 1 ν ν ν ν ν ν − − − = Υ (22)2.2
Parameter selection
In a patient specific model, as the one developed in this thesis, it is important to have accurate parameters in order to obtain reliable results. Obtaining correct parameters can be a hard task and much estimation needs to be done when there is limited data available either from material testing or from literature. Parameters needed for the FE element simulation of this thesis, to study the skulls biomechanical effects during an impact from a fall are obtained from literature and estimations. The biomechanical properties, elastic modulus, Poisson’s ratio and shear modulus are needed both for the skull and sutures. The additional parameters needed for the constitutive model used to predict failure (chapter 4.2.2) of the skull are longitudinal and transverse tensile strength, shear strength, transverse compressive strength, and normal and transverse shear strength.
The elastic modulus is a measure of the elastic deformation of the material when load is applied to it. It can be obtained from the slope of a stress-strain curve for each material (Figure 8). The shear modulus measures the resistance of a material to shear deformation. Higher value means that the material is more rigid to shear. The shear modulus is defined as the ratio of shear stress versus shear strain. The Poisson’s ratio is the ratio of lateral versus longitudinal strain in uniaxial tension. The value is dimensionless and ranges from -1 to 0.5 for stable materials and is equal to 0.5 for incompressible one. Tensile strength is the measure of materials strength in the form of how much tensile load it can withstand without failing. The tensile strength is also called ultimate strength and is defined as the maximum tensile force per unit cross sectional area of a specimen. A material deforms elastically under a slowly increasing tensile load up to a certain point and will return to its original size and shape if the load is removed. The material deforms permanently if the load is not removed at this point, which is called the yielding point. After the yielding point the material behavior depends on the material in question but finally it reaches its breaking point. The maximum of the stress-strain curve is known as the ultimate strength. Shear strength is the measure of the maximum shear stress which a material can withstand without failing. The shear stress acts parallel to the force direction and causes some portion of the material to shear with respect to another portion. Compressive strength is the measure of the maximum compressive load a material can take per cross sectional area without failing. For bones the compressive strength is usually greater than the tensile strength (Figure 9) [27] [28].
Figure 8: Stress-strain profile of bone samples taken at different orientations to the long axis of compact femoral shaft. Image is adapted from [29]
Figure 9: Failure strengths of anisotropic bone with sample orientation perpendicular to the long axis. Image is adapted from [29].
2.2.1 Parameters obtained from literature
Kriewall presented his data of elastic modulus, both parallel and perpendicular to the fiber direction, of parietal and frontal bone specimens obtained from a three point bending test [5] [8]. His specimens were taken from PMHS skulls of age from 24 to 40 weeks of gestation. He divided the results into two categories, preterm and term specimens. Data from 22 term specimens (age range 36 to 40 weeks of gestation) taken from parietal bones of 4 infant PMHS were presented. He found that the ratio of elastic modulus for specimens with the long axis parallel and perpendicular to the fiber direction was 4.2:1. Only data of 2 perpendicular specimens and 12 parallel specimens were presented for the frontal bones. The ratio parallel to perpendicular was found to be 1.8:1 [5].
Coats and Margulies tested the elastic modulus and ultimate stress for specimens with fiber direction perpendicular to the long axis and presented their data of parietal and occipital bone specimens from three point bending tests [12]. The specimens were taken from infant PMHS skulls of age from 21 weeks of gestation to 13 month old. They found significant influence of infant age on the elastic modulus. The data from Coats & Margulies was plotted and used as a basis to obtain age specific material parameters for the 2 month old infant skull model constructed in this thesis (Figure 10 to Figure 12). The ratio obtained by Kriewall was used to find the parallel elastic modulus for the skull bones (Table 1).
Figure 10: Elastic modulus of parietal bone vs. age of specimens with fiber direction perpendicular to the long axis. Data is obtained from [12]. Elastic modulus for parallel to the fiber direction is scaled by using ratio difference stated in [5].
Figure 11: Elastic modulus of occipital bone vs. age of specimens with fiber direction perpendicular to the long axis. Data is obtained from [12]. Elastic modulus for parallel fibers is scaled by using the ratio found in [5].
The plots presented in Figure 10 and Figure 11 show weak relationship between age and elastic modulus with few outliers which are highly influential on the linear regression line. The overall patterns indicate that the elastic modulus increases more before few weeks after birth than in the following months. None the less, due to limited data available, linear relationship is assumed. As Coats and Margulies did not obtain parameters for the frontal bone and Kriewall did not obtain parameters for the occipital bone, the difference of their values of termed specimens (36 - 40 weeks gestation) from parietal bone was measured. It was found that their values differed by a factor of 2 for specimens taken perpendicular to the fiber direction. This value was used to scale down the elastic modulus of the frontal bone found by Kriewall for its value to be in compliance with the values obtained from Coats and Margulies. Assuming linear relation between age and elastic modulus, shown in Figure 10 and Figure 11, the elastic modulus
y = 8.48x + 433.3 0 200 400 600 800 1000 1200 1400 -05 00 05 10 15 E la st ic m o d u lu s (M P a ) Age in months Parietal bone E_perpendicular y = 18.20x + 287.6 0 200 400 600 800 1000 1200 1400 -05 00 05 10 15 E la st ic m o d u lu s (M P a ) Age in months Occipital bone E_perpendicular
obtained for parietal, occipital and frontal bone, for the 2 month old infant, are presented in Table 1.
Table 1: Elastic modulus for parietal, occipital and frontal bone both for parallel and perpendicular fiber direction. Data is both obtained from [5] and [12] and estimated according to obtained ratios explained in the text. E_perpendicular (MPa) E_parallel (MPa) Parietal bone 450 1887 Occipital bone 324 1358 Frontal bone 1092 1929
Ultimate stress for parietal and occipital bone was obtained and plotted from data presented in [12]. From the graphs in Figure 12 linear relation may be assumed between elastic modulus and ultimate stress.
Figure 12: Ultimate stress vs. elastic modulus for parietal bone (left) and occipital bone (right). Data obtained from [12].
The ultimate stress values for the 2 month old infant skull bones are listed in Table 2. The ultimate stress for the frontal bone is estimated from the linear relationship presented for the parietal bone in Figure 12 (left).
Table 2: Ultimate stress for parietal, occipital and frontal bone both for parallel and perpendicular fiber direction. Data is both obtained from [12] and estimated according to obtained ratios explained in text
Ultimate stress_perp-Ultimate stress_parallel y = 0.06x + 3.82 0 10 20 30 40 50 60 70 80 90 0 500 1000 1500 U lt im a te s tr e ss M P a
Elastic modulus (MPa) Parietal bone E_perpendicular y = 0.029x + 5.27 0 10 20 30 40 50 0 500 1000 1500 U lt im a te s tr e ss M P a
Elastic Modulus (MPa) Occipital bone
The Poisson’s ratio is not known for infant skull bones and is therefore estimated to be the same as for an adult skull bone in radial compression with ν12=0.19 and in transverse compression
22 . 0
23=
ν [30].
The elastic modulus of the sutures is obtained from [12] and [14] as an average value for infants less the one years old, Esuture =8.1 MPa. The material is assumed to be nearly incompressible with Poisson’s ratio νsuture=0.49 and density equal to that of dura mater,
3 1130 m kg suture = ρ [14] [31]. 2.2.2 Parameters estimations
Some parameters of the skull bones cannot be obtained from literature and need to be estimated. The parameters E1,E2,ν12andν23are obtained from literature. Assumption is made thatE2=E3,
32
23 ν
ν = and ν21=ν31 due to the transverse isotropic symmetry plane. The Poisson’s ratio ν21
is found using equation (15) and the shear modulus G23 is found using equation (14). Huber
[32] proposed that the shear modulus for an in-plane orthotropic material could be predicted by equation (23).
(
12 21)
2 1 12 1 2 + ν ν = E E G (23)Craig and Summerscale [33] predicted from equation (23) that the bulk modulus could be estimated for a square symmetric material using equation (24).
(
3)
23 31 12 3 3 2 1 2 1 3 + ν ν ν = EE E K (24)For transversely isotropic material the elastic modulus is equal within the plane of symmetry and the Poisson’s ratio ν21=ν31, thus equation (24) becomes:
(
3)
23 21 12 3 2 2 1 2 1 3 + ν ν ν = EE K (25)The values for the shear strength are not available and therefore need to be estimated according to best knowledge. One way is to scale the shear strengths according to the ratio of the shear modulus for a certain direction to the ratio of the elastic modulus and strength for one of the uniaxial directions. The shear strength is assumed to be symmetric. This would result in the estimation equation (26). i ult ij ult ij ij E G G S ≈ ε ≈ σ (26)
The transverse compressive strength is estimated by scaling the transverse tensile strength according to ratio found between compact femoral shaft of transverse compressive strength and transverse tensile strength [29] (Figure 9). The bone should be stronger in compression which is
which is used to predict failure of the skull (Chapter 4.2.2), in the FE program LS-Dyna, are listed in Table 3.
Table 3: Parameters for constitutive model of skull bone fracture analysis for 2 month old infant
Abbrevi-ation Parietal bone: 2 month old Occipital bone: 2 month old Frontal bone: 2 month old Unit
References, formulas and estimations
1
E 1887 1358 1929 MPa
Elastic modulus parallel estimated from linear relation of parallel and perpendicular data from Kriewall and perpendicular data from Coats & Margulies 2006, [12], see Figure 10
2
E 450 324 1092 MPa Elastic modulus perpendicular. Coats &
Margulies 2006, [12], see Figure 10
12
ν
0.19 0.19 0.19 − Poisson’s ratio. McElhaney et al, [30], from adulthuman cranial bone
21
ν
0.045 0.045 0.11 − Equation (15)23
ν
0.22 0.22 0.22 − McElhaney et al, [30], from adult human cranialbone
23
G 184.4 132.8 447.5 MPa Shear modulus transverse to plane of symmetry,
equation (14)
12
G 421.6 303.5 634.9 MPa Shear modulus in plane of symmetry is estimated
using equation (23) , [34] [34].
ρ 2090 2090 2090 kg/ m3 Density, Kriewall, [8], from infant of 42 week of
gestational age
1
S 111.9 44.5 114.4 MPa Longitudinal tensile strength, Coats & Margulies
2006, [12]. See Table 2
2
S 29.0 14.6 66.4 MPa Transverse tensile strength, Coats & Margulies
2006, [12]. See Table 2
12
S 25.0 9.9 37.7 MPa Shear strength. Estimated using equation (26)
2
C 38.7 19.5 88.5 MPa
Transverse compressive strength, Scaling from “Frankel VH, Nordin M: Basic Biomechanics of the Skeletal System. Philadelphia, Lea & Febiger, 1980”
N
S 29.0 14.6 66.4 MPa Normal shear strength, assumed to be equal to
Chapter 3
Creating Finite Element mesh from
CT data
The infant head model, consisting of skull, suture, brain and scalp, is constructed from CT images of a 2 month old infant. For a realistic simulation to be carried out the scalp needs to be included as it may increase the contact area of head to plate [14]. The emphasis of the image processing is on obtaining relatively accurate and good representation of skull and sutures and approximate geometry of the scalp and brain. The chapter covers basic introduction to CT imaging and then goes through each step of the image segmentation process up to the point where a working model has been obtained. After a good representation of the head has been established, the model is meshed for it to be used in the FE simulation.
3.1
Introduction to computed tomography imaging
Computed tomography (CT) is an imaging technique which uses X-rays in the production of cross-sectional images of an object. The röntgen ray was discovered by Wilhelm Konrad Röntgen in 1895. A method to reconstruct the projection was developed by Johann Radon in 1917 and in 1972 the first CT scanner was developed by Godfrey N. Hounsfield and Dr. Allan M. Cormack [35]. The initial intensity of photons, I0,
generated by the thin X-ray beam, builds up each cross sectional image by scanning the total field of view (FOV) by a set of lines from a number of angels and distances from the center. This is done by using an x-ray tube which emits radiation transverse through the patient while rotating helically around him. The detector on the other side of the tube measures the residual radiation, It. The linear attenuation coefficient µ is a function of both the photon
energy used during the scanning and the specific material in each pixel. For a homogeneous material and an x-ray beam consisting of single energy photons the residual energy is given by
Figure 13: CT scanner rotating helically around the object
x t I e
I = 0 −µ∆ (27)
where ∆x is the length of the X-ray path through the material. The attenuation coefficient, µ, decreases, usually, with increased photon energy. The attenuation coefficient is used to calculate the Hounsfield unit (HU) of each pixel in the 2D slice of the FOV.
1000 2 2 ⋅ − = O H O H HU µ µ µ (28) where HO 2
µ is the attenuation coefficient of water. The HUs, representing the CT image, are in the range of −1000 to 3095 HU, where -1000 HU is air and 0 HU is water. The HUs are a gray scale values with a range too wide for the eye to distinguish between its values in one window. Thus when the image is viewed in an imaging program, only a certain window level is selected by the observer depending on the HUs of the object of interest. All HUs below and above this range are set to black and white respectively (Figure 14 and Figure 15). Table 4 lists some relevant HU values for this project.
Figure 14: Window level setting adjusted for visualization of bones
Figure 15: Window level setting adjusted for visualization of brain
3.2
Image segmentation
The CT image data set used in this project is of a 2 month old infant head with image resolution of 512x512x327 pixels. The size of the pixels are 0.4063 mm in x and y direction and slice thickness 0.5 mm. The images were taken by Siemens Definition scanner with attenuation strength 120 KeVand tube current 165 mAs. The image data is compiled in a DICOM format which stands for Digital Imaging and Communications in Medicine, and is a format used
programmed modules for image processing but also has the option of implementation of new methods using Python or C++ language [37]. 3D slicer needs to run on a computer that has large RAM (4 GB or more) and a fast graphic card [36]. The same applies when working with large datasets in DeVIDE.
It is important to find a proper segmentation method for the task at hand as it can be difficult to obtain an optimally segmented region with spatially accurate boundaries but still smooth enough to represent the correct morphology as noise and nearby objects may interfere. A range of segmentation methods exist such as pixel, edge and region based methods. Pixel-based segmentation is based on the gray value or HU scale which can be selectively chosen within some threshold interval. An edge-based method uses edges of an object detected by local discontinuities to separate the image into regions. A region-based method constructs homogeneous regions propagating from the origin of seed points placed by the user at regions of interest (ROI) [38].
3.2.1 3D Slicer
The dataset in DICOM format is imported into the 3D slicer where the sagittal, coronal and axial planes can be viewed (Figure 16). A pixel based method was chosen to obtain the model of the skull, sutures, brain and scalp using a threshold level to select the HU range which represents each part. The HU values used are based on values from literature (Table 4) as well as visual evaluation.
Figure 16: Window view in the 3D slicer environment. Axial, sagittal and coronal planes are shown in the three vertical images on the right.
Table 4: Table with HU of materials needed for segmentation of the infant head
Material HU Reference
calcium 100–300 [39]
bone 600 - 3071 [39]
Brain (white and gray matter)
15–35 (125keV), 30–40
[40] [39]
Sagittal sinus 50–68 (120keV) [41] [42]
Scalp (fat and soft tissue)
-200 – -30 -30 – 10
[43]
Collagen 250 [44]
The infant skull bone is compliant and contains both calcium and collagen. From literature, calcium has been found to be in the range of 100 to 300 HU, collagen around 250 HU and the lower threshold of an adult skull bone around 600 HU. Guided by the literature and visual evaluation, the infant skull bone was found to be within a threshold with a lower limit of 130 HU and an upper limit of 3071 HU, which is the maximum HU for this dataset (Figure 17). The lower jaw of the viscerocranium was separated from the skull by manual editing as it is not needed for the skull model (Figure 18).
Fat and soft tissue of a healthy adult head has been found in literature to be 100 to 30 HU and -30 to 10 HU respectively [43]. By visual evaluation, the range for an infant scalp, consisting of fat and soft tissue, was found to be in the range -100 to 15 HU (Figure 19).
Figure 17: Sagittal view with thresholded skull area indicated by color.
Figure 18: Coronal view of the jaw, in pink, was separated from the skull, yellow.
Figure 19: Fatty and soft tissue of an infant scalp is shown in brown.
(Figure 23). It was not possible to use thresholding or region growing methods for the sinus area as it is not only the sinus that is located in the gap but also partly dura mater and falx cerebri (Figure 6). These areas were therefore manually edited by going through every slice and selecting the regions missing for the completion of a solid brain part within the skull (Figure 24).
Figure 20: White matter of 15-35 HU
Figure 21: Gray matter of 35-40 HU Figure 22: Continuous area of 15-63 HU representing brain for the FE model
Figure 23: Axial view of dorsal (top) part of skull and brain showing also the superior sagittal sinus crossing the two hemispheres.
Figure 24: Axial view of solid brain representation after manual editing has been applied to include the sagittal sinus.
Figure 25: Brain before and after addition of sinus (left and right respectively).
The sutures were the most difficult part to segment. No previously published HU values were found for this membranous tissue built up of collagen fibers. The connection of the sutures to the skull is also in continuous contact with the dura mater and periosteum (Figure 6) making it more difficult to obtain only the geometry for the sutures when segmented. By knowing the location of the sutures and fontanelles, which were easily visible after the segmentation of the skull (Figure 4), it was assumed that the HU range was somewhere in between 63 and 129 HU (Figure 26). This range of HU also captures some of the dura mater and periosteum which surrounds the inner (concave) and outer (convex) sides of the skull respectively as well as cartilage of the ears (Figure 27). Narrowing this HU interval resulted in a large discontinuity of the supposed suture material. Manual editing was applied to disconnect some well-connected and unwanted parts to the sutures. Further processing was carried out in DeVIDE.
Mask images of each part were saved in compressed Meta format as .mhd file and .zraw file. The .mhd file is a header containing information such as dimension and element size of the image data. The .zraw file contains compressed binaries of the image.
Figure 26: Membranous tissue in the range of 63-129 HU. Left gap is part of metopic suture and right gap part of coronal suture.
Figure 27: Top view of skull and membranous tissue.
Figure 28: Sagittal view of skull, membranous tissue and scalp.
Figure 29: Detailed view from Figure 24 showing where membranous tissue covers the anterior fontanelle.
3.2.2 DeVIDE
The network in Figure 30 takes as input the geometry of the sutures and the network in Figure 31 takes as input the geometries of the brain, skull and scalp one at a time. All the image data is in compressed Meta format. A double threshold is used to isolate the masks of interest and set the values within the threshold range equal to 1 and the values outside the threshold equal to 0. The network in Figure 30 was constructed to continue the process of isolating the sutures from the dura mater and other tissues with the same HU values. The module fast marching is used to obtain the morphology of the sutures by selecting only pixels with continuous connection to the seed points place at ROI (Figure 33). Both networks use closing module with convolution kernel first with dilation followed by erosion. This closing filter is used to connect voxels of small gaps and discontinuities in the geometries. Figure 33 shows the geometry of the sutures obtained from the 3D slicer, Figure 34 shows the sutures after seed connectivity (fast marching) and
Figure 30: Network for isolation of sutures
Figure 31: Closing and filling of holes for skull, brain and scalp
Figure 32: Image mathematics performed to prevent penetration of objects
Figure 33: Sutures from 3D slicer Figure 34: Sutures after fast
marching
Figure 35: Sutures after closing
Figure 36: Original sutures are shown in green and added voxels in blue.
Figure 37: Brain (without sinus) containing holes in volume.
Figure 38: Brain (without sinus) after
imagefillholes.
To prevent intersection of two objects a mathematical operation is performed for the sutures, skull and brain. The network in Figure 32 starts by subtracting the sutures from the brain using the module imagemathematics. Both objects have set of mask values equal to 1 such that voxels of intersections are equal to 0, the relative complement of sutures in brain is equal to 1 and the relative complement of brain in sutures is equal to -1. All values equal to -1 are set to 0 with the second imagemathematics module and then the process is repeated for the brain and the skull.
the brain and the skull with reduced volume equivalent to the expansion of voxels of the sutures into areas of the other objects. The results of sutures, skull, brain and scalp are all written to separate uncompressed Meta format.
Figure 39: Final volume of the brain Figure 40: Final volume of the skull Figure 41: Final volumes of the skull, brain and sutures together
3.2.3 Matlab
The Meta images are imported into Matlab with the function metaImageRead.m from the Matlab program Slicer©. The function takes in datasets in Meta format and exports it as 3D arrays of voxel classification. Sutures and skull are added together to form one continuous geometry consisting of two parts, such that one part has voxels of value equal to 1 representing its volume and the second part with voxels of value equal to 2. The three datasets, skull with sutures, brain and scalp, now as 3D matrices, are saved in separate MAT files with the addition of an information file containing the size of the matrix and a second file containing voxel size. The three files within each MAT file are the inputs for the Matlab program Dicer (see Chapter 3.3.1) which is used for the meshing of the volumes.
For statistical purpose the segmented volumes are mapped to the original image to evaluate the number of voxels lying outside of its original HU region. The volumes for each object may have increased after the image processing of manual editing, closing and imagefillholes. A Matlab function, mask_vs_HU.m, was written which maps the mask to the original image to obtain the number of added voxels to the object.
Table 5: Voxel increase in geometries after image processing
Volume Number of voxels in Number of voxels in Percentage increase
From Table 5 it can be seen that volume increase of the brain and skull is very small and is probably mostly due to small holes and discontinuities which have been filled and fixed. The increase of the volume of the sutures is larger than expected. This could be a result from manual editing and inclusion of extra voxels during closing and imagefillholes filters. The sutures had discontinuities and contained many holes in the original image data due to its thin structure. The discontinuities and holes are mostly due to the resolution of the image data. The resolution may cause the average HU values on edges of two aligned anatomical structures to become skewed. Some small area belonging to the thin structure of sutures may therefore have been assigned to voxels of surrounding structures during the acquisition of the image data. The large increase of number of voxels after processing of the sutures data is party necessary to obtain continuous geometry. Manual editing might also have to be done with more care.
3.3
Hexahedral meshing
Many programs exist for meshing, but a wide range of them were developed for meshing of CAD (computer aided design) objects and are not suitable for meshing of complex geometries obtained from medical images (MRI, CT, etc.). A better image based approach is to create the mesh from the voxels of the image data. The mesh is then composed of brick elements with the geometric accuracy depending on the resolution of the image data [45].
3.3.1 Dicer
It is a great challenge to obtain a quality mesh for a complex geometry like the brain, skull, sutures and scalp for an application of FE method. To mesh the geometries, voxel hexahedral meshing was used. The segmented CT images are imported to the program “Dicer” [45] having been processed by 3D Slicer, DeVIDE and metaImageRead.m, to get the image data to a 3D array of voxel classification. The Dicer algorithm is based on voxel hexahedral meshing where every voxel of the image data becomes a hexahedral element.
The program, Dicer was chosen for the meshing of the head as it provides automatic hexahedral meshing for complex geometry in relatively short time. It has the option of merging planes to get larger elements, instead of resampling the image data, and it has error correction and smoothing algorithms [45]. The head geometry contains large number of voxels, it is constructed of more than one part with each part highly detailed and it has very thin morphology in some areas.
CT images have stair-stepped surfaces and the step size depending on the image resolution. The surface quality can be improved by applying smoothing. The term ‘smoothing’ refers to the moving of nodes of the mesh, maintaining a fixed connectivity, in order to improve the mesh quality [46]. When smoothing is applied to areas of only a few pixels in thickness some elements can become distorted [45]. The Dicer addresses this problem by having the option of merging planes together, such that two or more elements become one, and such that the
smoothened. Dicer has the option of choosing the smoothing coefficients for all boundary nodes and corner nodes of a value between 0 and 1, 0 resulting in no smoothing and 1 resulting in collapse of all nodes to a single point in space. For value in between 0 and 1, the new location of a boundary node is found from weighting the average of its original location and the centers of the final position of neighboring boundary nodes. This smoothing coefficient has to be selected with care to prevent element distortion [45]. To further prevent distortion of elements all internal nodes (nodes not lying on boundaries) are smoothened after all boundary nodes have been fixed in space [45].
Merging of planes is supposed to preserve the morphological details of the model better then resampling the image data to obtain lower resolution and less number of elements in a hexahedral mesh [45]. Dicer has the option of choosing how many planes should merge for each direction.
An estimate of the number of elements in the final mesh after merging was calculated from equation (29).
∏
= ≈ z y x i i original final n mesh mesh , , (29)With meshoriginal as the number of voxels in the original mesh and nias number of merging planes in directioni. This is only an estimation as the complexity of the geometry does not allow merging of some element planes and thus the number of final elements was found to be somewhat larger. For example for a cuboid with even number of elements in each direction, merging 2 planes together in all directions reduces the number of elements in the final mesh to 12.5 % of the original mesh.
The Dicer program is a non-commercial program created in Matlab. The code was provided to the department of STH from the author of the program. Some problems in the code were fixed in order to make the program work. The Dicer has the option of writing the output file in Abaqus and Patran format, but after some checking their codes were not functioning correctly. The Abaqus output writer was fixed as the Abaqus format can be opened by many programs, for example with Hypermesh which was used to view the results of meshing, while working on fixing the Dicer. A function was then written, lsdynawrite.m, as an addition to the Dicer which writes the output of the mesh to LS-Dyna .k format and can therefore be directly imported to LS-Dyna Prepost (see Chapter 4.).
3.3.2 Mapping
voxels as the skull and for this reason the mesh size of the brain is increased to reduce the number of elements in the model. Both fine and coarse mesh are wanted for the skull and sutures. The skull and sutures are therefore merged as one solid to obtain the coarse mesh and then afterwards, sutures are mapped to the solid.
Figure 42: The format of output file from centercoordinate.m
∑
= = = = 8 1 1n i ji j c n c c , j=x,y,z (30)The process of meshing and mapping is described by the flowchart in Figure 44. By meshing an object in Dicer with an element size corresponding to the image resolution, one obtains a detailed mesh with each part identity (pid) correctly assigned to the right element. The same object, but now with all parts as one, is meshed again with a chosen number of merged planes in each direction and smoothing values for boundary edges and corners. Next step is to map the correct pid from the fine mesh to the coarse mesh. A function has been constructed,
centercoordinate.m, as an extension to Dicer which finds the center node (equation (30)) of each
element in the coarse mesh, specified with an element identity (eid), and writes it to a .m file (Figure 42). Another function has been constructed, nodecoordinate.m, which finds the coordinate of each of the eight nodes belonging to each eid of the fine mesh and writes it to a .m file for each pid separately (Figure 43). The mapping is carried out by checking if a point from the coarse mesh is within the edges of an element from the fine mesh. This is done for all center nodes and elements in the coarse and fine mesh respectively.
A program, mapping.m, was constructed having three inputs: a fine mesh file for a specific part of the whole object, the whole coarse mesh file and a part identity (pid) to be assigned to the coarse mesh if a point is within an element. The Matlab function inhull.m is called within the
mapping.m which carries out the search of a point in an element. The reason for only taking in
one part of the object at a time is to make each run in a manageable time and observe the result. When mapping is completed each element of a column containing part number for the coarse mesh file should have a pid. In the case of more than two parts (p=2) one can map p-1 materials and assign the elements with no pid to pid = p. The .k file containing elements and nodes of the coarse mesh can now be written. The pid column is saved and loaded into the function
lsdynawrite.m and the meshing is run again with same parameters for plane merging and
smoothing but now the correct pid is given. It is important to comment out the drawing feature in the function dicer.m as the drawing takes extremely long time and is not in any case usable as the .k file can be opened and viewed very easily in LS-Dyna prepost or Hypermesh.
