Skeleton-based 3D Femur Shape Analysis System using Maximum-Minimum Centre Approach

Femoral shaft fractures are correlated with frequent morbidity and mortality. It is a major musculoskeletal disorder caused by tremendous force applied to the femur. One of the most common surgical treatments for fixation is intramedullary nailing (IM), which utilises a specially designed metal rod and screws to be implanted into the medullary canal. However, severe bowing of the femur can result in a mismatch between the IM nail and the alignment of the femur. Such mismatch is a risk factor for anterior cortical perforation off the distal femur with subtrochanteric fractures and leg length discrepancy with fractures of the femoral shaft. Therefore, accurate three-dimensional (3D) preoperative planning is mandatory to facilitate the implant s design based on the obtained geometric data, especially for the fractured bowed femur. This paper presents an automatic 3D femur shape analysis system (3D-FSA) based on the extracted skeleton of each individual patient to provide an accurate 3D preoperative simulation. The structure of the 3D femur was generated using a set of computed tomography (CT) images. By using the maximum-minimum centre approach for skeletonization, significant geometric and topological information was captured. The proposed approach can potentially assist in implant measurements.


I. INTRODUCTION
T HREE-DIMENSIONAL (3D) skeletonization is the process of generating a skeleton, sometimes called the curve skeleton. It provides an alternative to capture the inner structure of an overall complex 3D mesh. These computed skeletons consist of significant geometric and topological information that are used extensively to produce segmentation for various analyses and visualisations in medical imaging [1], robotics [2], reverse engineering [3] and video surveillance systems [4].
For decades, many research efforts have focused on finding the optimal approach to generate skeletons from these digital meshes. These approaches are mainly divided into three categories: (i) Medial Axis (MA) [5], [6], (ii) Reeb Graph [7] and (iii) Geometric Contraction [8]. This paper mainly focused on the medial axis, which is a set of curves defined as the locus of the centres of balls that are tangent to and maximally inscribed in the 3D mesh, that is, such balls contained in the mesh are tangential to three or more points of the mesh without crossing the boundaries. Hence, the medial axis transform (MAT) is the combination of the medial axis and the radius function that determines those maximally inscribed balls, and the medial axis skeleton is generated.
The concept of the medial axis was first presented by Blum [9] for biological shape analysis and is applied in this research. Since it is difficult to compute MAT accurately, there are several strategies used to construct a skeleton that approximates the medial axis for different purposes, such as Voronoi diagrams [6], [10], 3D voxel grids [11], [12] and others. These strategies are common and difficult to implement for skeleton representations.
In this paper, the maximum-minimum centre approach is proposed for the skeleton extraction. This method is an automatic and straightforward method that can be used, especially on the human femur. Although the development of the skeletonization is relatively well established in medical research, skeleton computation in 3D human femurs is relatively unexplored in the study of orthopaedics.
The human femur is the longest, heaviest and strongest bone in our human body. Different kinds of trauma with many forces can result in damage to this bone, such as in some motor vehicle accidents or motorcycle crashes. This can also happen in a lower-force accident, such as fall from slippery floor, ladder landing on foot among the older people due to their weaker bones or osteoporosis. There are usually two types of treatments for femoral fractures: nonsurgical and the surgical. Nonsurgical treatment requires only a cast to heal. This is mostly applies to very young children. In contrast, surgical treatment requires surgery to repair and heal the broken bones within 24 to 48 hours.
Open reduction and internal fixation (ORIF) is a surgical treatment that consists of two procedures under anaesthesia: open reduction and internal fixation [13], [14]. The first step is open reduction, where the orthopaedic surgeon will cut the skin and adjust the fractured bone to the normal position. The second part is internal fixation, where the bone fragments are held together with fixation devices such as plates, screws, stainless steel pins and wires until bone union [15], [16].
In addition, intramedullary nailing (IM) is the most common surgical treatment for certain diaphysial fractures of the femur [17]. During this procedure, a specially designed metal rod, also known as an intramedullary nail (IM nail), as illustrated in Fig. 1, is inserted into the medullary canal and passes across the fractured bone to keep it in position. Finally, screws are placed above and below the fracture to hold the femur in the correct alignment while the bone heals [18].
Nonetheless, orthopaedic surgeons often encounter severe problems such as iatrogenic fracture and complications during implant insertion due to exaggerated femoral bowing and narrowing of the diameter of the medullary canal [19]- [23]. Severe bowing of the femur can result due to a mismatch between the IM nail and the alignment of the femur. This mismatch is because most of the implants are designed for the femur that do not match the exaggerated bowing. Such mismatch is a risk factor for anterior cortical perforation off the distal femur with subtrochanteric fractures and leg length discrepancy (LLD) with fractures of the femoral shaft. LLD after IM nailing of femoral shaft fractures is common and is reported in 20 to 30 % of cases [24]. Several techniques are proposed in [25]- [28] to prevent LLD intraoperatively during femur fracture fixation. In addition, [29], [30] measured the LLD for a number of patients who had treated their femoral shaft fracture with IM nails for early intervention and repeat surgery.
Therefore, several morphological studies have been conducted on the differences in the magnitude of femoral bowing among the Asian population [31]- [34] and Western population [35]- [40]. In addition, the other significant parameters of the human femur, such as femoral length, diameter of medullary canal, radius of curvature (ROC), femur length FIGURE 1: An example of a specially designed metal rod, also known as an intramedullary nail (IM nail). and etc., were assessed using either roentgenogram [41]- [43], computed tomography (CT) [34], [44]- [46], 3D spaces [32], [47]- [51] or 3D printed models [21]. These obtained parameters from all walks of life, various races and ages are used as a reliable frame of reference during the preoperative templating and preoperative design of custom-made implant devices [22].
Apparently, femoral bowing can be determined by the measurement of the angulation between the proximal and distal quarters of the femoral diaphysis, and the cut-off value was defined between 5.25 • and 7 • [38], [52], [53]. Thus, a preoperative planning template in orthopaedic surgery is an essential prerequisite to estimate the correct nail diameter and length for the success of the procedures. Conventionally, the detailed surgery plan and measurement are written down as blueprint and performed on hardcopy radiographs using various methods. This practice has become less practical due to a lack of consistency in radiographic magnification [54] and the rapid development of computer-aided design (CAD) software in deriving 3D meshes from reconstructed lateral radiographs or images using CT images [33].
This paper presents an accurate and automatic 3D preoperative simulation called 3D-FSA, that produces a skeleton of the 3D femur using the maximum-minimum centre approach for the analysis of femoral shaft geometry, which focuses on the femoral bowing and femoral width [51]. These parameters are very useful for orthopaedics to accurately perform preoperative evaluation and engineers in developing the new intramedullary nailing system.

II. ALGORITHM OVERVIEW
There are four main modules for our proposed approach, as illustrated in Fig. 2. It begins with the preprocessing that takes individual snapshots of the individual human femur into a set of cross-sectional images, which is also known as CT imaging. These images were approved by the Institutional Review Board of the participating centre and saved in Digital Imaging and Communications in Medical format [51]. Then, they were imported into 3D modelling software (AVIEW Modeller; Coreline Soft, Seoul, South Korea [55]) to produce 3D samplings of anatomical elements for the human femur. With the use of reconstruction and parametrisation from such datasets, the structured data can be obtained to form a 3D mesh representing both the bone surface and the corticocancellous interface in the 3D-FSA, as shown in Fig. 3.
The structured data used in this research are in .obj format, which is explained in Section III. Once the 3D femur is con-structed, its skeleton is obtained in the module of 3D skeletonization (explained in Section IV) using the maximumminimum centre approach. Thus, the obtained skeleton is used to perform the necessary shape analysis due to its compact representation of the femur. More detailed descriptions for each module are presented in Section III and Section IV.

III. 3D FEMUR CONSTRUCTION
The 3D-FSA is a shape analysis system that is specialised in representing and manipulating triangulated surfaces for both the human bone and the corticocancellous interface using the maximum-minimum centre approach (Section IV). The system provides a support platform for 3D femur formation that handles meshes in .obj file format, as shown in Fig. 3.
The .obj file format is a geometry definition file format that is equipped with all the necessary data formats, namely, the position of each vertex, the respective normal value for each vertex that represents the reflection of lighting or Phong shading, and the list of vertices that define the faces for the 3D femur; thus, the analysis and experiment for this research become more accurate, consistent and robust. The illustrations for the 3D femur in vertices view can be seen in Fig. 3(a), and in the faces view can be seen in Fig. 3 According to [56], .obj is one of the crucial file formats in both 3D printing and 3D graphics applications. It was originally created by Wavefront Technologies for its Advanced Visualiser application to store geometric objects composed of lines, polygons, and freeform curves and surfaces. In addition, the colour and texture information can be stored in this file format too. Nevertheless, no scene information, such as light position or animations, can be stored in it.
The .obj file format is simple and open format with wide export and import support among Computer Aided Design (CAD) software. For example, the Autodesk 3ds Max is used to convert the rendered 3D mesh into .obj file format and read it using the Notepad in Windows or TextEdit in MacOS. Hence, .obj becomes famous as a file format for 3D meshes.

IV. 3D SKELETONIZATION
This is the most significant module in our algorithm after Section III, which calculates the MAT and produces the skeleton from the 3D femur. MAT was introduced by Blum [9] and is mainly composed of two properties: the medial axis (MA) and the radius function (refer to (1)).
Let Ω be a connected bounded domain in R with n−dimensions and B to be the loci of all maximal inscribed disks that meet two or more boundary points without crossing any of the boundaries in Ω. We define the MA, denoted by M A(Ω), to be the set of centres of the disks in B [57] as written in (2).
Each vertex p in the M A(Ω) is defined as a set of the pairs consisting of the centre and the radius r of the disks in B r (p) ∈ B, which forms the medial circle, and the volume is enclosed by the surface in Ω, which is exactly the union of these circles, as presented in (2).
Therefore, to obtain the skeleton of the 3D femur, we find the maximum radius of the circle that is well suited for each slice of the irregular inner 3D femur. Before further explanation, we review some basic facts about graphs in the context of 3D skeletonization. A bidirectional weighted graph G = {V, E} consists of a set of vertices V and a set of bidirectional edges E that connect them. The 3D femur is composed of a large number of faces F that are made of three connected edges each, or in other words, a triangle.
In this paper, 4 processes are involved, as depicted in Fig. 4, to compute a reliable and consistent skeleton from a 3D femur. A more detailed explanation is presented in the following subsections.

A. AXIS ADJUSTMENT AND FEATURE DETECTION
At the beginning of the program, this process involves detecting and recognising the distinctive features of the 3D femur. It examines every vertex to decide the position of the distal femur (the area of the leg just above the knee joint [18]) and proximal femur (includes the femoral head, neck and the region 5-cm distal to the lesser trochanter [58]).
With the obtained measurements, the dimensions, including the width, height and depth of the 3D femur, are determined. Moreover, the (x, y, z)-axes are modelled in the coordinate system of 3D-FSA using the appropriate geometric transformation, such as translation. Fig. 5 illustrates the outcome produced by this process.
In addition, two bounding boxes that are drawn as rectangular (pink colour), as shown in Fig. 5, are constructed. These bounding boxes are served as the indicator of the beginning (the distal femur) and ending (the proximal femur) bases for the 3D femur.
The green cylinder shape made of two separate circles with the calculated distance h, between these circles is defined as the sliding window (depicted in Fig. 6). The value h can be varied depending on the decision from the user.
The purpose of this sliding window is to slide over the whole 3D femur from the beginning to the end to capture and slice different portions of it.
Each slice that is captured by the sliding window is used to collect the adjacent faces within that portion, as shown in brown in Fig. 6, and perform some calculations in Section IV-B and Section IV-C.    h. During this process, the vertices v z within this selected portion from z to z + h are collected and added to V s (refer to (3)).

B. MAXIMUM-MINIMUM CENTRE APPROACH
Then, the adjacent faces to each vertex in V s are captured and stored in F s with the total number of m faces. We take the advantage of the Lindstrom approach [59], which introduced the simplex operators for fast and memory-efficient 3D mesh simplification. The equation is denoted in (4) to acquire these adjacent faces. Fig. 6 illustrates the sliding window (cylinder shape in green colour) and the adjacent faces within the selected portion (highlighted faces in brown colour).
where  With the acquired adjacent faces in F s , the area A(f i ) for each face i in F s and its centroidf i are calculated using the cross product and the formula in (5), respectively. This is followed by multiplying both calculated values, to derive the weighted area for each face i using (6). The weighted normal is also derived by multiplying the face area A(f i ) together with the face normal N (f i ), as defined in (7).
where N (f i ) = normal of a f ace f i Then, the weighted centre verticesF s as defined in (8) for all the faces in F s are derived by dividing the accumulation of the weighted areas in (6) and the accumulation of the areas A s . Moreover, the accumulation of the weighted normal for all the faces in F s is calculated to be used as the numerator of the fraction as denoted in (9). Both the values obtained from (8) and (9) are the average of the areas and the common orientation of the faces in F s , respectively.
where m = the total number of f aces in f i .
The weighted centre verticesF s define the medial axis of the 3D femur as plotted using small and green spheres in Fig.  9. It defines the loci of all maximal inscribed disks that meet the inner boundary without any overlap, namely, the skeleton, as portrayed in Fig. 15. In addition, F s is crucial in finding the thickness of the medullary canal in Subsection IV-C and the radius angle calculation in Subsection IV-D.
Whereas the weighted normalN s is used to find the residual error rates for the average of the weighted normal as in (10). These values are used to derive the equation for the least square plane fitting, which is then used to render the projection plane in the centre, as illustrated in Fig. 8 (peach  in colour).
During the sliding process, the sliding window serves a slice plane, which slices the 3D femur from the beginning until the end, with a thickness of h, as shown in Fig. 8 (the 2 green circles). The vertices V s within this sliding window are all projected onto the formed projection plane.
The blue vertices in the centre as shown in Fig. 8 and Fig.  9, represent the projected vertices from the sliding window during the sliding process. We call these projected vertices the inner slice polygon Ω s .
The inner slice polygon Ω s which designates the medullary canal of the 3D femur. The explanation of how to find this thickness is explained in Section IV-C. VOLUME 4, 2016 FIGURE 8: The sliding window (2 green circles), the slice plane in rectangle shape (grey in colour) and the projection plane in rectangle shape that located in the centre (peach in colour). Finally, the maximum-minimum centre is acquired from the results generated from the weighted centre verticesF s and the maximum-minimum thickness of the medullary canal (Section IV-C as illustrated in the maximal green circle B s (p, r) in Fig. 9.
The definition and notation used in this step are summarised as follows: 1) Slice plane : Calculate the weighted average of the face normal 2) Sliding window as base plane : Face normal average. 3) Inner slice polygon Ω s : Projection face centre to base plane. 4) Max-Min Centre : B s (p, r) is the maximal circle in Ω s .

C. MAXIMUM-MINIMUM THICKNESS CALCULATION
The anatomical parameters of the femur, such as the diameter (thickness) of the isthmus and the medullary canal, are significant in the preoperative measurements [19], [49], [60].
To obtain scaled measurements, the diameter (thickness) for both the isthmus and medullary canal was directly obtained by detecting the contours of the 3D femur in this step of the 3D-FSA.
A contour is the representation of a set of linking edges corresponding to a region boundary [61]. Contours are a useful tool for shape analysis and object detection and recognition. In 3D-FSA, contour N is represented as an ordered list of vertices that are used to construct the linking edges. Contour 1 and contour 2 of the medullary canal are defined as illustrated in Fig. 10 with the output produced in Section IV-B, such as the projected plane P , projected weighted centre vertices P (F s ) and the inner slice polygon Ω s .
As shown in Fig. 10, contour 1 is defined as the outer curve of the inner boundary of the medullary canal, and all the projected linking vertices a are saved in set A, whereas contour 2 is defined as the outer curve of the inner boundary of the medullary canal, and all the projected linking vertices b are saved in set B. The relations on these sets are depicted in (11).
The diameter (thickness) of the medullary canal is then measured by using the instructions written in Algorithm 1. Starting from vertex a i ∈ A of contour 1, compare it with all the vertices in set B of contour 2. In each comparison, only the minimum distance is acquired between these values and saved it in minDist, as illustrated in 11(a). This process is repeated until all the vertices a ∈ A are being compared. When a vertex a i ∈ A is changed in each execution, the maximum value maxDist is selected after comparison with the acquired minDist.
After all the vertices a ∈ A of contour 1 are examined, the algorithm continues with the vertices b ∈ B which go through the same comparison process, and the minimum distance minDist is collected. The minimum distance minDist is then compared with the maxDist to obtain the maximum distance of all, as depicted in Fig. 11(b).
There are three significant variables generated from this step: 1) c i , 2) c j , and 3) maxDist. These variables are used to determine the diameter of the medullary canal by finding the closest distance from contour 1 to contour 2, and the radius r is obtained using (12), which forms the maximal circle B s (p, r) of the 3D femur, as illustrated in Fig. 12.  Calculate dist = Distance(ax, ay, az, bx, by, bz) 9: if (dist < minDist)) then As defined in [62], any approximate circle's radius at any particular given point is called the radius of curvature (ROC) of the curve. As we move along the curve, the ROC changes. In differential geometry, the ROC is denoted as R. The maximalcircle B s (p, r) and the medullary curvatures of 3D femurs are produced in Section IV-C to determine the amount by which the curve derivates itself from being flat to a curve and from a curve back to a line. The ROC is the reciprocal of the curvature. This is another important preoperative parameter studied by many researchers [19], [21], [33], [34], [42], [46], [49], [51]- [53], [63]- [69] due to the difference between the curves of the femur and the contemporary femoral nails implicates the inadequacy of the implant design.
Therefore, 3D-FSA provides the computation for finding the curvature of the curve for the 3D Femur as written in Algorithm 2. The basic concept behind this computation is as follows: "The circle defined by centre (a, b) and radius r will yield the least mean square value for the expression (x − a) 2 + (y − b) 2 − r 2 and the circle's curvature will be 1 r ". Fig. 13 demonstrates two examples for the involved range of this calculation with various starting indices but similar ending indices. Additionally, Fig. 14 visualises the projected vertices used for the calculation of the minimum radius of curvature, together with the projected angle and projected plane in faces view and vertices view.

V. RESULTS AND DISCUSSION
The results generated by 3D-FSA can be viewed in Fig. 15 and utilised by the authors [51] in their morphological study of the femoral geometry that focuses on the bowing and width among the Korean ethnicity for the age range between 20 and 89 years old [51], without the concern of any implant, deformity, or surgical history of the femur. There were a total of 1400 participants, of whom 2800 femurs were enrolled in the study.
The collected CT images of both femurs are fully scanned, and the conversion to 3D file formats is completed with 3D modelling software. The 3D-FSA then reads the file and constructs a 3D femur. Next, 3D skeletonization is applied to obtain a compact representation of the femur, i.e., the skeleton of the femur. Both the 3D femur and the skeleton are used to obtain the following parameters (together with the appropriate defined location) for the analysis of the femur shaft: The main three parameters are depicted in Fig. 16. To validate the measurement of 3D skeletonization, 50 participants who were evenly distributed by age and sex were selected to repeat the measurement after 2 months. Moreover,  the intraclass correlation coefficient of each parameter was calculated, and the results showed that the 3D-FSA produced an accurate measurement compared to the two-dimensional measurement. reproducible results [51]. This study also revealed that the age-and sex-related factors affect the femoral bowing and the diameter of the medullary canal. The results showed that women had more bowed femurs than men. Moreover, as the age increases, the femur becomes more curved, and this tendency is more apparent in women. As a result, these patients are more vulnerable to insufficiency fractures. Several studies [19], [33], [34], [37], [70], [71] have shown that Asians have a smaller ROC than Caucasians, especially in women, and that femur bowing is correlated with race. [21], [51] presumed that if 700 mm is the minimum ROC possible with the current nails, the current intramedullary nail would result in a mismatch problem in 11.5% of the Korean population because of their severely bowed femur.
Their results demonstrated that during the ageing process, both the outer diameter and the diameter of the internal medullary canal were increased among women older than 50 years. Moreover, the linear regression analysis showed that femoral bowing was not related to medullary canal widening but was related to ageing rather than osteoporosis.
However, there are several limitations in the study as follows: 1) The collected data are not representative of each individual change but rather a representative based on their ages.
• The average height of men increased from 169 cm for those who were born after the Korean War in 1954 to 175 cm for those who were born in the year of 1983 [72]. • [73] presented that there was a good correspondence between the rate of bone loss from a sample in a cross-sectional study and from longitudinal data over a 16-year of observational period.
2) The measured data did not have the resolution to assess cortical porosity or change because the images used were conventional CT images rather than micro-CT data.
To enhance the efficiency and accuracy of the 3D-FSA, 3D mesh segmentation and artificial intelligence (AI) techniques should be deployed with a sufficient amount of data. As created AI has globally-attributed to mostly one goal, it obtains a higher level of accuracy and surpasses the existing benchmarks. More data will result in better and more accurate models.
The 3D-FSA would provide a great education platform for better understanding of the femur geometry. In addition, it provides a detailed guidance in preparing and developing a new intramedullary nailing system.

VI. CONCLUSION
An automatic 3D-FSA based on the extracted skeleton of each individual patient is developed and used to provide an accurate 3D preoperative simulation. Moreover, 3D-FSA has provided an analysis platform for the geometry of the femoral shaft in assisting the morphological study that focuses on the bowing and width among the Korean ethnicity for the age range between 20 and 89 years old [51]. A total of 2800 femurs were enrolled in the study, and it was concluded that femoral bowing and the width of the medullary canal gradually increase over time, especially among women. The results of these 3D analyses generated by 3D-FSAs are accurate and reproducible. Therefore, it is the best reference for surgical preparation and implant designs.