An open-access hip joint model repository suitable for finite element method

Background and objective: population-based ﬁnite element analysis of hip joints allows us to understand the effect of inter-subject variability on simulation results. Developing large subject-speciﬁc population models is challenging and requires extensive manual effort. Thus, the anatomical representations are often subjected to simpliﬁcation. The discretized geometries do not guarantee conformity in shared interfaces, leading to complications in setting up simulations. Additionally, these models are not openly accessible, challenging reproducibility. Our work provides multiple subject-speciﬁc hip joint ﬁnite element models and a novel semi-automated modeling workﬂow. Methods: we reconstruct 11 healthy subject-speciﬁc models, including the sacrum, the paired pelvic bones, the paired proximal femurs, the paired hip joints, the paired sacroiliac joints, and the pubic symphysis. The bones are derived from CT scans, and the cartilages are generated from the bone geometries. We generate the whole complex’s volume mesh with conforming interfaces. Our models are evaluated using both mesh quality metrics and simulation experiments. Results: the geometry of all the models are inspected by our clinical expert and show high-quality discretization with accurate geometries. The simulations produce smooth stress patterns, and the variance among the subjects highlights the effect of inter-subject variability and asymmetry in the predicted re-sults. Conclusions: our work is one of the largest model repositories with respect to the number of subjects and regions of interest in the hip joint area. Our detailed research data, including the clinical images, the segmentation label maps, the ﬁnite element models, and software tools, are openly accessible on GitHub and the link is provided in Moshfeghifar et al.(2022)[1]. Our aim is to empower clinical researchers to have free access to veriﬁed and reproducible models. In future work, we aim to add additional structures to our models. © 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The finite element (FE) method brings hip joint models to life by simulating their behavior under given conditions.This method is a numerical approach to solving partial differential equations by discretizing the hip joint domain into a finite mesh [2] .The ap-Fig. 1.The LibHip modeling flowchart; The workflow is aimed at minimum user interaction and proper FE model features.The blue box presents the shape-related process, and the green boxes indicate discretization-related processes.Our approach starts from the 3D representation of bones and generates the volume and surface mesh of the hip joint area.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)Fig. 2. The overview of our modeling features; Middle: the regions of interest cover the sacrum, the paired pelvic bones, the paired proximal femurs, and the three joints connecting these bones together; Left : all the bone-cartilage interfaces have conforming meshes.Right : congruent interfaces between the articular cartilages in the hip joint.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Problem statement
Despite the advantages of FE analyses, most studies assessing hip joint behavior are based on models from a single or a few subjects.This approach does not account for the effect of the intersubject variability and leaves uncertainties regarding the generality of the predicted results [13] .On the other hand, populationbased studies are limited by clinical resources and the significant amount of time and manual effort to prepare large-scale hip joint models.The area around the hip joint exhibits multiple structures in contact, and each structure is typically segmented from clinical images, transformed into surface meshes, and filled with volume mesh elements.This process is time-consuming and usually requires extensive manual work by FE specialists.
A common requirement in generating a volume mesh is an explicit definition of the surface mesh, meaning the surface has to be closed, non-intersecting, and manifold.While these requirements are necessary for identifying the interior part of the surface mesh, they are usually not satisfied by the quality of the segmentation results and may require extensive manual post-processing.Further, most FE formulations require volume meshes to share faces on adjacent surfaces in the hip joint area, a property that is challenging to enforce for many volume meshing tools.If this property is not satisfied, additional and often insurmountable boundary conditions appear in the contact formulation, which might lead to invalid simulation results.

Contributions
We present a hip joint model repository designed for FE simulation studies.We are combining cutting-edge geometry processing algorithms to reduce the amount of user interaction in developing FE models and tackle the shortcomings of conventional methods.Fig. 1 summarizes our FE modeling workflow.The significant technical contributions and unique traits of these models are outlined below and illustrated in Fig. 2 .
• Open-access data aligned with the FAIR guiding principles [14] : our population study includes 11 male subjects with no diagnosed disease related to the hip joint.All the models, including the clinical images, the segmentation label maps, and the FE models, together with the algorithms to prepare them, are available on Github with the link provided in Moshfeghifar et al. [1] .The vision of this approach is to empower clinical researchers to have easy access to verified and reproducible FE models.A senior consultant radiologist has clinically validated the geometry of all the models.The FE model's simulation performance is verified by running FE analyses and analyzing the estimated results.• A full pelvic girdle with both the hip joints: all the bones in the hip area encompassing the sacrum, the paired pelvic bones, and the paired proximal femurs are segmented directly from CT scans.Having the bone models from both sides of the body preserves bilateral anatomical information in each subject and considers the effect of the pelvic girdle in the hip joint behavior.
• Subject-specific cartilages: we reconstruct the sacroiliac joints, the pubic symphysis, and the articulating hip joint cartilages semi-automatically using the bone geometries [15] .This approach provides fast, subject-specific cartilages with nonuniform thickness.Modeling cartilage in the inter-bone cavities ensures inter-connectivity of the weight-bearing path in the pelvic girdle.Our joint geometries maintain conforming boundaries in their bone-attached interfaces.Further, we provide congruent interfaces between the two articular cartilages in the hip joint.The high congruence level prevents potential spikes and peak stresses in the articulating surfaces, resulting in continuous stress distribution.• Highly accurate discretization: we employ a multi-body volume mesh generation that simultaneously generates the whole complex's volume mesh, ensuring neither overlapping nor gaps in the interfaces.This method welds the interface nodes in the meshing step, avoiding further contact definitions in the simulation setup.We use the Hausdorff distance to measure the geometrical accuracy in our models.

Background
Assessing the hip joint behavior through FE analysis requires three main steps: the first is to generate a proper approximation of the anatomical structures; Second, this anatomical information is transformed into surface and volume representations to build a FE model; And lastly, the behavior of this model is studied under a chosen analysis scenario.To gain realistic predictions, we need detailed representations of the hip joint area, including soft and hard tissue, and a high-quality discretization to reduce the prediction errors [16] .

Anatomical terminology
The skeletal structure in the hip joint area is formed by the pelvic girdle reaching the femoral bones.The pelvic girdle consists of the paired pelvic bones and the sacrum in the middle [17] .The sacrum reaches the pelvic bones in the sacroiliac joint , and the two pelvic bones connect at the pubic symphysis .The acetabular part of the pelvic bone articulates against the femoral head to form the hip joint .The lunate surface of the acetabulum and the femoral head, except the fovea pit, are covered by the hip articular cartilages [17] .

Hip joint modeling approaches
Computer-aided design models are simple solutions to capture the overall shape of the hip joint area, e.g., represent the

Table 1
The comparison between some of the most related population-based studies and our work; Most of these works only cover one side of the body and increase the subject numbers by manipulating the existing models.Our work is one of the largest model repositories regarding the number of subject-specific data and the region of interest in the hip joint area.hip joint as a perfect ball-in-socket joint.Even though convenient, such simplifications overestimate the stress predictions in the hip joint [18,19] .Subject-specific modeling is another popular choice that solves some of the issues arising from computeraided design models, thanks to having detailed anatomical information.The skeletal structures are mainly reconstructed based on clinical image modalities such as computed tomography (CT) and magnetic resonance imaging (MRI) [4,[20][21][22] , using manual [18,23] , semi-automatic [24,25] , or fully-automated segmentation approaches [26] .Further, subject-specific modeling approach accounts for articulating cartilages.The cartilage thickness and interfaces are often either delineated from the images, assigned uniformly, or approximated by a radius representing the smoothed joint space midline [18,23,25] .Given the tight space in the hip joint, it is difficult to identify the opposing cartilage interfaces [27] .CT arthrography and manual traction increase the joint space and enhance the cartilage visibility at the cost of being an invasive intervention [23,24] .In this study, we semi-automatically segment the bone geometries from CT scans and use a fast geometry processing method to create subject-specific cartilages from bone geometries.This method is based on our previous work to create cartilages where cartilage segmentation is not feasible [15] .FE analysis over a single hip joint provides valuable information regarding that specific model [20,22,28] .However, the hip joint shape and the relative alignment of structures vary among individuals, and adding more models in the simulations allows us to account for this variability.It is suggested that there is a significant variation in the hip joint contact mechanics due to the geometrical differences [19,29] .Currently, a limited number of works use multiple hip joint models to study the effect of inter-subject variability, and in most cases, the models are not freely available.Table 1 compares the most related studies to our work.

Study
Preparing population-based models is challenging: a limited number of models or clinical images of healthy hip joints are available, and the ethical issues regarding radiation exposure prevent custom image acquisition.Unlike Henak and Harris et al., who have performed custom image acquisition on both healthy and pathological subjects [19] , most studies only have access to either healthy or diseased images.Thus, they generate more models by post-processing the existing models.For instance, Zhao et al. made multiple diseased models by changing the morphological parameters of a single healthy model [4] .Russel et al. used 11 diseased subjects and built an age-matched healthy hip model by deforming the Visible Human Project model [31,33] .These approaches fail to provide geometries arising directly from subject-specific images.Our work uses 11 healthy images from The Cancer Imaging Archive.
We have chosen open-access images to ensure reproducibility and exterior validity in our method [34] .
The region of interest is another challenging issue in population-based modeling.Most modeling approaches do not necessarily include the whole anatomical structures in the hip area.For example, Henak et al. and Harris et al. have chosen to model the hip joint from only one side of the body.The geometry of the rest of the pelvic girdle is ignored, and their effect is presented by boundary conditions [19,29] .To make bilateral models, Zhao et al. assumed the hip area to be symmetrical along the sagittal plane and flipped the model of one side to generate a symmetrical model [4] .Although these studies provide fundamental predictions of the hip joint behavior, their modeling assumptions are simplified with no consideration of the whole hip area and the bilateral variation, which may affect the hip joint overall behavior.In our work, we provide models from both sides of the body which contain bilateral geometrical information in each subject.Additionally, we present the geometrical variation within our population study using well-known anatomical metrics.

Surface and volume discretization
Segmented anatomical structures are transferred to surface meshes using the marching cube algorithm, and ideally, they should be water-tight two manifold surfaces.Unfortunately, this is rarely the case for clinical data: errors can occur in the segmentation due to insufficient input data, the presence of small features, or numerical issues in the contouring methods.Traditionally, these imperfections need to be cleaned up before volumetric meshing, a laborious task, even when using state-of-the-art remeshing tools [35] .
When all surfaces are clean, a volume mesh is created by filling the interior part of each surface.Tetrahedral elements are the most commonly used types of finite elements due to their good performance and ease of generation [36][37][38] .Various approaches, such as Delaunay meshing, advancing front, or grid stuffing, are used for volumetric meshing.Due to the high failure rate of these methods, the surface cleaning must be iteratively repeated until the meshing algorithm succeeds on all the sub-volumes [37] .Our work adopts a radically different approach to remove this manual and time-consuming step: using fTetWild, we create a single discretization for all surfaces, including their imperfections, and we then extract a consistent volumetric partitioning only after the volumetric meshing is performed [38] .This approach is robust against bad input surface triangles and automatically deals with holes, selfintersections, or inconsistent inside/outside orientation.Another property of fTetWild is that the output mesh resolution is decoupled from the input resolution and is controlled by an explicit parameter measuring the maximal tolerated geometric error [38] .This is ideal in our setting, as it allows us to generate coarse meshes starting from the fine surface meshes extracted from clinical images.In contrast, traditional approaches have limited control over mesh resolution, as they cannot easily coarsen the surface mesh [39] .
Further, conventional volumetric meshing tools mesh the different domains of the hip joint multi-body separately in a one-by-one approach.In contrast, we process all the surfaces of our multibody model simultaneously , considering them as faces in a single, conforming tetrahedral mesh, and thus allowing them to overlap, have gaps, or other imperfections.We then perform inside/out tests and constructive solid geometry operations as a filtering operation on the resulting tetrahedral mesh, using signed distance as a geometrically robust filter.This approach benefits from naturally guaranteeing conformity between adjacent surfaces in a multibody model.We refer to [37,38] for details on the fTetWild meshing algorithm.Fig. 3 compares the bone-cartilage interface between the one-by-one and simultaneous volume mesh generation approaches.
By conforming interface, we mean that the bone-cartilage boundaries have the same discretization in the shared interfaces.If the discretization connectivity differs in the shared interfaces, but the interface shape is still the same, we have a congruent interface.Fig. 4 illustrates the conforming and congruent interfaces in the hip joint.The bone-cartilage interfaces conform, while the articular cartilages have congruent interfaces.Conforming interfaces allow us to merge the discretization of the two sub-domains.This is pre-ferred in simulation methods as it reduces the number of variables and the need for adding extra constraints to keep vertices together during simulation.Congruent interfaces in the hip joint allow a smooth load transition between the articular cartilages.Joint congruency evaluates the encasing of the articular cartilages in the hip joint and is a way to compare healthy and pathological joints [19] .

Finite element method solvers
Various research libraries, commercial codes, and software are available to solve elastic problems with the FE method.However, they are sensitive to user-defined solver settings, unrelated to the physical quantities of interest.For example, contact parameters and the time step size can significantly affect the stability and accuracy of the estimated solutions and even lead to inverted elements.Additionally, in some contact models, such as the ones used in Abaqus and FEBio, the user must assign pre-defined contact surfaces.These two surfaces must have specific and different mesh resolutions [40,41] .FEBio also requires an initial slight penetration between these surfaces to detect the contact surfaces robustly [41] .
PolyFEM, in contrast, adopts the Incremental Potential Contact (IPC) formulation [42,43] .This algorithm requires no manual selection of the contact surfaces and no prescribed mesh resolution.This solver is robust to large deformations and ensures that no elements are inverted using explicit line-search checks [42] .The downside is that it requires an initial configuration free of penetrations (and it preserves this condition throughout the simulation [44] ).Further, PolyFEM allows adaptive p-refinement on Lagrangian tetrahedral elements, which we use to increase accuracy by employing a higher-order basis in thin cartilage layers [45] .In this work, we want to generate models compatible with both strategies.We thus generate two sets of FE models for each subject: with and without a gap between the articular cartilages.We show FE solutions computed using both approaches in Section 4.3 .

Material and methods
In this section, we explain our modeling workflow, illustrated in Fig. 1 .We start from the data which we use to reconstruct the anatomical structures.Next, we describe the volume mesh generation step and explain the steps to ensure high-quality meshes.Finally, we describe the FE setup we employ to test the quality of our models.

Image data
The input to our modeling workflow is the surface mesh of the bony structures.We have chosen to create our input models based on CT scans as this modality is suitable for defining the location and topography of the bones in the hip joint area [46] .To ensure reproducibility, we obtain the scans in DICOM format from the open-access Cancer Imaging Archive [34] .
We select our subjects from the CT Colonography, the Lymph Nodes, and the Cancer Genome Atlas Urothelial Bladder Carcinoma (TCGA-BLCA) collections [47][48][49][79][80][81].The scanning procedure in these collections covers the hip joint area explained in Section 2.1 .The subjects are in a supine position during the image acquisition, which is the closest to an unloaded joint state.We chose 11 subjects of the same gender and age range, with no reported disease related to the hip joint, the slightest rotation in the body, high image resolution, and minimum image noise.We crop the CT scans to the hip joint area to minimize the computational load during segmentation.Next, the cropped images are stored in NIFTI format to store each subject as a single file in our repository while preserving all the essential metadata.Table 2 outlines the image properties and traits for each subject.

Table 2
The specification of image data used to create the input models; All images are selected from the cancer imaging archive.These images consists of adult male subjects with no reported disease related to the hip joint area.These data are provided to ensure complete transparency in our work.

Bone geometry reconstruction
We obtain an explicit surface representation of all the bones using a semi-automated approach implemented in the 3D slicer software package and a previous work by Xu et al. [ 50,51,78 ].This method entails initial region labeling, contouring [52] , followed by manual refinements to ensure accurate 3D approximations with no rough surfaces, holes, and irrelevant connected tissues.
These bone segmentation are verified by a senior consultant radiologist.The clinical expert initially scrolls through all the segmented slices in each subject and verifies the bone contours and the existing gaps in the inter-bone cavities.Then, he verifies the anatomical shape and smoothness of the reconstructed 3D surfaces.

Bone anatomical measurement
The shape of the hip joint area varies among our subjects.We provide the anatomical measurements of the bones and cartilages to quantify this variation and characterize each subject's size and shape.Most of the anatomical measurements in the literature are based on a 2D assessment of the location of measurement [53] or are based on cadaver skeletal measurements [54] .This procedure may not be sufficient for subject-specific 3D geometry analysis [55] .Therefore, we have chosen to obtain the subject-specific anatomical measurements from the 3D surface mesh.These parameters are illustrated in Fig. 5 .
We fit a sphere to the femoral head and choose the center and radius of this sphere as the hip joint center (HJC) and the simplified femoral head radius (SR) , respectively.To minimize the bias from manual fitting, we use a least-squares method for spherical fit [56] .We further measure the distance from the HJC to the actual femoral head surface mesh (AR) to see how the geometry of the femoral head deviates from a sphere.We define the visible femur length (VFL) as the length of the line connecting the most proximal point of the femur to the mid-point of the most distal part of the same bone; The femoral neck-shaft angle (NSA) measures the angle between the neck and shaft axes.We first extract the centerline of the femoral bone geometry using the VMTK extension in the 3D slicer software package [57] .Then, we apply a least-squares Linear Regression method to find the best fitting lines to the neck and shaft part of the centerline [58] .The angle between these two lines shows the neck-shaft angle.The width of the pelvis (PW) is explained as the distance between the most lateral point of the pelvic bones to the femoral head center; The i nter-hip separation (IHS) is assumed as the distance between the paired hip joint centers; The height of ilium (IH) is the vertical distance between the most supe-Fig.5.The anatomical measurements of the hip joint area; These parameters quantify the geometry of our models and measure the inter-subject and bilateral variability among subjects.Left: a red sphere is fitted to the femoral head.The center of this sphere is the hip joint center (HJC), and the radius (SR) is defined as the simplified femoral head radius.We define the distance between the HJC and the femoral head boundary as the actual femoral head radius (AR).The SR and AR comparison shows the lost subject-specific geometry after simplification.The visible femoral length (VFL) is the most extended visible length in the femur bone.The neck-shaft angle (NSA) shows the angle between the femoral neck and shaft axes.Middle and Right: observe the width of the pelvis (PW), the inter-hip separation (IHS), and the height and width of the ilium (IH, IW).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)Fig. 6.The comparison between the surface mesh output from 3D Slicer and remeshed surface mesh in fTetWild; The fTetWild result shows a clean mesh with no dense or erroneous triangles.The coarser mesh on the pelvic side guarantees robust contact definition in the FE analysis.
rior part of the pelvic bones and the hip joint center; The width of ilium (IW) is defined by the horizontal distance between the hip joint center and the most lateral point of the pelvic bones.We measure the left and right HJC, SR, AR, VFL, IH, and IW separately to capture the bilateral variation with respect to the sagittal plane.

Bone surface mesh re-meshing
As explained in Section 2.3 and illustrated in Fig. 6 , the bone surface mesh extracted from the CT scans are typically dense and may have poor quality, and other imperfections.We need to improve the quality and resize the triangles for two particular rea-Fig.7. The cartilage generation procedure for the sacroiliac joint and the pubic symphysis; The primary and secondary bones are shown with V P , F P and V S , F S , respectively.The cartilage base starts from the primary bone and extrudes to the secondary bone.The blue and green surfaces indicate the primary ( F D C ) and secondary interface estimations ( F E C ), respectively.These subsets are selected based on the distance to the opposite bone ( Eq. ( 1) ).The pink ring indicates the connecting mesh between the two interfaces ( F R C ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)sons: 1) our cartilage generation algorithm in Section 3.5 uses the bone surface mesh to generate the base of the cartilages [15] .Thus, the cartilage's surface mesh quality depends on the quality of its underlying bone mesh quality; 2) Some FE solvers enforce specific requirements on the mesh resolution for robustly detecting contact points [40,41] .
We use fTetWild to mesh the volume of the surface boundary and re-extract the surface of the volume mesh as the resulting remeshed surface [37] .This meshing method allows vertices of the boundary of the volumetric mesh to move within an envelope size of the winding number field level-set that defines the surface geometry [37] .We tune the envelope (ε) and the ideal edge length (l) parameters in fTetWild to provide different mesh densities between the articular joints in the hip joint.This guarantees coarser meshes on the acetabular cartilage.

Cartilage geometry reconstruction
We apply a specialized geometry processing method to generate subject-specific cartilages for the hip joints, the sacroiliac joints, and the pubic symphysis.This method was initially introduced by Moshfeghifar et al. [15] to generate subject-specific hip joint cartilages with conforming bone-cartilage interfaces and nonuniform thickness.Our work adds new ideas to this algorithm to improve the hip joint results and extend this method to the other two joints.
We model the sacroiliac joint and the pubic symphysis as single-piece cartilages, filling the inter-bone cavity.The hip joint is modeled as double-piece cartilage.The joint space is divided between the acetabular and the femoral layers, allocating roughly half of the joint space to each cartilage's thickness.

Basic modeling concept
Each cartilage is generated based on the shape and distance of the involved bones in the joint.The geometry of each bone is defined as (V, F ) , where V ∈ R N×3 is a set of mesh vertices, and F ∈ Z K×3 is a set of mesh faces; N and K denote the number of vertices and faces, respectively.The cartilage reconstruction starts from one of the bones, referred to as the primary bone (V P , F P ) , and grows toward the second bone, referred to as the the secondary bone (V S , F S ) .The steps are summarized below and illustrated in Figs.7 and 8 .C .This subset is selected based on the distance to the opposite bone.The initial primary estimation in the femur bone (yellow) does not comprehensively cover the femoral head.Thus, we apply a curvature-based region filling to grow the yellow cartilage to the correct portion on the femoral head (blue).We make a copy of the primary interface to create the top surface of the cartilage F E C .We initially extrude a part of this surface to a maximum of half the distance to the second bone (green) and then connect the extruded subsets via a smooth blend (pink).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Primary interface estimation: the main idea in this step is that only the sufficiently close parts of each bone have cartilages attached to them.Hence, we select an initial subset of the primary bone as the base of the cartilage based on its distance to the secondary bone.Our choice is based on the signed distance between the face barycenters of the primary bone and the secondary bone vertices.Provided the distance filter parameter, δ, we assign this set of faces as: where BC ( f ) is the barycenter of face f .We select the value of δ from literature providing the expected distance in the joint cavity [59] .The accuracy of the surface estimation highly depends on this value.We denote this subset of faces as F D C .We further trim the outer boundary and discard the faces with two boundary edges to provide additional robustness to the initial guess.The trimming helps ensure the primary interface does not cross the natural ridges on the bones.Secondary interface definition: In the next step, we define the top surface of the cartilages and denote the subset as F E C .In single-piece cartilage joints this is a subset of the secondary bone, while in the double-piece cartilages, we extrude the primary interface half-way towards the secondary bone and assign it as the top surface.Closed surface: we interpolate between the boundary of the primary and secondary interfaces to create a single closed surface (V C , F C ) .

Single-piece cartilage joints
Fig. 7 illustrates the cartilage generation steps for the sacroiliac joint and pubic symphysis.In the sacroiliac joint, we refer to the sacrum bone as the primary bone and the pelvic bone as the secondary bone.In the pubic symphysis, the left pelvis is considered as the primary bone and the right pelvis is the secondary bone.
Primary and Secondary interface estimation The geometry of the bones in the sacroiliac joint and pubic symphysis have welldefined ridges nearby the joint space.Thanks to this structure, these bones are sensitive to the distance filtering parameter and allow an accurate estimation of the primary and secondary inter- ) into a single mesh.

Double-piece cartilage joints
Fig. 8 illustrates the cartilage generation procedure in the hip joint cavity.We generate one cartilage attached to the acetabulum bone and one attached to the femoral head bone.The acetabular cartilage is in line with the lunate surface of the acetabulum, and the femoral cartilage covers most of the femoral head, excluding the fovea pit [17] .These two cartilages follow the same curvature and thin out as they get closer to their borders.Focusing on the femoral cartilage, we describe the femur as the primary bone (V P , F P ) and the pelvic as the secondary bone (V S , F S ) .On the acetabular side, the primary and secondary bones are swapped, assuming the pelvic as the primary and the proximal femur as the secondary bones.
Primary interface estimation The lunate surface forms a plateau in the acetabulum bone; Thanks to this structure, the pelvic side is sensitive to the distance filtering parameter ( δ), allowing an accurate estimation of the primary interface.The initial estimation on the femoral side does not yet comprehensively cover the femoral head.We observe the femoral cartilage border as a change in the curvature between the femoral head and the femoral neck.Thus, we further apply a curvature-based region filling approach to grow the initial guess to the correct portion on the femoral head.The details of this step can be found in Moshfeghifar et al. [15] .Fig. 8 illustrates the initial and the final primary interfaces on the femur in yellow and blue color, respectively.
Secondary interface definition We assign a thickness profile to the primary interfaces to create the top surface of the cartilage.We make a copy of F D C and denote it as the extruding surface (F E C ) .
Initially, we select a subset of F E C and extrude it towards the secondary bone as: Where n is the unit outward normal direction of vertex v in F D C .The number 1 / 2 describes the extrusion height which is equal to the mid-distance between the two bones.This guarantees congruency between the cartilage-cartilage interface.To ensure an excellent blending profile in the corners, we extrude the border meshes based on a sin function.These two surfaces are shown in green in Fig. 8 .
We connect the two extruded subsets via a smooth blend.We apply a biharmonic weighting scheme to compute a blended extrusion on the remaining of F E C which we did not initially select for extrusion.This subset is shown in pink color in Fig. 8 .A detailed explanation of this blending method is available at [61] .
Closed surface We invert the triangles in the primary interface and combine it with the secondary interface into a single mesh.Observe the closed surfaces in Fig. 9 .As explained in Section 2.3.1 , we provide two versions of hip joint models for each subject: with and without a gap between the articular cartilages.This can be done by controlling the extrusion height when defining the top surface.We found that reducing the extrusion height from the midpoint 0.5 to 0.45 leaves a small gap between the two articular cartilages.In this version, the cartilages will come into contact after applying a load in the FE simulation.

Cartilage anatomical measurement
The average cartilage thickness and the bone coverage area are measured for the paired hip joints, the paired sacroiliac joints, and the pubic symphysis.The bone coverage area for the single-piece cartilages is measured for both the primary and secondary bones.

Volume mesh generation
We employ the same meshing tool introduced in Section 3.4 to generate volume meshes for three reasons: 1) we want to have control over mesh size in our models and have the option to generate coarser meshes than the input mesh size; 2) Even though all the shared interfaces in our surface mesh models are conforming by design, we still need to ensure these properties are preserved after the volume mesh generation; 3) We want high geometrical accuracy after volume mesh generation.
Using fTetWild, we create volume mesh for all the sub-domains simultaneously rather than building our discretized model one-byone [37] .We apply a union operation on the input surface meshes and calibrate the ideal edge length (l) parameter in fTetwild to control the element size [37] .We select a small envelope ( ) value to preserve the anatomical details after volume mesh generation.The optimal mesh resolution and meshing parameters for subject m1 are obtained by performing a mesh convergence study on seven resolutions [1] .The meshing parameters for each subject are then calibrated to produce similar mesh properties as the optimal mesh settings.The meshing parameters and the number of elements are provided in Section 4.2 for all the models.fTetWild constructs volume meshes inside and outside the model, filling a bounding box around the model.Fig. 10 shows the raw output of fTetWild.These tetrahedrons still have no inside/out classification.We apply a post-processing step to extract the interior volume of each domain and filter out the elements that do not belong to any of the domains.This procedure is demonstrated in Fig. 10 .We calculate the signed distance of the centroid of each tetrahedron toward the surface of each domain.If the signed distance is negative, we consider that tetrahedron inside that particular domain.The tetrahedrons which have no negative signed distance towards any domains in our model, are filtered out before exporting the final mesh.

Mesh quality
We ensure that the final volume mesh elements are in good quality with no flat elements.The most convenient approach involves the dihedral angles, which is the angle between adjacent facets in a tetrahedral.We expect to see no dihedral angle close to where θ min and θ max denote the minimal and maximal dihedral angle in our models and θ is the dihedral angle of each element.
To give quantitative information about the number of poor elements we provide the percentage of dihedral angles falling out of this range.
Further, we test the robustness of our discretization results by checking other quality metrics suggested in literature.These include the volume-edge ratio (Q v l ) [63] , the radius ratio (Q rl ) [64] , and the volume-area ratio (Q rr ) [65] .

Finite element simulation
In this section, we demonstratethe performance of our models in different simulation setups and show that our models are compatible with different FE solvers.A pseudo-stance scenario under dynamic structural mechanics analysis is set up in the FEBio and PolyFEM solvers.As mentioned in Section 3.5.3, we provide 11 FE model with two hip joint versions: with and without a gap between the articular cartilages.Since FEBio requires an initial slight penetration between the contact surfaces, we use the model versions with no gap in the hip joints.PolyFEM, in contrast, requires an initial configuration free of penetrations; thus, we use the model versions with a small gap between the articular cartilage layers.
Note that the proposed choice of material model, load, and boundary conditions can be replaced or adjusted depending on different applications.We use the most simple material choices that are easily accessible from literature.Since the deformation of the bone is negligible under a pseudo-stance position, we consider the bone to be homogeneous, and the bones and cartilages are presented by an unconstrained Neo-Hookean material model.This material has non-linear stress-strain behavior but reduces to the classical linear elasticity model for small strains and small rotations [66] .The material properties for both experiments are from literature: the Elastic modulus (E) for the bones and cartilages is 17 GPa and 12 MPa, respectively, and the Poisson's ratio (ν) is equal to 0.30 and 0.45, respectively [6] .

Experiment A: simulation performance of the FE models using FEBio
This experiment tests the simulation performance of all 11 models.Each model consists of 12 deformable bones and cartilage The sacrum is fixed in all directions, and the distal end of the femurs are pushed towards the pelvic girdle, using a force-controlled simulation in FEBio.Right: Half of the body is removed, and the effect is applied as fixed boundary conditions in the sacroiliac and pubic symphysis attachment areas.The distal part of the femur is moved towards the pelvic girdle using a displacement-controlled simulation in PolyFEM.
parts.In Fig. 11 , we fix the pelvic girdle by restricting the sacrum's displacement and rotation in the x, y , and z -direction.The distal end of each femoral bone is tied to a rigid body.This rigid body has a force applied in the z -direction and is restricted in the other directions.The rigid force starts from zero and increases linearly to 430N on each femur.The articular interfaces in the hip joints are selected as the contact surfaces, and an augmented surface contact algorithm with friction-less tangential interaction is applied between them.

Experiment B: simulation performance of selective domains using PolyFEM
As explained in Section 2 , most of the hip joint populationbased studies use the geometry of only one side of the body and apply the effect of the other side as boundary conditions [19,29] .This experiment shows that one can use selective domains from our multi-body models for running simulations.We remove the geometry of half of the body and immobilize the pelvic bone by applying fixed boundary conditions to the pubic symphysis and the sacroiliac attachment surfaces in the x, y , and z -direction [67][68][69][70] .The distal ends of the two femurs are restricted in x and y directions, and a prescribed z -displacement equal to 1 mm is applied Fig. 12. Subject-specific finite element models of 11 male subjects with no diagnosed diseases related to the hip joint.Each model consists of the sacrum, the pelvic bones, the proximal femurs, the sacroiliac joints, the pubic symphysis, and the hip joints.The bones are directly derived from the CT scans, and the cartilages are generated subject-specifically using the bone geometries.Comparing the overall shape and the anatomical measurements among these 11 models indicates a considerable inter-subject variability in our study group.Besides, the difference between the left and right sides of the body shows bilateral variation and asymmetry in the hip joint area.directly to this area.As PolyFEM uses continuous collision detection to dynamically detect the contact surfaces, there is no need to manually specify the contact surfaces beforehand.Moreover, we benefit from the adaptive p-refinement option in PolyFEM and assign Tet10 elements to the cartilage layers to increase the simulations' accuracy.

Results and discussion
In this section, we present and discuss our results, starting from the accuracy of the anatomical approximations.Next, we show the mesh quality measurement results and describe our findings from the FE simulations.

Anatomical structures
The FE models of the 11 subjects are illustrated in Fig. 12 .The anatomical measurements for the bones and cartilages are measured in Tables 3 and 4 , respectively.The measurements are used to characterize each subject's size and shape and quantify our subjects' geometrical variation.
These measurements agree with the range of values reported in the literature, showing proper approximations of both bony and cartilage structures.For example, our models' mean femoral and acetabular cartilage thickness is equal to 1.56 mm and 1.38 mm, respectively.These numbers fall in the same range of 1.15 mm to 1.60 mm for the acetabular cartilage, and 1.18 mm to 1.78 mm Fig. 13.The comparison between the spherical and detailed design of the femoral head; We fit a sphere to the femoral head, and the sphere's diameter is defined as the simplified femoral head radius (SR).We define the distance between the sphere's center (HJC) and the femoral head boundary as the actual femoral head radius (AR).The variation of AR in subject m1 is shown as a probability density, and SR is shown with a red line.The red line and the beige distribution indicate the amount of subjectspecific data lost when simplifying the femoral head with a sphere.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3
The bone anatomical measurements; These measurements agree with the range of values reported in literature [53,55,74] , showing proper approximations of the bony structures.Comparing the anatomical measurements among the 11 models indicates a considerable inter-subject variability in our study group.This differences em phasize the im portance of subject-specific modeling and the role of population-based studies.The difference between the left and right side shows the bilateral variation and asymmetry in this area.

Pelvic girdle
Femur bones

Table 4
The joint anatomical measurements; These measurements agree with the range of values reported in literature [18,67,71,72,75] , showing proper approximations of the cartilage structures.Comparing the anatomical measurements among the 11 models indicates a considerable inter-subject variability in our study group.These differences emphasize the importance of subject-specific modeling and the role of population-based studies.The difference between the left and right sides shows the bilateral variation and asymmetry in this area.Comparing the overall shape and the anatomical measurements among these 11 models indicates a considerable inter-subject variability in our study group.These differences emphasize the importance of subject-specific modeling and geometrical variation in FE studies.Further, the difference between the left and right measurements shows asymmetry in this area which has been idealized in some studies [4] .In Fig. 13 , we compare a spherical head shape with a subject-specific femoral head.The difference between the radius of these two demonstrates how the geometry of the femoral head deviates from a sphere, and the amount of subject-specific data one loses when simplifying the femoral bone.

Model
Our modeling workflow is implemented in Python, mainly using the Libigl library [76] .To quantify the speed of the cartilage generation for each subject, we time all the steps described in Section 3.5 to generate the two hip joints, the two sacroiliac joints, and the pubic symphysis.The geometry processing is done using a MacBook Pro-2018 with a 2.7GHz quad-core Intel i7.We observe that, on average, the cartilage reconstruction takes 4 minutes for each subject.While this time is significantly short, generating good cartilages can take longer as we need to calibrate the free parameters in the cartilage generation pipeline [15] .These parameters did not deviate significantly among our subjects, so minimal time was needed for parameter tuning.
It should be noted that our cartilage generation algorithm is sensitive to the geometry of the input bone mesh.Any bone abnormalities or segmentation errors affect the cartilages' final shape.For example, even though there was no reported information regarding diseased hip joints, we found that some subjects, such as m4 have a hip cam impingement disorder which is a bump close to the femoral head.As we generate the femoral cartilage based on the femur bone curvatures, this bone abnormality affects the shape of the cartilage boundary.

Mesh quality
Table 5 summarizes the volume mesh properties of each subject.The Hausdorff distance between the surface boundary before and after volume mesh generation is, on average, 0.22 mm.To better understand the Hausdorff distance size, we additionally present it in percentage terms.We divide this distance by the diagonal of the bounding box in the initial surface and denote it by % .The average percentage of 0.05% indicates highly-accurate geometries.
The mean dihedral angle of 70.29 • in Table 5 indicates highquality elements.On average, the percentage of the low-quality elements is 0.21%.This number is less than 1%, which is negligible compared to the total number of meshes in the model.The normalized quality histogram of the other metrics is presented in Fig. 14 .The histograms are tilted towards one, with an insignificant number of elements close to zero.The number of elements starts growing around 0.40, indicating high-quality elements based on these measures.

Finite element simulation
Fig. 16 and 15 illustrate the von Mises stress results of experiment A and experiment B , respectively.We use a MacBook Pro 2018 with a 2.7 GHz quad-core Intel i7 to run the simulations of experiment A. Each simulation takes around 8 minutes to converge.For experiment B, we use a workstation with an AMD Ryzen Threadripper PRO 3995WX CPU.We use a max of 16 threads, which leads to a runtime of around 130 minutes for each simulation.In both experiments, we experienced no convergence issues related to the discretized geometries.
Fig. 16 shows that even though all the subjects have the same FE simulation setup, the differences in the bone and cartilage geometries affect the stress distribution pattern and values across the models.The stress pattern in the hip joints shows gradual changes in the joint cartilages with no spurious stress peaks.The stress ranges between 0 MPa to 0.58 MPa in the hip joint cartilages, which is lower than the reported values in the literature [6] .The results reported in our work only serve as a quality test of the simulation properties of our models; Thus, to get closer results, we need a more advanced simulation setting to model a real-stance position, such as proper material properties and boundary conditions.

Table 5
The volume mesh properties of our models; The fTetWild parameters allow preserving the anatomical details and controlling the element size.We use the dihedral angle θ to quantify the quality of the elements.To better understand the amount of low-quality elements in each model, we provide the percentage of dihedral angles falling out of the accepted range.The Hausdorff distance (HD) shows how the geometry is changed after the volume mesh generation.We also present it in percentage terms to better understand the Hausdorff distance size; We divide this distance by the diagonal of the bounding box in the initial surface and denote it by % .The average of the results is shown in Bold , indicating high-quality discretization with accurate geometries.: envelop of size epsilon; l: ideal edge length for the pelvic girdle and the legs; Mean: the mean dihedral angle; SD: standard daviation; Q θ : the dihedral angle mesh quality metric; % θ < and % θ > : percentage of elements which the dihedral angle falls below or above the range of good-quality tetrahedral elements, respectively; HD: Hausdorff Distance; % : the percentage of the HD with respect to the bounding box.

Conclusion and future work
Our work aimed to provide multiple subject-specific finite element models of the hip joint for simulation purposes.We provide verified bilateral models covering the bones and cartilages in the hip joint area.Our models allow scientists to account for the intersubject and bilateral variability when studying the hip joint behavior.The FE analysis experiments in this work qualify the usability of our models in FE simulations.We have made different versions of each subject to ensure the models are compatible with different FE method solvers and scenarios.We suggest that it is not suitable to only use a single subject in the FE analyses.
Our modeling pipeline and geometry processing codes are available on GitHub, located at [1] .The codes cover the input bone remeshing process, the cartilage generation process, the multi-body volume generation process, the mesh quality assessment codes, and a sample FEBio model file generator.Even though the FE models are the main output in our work, one can use each section's initial input and output separately for other purposes.For example, the CT scan segmentation label maps can be used as training data to develop automated bone segmentation models.
Additionally, the pipeline is compatible with adding more data or replacing each section.For example, additional bone segmentation can be added to the pipeline's input to develop more FE models.One can also directly segment the cartilages from image modalities, replace them with the generated cartilages in our method, and still generate high-quality FE models.The concept behind our pipeline is not limited to the hip joint only, and similar concepts can be applied to the jaw dataset or dental applications [77] .
The proposed cartilage generation method in our modeling workflow is fast and generates proper cartilages.However, we must note that our cartilage generation method is sensitive to the input bone geometries, and any abnormalities or segmentation errors may lead to less accurate results.The surface and volume discretizations have been proved to have high-quality elements using several mesh quality metrics.A feature of these models is that we have conforming interfaces in our multi-body models.Thus, one does not need to manually define fixed contact definitions in the bone-cartilage interfaces when running simulations.Further, the high congruency level between the articular cartilages in the hip joint allows a smooth transition of force in the hip joint, leaving no stress peaks due to improper cartilage geometries.
It is important to note that all the subjects in our work are in the supine position during the image acquisition.They are lying down on the scanning tray, and the relative location of the femur and pelvic bones are affected by this position.One solution is to relocate the bone geometries to a neutral position before setting up any custom simulations setup.We leave this for future work.Future work could extend to using these FE models to estimate the stress distribution in a complete gait cycle.Our models currently lack other soft tissues.Future work can add more structures in the hip joint area, such as the hip joint labrum, capsule, and ligaments.

Fig. 3 .
Fig.3.The comparison between simultaneous and one-by-one meshing approaches in the sacroiliac joint; The mesh in our work (a) is generated using the multi-body meshing approach.This method guarantees conforming interfaces with the hip and sacrum bones.Figure (b) shows the result of the one-by-one mesh generation approach, resulting in overlapping (c) or separated (d) interfaces.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig.3.The comparison between simultaneous and one-by-one meshing approaches in the sacroiliac joint; The mesh in our work (a) is generated using the multi-body meshing approach.This method guarantees conforming interfaces with the hip and sacrum bones.Figure (b) shows the result of the one-by-one mesh generation approach, resulting in overlapping (c) or separated (d) interfaces.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8.The acetabular and femoral cartilage generation in the hip joint; The blue surface indicates the primary interface estimation on the bones F DC .This subset is selected based on the distance to the opposite bone.The initial primary estimation in the femur bone (yellow) does not comprehensively cover the femoral head.Thus, we apply a curvature-based region filling to grow the yellow cartilage to the correct portion on the femoral head (blue).We make a copy of the primary interface to create the top surface of the cartilage F E C .We initially extrude a part of this surface to a maximum of half the distance to the second bone (green) and then connect the extruded subsets via a smooth blend (pink).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9.The reconstructed hip joint cartilages on the acetabular (orange) and femoral (yellow) sides; The acetabular cartilage is in-line with the lunate surface of the acetabulum, and the femoral cartilage covers most of the femoral head, excluding the fovea pit [17] .The cartilages have congruent interfaces and thin out as they get closer to their borders.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Overview of the fTetWild volume mesh generation; The raw output is a tetrahedral mesh of the bounding box.This volume yet has no inside/outside classification.We apply a post-processing step to filter out the exterior volume and keep the interior elements in each sub-domain.The volume mesh is extracted using the signed distance measure w.r.t the input surface.Negative and positive distances indicate interior and exterior elements, respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)0 • or 180 • [62] .We employ a method proposed by Freitag et al. to define the lower and upper bound of the dihedral angle as:

Fig. 11 .
Fig. 11.Overview of the pseudo-stance simulation scenario in experiment A (left) and experiment B (right).Left:The sacrum is fixed in all directions, and the distal end of the femurs are pushed towards the pelvic girdle, using a force-controlled simulation in FEBio.Right: Half of the body is removed, and the effect is applied as fixed boundary conditions in the sacroiliac and pubic symphysis attachment areas.The distal part of the femur is moved towards the pelvic girdle using a displacement-controlled simulation in PolyFEM.
IHS: inter-hip separation in mm; L-, R-IH: height of the left and right ilium in mm; L-, R-IW: width of the left and right ilium in mm; PW: the width of the pelvis (PW) in mm; L-, R-FHSR: the left and right femoral head sphere radius in mm; L-, R-FHAR: the left and right femoral head actual radius in mm; L-, R-VFL: the left and right visible femur length in mm; L-, R-NSA: the left and right femoral neck-shaft angle in degrees.

Fig. 14 .Fig. 15 .
Fig. 14.Normalized mesh quality histogram for m1 ; The quality metrics include the volume-edge ratio(Q v l ) in blue[63] , the radius ratio (Q rl ) in green[64] , and the volume-area ratio (Q rr ) in red[65] .The closer the distribution to one, the higher the mesh quality we have.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .
Fig. 16.The FE analysis results of the hip joint in experiment A , using FEBio; The color bar shows the normalized von Mises stress in the hip joint cartilages.Even though all the subjects have the same FE simulation setup, we observe that the differences in the bone and cartilage geometries affect the stress distribution pattern and values across the subjects.The stress distribution pattern in these joints shows no spurious stress peaks and gradual changes in the joint cartilages.The overall stress values range from 0 MPa to 0.58 MPa.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) [33]subject-specific geometry from images; CAD: computer aided designed models; SS+mirror: half of the model is derived from images and it is flipped to generate the other side; SS+AM: the diseased models are derived from images and the healthy model is generated by deforming the Visible Human Project Model[33]; PB: pelvic bone, PF: proximal femur, SB: sacrum bone, HJ: hip joint, HJ-L: hip joint labrum, HJ-C: hip joint capsule, SIJ: sacroiliac joint, PS: pubic symphysis; I: image availability, FEM: availability of the finite element models.