Quantitatively Structuring A Trabecular Network Using Computational Geometry Techniques

Osteoporosis or ‘porous bones’ is the major cause of bone fracture in elderly post-menopausal women. As such, predicting fracture occurrence can significantly reduce the associated morbidity and reduce health-care expenditures. Examining local micro-structure of trabecular bone is becoming increasingly important for investigating bone fracture behavior. Conventional voxel based measures of bone architecture cannot provide morphology information regarding the fractures occurring at individual trabeculae. In this paper, the Delaunay triangulation, a computational geometry technique is used for assessing the effect of trabeculae structure on local micro-fracture, from sequential 2D images of the cross section of a hamster femur surface undergone a multistage loading until fracture occurrence. The most significant structural parameters to characterizethe strength of cancellous bone and the change in morphology of individual trabeculae during the fracture process are investigated. The results from different regions were compared to describe the possible architectural factors affecting the local micro-damages.


Introduction
Osteoporosis is a major public health problem which will worsen as the global population ages.Accurate diagnosis of osteoporosis can reduce the occurrence of fractures and the associated mortality rates and healthcare costs.The current standard in bone health assessment is to measure Bone Mineral Density (BMD) using dual x-ray absorptiometry (DXA) [1].However, BMD is not a comprehensive indicator of bone strength, which is also highly dependent on the structure of trabecular bone [2].In basic bone research, structural parameters are key measures in characterising the structure of cancellous bone [3].The effect of morphological heterogeneity has been studied using these parameters.Several studies have been performed to identify the contribution of these indices to bone strength and on how they are altered with age and by pharmaceutical intervention [4].This was based on the intuitive idea that the mechanical properties of trabecular bone depend strongly on its apparent density or bone volume [5,6].The most widely used structural indices are trabecular width (Tb.Wi), trabecular separation (Tb.Sp), trabecular number (Tb.N) and mean lntercept length.The calculation of these indices is mainly voxel based, only the means of physical features of a bone structure can be calculated [7,8].Meanwhile, Qin et al [9][10][11][12] presented a multifield theory for analysing bone remodelling and porosity which can be used for analyzing bone structure and density behaviours due to both mechanical and electric fields changes.According to its definition and theory, bone strength depends not only on the quantity of bone tissue and density but also on the quality, which is characterized by the bone architecture, mineralization, turnover, and micro-fracture accumulation.Besides, micromechanics approaches [13][14][15][16] were also used for evaluating effective bone or dental properties [17,18] and then evaluating bone strength.
As we know, in a trabecular structure the longitudinally oriented trabeculae are longer and thinner than trabeculae in transverse direction, and width of individual trabeculae is not uniformly distributed.In individuals with osteoporosis the longitudinal trabeculae are typically maintained, and in fact, become thicker to withstand the loads typically applied to them; while the transverse trabeculae are generally lost to a greater extent compared with controls without osteoporosis [19,20].As the cross-struts between longitudinally oriented trabeculae become discontinuous, the remaining trabeculae become relatively longer.Common voxel based measures of bone architecture do not provide morphology information regarding the fracture occurring at individual trabeculae, thus cannot be used for accurate fracture prediction and monitoring response of bone architecture to drug treatments [21][22][23].
In this paper, a new geometry computational approach is developed for measuring structural histomorphometric parameters of individual trabeculae.The most significant structural parameters to characterizethe strength of cancellous bone and the change in morphology of individual trabeculae during the fracture process are investigated.

Trabecular Bone
Trabecular bone is a form of skeletal tissue found at the ends of long bones, near joints and inside vertebrae.Though its composition of hydroxyapatite, collagen and water is similar to dense cortical bone, there is significant difference in structure.Trabecular bone is highly porous, made up of elements known as trabeculae.The pores of the trabecular structure are filled with bone marrow.The complex structure of a trabecular bone gives it unique mechanical properties [24].In order to characterize properties of individual trabeculae, the topological structure that exists in trabecular bone is presented here.In the skeleton of 2D micro-architectural images, individual trabeculae are joined together at nodes, and stop at free ends and the image border.These trabeculae are usually classified into 7 different groups (see Table 1 and Figure 1 for details).

Bone Segmentation
Delaunay triangulation geometry has been shown efficiency in dividing the image into simpler segments [25].The proposed method makes use of this geometrical computation technique to decompose the image into sets of triangles for representing different trabecular structures.To construct the image partition, boundary points of individual trabeculae are taken as the input of Delaunay triangulation construction.These points are found by a contour extraction method [26] (Figure 2).The result is a set of segments which can be divided into two distinct groups.Using the gray-scale value of the centre of each triangle, it is identified as a bone segment or a marrow segment.Bone segments are divided into four different groups (See Table 2 and Figure 3) where neighbouring triangles are defined as triangles sharing a common edge [27]: • Isolated triangle: a bone segment with no other bone segment(s) in its neighbourhood.

•
End triangle: a bone segment with only one other bone segment in its neighbourhood.

•
Normal triangle: a bone segment with two other bone segments in its neighbourhood.

•
Junction triangle: a bone segment with three other bone segments in its neighbourhood.

Figure 3. Different segment types
A node and a free end in a trabecular structure correspond to a J-T and an E-T respectively.N-Ts form the main body of a trabecula.Individual trabeculae are extracted by tracing adjacent N-Ts.In the partitioned image, Quantitatively Structuring A Trabecular Network Using Computational Geometry Techniques each E-T is assigned to only one trabecula.Therefore, E-T can be used as a start point in tracing the N-Ts.The tracing process is terminated by approaching a junction (for Nd./Fe.trabeculae) or unprocessed end triangle (for Fe./Fe.Tr trabeculae abeculae).Applying this method to all identified E-Ts in trabecular structure and tracing N-Ts to find the associated branch, results in extraction the entire structure (Figure 4).

The skeleton of a Trabecular Structure
The skeleton of a trabecular structure is represented by characteristic local symmetry points.Two different characteristic points exist in bone segments: • The characteristic skeleton point of an N-T is the mid-point of a straight-line segment connecting the mid-points of two internal edges of the triangle [28].
An edge of an internal triangle is an external edge if it is a contour segment.Otherwise, it is an internal edge.

•
The characteristic skeleton point of an I-T, E-T, or J-T is the centroid of the triangle.Since the skeleton obtained from this process contains many unwanted short line segments, it is called the primary skeleton of the shape (Figure 7).

Artifact Removal
Initial skeletonisation of a shape generated via Delaunay Universal Journal of Engineering Science 2(1): 30-37, 2014 33 triangulation may not always conform to human perceptions of the shape.For instance, a J-T may be found in a place which should not be an intersection of the shape (Figure 8 (A)).Consequently, the primary skeleton generated using Delaunay triangulation contains many unwanted skeletal elements which may induce error in measured properties and increase the time and computation required to evaluate the entire trabecular structure.Skeletonisation artifacts include two different types of artifacts: • • Periphery artifacts correspond to shape protrusions which are perceptually insignificant.

•
• Intersection artifacts refer to extra skeleton branches between split skeleton junction points after skeletonisation.To reduce periphery artifacts, the initial skeletonisation of the shape needs modification.Two criteria were used in artifact identification: 1.It is noted that the length of a branch should be greater than or equal to its width.Thus, a protrusion of a shape is considered as being perceptually insignificant if its length is shorter than its maximum local width.The shortest trabecula is one whose length is equal to its maximum local width (Figure 9).
2. If the length of a protrusion skeleton laying outside the junction segment is less than its length inside the junction segment, it is considered insignificant and hence, can be removed (Figure 9).In the proposed method, these lengths are approximated using first skeletal line segment in each branch.
The characteristic point of a junction polygon which is comprised of two or more adjacent J-Ts has to be obtained via a different approach.The good continuation principle of human perceptual organisation [29] suggests that visual elements that require the fewest number of changes will be grouped into wholes.Thus, characteristic point of a junction polygon can be approximated by assuming that trabeculae do not change their orientation after entering the intersection.As such, the new characteristic point of the junction polygon is defined as the mean value of the center of junction segments weighted with the number of branches entering each junction segment.For several images, this procedure was repeated consecutively and terminated in the case that no artifact was found in the skeletonised image.The removal of all artifacts was validated for all loading stages after 3 iterative operations using visual inspection (Figure 8 (B)).

Trabecular Parameters Calculation
Bone surface fraction (BS/TS), trabecular width (Tb.Wi), trabecular spacing (Tb.Sp) and trabecular number are the most commonly used parameters to quantify structure of trabecular bone.In our method, they are measured as follows.
• BS/TS BS/TS may be calculated in two ways.The first is a mesh-based BT/TS which is done by dividing the total area occupied by bone segments by the total area of the 2D projection of trabecular structure: The elastic modulus and strength of trabecular bone are primarily determined by volume fraction, the variations of which cause the heterogeneity of trabecular bone [24].This index is probably the most known and the most used output of a micro-CT analysis performed on bone.It indicates the fraction of a given area of interest (i.e. the total area) that is occupied by mineralised bone (Bone area (BS)).
Quantitatively Structuring A Trabecular Network Using Computational Geometry Techniques

•
Tb.Wi Tb.Wi is calculated as the area weighted mean width using the area (A i ) and width (w i ) of individual trabeculae. .
Tb.Wi represents the mean thickness value of all the beams of trabecular bone of the 2D structure.
• Tb.Sp Area weighted mean can also be used to estimate the average separation of the struts when applied to the inverse of the structure.Tb.Sp here is defined as the mean width of the inter-trabecular spaces or pores.
• Tb.N Trabecular number, Tb.N is taken as the inverse of the mean distance between the mid-axes of the structure to be examined.

BS TS Tb N TbWi =
Tb.N gives an idea of the percentage of trabecular bone present in the volume against the mean trabecular thickness.

Numerical Results in Fracture Process
To study the failure behaviour of trabecular architecture, time-lapsed 2D images of a cancellous at 6 stepwise micro-compression loading were acquired [30].In the uni-axial compression, the mounted specimen is subjected to a set of progressive mono-cyclic nominal strains at a constant strain rate.The specimen is allowed to relax after each incremental strain.This process is continued until fracture occurrence at the fourth stage.
The BS/TS of trabecular structure for different stages of loading was measured and summarised in Table 3.It can be seen that the incremental load applied to the specimen increased the value of bone surface fraction for the entire structure before fracturing.This is mainly due to reduction in the area occupied by pores in trabecular structure.However, it can also be justified by 3D nature of trabecular structure.Applying the compression load to this structure results in displacement of trabeculae in three different dimensions.As such, movements along the dimension perpendicular to the imaging plane results in variation of BS/TS at different stages of loading.From Table 3 we see that our vector based calculation has the same trend with the voxel based method.
Trabecular width for different stages of loading was measured and along with surface fraction was used to calculate trabecular number (Table 4).During fracture process, the value of Tb.Wi increased, which can also be explained by 3D displacements of individual trabeculae.Despite the same trend in values of Tb.Wi and BS/TS, increase in the value of Tb.N indicates a higher load sensitivity of bone surface fraction, and hence, reduction of average trabecular width when compared against BS/TS.To study the effect of heterogeneity in fracture mechanism of trabecular structure, the resulting images before loading (A) and after fracture occurrence (B) were shown in Figure 10.The initial structure of trabeculae was divided into sub-regions perpendicular to the direction of loading.The size of the sub-region was chosen based on minimization of the least subregion BS/TS in trabecular structure.The deformation of the subregion with the minimum BS/TS relative to the other subregions comprising the specimen (including the deformation of individual trabeculae) was tracked during sequential compression steps to identify the weakest (least rigid) subregion of the specimen.The first and the last sub-region were excluded from the measurements in order to avoid the effects of end attachments in the calculations of micro-structural indices.The BS/TS for different regions are shown in Table 5. BS/TS measurement for different sites demonstrated that the least rigid segment of the specimen where the most deformation occurred had the minimum value for BS/TS.The same results were reported by previous studies [31].The change in geometrical parameters of individual trabeculae during fracture process for different subregions and orientations was computed and summarised in Table 6.As can be seen, longitudinal trabecula were greater in number and thinner than trabeculae in transverse direction.No significant change in the length of longitudinal trabeculae was observed.However, majority of these trabeculae thickened during the fracture process.In oblique direction, 66% of total trabeculae oriented in this direction, were stretched as a result of incremental load.The main changes in trabeculae geometry were observed in transverse direction.All the trabeculae oriented in this direction were stretched with 17% average nominal strain.Half of these trabeculae experienced 20% reduction in average width, while the rest thickened with 7% mean nominal strain.

Conclusion
Many researchers have examined the effect of changing bone density on the structural properties of whole trabecular structure, but few have evaluated the biomechanical consequences of specific changes in trabecular morphology due to compression.
In this paper, structural histomorphometric parameters for morphology measurements are obtained in two ways: mesh based measurement at multiple locations by Delaunay Triangulation or voxel based calculation.Changes in dimensions of individual trabeculae subjected to incremental load have been studied, especially, the values of the most significant structural parameters of the cancellous bone, BS/TS, Tb.Wi, and Tb.N in different stages of incremental loading were calculated.The values of Tb.Wi obtained from traditional voxel based approach were on average 20% larger than the results from mesh based method.The Tb.Wi obtained from voxel based method can be a good measure for comparing the thickness of trabecular structures.However, as this is a scalar this measure may not be able to describe all structural changes.This calculation does not differentiate between length and width, and is statistically equivalent to the mean of all intercept length measures.Hence, due to larger number of thinner horizontal trabeculae as compared with the number of thicker vertically oriented trabeculae, Tb.Wi falls between the mean true width computed for the horizontal and vertical trabeculae.Thus, it could be useful to rather compare the histograms, where the two cases could clearly be distinguished.While the width obtained from the proposed method is based on true width estimation of individual trabeculae.
It was also observed that least rigid segment of the specimen where the most deformation occurred had the minimum value for BS/TS.Moreover, thinning of trabeculae did not occur uniformly during overloading of trabecular bone.Trabeculae oriented along transverse direction tended to be perforated, become thinner, or eventually disappear while those oriented longitudinally tended to retain their thickness.

Figure 2 .
Figure 2. a) Trabecular bone image obtained from microCT; b) Binary representation; c) boundaries of individual trabeculae; d) Bone and marrow segments

Figure 4 .
Figure 4. Different branches extracted are shown in different colours

Figure 5 .
Figure 5. Local length and width of individual trabeculae Total area of each trabecula is calculated form summation of rectangular areas (Figure 6) obtained via local measurements of each skeletal line segment.

Figure 6 .
Figure 6.Local rectangular sub-region of individual trabeculae

Figure 7 .
Figure 7. Characteristic points of bone segments shown in blue colour

Figure 8 .
Figure 8. Medial axis of trabecular structure before (A) and after (B) artifact removal

Figure 9 .
Figure 9. Maximum local width of a protrusion (bold black) and its external (blue) and internal (red) medial axis

Figure 10 .Figure 11 .
Figure 10.Trabecular structure obtained from microCT and binary representation before loading (A) and after fracture occurrence (B)

Figure 12 .
Figure 12.Orientation profile of a sample trabecula shown in Figure 10

Table 2 .
Segment typesTriangle typeNo of adjacent bone segment

Table 3 .
Nominal strain and BS/TS of different stages of loading

Table 4 .
Tb.Wi and Tb.N at different stages of loading

Table 5 .
Nominal strain and BS/TS of different sub-regions in trabecular structure