Computer Methods and Programs in Biomedicine

Background and objective: Pre-operative surgical planning using computer simulation is increasingly stan- dard practice before Total Hip Arthroplasty (THA), in order to determine the optimal implant positions, and thereby minimise post-operative complications such as dislocation, wear and leg length discrepancy. One of the limitations of current methods, however, is the lack of information on the subject-speciﬁc reference range of motion (ROM) that could be used as targets for surgical planning. Only a limited num- ber of hip motions are considered, which are neither subject-speciﬁc, nor representative of all the hip motions associated with all the activities of daily livings (ADLs). In this paper, therefore, a method was developed to calculate subject-speciﬁc representative bony range of motion (B-ROM) that would cover all the possible joint motions and presented in terms of pure joint motions. Methods: Only 3D bone geometries of femur and pelvis, constructed from personalised CT scan, were used as inputs for healthy hip joint whereas implant geometries and their positions on native bone ge- ometries were required for planned treatment side or replaced side. Hip joint motion simulation was carried out using six different Tait-Bryan intrinsic rotation sequences of three pure joint motions - ﬂexion- extension, abduction-adduction and internal-external rotation, and B-ROM was then identiﬁed for any of these six different sequences which caused earliest feasible impingement. The B-ROM could be used as a list of ROM data points or visualised as multiple 2D surface plots or a 3D envelop. Using the developed method, the B-ROM of a contralateral healthy hip joint of a patient can be used to deﬁne the subject-speciﬁc target ROM values to inform the surgical planning of the arthritic hip side so that the patient’s natural ROM could be restored as closely as possible by the planned implant placements. This was demonstrated with a clinical veriﬁcation study using ‘non-dislocating’ and ‘dislocating’ THA patients. Results: The results supported the study hypothesis that the percentage of intersected volume of the healthy and replaced side B-ROM was higher for the ‘Non-Dislocator’ patient (95%) compared to ‘Dislocator’ (78%). Also, the results showed that the only one sequence (ﬁrst ﬂexion-extension, then abduction- adduction and ﬁnally internal-external rotation) was not adequate to identify all the possible limiting B-ROM, and therefore, all the six rotation sequences should be considered. Conclusions: The method encompasses every potential ADL, and as a result, more comprehensive surgical planning is possible, as the implant positions can be optimised in order to maximise impingement-free ROM, and consequently minimise clinical complications.


Introduction
Total Hip Arthroplasty (THA) is a successful technique to restore lost or restricted mobility of patients suffering from osteoarthritis (OA), rheumatoid arthritis (RA), and acute trauma [1][2][3] . Recurrent dislocation after THA is the primary cause of revision surgery [ 4 , 5 ]. Although the overall incidence of dislocation has decreased over the past two decades [ 6 , 7 ], a significant number of THA patients continue to experience repeated episodes of dislocation. Such multiple dislocations, defined as two or more occurrences [8] , occur in over 60% of patients at a minimum follow-up of one year after the first dislocation [ 7 , 8 ] and over half of these patients require revision surgery [ 7 , 8 ]. Impingement is a major cause of inadequate range of motion (ROM), and subsequently, post-THA dislocation [ 9 , 10 ]. This impingement can be broadly classified into prosthetic impingement (PI) and bony impingement (BI) [11] . PI occurs when the prosthetic femoral neck comes in contact with the rim of the liner/cup. BI or bone-to-bone impingement (BTBI) occurs between the osseous femur and the osseous pelvis. Also, there is implant-to-bone impingement (ITBI) which occurs when the prosthetic femoral stem comes in contact with pelvis, or the bony femur comes into contact with the rim of the liner/cup. The risk of these impingements and recurrent dislocations is related to a number of factors, such as implant design, orientation and positioning within the bone, related surgical procedures, bony structures around the hip and other patient related aspects such as gender, age, and history of previous hip surgery [ 8 , 12-15 ].
Surgeons and health care providers have therefore employed various mitigating strategies, such as careful pre-operative planning, thorough intra-operative assessments, and meticulous postoperative care [7] . Pre-operative surgical planning is now regarded as standard practice before THA, with some systems incorporating a computer simulation to identify optimal implant positions to avoid post-THA PI [16][17][18][19][20] and BTBI and ITBI [21][22][23] . In such state-of-the-art planning, a limited number of hip joint motions or activities of daily living (ADLs) are simulated using indicative ROM values collected from the literature (such as flexion of 120 0 , extension of 20 0 etc) to check for impingement (PI or BTBI or ITBI). There are, however, two major limitations associated with this approach. (a) Subject-specific reference values of ROM for a particular activity (or hip motion) are unknown, requiring therefore use of an arbitrary reference value [ 16 , 17 , 24 ]. For example, a reference value of hip flexion of 120 0 may be used to detect impingement during surgical planning, but the actual value required could be 130 0 for a particular patient. Consequently, any implant positioning based on an analysis of 120 0 of flexion might not be suitable for patients that require 130 0 . (b) Only a limited number of hip joint motions can be simulated, which do not represent all potential ADLs. Therefore, in order to get best outcome from a THA planning tool, it is important to know (a) how many potential hip motions should be checked for impingement analysis, and (b) what the subjectspecific limiting values for those hip motions should be. Turley, et al. [25] developed a benchmark ROM for hip joint movement by tracing the location of the knee centre which was considered as a position vector from hip joint centre. As a result, the effect of internal and external hip rotation was not directly included in the method. Also, the method used generalised limiting values of pure joint motions from the literature to construct the benchmark ROM. One approach would be the use of wearable motion sensors, to measure subject-specific ROM for all ADLs. However, this would require significant experimental set-up and long execution time and is unlikely to be clinically practical. Also, there would be additional challenges such as registration of motion sensor measurement to the CT data as the planning is generally performed using the pelvis and femur geometries, constructed from CT images at supine position. Furthermore, patients generally come to the hospital when they already have hip problems. Hence, measuring the hip joint motion in that stage would provide ROM for an arthritic or diseased hip instead of healthy hip ROM, and therefore, cannot be used as reference ROM.
The purpose of this study was, therefore, two-fold. Firstly, a method was developed to calculate a subject-specific bony range of motion (B-ROM) boundary that would cover all the possible joint motions in 3D and presented in terms of pure joint motions i.e., amount of flexion-extension, abduction-adduction, internalexternal rotation and their combinations. Secondly, the method was used to determine the B-ROM for both the contralateral hip and the planned treatment side of a patient, where the contralateral hip is healthy. The B-ROM of this healthy side could be used as subject-specific reference values for surgical planning, and thereby enable planning that more reliably maximises the impingementfree ROM. This would therefore provide patient-specific and comprehensive ROM values for the planning of implant positions, in contrast to existing methods that use limited hip motions without subject-specific reference values.
The paper is organised as follows. The conceptual novelty of constructing a subject-specific B-ROM was described in the first part of the Materials and Methods section followed by the description of the developed method. Thereafter, a clinical verification study was included to highlight the potential applications of the method. The rest of the paper described the results followed by discussion.

Conceptual novelty: bony range of motions (B-ROM)
In the state-of-the-art THA pre-planning, a small number of hip joint motions are used along with the generalised reference values collected from the literature [ 16 , 17 ]. This has two major limitations. Firstly, it is not a well representation of all the possible ADLs that the patient would be able to perform. Secondly, the reference values of hip motions are not subject-specific. In this work, a novel method was developed to calculate subject-specific B-ROM that was a well representation of all the possible hip joint motions associated with all the ADLs for a particular subject. ADLs are functional hip motions that are constrained by bony structures, soft tissues and ligaments. The developed B-ROM is a theoretical representation of all the hip motions permitted by the bony geometry without directly considering the restrictions imposed by the soft tissue and ligaments. Therefore, it was a reasonable assumption to consider that the B-ROM would cover the entire spectrum of ADLs for a particular subject as the inclusion of soft tissues and ligaments could only reduce the B-ROM further. The B-ROM was represented through pure joint motions (i.e., flexion-extension, abduction-adduction and internal-external rotation) and their combinations. It was presented as a list of various hip motions or visualised using 3D envelope as well as multiple 2D surface plots.
In the ideal situation, the pre-arthritic healthy B-ROM of the planned treatment side would be known, and this could be used to assess whether patient's natural B-ROM would be optimally restored by the planned THA implants. However, directly measuring pre-arthritic healthy B-ROM of the treatment side is not generally possible because patients generally come to the hospital when they already have hip arthritis, and their normal anatomy has already changed. In such cases, when a patient's contralateral hip is healthy, the B-ROM of the healthy side could be used as a subjectspecific benchmark. It was based on the assumption that the B-ROM for both healthy hip joints would be similar, if not exactly same. Therefore, B-ROM of the contralateral healthy side could be used to determine the best possible implant positioning on native bone geometries to best restore the patient's native and natu-ral impingement-free hip motions. In such cases, when a patient's contralateral hip is healthy, the B-ROM of the healthy side could be used as a subject-specific benchmark. Although it is recognized that the contralateral healthy BROM may not be exactly same as the ideal pre-arthritic healthy B-ROM of the planned treatment side, it is likely to be a reasonable approximation.
In this paper, a healthy hip was defined as a radiologically normal hip, with no evidence of osteoarthritis (grade 0, as defined by Kellgren and Lawrence [26] , and later accepted by the World Health Organisation (WHO) [ 27 , 28 ]).

Development of the method
The construction of subject-specific B-ROM was implemented through five steps as explained in the following sections.
Inputs: 3D geometries of pelvis and femur were required along with pelvis and femur landmarks which were obtained using the following steps: (a) CT images of a patient requiring THA surgery, (b) construction of subject-specific bone geometries from CT images, (c) identification of bony landmarks. The CT images used in this study were available on Corin's image library. They were scanned for pure clinical reasons before hip replacement according to the 'OPS Dynamic Hip Analysis' low dose CT protocol (Optimized Ortho, Sydney, Australia) [29] that was defined based on the work of Huppertz, et al. [30] . The CT slice thickness was 2.0mm with in-plane resolution of 0.4mm. The CT scans, used in this study, were fully anonymised in accordance with the institution's ethical guidelines, and research ethics approval by Bellberry Human Research Ethics Committee (BHREC) (2012-03-710). Construction of bony geometries from CT images and identification of landmarks was performed using Simpleware TM ScanIP software (Synopsys, Inc., Mountain View, USA). The DICOM CT images were imported in ScanIP followed by cubic/isotropic resampling and cropping the region of interest. The pelvis and femur were then semi automatically segmented using lower and upper greyscale threshold values, and subsequently, morphological close operation was performed to close the holes within the segmented masks. Finally, the 3D geometry was constructed from the mask, and subsequently, exported as STL files for the simulation. Four landmarks were identified for pelvis as follows ( Fig. 1 a,  In the case of the planned treatment side, two additional inputs were required: (a) CAD model of the planned implants (specifically, the cup, liner and stem) in STL format, and (b) planned implant positioning on native bone geometries (e.g., cup/liner inclination and anteversion angle, stem offset, surgical neck cut of femur etc.). With these inputs, the implant geometries could be positioned onto the native bone geometry according to the surgical plan. Implant selection, positioning and identification of all bony landmarks was carried out by experienced engineers in Corin Ltd along with the help of surgeon feedback. This is a regular process for the Corin OPS THA pre-planning service [29] .
Step 1: Identification of hip joint centre The centre of a sphere, which was best fitted on the native femoral head, was defined as the Hip Joint Centre (HJC). This entire operation was implemented in Matlab (Version 2018b, The MathWorks, Inc., Natick, Massachusetts, United States). Firstly, the data cursor mode was turned on using 'datacursormode' function, and subsequently, multiple points on the spherical articular surface of femoral head were selected manually ( Fig. 2 a). It was observed that more than 20 points was sufficient to get a reliable HJC. Thereafter, the selected points were used to create a best fit sphere, and subsequently, the centre of this sphere was defined as HJC ( Fig. 2 b). Similar operation could be performed in any other reverse engineering software such as Geomagic (3D System, Morrisville, North Carolina), and thereafter, the HJC could be exported to use in the simulation. For planned treatment side, centre of the head implant was considered as HJC.
Step 2: Construction of pelvic and femur coordinate system and alignment to the global coordinate system: The Pelvic Coordinate System (PCS) was constructed using the four pelvic landmarks ( Fig. 1 a). First, the midpoint (PTUB Mid ) of PTUB Right and PTUB Left was calculated. Anterior Pelvic Plane (APP) was then defined using three landmark points -ASIS Right, ASIS Left, and PTUB Mid ( Fig. 1 a). The medial-lateral x-axis ( X P ) was defined by a line connecting from ASIS Left to the ASIS Right ( Eq. (1 )). The normal to the APP represented the anterior-posterior y-axis ( Y P ), and its positive direction was from posterior to anterior side ( Eq. (2 )). Finally, the superior-inferior z-axis ( Z P ) was orthogonal to the other two axes ( Eq. (3 )). The PCS is represented by the rotation matrix R W P in Eq. (4 ) where the superscripts represented i, j , and k component of each unit vector X P , Y P , and Z P . The origin of the PCS was the HJC. (2) Similar to the PCS, the origin of the Femur Coordinate System (FCS) ( Fig. 1 b) was the HJC. The knee centre (KC) was defined by the mid-point of EPI Lateral and EPI Madial . The superior-inferior z-axis ( Z F ) was a line running in the positive direction from the KC to the HJC ( Eq. (5 )). The anterior-posterior y-axis ( Y F ) was constructed perpendicular to a plane which was defined by three points -HJC, EPI Lateral and EPI Madial ( Eq. (6 )). The medial-lateral x-axis ( X F ) was orthogonal to the other two axes ( Eq. (7 )). The FCS is therefore represented by the rotation matrix R W F in Eq. (8 ) where the superscript represent i, j , and k component of each unit vector X F , Y F , and Z F .  The PCS and FCS were not necessarily aligned together when they were constructed using the input landmarks ( Fig. 1 c). Therefore, the pelvis and femur geometries were first translated in such a way that the HJC coincided with the origin of the World Coordinate System (WCS). Thereafter, the pelvis geometry was rotated using the rotation matrix ( R W P ) T to align PCS with WCS so that X P and Y P coincided with X W and Y W respectively ( Fig. 1 d). Similarly, FCS and WCS were aligned by rotating the femur geometry using the rotation matrix ( R W F ) T so that X F and Y F coincided with X W and Y W respectively ( Fig. 1 d). This position was termed as neutral position where PCS, FCS and WCS were coincident ( Fig. 1 d).
Step 3: Removal of femur head: The radius of the best fit sphere, which was created to determine HJC, was increased with very small step size (0.1mm) until the expanded sphere covered the entire femoral head, and consequently, the expanded sphere could intersect the femur in exactly two regions -(a) femoral head and (b) rest of the femur geometry. Due to this criteria and use of very small step size (0.1mm), the variation in expanded sphere dimension was negligible. Subsequently, the entire pelvis and femur with resected head were used as two distinct geometries for impingement identification ( Fig. 3 a) for contralateral healthy side. Original best sphere could not be used for the resection as it would not be possible to resect the entire femur geometry into exactly two regions. It was observed that the average differences in diameter between best fit and enlarged sphere was 1.5 to 3 mm. The purpose of the simulation was to identify any potential collision/impingent between femur and pelvis of a healthy hip joint so that this could be used as a personalised reference values of ROM during pre-planning. As the spherical articular surface of femoral head could not impinge with acetabulum during hip joint motion of a healthy hip, this part was removed from the analysis ( Fig. 3 a). If surgical neck cut was performed for a healthy hip, more bone would have been removed specially from the neck area. As a result, the analysis would potentially suggest a larger range of movement than would actually be possible when the femoral neck bone is present.
For same reason, head implant was not considered while calculating limiting ROM for planned treatment side. The pelvis, cup and liner geometries were assembled together to be considered as a single geometry whereas stem with surgically neck cut femur were assembled together in a single geometry for impingement identification.
Step 4: Simulation to identify B-ROM and visualisation: The pure joint motion was defined as follows: (a) Rx represented the rotation with respect to global x-axis ( X W ) to generate flexion-extension movement of angle α, (b) Ry depicted the rotation around global y-axis ( Y W ) to generate abduction-adduction of angle β, and (c) Rz created the rotation around global z-axis ( Z W ) to generate internal-external rotation of angle γ .
Any combined hip joint motion could be represented as a multiplication of these matrices in a given sequence. For example, xyz sequence for intrinsic rotation is represented as Rx * Ry * Rz. As matrix multiplication is non-commutative, it leads to six different sequences. As a result, the femur could be in six different positions due to six different rotation sequences for a given set of α, β, and γ . In this paper, 'Tait-Bryan' intrinsic rotation definition was used to define six rotation sequences as follows: FE-AA-IRER (Rx * Ry * Rz), FE-IRER-AA (Rx * Rz * Ry), AA-FE-IRER (Ry * Rx * Rz), AA-IRER-FE (Ry * Rz * Rx), IRER-FE-AA (Rz * Rx * Ry), and IRER-AA-FE (Rz * Ry * Rx) where FE, AA and IRER represented Flexion-Extension, Abduction-Adduction and Internal Rotation-External Rotation respectively.
The simulation was developed in such a way that the following combinations were simulated for a particular IR or ER: (a) Abduction-Flexion (Quadrant I), (b) Adduction-Flexion (Quadrant II), (c) Adduction-Extension (Quadrant III), and (d) Abduction-Extension (Quadrant IV) ( Fig. 4 ). Therefore, 2D B-ROM would represent the limiting hip joint motion for a particular IR or ER through various combinations of Flexion-Extension and Abduction-Adduction. When all the ER and IR were considered, it would produce a stack of 2D B-ROM that was termed as 3D B-ROM. This is represented through a pseudo code in Fig. 4 .
The outermost loop was termed as 'Loop γ ' that was used to generate the amount of IR/ER with variable γ and step size γ . Within the 'Loop γ ', there were four main blocks that were used to generate data for Quadrant I to IV, and detailed pseudo code of Quadrant I is only shown in Fig. 4 . Each block consisted of two loops, termed as 'Loop α' with variable α and step size α and 'Loop β' with variable β and step size β. The step size γ was used as 5 0 whereas step size α and β were used  of α, β, and γ using six rotation sequences, and impingement was checked for each femoral position. In this work, a Matlab based function 'fastMesh2Mesh', developed by Seers [32] based on Möller-Trumbore 'ray triangle intersection' (RTI) algorithm [33] , was used for collision detection. Alternatively, a collision check function within Matlab using the MEX API interface was used. The function was developed based on the 'collision detection' algorithm available in the Proximity Query Package (PQP) library. The output (intersection points) from the function were compared with the intersection points, computed by Mimics 3-matic software when same STL geometries were used [21] . It was observed that the output points were similar, and therefore, the function was used thereafter for the ease of automation. When there were impingements, the corresponding rotation sequences were stored as 'Se-qIds' along with the limiting values of α, β and γ as α Lim , β Lim , and γ Lim respectively ( Fig. 4 ) before exiting from 'Loop β'. When there was a limiting value with α Lim = 0, the 'Loop α' stopped, and subsequently, it exited from one block and entered to the next block. This process continued until all the four blocks were executed to generate the data points for Quadrant I to IV for a given value of γ . The maximum and minimum values of α, β, and γ were used as ±180 0 , and therefore, no assumptions or literature values were used in selecting the maximum values of Flexion-Extension, Abduction-Adduction and IR-ER. If the B-ROM was de-veloped for contralateral healthy side, only bone-to-bone impingement was considered whereas for planned treatment side, implantto-implant, bone-to-bone and implant-to-bone impingements were considered to get the limiting ROM. Only bony geometries were considered for B-ROM calculation without considering soft-tissue and ligament constraint. Therefore, the calculated B-ROM would be larger compared to the functional ROM, and as a result, some of the impingement points would not be feasible in reality. Therefore, in order to identify the feasible impingement area, only the impingement points which were within a 'feasible sphere' with centre at HJC and radius from HJC to Lesser Trocheneter (LT) were considered. Consequently, the B-ROM that produced impingement within this 'feasible sphere' were finally considered. It was assumed that this process partially included the effect of soft tissue and ligament in constraining the ROM. For example, the B-ROM points were removed that caused impingement in the pubis area, and this was not possible in reality due to soft tissue and ligament constraints. However, it was not claimed that this would entirely include the effect of soft tissue and ligaments constraints on restricting hip ROM.
The analysis consisted of mainly two steps -(a) producing 'Inputs' that required manual operation such as construction of bone geometries from CT images and identification of bony landmarks; (b) identifying BROM that was fully automated procedure implemented in Matlab. The first step was generally completed within 30-40 mins by experienced engineer. The simulation (i.e., the second step) time was approximately 3-4 h for contralateral healthy side and 5-6 h for planned treatment side when it was executed in a computer with Intel(R) Xeon(R) CPU E3-1535M v6, 3.10 GHz processor, 64 GB RAM, 4 cores (8 logical processor) and Windows 10 operating system and Matlab 2018b version.

Output:
The immediate output from the model was the list of four set of values: α Lim , β Lim , γ Lim , SeqIds ( Fig. 4 ). The first three values represented the limiting values of flexion-extension, abductionadduction and internal external rotation respectively for a combined motion that generated an impingement. The fourth term, 'SeqIds' defined the combined motions by describing the particular 'Tait-Bryan' rotation sequences. This output was visualised using 2D and 3D plots of point clouds, and subsequently, 3D B-ROM envelop was generated, and saved as STL file for future use. An additional information was also stored -the actual impingement location on native pelvis and femur geometries, which was associated with the calculated B-ROM.

Case studies a) Case Study 1: All sequences vs XYZ clinical sequence
This case study was performed to justify that only one particular 'Tait-Bryan angle' rotation sequence was not enough and not well representative of all the hip joint motions associated with all the ADLs. Generally, FE-AA-IRER was the most commonly used rotation sequence to define femoral movement in the previous studies [ 18 , 25 , 34 ]. Therefore, a comparative analysis was performed based on the calculated B-ROM and related impingement area when all the six rotation sequences were used instead of using FE-AA-IRER sequence only for a subject with healthy left hip joint.

b) Case Study 2: Clinical study of a 'Dislocator' and a 'Non-Dislocator' patient
A clinical verification study was performed to demonstrate the potential use of the developed B-ROM. Two patients were selected based on the outcome of their THA after 2 years: (a) a 'Non-Dislocator' patient where there had been no postoperative episodes of hip dislocation; (b) a 'Dislocator' patient where there had been at least one clinical episode of dislocation of the THA. Requirements of the 3D models of the implants, positioned onto the native bone geometries, and bony landmarks were detailed in the 'Input' section of the method. Table 1 summarises the patient characteristics and related surgical information. The cup inclination and anteversion angle were within the safe zone defined by Lewinnek et al. [35] . The B-ROM of the contralateral healthy side of each patient was calculated. This contralateral healthy B-ROM was then used as target ROM to evaluate the B-ROM of their replaced side as discussed in the Section 2.1 . It was hypothesised that the intersected volume of the healthy and replaced B-ROM would be higher for the 'Non-Dislocator' patient compared to the 'Dislocator'. Higher intersected volume would indicate that the B-ROM of the replaced side more closely followed the contralateral native healthy B-ROM, and therefore, chances of impingement would be less. This would not imply that the impingement was the reason for dislocation for the 'Dislocator' patient, it would be however consistent with the association between impingement and dislocation observed in clinical practice.

Output from the developed method
The immediate output from the model was the list of four set of values: α Lim , β Lim , γ Lim , and SeqIds ( Fig. 4 ). This was used further to create various visual representations ( Fig. 5 ). 2D surface plot was created to depict various combinations of FE and AA for a particular IR or ER ( Fig. 5 a). The SeqIds was not shown explicitly. When multiple 2D surface plots were stacked together along the IRER axis (or Z-axis), 3D point cloud was created ( Fig. 5 b) that was subsequently converted into a 3D envelop and saved as STL file format ( Fig. 5 c). The 3D B-ROM, shown in Fig 5 b and d, was associated with the impingement area that were indicated by the green area on the pelvis geometry ( Fig. 5 d). The red impingement area was considered not to be feasible as it was outside of a 'feasible sphere'. The B-ROM associated with the red impingement areas were not included in 2D or 3D B-ROM plot ( Fig. 5 a-c).   using only one particular sequence such as FE-AA-IRER as performed in the previous studies [25] . Three 2D B-ROM plots for ER of 15 0 , no rotation (i.e., 0 0 ), and IR of 10 0 are shown in Fig. 6 a to highlight the differences in B-ROM boundaries due to use of six sequences and only one sequence. There were mainly two differences as follows. (a) The B-ROM boundary generated using all of the six 'Tait-Bryan angle' rotation sequences was much smaller (blue boundary) than the B-ROM boundary calculated using only FE-AA-IRER sequence (black boundary) ( Fig. 6 a). This meant that there were some combinations of FE-AA for a particular IR or ER that was not captured by FE-AA-IRER sequence alone. (b) The B-ROM for all the six sequences provided a closed boundary for Abd-Flex (Quadrant I) and Abd-Extn (Quadrant IV) combinations as compared to the B-ROM results for FE-AA-IRER sequence. Therefore, it appeared that there would not be any impingement where the boundary was not closed for FE-AA-IRER sequence (Quadrant I and IV). In reality, impingement might still occur as shown by the closed blue boundary specially in those quadrants. Fig. 6 b shows that the 3D B-ROM associated with FE-AA-IRER sequence is much larger than the proposed B-ROM in this study for all the IR and ER. Fig. 6 c and 6 d show the associated feasible impingement areas on pelvis and femur for using FE-AA-IRER sequence and all of the six sequences. It was observed that the impingement area on the acetabulum area was similar for both cases although the pelvic ischium area showed larger impingement area for considering all the six rotation sequences. On the other hand, the impingement area was larger in illum region of the pelvis for FE-AA-IRER sequence as the amount of abduction was more (Quadrant I in Fig. 6 b) that were discarded while considering all sequences. However, the impingement area on femur was much larger for all six rotation sequences in comparison to FE-AA-IRER sequence. It highlighted the importance of considering all the six sequences as it would assess all the various femoral position for impingement checking in comparison to using only FE-AA-IRER sequence.

Case study 2: clinical study of a 'dislocator' and 'non-dislocator' patient
The B-ROM for contralateral healthy side and replaced side were calculated for both 'Non-Dislocator' and 'Dislocator' patient. When contralateral healthy B-ROM was compared to the B-ROM of the replaced side through 3D B-ROM envelop, it was observed that the intersected volume was higher for 'Non-Dislocator' patient compared to the 'Dislocator' subject ( Fig. 7 a and 7 b). Higher intersected envelope demonstrated that the replaced B-ROM more closely followed the reference contralateral healthy B-ROM.
The percentage of intersected volume between the B-ROM of healthy and replaced side was measured by calculating the intersected area for each 2D B-ROM using 'polyshape' and 'intersect' functions in Matlab (Version 2018b, The MathWorks, Inc., Natick, Massachusetts, United States). Thereafter, summation of all the intersected area was divided by the summation of all the healthy B-ROM area to get the percentage of intersection. For 'Non-Dislocator' subject, the percentage of intersection was 95% whereas for 'Dislocator' subject, it was 78%. The result was perfectly aligned with the study hypothesis that the intersected volume between the B-ROM of healthy and replaced side would be higher for the 'Non-Dislocator' patient compared to the 'Dislocator'. Although it was not demonstrated that impingement was the reason for dislocation in this particular case, it was however consistent with the association between impingement and dislocation observed in clinical practice.

Discussion
In this paper, subject-specific B-ROM was developed using only 3D bony geometries. The novel aspects of the proposed method were as follows. (a) A method was developed to calculate theoret-ical B-ROM that would be a close representation of all the possible hip joint motions that a subject can perform during all possible ADLs. The B-ROM was represented though pure joint motions i.e., amount of flexion-extension, abduction-adduction, internalexternal rotation and their combinations, and subsequently, visualised using multiple 2D plots (each 2D plot was for a particular internal/external rotation with FE-AA combinations) or 3D envelope (for all combinations of FE, AA and IRER) or could be simply exported as list of data points for user specific use. (b) The B-ROM of a contralateral healthy hip could be used to define subject-specific reference values to determine best possible B-ROM of the planned treatment side. This subject-specific and comprehensive analysis would therefore ultimately provide greater confidence when it comes to pre-planning. (c) The developed B-ROM only used CT developed bone geometries (pelvis and femur) for healthy hip joint side. For a replaced/planned hip, additional data such as implant geometries along with their positions on native geometries are required. No additional data was required to develop subject-specific B-ROM. This method could therefore be implemented in any existing pre-planning tool.
There are six different sequences to define the 'combined motion' according to the 'Tait-Bryan' angle sequence with intrinsic rotation definition: FE-AA-IRER, FE-IRER-AA, AA-IRER-FE, AA-FE-IRER, IRER-FE-AA and IRER-AA-FE. Therefore, with a given set of values of angular movement ( α, β, and γ ), the femoral positions associated with the six different sequences of combined motions were simulated, and subsequently, impingement checking was performed to get limiting B-ROM. It was observed that only a particular rotation sequence (such as FE-AA-IRER) was not enough to describe all of the possible femoral motions. In the previous work [ 18 , 25 , 34 ], generally FE-AA-IRER was used to define the combined motion of femoral movement. However, as shown through Case study I, only FE-AA-IRER sequence would produce a B-ROM which might be larger compared to the B-ROM constructed using all the six sequences ( Fig. 6 a). Thus, any planning using the developed B-ROM would provide additional information regarding all of the limiting femoral positions which could not be captured using one sequence alone as performed in the work of Turley et al. [25] . In addition, the Turley et al. [25] developed a benchmark ROM for ADLs but it used generalised pure joint motion values from the literature. Moreover, it traced the knee centre which was considered as a position vector from hip joint centre, and therefore the effect of internal and external hip rotation was not directly included in the method. In the developed method, both of the limitations were mitigated by representing the amount of rotations instead of tracking any point on femur.
In majority of the cases, when resection of femur was performed with best fit sphere (instead of enlarged sphere), there were more than two regions or continues connections between femoral head and neck area after intersections, and therefore, this was not suitable for simulation. As a result, a direct comparison of the BROMs when best fit and enlarged spheres were used for resection was not always possible. However, a sensitivity study was performed with three set of sphere diameters -(a) Sphere A with ideal enlarged diameter identified using abovementioned criteria, (b) Sphere B that was produced by increasing the diameter of the Sphere A by 1mm, (c) Sphere C that was produced by reducing the Sphere A diameter by 1 mm (in this case, the additional regions created after intersection was manually deleted so that it was useful for simulation). Various hip joint motions were simulated for each of the three cases to identify differences in BROM. It was observed that the average variation in BROM was ±2 0 due to 1mm change (increase or decrease) in sphere diameter from the ideal enlarged sphere diameter.
In clinical context, the developed method can be used as follows. The B-ROM can be generated for both contralateral hip side and planned treatment side of a patient. When a patient's contralateral hip is healthy, the B-ROM of this healthy side can be used to describe subject-specific target values to determine best possible B-ROM of the planned treatment side. There are various approaches to utilise the reference B-ROM. (a) The list of B-ROM data points ( α Lim , β Lim , γ Lim , SeqIds) for contralateral healthy side can be directly used as 'targets' to evaluate the chance of postsurgical impingement for the planned treatment side after device placement. Implant positions that successfully avoid impingements for most of the 'targets' B-ROM can be selected, provided they follow other structural and physiological constraints. (b) Various 3D B-ROM envelopes can be created for alternative planned implant positions. The 3D planned B-ROM are compared with contralateral healthy B-ROM (as demonstrated by Case study 2), and the best plan is selected where the planned B-ROM cover the entire healthy B-ROM, or the intersection of the two B-ROMs is greater. However, both the intersection volume and location of the intersection are important when they were partially intersected. Although this approach provides a better visual representation as the surgeon can check how the planned B-ROM is changing by altering the implant positions, this approach is computationally expensive. (c) If the surgeons are more interested in some particular motion rather the entire 3D B-ROM, this can be easily extracted from the 3D B-ROM enveloped and subsequently, can be used in various other representation styles. (d) A computer animation can be used to check how the 3D B-ROM is directly related to the impingement area. This would help surgeons to decide which hip joint motions are more critical for surgical planning. (e) Although the B-ROM was developed to use in pre-planning tool, this could also be used to plan revision surgery and to check the performance of an already replaced hip joint. This was demonstrated through case study 2 where the patient was selected in such a way that their contralateral side of the replaced hip was healthy. In order to differentiate between bony and prosthetic impingement for planned treatment side, the following approaches could be employed. In the first approach, two simulations would be executed with two different geometries. The first simulation run would be exactly same as explained in the method section (i.e., the pelvis, cup and liner geometries are assembled together to be considered as a single geometry). As a result, the calculated B-ROMs would include impingements due to both osteophytes and cup/liner position. The second simulation run would be performed using the pelvis geometry only (i.e., cup and liner geometries would not be assembled with pelvis). The B-ROMs from the first simulation would then be used for motion simulation, and subsequently, impingement detection. The set of B-ROMs that would produce impingement in second simulation would be due to the osteophytes. Otherwise, it would be due to the cup/liner geometry. The second approach would be to identify the location of the impingement points for each of the limiting B-ROMs. If the coordinates of the impingement points are within the native pelvis geometry, the associated BROM would be due to the osteophytes. The other B-ROMs would be, therefore, resulted due to the cup/liner position.
There were few limitations in the developed method. Firstly, B-ROM was developed by considering only bone geometries without directly including the restriction imposed by the soft tissues and ligaments. Therefore, the B-ROM was not actually a functional ROM (F-ROM) of a hip joint. As a result, some areas (or data points) of B-ROM might have overestimated the ROM in comparison with the functional ROM (F-ROM) boundary. This can be addressed by performing F-ROM measurement of a healthy subject or cadaver hip using optical or wearable systems and correlate the F-ROM with B-ROM. This is the immediate planned future work. Secondly, the method was primarily developed for the unilateral osteoarthritic patient where the contralateral side is healthy. For bilateral osteoarthritis, this method could not be used directly as there would be no contralateral healthy B-ROM. Therefore, the challenge would be to define a B-ROM that could be used as targets for bilateral osteoarthritis sides. This could be possible by extending the method as follows. The pre-arthritic healthy geometries of pelvis and femur of the bilateral osteoarthritic patient could be predicted using statistical shape modelling (SSM). Thereafter, the predicted prearthritic healthy geometries would be used to estimate the healthy B-ROM, and subsequently used it as targets for planning. Another approach would be to define generalised healthy B-ROMs for various group of populations who have common characteristics. Thereafter, if a new patient comes in future, the patient would be associated with a particular classified group based on the common characteristics (e.g., femoral offset, femoral version, side, gender etc.), and the healthy generalised B-ROM for that group would be used as a target B-ROM for surgical planning. Thirdly, variation in geometry construction due to different image resolution were not considered in the model. Future work will be performed to analysis the uncertainty of the model. Finally, this paper has only demonstrated the developed method with its potential use through a verification study of two patients. A larger clinical trial using larger population of patients is required to explore its true clinical potential.

Conclusions
A method was developed to define subject-specific B-ROM which was represented in terms of pure joint motions -flexionextension, abduction-adduction, internal-external rotation and their combinations. The B-ROM was constructed purely from CT developed bone geometries (pelvis and femur) for healthy hip joint, and additional implant geometries (liner, cup and stem) were employed for planned treatment side. The method could be readily included into any exiting 3-dimensional pre-planning tool for hip replacement. The novel aspects of this work were that the developed B-ROM produced a close representation of all the possible ADLs for a particular patient, and the contralateral healthy B-ROM provided better personalised target values with which the ROM of the planned treatment side can be compared. This may provide greater confidence in surgical planning, compared with current state-of-the-art processes that consider only a limited number of selected ADLs with generalised ROM collected from the literature.

Declaration of Competing Interest
The authors declare no conflict of interest.