Reconstruction of Three-Dimensional Tibiofemoral Kinematics Using Single-Plane Fluoroscopy and a Personalized Kinematic Model

: Model-based 3D/2D image registration using single-plane ﬂuoroscopy is a common setup to determine knee joint kinematics, owing to its markerless aspect. However, the approach was subjected to lower accuracies in the determination of out-of-plane motion components. Introducing additional kinematic constraints with an appropriate anatomical representation may help ameliorate the reduced accuracy of single-plane image registration. Therefore, this study aimed to develop and evaluate a multibody model-based tracking (MbMBT) scheme, embedding a personalized kinematic model of the tibiofemoral joint for the measurement of tibiofemoral kinematics. The kinematic model was consisted of three ligaments and an articular contact mechanism. The knee joint activities in six volunteers during isolated knee ﬂexion, lunging, and sit-to-stand motions were recorded with a biplane X-ray imaging system. The tibiofemoral kinematics determined with the MbMBT and mediolateral view ﬂuoroscopic images were compared against those determined using biplane ﬂuoroscopic images. The MbMBT was demonstrated to yield tibiofemoral kinematics with precision values in the range from 0.1 mm to 1.1 mm for translations and from 0.2 ◦ to 1.3 ◦ for rotations. The constraints provided by the kinematic model were shown to effectively amend the nonphysiological tibiofemoral motion and not compromise the image registration accuracy with the proposed MbMBT scheme. using a random forest (RF) model. In the current study, the ratio, φ 3 D 2 D , between the 3D ligament length ( L 3 D ) and its projected 2D ligament length ( L 2 D ) on the mid-sagittal plane was assumed to be associated with tibiofemoral motion and was subject-speciﬁc. To provide personalized ligament length variation during entire activity, for each ligament, an RF model was trained to predict this ratio φ 3 D 2 D at each instant. The input feature vector was composed of the ratio obtained from the non-weight-bearing extended tibiofemoral pose (i.e., during CT scan), φ CT , and the deviations of ﬂexion/extension (FE), adduction/abduction (AA), internal/external rotation (IER), anterior/posterior (AP) translation and proximal/distal (PD) translation of tibiofemoral joint with respect to values at fully extended tibiofemoral pose. The subject- and task-speciﬁc lengths of the ACL, PCL and MCL at each instant were thereafter predicted following a leave-one-out cross-validation scheme, in which all the kinematic data reconstructed using the validated MBT approach [20] from the other ﬁve subjects (excluding the subject under analysis) were used to train the RF model.


Introduction
Accurate measurement of the three-dimensional motion and arthrokinematics of the knee is essential for the functional assessment of the joint. To this end, skin marker-based stereophotogrammetry is commonly used to measure the spatial poses of the thigh and shank. However, the soft tissue artifacts associated with skin markers inevitably reduced the accuracies of the measured tibiofemoral kinematics [1,2]. Appropriate kinematic models, either at the segment level [3] or at the multibody level [4], can be applied with the measured skin markers to better reproduce the spatial poses of the body segment of interest and thereafter more accurately estimate joint angles. Multibody kinematics optimization techniques require the kinematic model of the knee to strictly or loosely link the adjacent segments while providing constraints to unphysiological motions [5].
Various knee joint kinematics models were introduced and embedded into the multibody model with the intent of generating more accurate kinematics of the knee [6]. Among them, the spherical [4,7] and hinge models [8,9], the most simplified and commonly utilized in multibody kinematics optimization, were demonstrated to not enable physiological knee kinematics [6]. While the elastic knee model with a joint stiffness matrix showed the potential to better reproduce physiological kinematics [10], subject-specific customization of such stiffness coefficients should be further addressed. The spatial parallel mechanism represented by rigid links connecting adjacent joint articular surfaces is one of the most popular models in the community [6,[11][12][13][14]. The approach was featured by its direct connection with anatomical structures (i.e., the ligaments and articular contacts) and capability to customize to individual subjects and replicate joint passive motion [15]. Recent studies have shown that the knee model based on a spatial parallel mechanism can be personalized by using personal medical images and 3D kinematics data [16,17]. Accurate representation of knee articular features has shown its potential to drive accurate 3D knee joint kinematics during non-weight-bearing knee flexion and extension motions with submillimeter and subdegree errors [17]. However, the feasibility of replicating accurate knee kinematics during various weight-bearing activities using the precise knee model remains unclear.
Model-based tracking (MBT) techniques integrating the dynamic X-ray imaging system, computed tomography (CT), and 3D/2D image registration methods [18][19][20][21] enabled direct measurement of 3D skeletal kinematics while eliminating the use of skin markers. The 3D pose of a skeleton can be determined by registering its CT model-projected digitally reconstructed radiographs (DRRs) to the X-ray images. The biplane imaging configuration is currently the most commonly used setup to capture the 3D movement of a single joint during particular tasks [22][23][24], as it can provide stereo X-ray images for 3D/2D image registration. The limited intersection volume of measurement and the higher ionizing radiation dose were considered the primary concerns of biplane configuration. While the application of a single-plane imaging system can enlarge the area of measurement (given the same size of the X-ray detector) and lead to lower radiation exposure to subjects, limited measurement accuracy was observed, especially in out-of-plane motion components [25][26][27][28]. This was associated with both the substantial differences between the DRRs and the X-ray images [29] and the perspective projection of the imaging system. The out-ofplane motion of a 3D bone model would have less influence on the appearance of the DRR and thus the objective function values than those by the in-plane motion components [28].
Introducing additional information other than the images during 3D/2D image registration may help ameliorate the reduced accuracy of single-plane image registration. For the analysis of the total knee replacement kinematics, the anti-collision maneuver was introduced to detect the unrealistic penetration between the metal and insert components, which was shown to yield more accurate relative poses between the femoral and tibial components in the out-of-plane translation component [30]. A similar maneuver was also applied on the vertebrae and shown to successfully increase the accuracy of intervertebral motion measurement [31]. However, anticollision detection may not be applicable to the natural knee, as the femur and tibia are not compactly close to each other. Instead, a personalized model with an appropriate representation of anatomical structures may provide useful information that helps guide the motion pattern of the tibiofemoral joint. Personalization of the knee model in vivo is a challenge, as the model parameters were demonstrated to have substantial influences on the model-determined kinematics [32]. Taking magnetic resonance (MR) images is a feasible approach to customize the model parameters [17], but it requires subjects to undergo additional MR scanning and be subjected to intra-and interobserver variability in identifying ligament attachment sites [33]. A feasible alternative is to use a generic digital model embedding known model parameters [34,35]. Personalization of model parameters can be thereafter obtained by proper nonrigid shape transform of the generic model to best-fit available anatomical features of the subject under analysis. The approach can be particularly suitable for the scenario of MBT, as 3D surface bone models are normally obtained from the subject's CT. Since the knee model has the ability to guide the particular motion pattern given the subject-specific model parameters, properly integrating it with the image registration procedure may be beneficial to improve the accuracy of out-of-plane motion components.
The purpose of the present study was to develop a multibody model-based tracking (MbMBT) scheme based on simultaneous 3D/2D image registration of the femoral and tibial bones while considering the kinematic constraints provided by the personalized kinematic model of the tibiofemoral joint. To evaluate the performance of the proposed method, the reconstructed tibiofemoral kinematics of six young volunteers during functional tasks with the MbMBT approach were compared against those determined using a validated biplane MBT technique.

Imaging Acquisition
Six healthy male volunteers (age: 23.2 ± 1.7 years; mass: 74.2 ± 15.4 kg; height: 175.2 ± 4.2 cm) participated in this study. Prior to data acquisition, all subjects provided written informed consent, and the study was approved by the China Medical University and Hospital Research Ethics Committee. Computed tomography image stacks of the knee of all subjects (voxel size: 0.7 mm × 0.7 mm × 0.6 mm) were acquired from CT scans (Optima CT 660; GE Healthcare, Chicago, IL, USA). Semimanual segmentation with a region-growing algorithm was employed to isolate regions of the right femur and tibia/fibula, which were further used to create volumetric and surface bone models [21]. A clinical biplane X-ray imaging system (Allura Xper FD20/20, Philips Medical Systems, Amsterdam, The Netherlands) was set up with two X-ray units placed orthogonal to each other, from which medial-lateral (ML) and anterior-posterior (AP) view fluoroscopic images of an individual subject's right knee during tested activities were acquired (resolution: 512 × 512; grayscale: 8 bit). The imaging system was operated at tube voltages ranging from 50 to 59 kVp, a tube current of 5 mA and a frame rate of 30 fps per X-ray unit. The tested activities performed at self-selected speeds included isolated knee flexion from the standing position, knee bending in the lunge stance and sit-to-stand. Note that while the biplane images were collected, only ML view fluoroscopic images were employed for the proposed MbMBT analysis to mimic the condition that the single-plane X-ray imaging system was used for data acquisition.

Personalize Tibiofemoral Kinematic Model
Anatomical reference frames of the femur and tibia/fibula were defined on the surface bone models ( Figure 1A) in accordance with an automated determination method [36]. The personalized tibiofemoral kinematic model was modified from the parallel mechanism of the knee model [14,17]), in which three linear links connected endpoints of the anterior cruciate ligament (ACL), posterior cruciate ligament (PCL) and deep bundle of the medial collateral ligament (MCL) ( Figure 1B). A torus-on-curved surface contact mechanism modeled the congruent contacts between femoral condyles and the tibial plateau ( Figure 1C). The personalized ligament endpoint locations were estimated using a shape transform of generic bone model imbedding known endpoints to match the subject-specific surface bone models ( Figure 1B), as described in the following subsection. The femoral condyle surface was represented by a geometric torus through least-squares fitting of a torus model to vertices in the posterior part of the femoral condyle ( Figure 1C). Given different tibiofemoral poses, a "contact sphere" inside the revolution of the torus and closest to the tibial plateau was identified ( Figure 1C). The tibial plateau surface of each compartment was represented as the curved surface of a tube determined by least-squares fitting of a tube model to the vertex points in the medial facet of the tibial plateau ( Figure 1C). The centroid of the tube was set to coincide with the contact sphere centroid of the corresponding femoral compartment at the extended pose of the tibiofemoral joint (i.e., the pose during the CT scan). All model geometric parameters are described with respect to the corresponding reference frames. each compartment was represented as the curved surface of a tube determined by leastsquares fitting of a tube model to the vertex points in the medial facet of the tibial plateau ( Figure 1C). The centroid of the tube was set to coincide with the contact sphere centroid of the corresponding femoral compartment at the extended pose of the tibiofemoral joint (i.e., the pose during the CT scan). All model geometric parameters are described with respect to the corresponding reference frames. The personalized endpoint positions of the knee ligaments were estimated by besting-fitting the statistical shape model (SSM, colored bone model) embedding known endpoint coordinates to the subject's bone geometry (semitransparent model). The ligament endpoints from the deformed SSM were then projected onto the nearest face of the subject-specific model, determining the personalized ligament endpoint positions. (C) Congruent contact of the tibiofemoral joint was represented by the best-fitted bilateral torus (black dots) and the curved surfaces of tube geometry (red). Contact spheres (blue) nearest to the tibial plateau would also be identified given a specified tibiofemoral 3D pose.

Ligament Endpoints
The generic model of the femur and tibia imbedding ligament endpoints was built based on the statistical shape model (SSM) technique. The SSMs of the femur and tibia were generated using an SSM generation procedure as suggested in [37,38], and a set of CT-derived surface bone models was obtained from 30 healthy male subjects gathered from our former research projects (Ministry of Science and Technology in Taiwan, ROC). Overall, 60 knee joints, consisting of 30 right knees and 30 mirrored left knees, were recruited for building the SSM. The major steps to build an SSM are briefly described below: (i) The CT-derived surface bone models were converted to point distribution models (PDMs) possessing the same number of nodes at corresponding locations using the coherent point drift algorithm [39] (number of nodes: femur = 3051; tibia = 3032); (ii) generalized Procrustes analysis [40] was performed to best align all PDMs under the same reference frame; and (iii) principal component analysis was applied on aligned PDMs to establish an SSM consisting of a mean shape and modes of shape variations. A new shape (a new point cloud) can be thereafter generated by superposing deformations, controlled by parameters of shape variation modes, onto the mean shape.
Bone templates carrying ligament endpoint coordinates are available in the literature [35], which enables reasonably accurate estimation of personalized ligament endpoint locations, free from the requirement of MR imaging. To attach ligament endpoints onto the SSM, the SSMs of the femur and tibia were deformed and rigidly transformed to match the bone templates by minimizing the sum of squared point-to-face distances. The ligament endpoints were projected onto the nearest triangular face of the deformed SSM and expressed in a barycentric coordinate system. Through this procedure, the SSM of the femur and tibia imbedding the corresponding ligament endpoint were obtained. It can thereafter be used to best fit the CT-derived subject-specific bone models of the subject under analysis through the previously mentioned SSM deformation and rigid transformation procedure ( Figure 1B). The barycentric coordinates of ligament endpoints from the deformed SSM were converted to Cartesian coordinates and were subsequently projected onto the nearest face of the subject-specific bone model and expressed in the

Ligament Endpoints
The generic model of the femur and tibia imbedding ligament endpoints was built based on the statistical shape model (SSM) technique. The SSMs of the femur and tibia were generated using an SSM generation procedure as suggested in [37,38], and a set of CT-derived surface bone models was obtained from 30 healthy male subjects gathered from our former research projects (Ministry of Science and Technology in Taiwan, ROC). Overall, 60 knee joints, consisting of 30 right knees and 30 mirrored left knees, were recruited for building the SSM. The major steps to build an SSM are briefly described below: (i) The CT-derived surface bone models were converted to point distribution models (PDMs) possessing the same number of nodes at corresponding locations using the coherent point drift algorithm [39] (number of nodes: femur = 3051; tibia = 3032); (ii) generalized Procrustes analysis [40] was performed to best align all PDMs under the same reference frame; and (iii) principal component analysis was applied on aligned PDMs to establish an SSM consisting of a mean shape and modes of shape variations. A new shape (a new point cloud) can be thereafter generated by superposing deformations, controlled by parameters of shape variation modes, onto the mean shape.
Bone templates carrying ligament endpoint coordinates are available in the literature [35], which enables reasonably accurate estimation of personalized ligament endpoint locations, free from the requirement of MR imaging. To attach ligament endpoints onto the SSM, the SSMs of the femur and tibia were deformed and rigidly transformed to match the bone templates by minimizing the sum of squared point-to-face distances. The ligament endpoints were projected onto the nearest triangular face of the deformed SSM and expressed in a barycentric coordinate system. Through this procedure, the SSM of the femur and tibia imbedding the corresponding ligament endpoint were obtained. It can thereafter be used to best fit the CT-derived subject-specific bone models of the subject under analysis through the previously mentioned SSM deformation and rigid transformation procedure ( Figure 1B). The barycentric coordinates of ligament endpoints from the deformed SSM were converted to Cartesian coordinates and were subsequently projected onto the nearest face of the subject-specific bone model and expressed in the corresponding anatomical reference frame. In this way, the estimated personalized ligament endpoint locations were obtained for each subject ( Figure 1B).

Ligament Length
The ligament length was simply defined as the linear distance of respective endpoint positions (Figure 2A), and the length variation during tibiofemoral motion was predicted Appl. Sci. 2021, 11, 9415 5 of 16 using a random forest (RF) model. In the current study, the ratio, φ 3D 2D , between the 3D ligament length (L 3D ) and its projected 2D ligament length (L 2D ) on the mid-sagittal plane was assumed to be associated with tibiofemoral motion and was subject-specific. To provide personalized ligament length variation during entire activity, for each ligament, an RF model was trained to predict this ratio φ 3D 2D at each instant. The input feature vector was composed of the ratio obtained from the non-weight-bearing extended tibiofemoral pose (i.e., during CT scan), φ CT , and the deviations of flexion/extension (FE), adduction/abduction (AA), internal/external rotation (IER), anterior/posterior (AP) translation and proximal/distal (PD) translation of tibiofemoral joint with respect to values at fully extended tibiofemoral pose. The subject-and task-specific lengths of the ACL, PCL and MCL at each instant were thereafter predicted following a leave-one-out cross-validation scheme, in which all the kinematic data reconstructed using the validated MBT approach [20] from the other five subjects (excluding the subject under analysis) were used to train the RF model. using a random forest (RF) model. In the current study, the ratio, 2 D φ , between the 3D ligament length ( 3D L ) and its projected 2D ligament length ( 2D L ) on the mid-sagittal plane was assumed to be associated with tibiofemoral motion and was subject-specific. To provide personalized ligament length variation during entire activity, for each ligament, an RF model was trained to predict this ratio 3 2 D D φ at each instant. The input feature vector was composed of the ratio obtained from the non-weight-bearing extended tibiofemoral pose (i.e., during CT scan), CT φ , and the deviations of flexion/extension (FE), adduction/abduction (AA), internal/external rotation (IER), anterior/posterior (AP) translation and proximal/distal (PD) translation of tibiofemoral joint with respect to values at fully extended tibiofemoral pose. The subject-and task-specific lengths of the ACL, PCL and MCL at each instant were thereafter predicted following a leave-one-out cross-validation scheme, in which all the kinematic data reconstructed using the validated MBT approach [20] from the other five subjects (excluding the subject under analysis) were used to train the RF model.
Errors of the estimated ligament length were determined by the differences between model-determined 3D L and those derived from the tibiofemoral kinematics reconstructed using the biplane MBT. Root-mean-squared (RMS) errors were calculated over all image frames for each subject. The efficacy of the approach was evaluated by comparing the 3D L estimated using the model-predicted

Multibody Model-Based Tracking Scheme
The proposed MbMBT scheme mainly consisted of three steps ( Figure 3): (i) singlebody 3D/2D image registration; (ii) multibody image registration (MIR) with two-level optimization designated to correct the joint dislocation; and (iii) MIR with a hybrid objective function and refined kinematics model. Errors of the estimated ligament length were determined by the differences between model-determined L 3D and those derived from the tibiofemoral kinematics reconstructed using the biplane MBT. Root-mean-squared (RMS) errors were calculated over all image frames for each subject. The efficacy of the approach was evaluated by comparing the L 3D estimated using the model-predicted φ 3D 2D and constant φ CT . The Wilcoxon signed rank test was employed for statistical comparison at a significance level of 0.05.

Multibody Model-Based Tracking Scheme
The proposed MbMBT scheme mainly consisted of three steps ( Figure 3): (i) singlebody 3D/2D image registration; (ii) multibody image registration (MIR) with two-level optimization designated to correct the joint dislocation; and (iii) MIR with a hybrid objective function and refined kinematics model.

Single-Body 3D/2D Image Registration
The first estimation of pose parameters of each bone model (expressed via rigid transformation, T k , k = femur and tibia) was determined through an existing single-body 3D/2D image registration approach [21]. To reproduce the projection geometry of the fluoroscopy system, a perspective projection model featuring the focal point, focal length and virtual image plane was reconstructed using a custom-made calibration object [21]. A global reference frame (S G ) was defined originating at the bottom left corner of the image plane with the x-and y-axes directed parallel to the image plane and the z-axis normal to the plane. Within the field of view of the projection model, transformation of the volumetric bone model with respect to S G was estimated by best fitting its DRR to the fluoroscopic image. Normalized cross-correlation of image gradient (GC) was used to quantify the similarity between the DRRs and fluoroscopic images [29]. A genetic algorithm (GA; population size = 50; no. of generations = 60) was employed as the optimizer to search for optimal transformation (design variables: 6 degrees of freedom (6-DOF) of the bone) such that the GC values were maximized. femur and tibia k k = T ) was determined through an existing single-body 3D/2D image registration approach [21]. To reproduce the projection geometry of the fluoroscopy system, a perspective projection model featuring the focal point, focal length and virtual image plane was reconstructed using a custom-made calibration object [21]. A global reference frame ( S G ) was defined originating at the bottom left corner of the image plane with the x-and y-axes directed parallel to the image plane and the z-axis normal to the plane. Within the field of view of the projection model, transformation of the volumetric bone model with respect to S G was estimated by best fitting its DRR to the fluoroscopic image. Normalized cross-correlation of image gradient (GC) was used to quantify the similarity between the DRRs and fluoroscopic images [29]. A genetic algorithm (GA; population size = 50; no. of generations = 60) was employed as the optimizer to search for optimal transformation (design variables: 6 degrees of freedom (6-DOF) of the bone) such that the GC values were maximized.
For each image frame, given the registered transformation, deviations of FE, AA, IER and AP and PD translations of the tibiofemoral joint with respect to those at the extended tibiofemoral pose were calculated. These kinematic deviations together with CT φ were fed into the trained RF model to predict

Multibody Image Registration with Two-Level Optimization
In the second step of the MbMBT, a two-level optimization scheme [21] with MIR was employed to correct artifactual tibiofemoral dislocation. To this end, the femur and tibia/fibula models were simultaneously registered to the fluoroscopic image ( Figure 4A) while minimizing the deformations of the personalized tibiofemoral kinematic model. In the outer-level optimization, an overall 11 parameters describing 6-DOF of the femur and For each image frame, given the registered transformation, deviations of FE, AA, IER and AP and PD translations of the tibiofemoral joint with respect to those at the extended tibiofemoral pose were calculated. These kinematic deviations together with φ CT were fed into the trained RF model to predict φ 3D 2D for each ligament. At the same time, the registered transformations were used to compute L 2D . Multiplication of φ 3D 2D and L 2D would yield an estimation of L 3D used for the subsequent multibody image registration.

Multibody Image Registration with Two-Level Optimization
In the second step of the MbMBT, a two-level optimization scheme [21] with MIR was employed to correct artifactual tibiofemoral dislocation. To this end, the femur and tibia/fibula models were simultaneously registered to the fluoroscopic image ( Figure 4A) while minimizing the deformations of the personalized tibiofemoral kinematic model. In the outer-level optimization, an overall 11 parameters describing 6-DOF of the femur and 5-DOF of the tibia/fibula (excluding z-axis translation) were taken as the design variables, whereas the remaining parameter was the design variable in the inner-level optimization. In the scheme of the two-level optimization, given each set of outer-level 11-parameters, an inner-level subproblem is solved that gives an optimal solution for tibial z-axis translation. The outer-and optimally inner-level parameters were then joined together to spatially transform the volumetric models of the femur and tibia/fibula, from which the combined DRR was generated ( Figure 4B). The outer-level optimization thereafter took the GC values between the combined DRRs and the fluoroscopic image as the objective function value.
In the inner-level subproblem, the outer-level 11-parameters served as the fixed parameters, the tibial z-axis translation (Z tibia ) was determined as the resultant tibiofemoral pose led to minimize ligament elongation relative to L 3D and the shift of contact spheres of the bilateral condyles (S c ) with respect to the tube central axis in the corresponding tibial Appl. Sci. 2021, 11, 9415 7 of 16 plateau ( Figure 2B). To this end, Brent's algorithm [41] was employed to solve the following equation: where L * 3D indicates the ligament length obtained by the resultant tibiofemoral pose, and c represent different ligaments and joint compartments, respectively. an inner-level subproblem is solved that gives an optimal solution for tibial z-axis translation. The outer-and optimally inner-level parameters were then joined together to spatially transform the volumetric models of the femur and tibia/fibula, from which the combined DRR was generated ( Figure 4B). The outer-level optimization thereafter took the GC values between the combined DRRs and the fluoroscopic image as the objective function value. In the inner-level subproblem, the outer-level 11-parameters served as the fixed parameters, the tibial z-axis translation ( tibia Z ) was determined as the resultant tibiofemoral pose led to minimize ligament elongation relative to 3D L and the shift of contact spheres of the bilateral condyles ( c S ) with respect to the tube central axis in the corresponding tibial plateau ( Figure 2B). To this end, Brent's algorithm [41] was employed to solve the following equation: ( )  After completion of image registration over entire images, the registered transformations of the femur and tibia/fibula were used to derive the length of the ligament at each instant. Meanwhile, the magnitudes of PD compression of the joint, as characterized by the vertical shift of the contact spheres with respect to the tube central axis, S V , were obtained ( Figure 2B). A third-order polynomial equation was fitted to each of the abovementioned parameters against the flexion angles (θ F ), given parametric models of ligament length, F (θ F ), and the models of joint compression, F c (θ F )( Figure 2C). These parametric models were intended to provide a simplified description of subject-specific, motion-related and smoothly continuous soft tissue deformations of the tibiofemoral joint.

Multibody Image Registration with Hybrid Similarity Measure
The third step of MbMBT aimed to find compromises T f emur and T tibia , such that the combined DRR best matched the fluoroscopic image ( Figure 4C) while complying with the geometric constraints provided by the joint kinematic model. To achieve this goal, a hybrid objective function consisted of an image similarity term defined as the GC values between the DRR (I D ) and the fluoroscopic image (I X ) and a kinematic model term. The 12-DOF of the femur and tibia are design variables of the optimization. The objective function was formulated as follows, where || denotes the number of elements, δ = 1 is an empirical parameter controlling the weight of the kinematic model term and S * V,c indicates the PD compressions in bilateral joint compartments obtained by the resultant tibiofemoral pose.

Evaluation of MbMBT
For each subject, the subject-specific volumetric model of the femur and tibia/fibula were registered to ML view fluoroscopic images in each motion trial using the proposed MbMBT procedure, giving MbMBT-determined bone pose parameters. The standard reference poses of the bone were also determined with corresponding biplane images using a validated biplane MBT method [20]. Given the measured bone pose parameters, tibiofemoral kinematics were calculated following the convention proposed by Grood and Suntay [42]. The degree of agreement between MbMBT-determined kinematics and reference kinematics was evaluated via a coefficient of determination (R 2 ). Differences between MbMBT-determined and reference kinematic variables were computed. For each motion trial, mean absolute differences (MAD), mean differences (denoted as bias) and standard deviation (SD) of differences (denoted as precision) over all image frames were quantified. The means and SDs of the MAD, bias and precision across all subjects were computed. The Wilcoxon signed-rank test was applied to perform between-method comparisons of kinematic variables. The statistical analyses were conducted using MATLAB.

Results
The lengths of the ACL, PCL, and MCL predicted with the random forest regression model were compared with standard reference values derived using the biplane MBTdetermined tibiofemoral kinematics and the predicted values using a uniform scaling factor φ CT (Figure 5A-C) and were shown to have average RMS errors of 0.55 mm, 0.21 mm and 0.14 mm, respectively ( Figure 5D). The model-predicted PCL length was significantly more accurate than that obtained using a uniform scaling factor ( Figure 5D). While the model-predicted MCL length was not significantly more accurate, the averaged RMS error was nearly half the values obtained using a uniform scaling factor (0.14 mm vs. 0.26 mm, Figure 5D).
The MADs of 6-DOF of individual bone obtained using single-body model-based tracking (SbMBT, equivalent to step-1 of MbMBT) and the proposed MbMBT are shown in Figure 6, in which the median values of MAD for in-plane translations (X-axis, Y-axis) and rotation (Z-axis) were all below 0.3 mm and 0.2 • , respectively. However, the median MADs of the SbMBT-and MbMBT-determined out-of-plane translations (Z-axis) were 3.5 mm and 3.0 mm, respectively. The median values of MAD in out-of-plane rotation components were in the range of 0.3 • to 0.8 • for the two methods. No significant differences were found in any motion component.
For the description of tibiofemoral kinematics during isolated knee flexion, MbMBTdetermined joint center translations were shown to have precision values less than 1.1 mm and significantly more precise than those obtained with SbMBT (Table 1). A similar tendency can also be found in the MAD metrics (Table 1). No considerable differences were found in rotation components between the two methods, except for Add/Abd bias, for which the MbMBT yielded a slightly lower bias (0.5 ± 1.2 • ) ( Table 1). The MbMBT enabled the reproduction of Flex/Ext, I/E rotation and A/P translation highly representative of standard reference kinematics (R 2 > 0.95), moderately representative of joint translation in P/D and L/M components (R 2 : 0.59-0.70) and weakly representative of Add/Abd (Table 1).
For the measurement of lunge activity, precision values were less than 0.6 mm in MbMBT-determined translations and 0.4 • in rotations, in which three translation components and Add/Abd were significantly more precise in the MbMBT method (Table 2).
During the sit-to-stand motor task, the MbMBT gave precision values less than 1.0 mm in translations and 1.0 • in rotations ( Table 3). The use of MbMBT yielded significantly more precise 6-DOF kinematics of the tibiofemoral joint than the SbMBT method ( Table 3).
The MADs in the rotation components were also significantly smaller with the proposed MbMBT method (Table 3). The MADs of 6-DOF of individual bone obtained using single-body model-based tracking (SbMBT, equivalent to step-1 of MbMBT) and the proposed MbMBT are shown in Figure 6, in which the median values of MAD for in-plane translations (X-axis, Y-axis) and rotation (Z-axis) were all below 0.3 mm and 0.2°, respectively. However, the median MADs of the SbMBT-and MbMBT-determined out-of-plane translations (Z-axis) were 3.5 mm and 3.0 mm, respectively. The median values of MAD in out-of-plane rotation components were in the range of 0.3° to 0.8° for the two methods. No significant differences were found in any motion component.   The MADs of 6-DOF of individual bone obtained using single-body model-based tracking (SbMBT, equivalent to step-1 of MbMBT) and the proposed MbMBT are shown in Figure 6, in which the median values of MAD for in-plane translations (X-axis, Y-axis) and rotation (Z-axis) were all below 0.3 mm and 0.2°, respectively. However, the median MADs of the SbMBT-and MbMBT-determined out-of-plane translations (Z-axis) were 3.5 mm and 3.0 mm, respectively. The median values of MAD in out-of-plane rotation components were in the range of 0.3° to 0.8° for the two methods. No significant differences were found in any motion component.   Table 1. The means ± standard deviations of the mean absolute differences (MAD), bias and precision values over all subjects for six degree-of-freedom kinematics of the knee joint during isolated flexion. The coefficient of determination values (R 2 ) were obtained from linear regression analysis on the knee kinematics derived from the SbMBT and MbMBT methods against reference values. Asterisks indicate significant differences in the variables between the MbMBT and SbMBT.   Table 3. The means ± standard deviations of the mean absolute differences (MAD), bias, and precision values over all subjects for six degree-of-freedom kinematics of the knee joint during sit-to-stand. The coefficient of determination values (R 2 ) were obtained from linear regression analysis on the knee kinematics derived from the SbMBT and MbMBT methods against the reference values. Asterisks indicate significant differences in the variables between the MbMBT and SbMBT. The degree of agreement between the MbMBT-determined kinematics and the reference values during the two weight-bearing activities were greater than 0.94 and 0.92, respectively, except for the LM translation (Tables 2 and 3). While the LM translation was still the least accurate component, improvements in MAD, precision and degree of agreement over the SbMBT were observed (Tables 2 and 3).

Discussion
A new MbMBT scheme considering the kinematic constraints provided by the personalized model of the tibiofemoral joint was developed in the current study. In general, the MbMBT approach could more precisely determine the 3D tibiofemoral kinematics than the conventional SbMBT, especially in the translation components (Tables 1-3). This improvement can be attributed to the proposed tibiofemoral joint models, as they provide physical constraints to prevent the nonphysiological motion of the tibia relative to the femur. Simultaneous image registration of the two bones also helps ensure that the originally accurate kinematic components from the SbMBT would not be compromised by joint model imperfection and motion-related kinematic variability. However, the personalized kinematic model was shown to not effectively diminish the 6-DOF errors of individual bones ( Figure 6), similar to the strategies of the anti-collision maneuver at the joint [30,31].
To better reflect the physiological deformation of the ligament in the joint model, the present study first estimated the ligament length in 3D space (i.e., L 3d ) from its sagittal plane projection length (i.e., L 2d ) using the trained RF model with the prior tibiofemoral kinematics obtained from single-body 3D/2D image registration (step-1). The approach was developed to account for the fact that the ligament length changes along with the tibiofemoral motion [43] and presents individual differences. Since the fluoroscopic detector is often configured to capture ML view images of the knee joint, the out-of-plane translation error of the bone should have a minor influence on the estimation of the L 2d . Thus, the relatively accurate L 2d was taken to serve as the baseline length of the ligament, and its scaling factors to L 3d were predicted with the trained model. The results indicated that the predicted scaling could yield significantly more accurate PCL lengths than those obtained with uniform scaling ( Figure 5D). The length of the ACL predicted with RF regression showed the least accurate result among the three ligaments, in which the RMS error was near 0.6 mm close to that obtained with uniform scaling ( Figure 5D). This may indicate that the scaling factors between the sagittal and 3D lengths of the ACL were not sufficiently explained by tibiofemoral motions. Further exploration is required to better predict the ACL length variation.
The presented three-step MbMBT scheme originated from the observations that the SbMBT using single-plane images could normally determine accurate in-plane translations and rotations, moderate accuracies in out-of-plane rotations and poor out-of-plane translations [21,25,28], which often propagated to result in artifactual tibiofemoral dislocation. Thus, the MIR with two-level optimization in step-two was designed to solely amend the z-axis translation of the tibia by the tibiofemoral kinematic model, while the other DOFs were merely fine-tuned to ensure image matching between the DRR and the fluoroscopic image. With the strategy, the outer-level optimization left the control of the least accurate z-axis translation to the inner-level subproblem, where the tibiofemoral model dominates the relative z-axis translation of two bones. After completion of the step-two analysis, the model parameters, including variations in ligament length and joint compression, were updated according to the tibiofemoral kinematics to better represent the "subject-specific" and "motion-related" soft tissue deformation during a particular task. Polynomial fitting can also ensure the smoothness of the tissue deformation against knee flexion, which enables improvement of the smoothness of tibiofemoral motion across multiple image frames under the guidance of the joint model in step-three analysis.
The kinematic constraints provided by the personalized tibiofemoral model were demonstrated to effectively ameliorate three-axis translational precision of the joint during all the tested functional activities (Tables 1-3). In particular, the L/M translation, typically the least accurate component in SbMBT, obtained an over 65% improvement in precision (average precision values: from 3.6 mm to 1.1 mm during isolated flexion, from 3.2 mm to 0.6 mm during lunge, and from 2.9 mm to 1.0 mm during sit-to-stand; Tables 1-3). This improvement might be attributed to the step-two analysis, where the two-level optimization effectively separates the guidance of the joint model on the out-of-plane translation (z-axis) of the tibia. The other 11-DOFs were thereafter fine-tuned given the model-determined tibial z-axis translation to maximize the image registration criteria. This compromise strategy between the kinematic model and the image matching was shown to not notably influence the accuracies in individual bone 6-DOFs ( Figure 6) but slightly adjusted the relative poses of the two bones to increase the precision in determining the three translation components (Tables 1-3).
In the present study, the contribution of the personalized kinematic model to the accuracy of tibiofemoral rotations was relatively minor (Tables 1-3). The reasons may be twofold. First, the SbMBT-determined Flex/Ext and I/E rotation during the three activities were sufficiently representative of standard reference kinematics, as revealed in agreement analysis where the R 2 values were above 0.98. The capability of the model is not adequate to further diminish the errors in these two components owing to the model imperfections and uncertainties in the model parameters. Second, our MbMBT strategy was not intended to strictly guide the joint rotations by the kinematic model, as the most suitable model parameters for a particular motion task were unknown. In other words, discrepant tibiofemoral rotation patterns are sure to result in different behaviors in ligament elongation [43,44] and articular contacts [45] and should have different settings for the model parameters. Therefore, such subject-specific and motion-related model parameters were represented by polynomial equations determined using tibiofemoral kinematics after the completion of step-two analysis. The role of the kinematic model in the determination of joint rotations was to provide physical constraints to the bone pose during image registration, which enabled eliminating largely deviated joint angles and improving the angle smoothness throughout the motion task. The effectiveness of the model could be observed in the Add/Abd component (Tables 1-3), in which the R 2 values were slightly increased with the MbMBT approach. The precision values were also statistically improved in the two weight-bearing tasks (Tables 2 and 3).
The model inaccuracy is believed to have an impact on the performance of the MbMBT. In the current study, subject-specific ligament insertion points were not available, as MR images were absent. Instead, bone templates embedding known ligament endpoints were deformed to predict the subject's ligament insertion sites. However, the subjects enrolled in the present study were young Chinese people, but the reference ligament endpoint coordinates were not obtained from similar groups [35]. Potential morphological differences may be present in the bones and soft tissues among different ethnic groups and ages. Some studies have shown morphological differences in the knee [46] and ankle [47] joints between Caucasian and Chinese populations. A validation study is also required to assess the endpoint prediction accuracy and its influences on the performance of the MbMBT in the determination of tibiofemoral kinematics.

Conclusions
The present study proposed a new framework to integrate single-plane fluoroscopy and a kinematic model of the tibiofemoral joint for the measurement of 3D tibiofemoral kinematics. The new MbMBT scheme considering the kinematic constraints provided by the personalized kinematic model has been evaluated in vivo. The kinematic model, which consisted of three ligaments and an articular contact mechanism, was shown to help amend the nonphysiological motion of the tibia relative to the femur. The MbMBT was demonstrated to provide a more precise determination of the tibiofemoral kinematics in comparison to the SbMBT in all the tested motion tasks. In further studies, if subject-specific parameters regarding the ligament endpoints and elongation patterns are available, the presented kinematic model may be beneficial to further improve the performance of the MbMBT.