An Automated Three-Dimensional Bone Pose Tracking Method Using Clinical Interleaved Biplane Fluoroscopy Systems: Application to the Knee

Featured Application

A new fully automated model-based interleaved biplane image tracking scheme in conjunction with clinical asynchronous biplane fluoroscopy was proposed and evaluated for its performance on dynamic tibiofemoral kinematics measurements. The approach may encourage the further use of clinical imaging systems for the noninvasive and precise examination of three-dimensional, dynamic joint functions in clinical practice and extend the application of model-based tracking techniques for orthopedic biomechanical investigations in a budge-efficient way.

Abstract

Model-based tracking of the movement of the tibiofemoral joint via a biplane X-ray imaging system has been commonly used to reproduce its accurate, three-dimensional kinematics. To accommodate the approaches to existing clinical asynchronous biplane fluoroscopy systems and achieve comparable accuracy, this study proposed an automated model-based interleaved biplane fluoroscopy image tracking scheme (MIBFT) by incorporating information of adjacent image frames. The MIBFT was evaluated with a cadaveric study conducted on a knee specimen. The MIBFT reproduced skeletal poses and tibiofemoral kinematics that were in good agreement with the standard reference kinematics provided by an optical motion capture system, in which the root-mean-squared (Rms) errors of the skeletal pose parameters ranged from 0.11 to 0.35 mm in translation and 0.18 to 0.49° in rotation. The influences of rotation speed on the pose errors were below 0.23 mm and 0.26°. The MIBFT-determined bias, precision, and Rms error were comparable to those of the reported model-based tracking techniques using custom-made synchronous biplane fluoroscopy. The results suggested that the further use of the clinical imaging system is feasible for the noninvasive and precise examination of dynamic joint functions and kinematics in clinical practice and biomechanical research.

1. Introduction

Accurate and noninvasive measurement of the three-dimensional (3D) kinematics of human limbs is fundamental for exploring their biomechanical characteristics [1]. With quantified data of rigid-body kinematics and cartilage contact arthrokinematics of the tibiofemoral joint, both subject- and task-specific biomechanics and pathomechanics can be explored further [2,3], leading to a more objective and precise evaluation of functional performance in normally healthy [4,5,6], injured [7], and treated knee joints [8,9]. Direct measurement of skeletal motion has also been increasingly employed for the assessment of other motion capture instruments affected by soft tissue artifacts [10,11,12].

To gain insight into the internal structures of the joint, medical imaging facilities are typically employed to acquire cross-sectional slices or projection images of the target organ. Among these, X-ray fluoroscopy [13], four-dimensional computed tomography (4D CT) [14], and dynamic magnetic resonance (MR) image technologies [15,16] enable the recording of continuous two-dimensional (2D) or 3D images of joint motion at a limited temporal or spatial resolution. The registration of volumetric bone models obtained from high-resolution CT or MR images to a series of continuous images of dynamic joint movement, which is commonly called model-based tracking, can take advantage of the functionalities of different image modalities, yielding continuous six-degree-of-freedom (6-dof) data of skeletal and joint motion. Model-based tracking can be conducted fully based on MR images by registering MR-derived bone models to each frame of a relatively lower-resolution set of dynamic MR image slices [17,18,19] or volumetric images [20]. However, closed-bore MR image scanners restrict subjects from performing limited motions and are also not suitable for patients with metal prostheses.

The application of standard X-ray fluoroscopy or newly developed mobile fluoroscopic imaging systems enable more dynamic kinematic measurements during functional tasks [21,22]. By numerically registering the subject-specific model of the prosthesis or bones to the single-plane fluoroscopic image, model-based tracking routines reproduce 6-dof kinematics of natural knee [23,24] and replacement joints [25,26]. While experimental evaluations showed that the methods enabled the determination of accurate planar translations and rotations, the out-of-plane motion components were found to be substantially less accurate [27,28]. Introducing an additional synchronous X-ray unit yields a framework of stereo X-ray imaging [29], which has been widely demonstrated to be useful for determining 3D tibiofemoral kinematics in submillimeter and subdegree accuracy with biplane 3D to 2D image registration [30,31,32].

Synchronous biplane X-ray imaging of dynamic knee motion is currently achieved by applying custom-made high-speed stereo radiography [33,34] or integrating two C-arm fluoroscopies [35,36], but it is infeasible in clinical biplane X-ray imaging systems [37]. To minimize cross-scatter contamination [38], the clinical imaging system is normally operated with an asynchronous X-ray pulse and acquisition scheme between two X-ray units, generating interleaved biplane fluoroscopic images with a constant temporal offset. The interleaved images have adverse effects on the biplane model-based tracking of the joint motion, however, leading to erroneous skeletal kinematics estimation [39,40]. To address the limitations of the interleaved biplane images, new model-based tracking methods were introduced with strategies related to kinematic interpolation [39,40]. The method termed motion component partition and interpolation (MCPI) [40] was assessed for its accuracy in determining tibiofemoral kinematics in our previous study. While the MCPI improved translation accuracy over single-plane model-based tracking, rotational errors were not significantly improved [40]. Valid kinematic interpolation also relies on precisely registered bone poses in neighboring frames obtained by single-plane image registration [39,40], in which a failed determination of bone poses yields adverse effects on the outcome of kinematic interpolation and therefore the accuracy of subsequent model-based tracking. While precisely adjusting the initial pose of the bone may preclude the issue, minimization of manual intervention is preferable in an automated model-based tracking routine.

To accommodate the interleaved biplane fluoroscopic images obtained from clinical fluoroscopy and eliminate reliance on kinematic interpolation, the present work aimed to develop a new automated model-based tracking scheme for the noninvasive measurement of the 3D kinematics of tibiofemoral motion with an accuracy comparable to those using custom-made synchronous biplane X-ray imaging system. A motion experiment with a knee specimen was performed to evaluate the accuracy of the new model-based tracking method on the determination of 3D bone poses and joint kinematics during dynamic tibiofemoral motions. The influence of the joint rotation speeds on the measurement errors was also quantified.

2. Materials and Methods

2.1. Overview

The dynamic tibiofemoral motions were recorded using a clinical asynchronous biplane fluoroscopy imaging system. The CT data of the same joint were also acquired and used to create volumetric models of the femur and tibia. The 3D pose parameters of the bone in each fluoroscopic frame are determined using the new model-based tracking scheme composed of three major stages, namely (i) 2D/2D template registration; (ii) motion-compensated frame interpolation; (iii) biplane 3D/2D image registration. The registered bone poses for the entire image sequence were then used to compute 3D kinematics of the tibiofemoral joint. A motion experiment using a cadaveric knee joint was performed to evaluate the measurement accuracy of the proposed method. Detailed illustrations of the model-based tracking and motion experiment are provided in the following subsections.

2.2. Asynchronous Biplane Fluoroscopy Imaging

A clinical asynchronous biplane fluoroscopy system was utilized to acquire interleaved biplane 2D fluoroscopic images of the knee in movement from two different perspectives with a 1/2 frame time offset (Figure 1A). Two X-ray units (consisting of an X-ray source and flat-panel detector) positioned orthogonally enabled the acquisition of two perpendicular projections of the knee. Each X-ray unit was represented by an X-ray point source and an image plane, and the overall configuration of the two X-ray units could be virtually reconstructed with known intrinsic and extrinsic parameters. A custom-built calibration box consisting of two parallel acrylic sheets embedded with lead beads at known positions was imaged by fluoroscopy to determine intrinsic parameters, $Q_j$, associated with the $j\mathrm{th}$ X-ray units ($j=1\mathrm{or}2$) [41]. Each image plane was attached with a coordinate system, $F_{xyz}^j$, with the $x_j$ and $y_j$ axes parallel to the image plane and the $z_j$ axis normal to the image plane (Figure 1C). Transformation (or extrinsic parameters) between the $F_{xyz}^j$ of the two X-ray units was determined following an established procedure [42] using another custom-built calibration object [40]. A global coordinate system $G_{XYZ}$ was also defined to coincide with the $F_{xyz}^j$ of the 1st X-ray unit.

2.3. CT-Based Bone Model

The knee joint was CT scanned to acquire the volumetric image of the CT radiodensities of the knee. Regions of the femur and tibia were segmented using a region-growing method [43]. A volumetric bone model was created by extracting a subvolume from CT images containing the segmented region of the bone of interest. The Hounsfield unit values of the voxels belonging to the region of the bones were transformed to attenuation coefficients [44], while other voxel values were directly set to zero (≈ attenuation coefficient of air). Surface bone models were constructed with a marching cubes algorithm [45] and segmented bone regions. Anatomical frames were then created using an automatic determination method [46].

2.4. Model-Based Interleaved Biplane Fluoroscopy Image Tracking

The new model-based tracking scheme for automated determination of 3D bone poses, termed model-based interleaved biplane fluoroscopy image tracking (MIBFT), was accomplished in three stages: (1) 2D/2D template registration, (2) frame interpolation, and (3) biplane 3D/2D image registration. The 6-dof of the bone at frame k is expressed by a transformation matrix $T({\overrightarrow{\mathrm{P}}}^k)$, parameterized by pose parameters ${\overrightarrow{\mathrm{P}}}^k={\left\{T_X^{},T_Y^{},T_Z^{},{\theta }_X^{},{\theta }_Y^{},{\theta }_Z^{}\right\}}^T$. The three translation components $T_X^{},T_Y^{},T_Z^{}$ and rotation components ${\theta }_X^{},{\theta }_Y^{},{\theta }_Z^{}$ are represented with respect to $G_{XYZ}$.

2.4.1. Two-Dimensional/2D Template Registration

A 2D/2D template registration procedure was devised to provide an initial estimation of the transformation of the bone model for each image frame, $T_{treg}=T(R_{treg}^k,v_{treg}^k)$, consisting of a rotation matrix $R_{treg}^k({\theta }_X^{},{\theta }_Y^{},{\theta }_Z^{})$ and a 3D position vector $v_{treg}^k(T_X^{},T_Y^{},T_Z^{})$. For the first image frame being analyzed, $I_{}^k$, a 2D template of the bone of interest, ${\mathrm{M}}_{j,j=1\mathrm{or}2}$, was initialized by extracting a bounding image region for the bone from the preceding image frame $I_{}^{k-1}$ (Figure 2B). The border of the template was determined by the outermost projection boundary of the registered surface bone model (Figure 2A). The template ${\mathrm{M}}_j$ was allowed to approximately match the succeeding image frame $I_{}^{k+1}$ by adjusting its in-plane positions and rotation ${\overrightarrow{p}}^j={\{t_{x_j}^{},t_{y_j}^{},r_{z_j}^{}\}}^T$, as $I_{}^{k-1}$ and $I_{}^{k+1}$ were both acquired from the $j\mathrm{th}$ X-ray unit (Figure 2B).

A sampling importance resampling particle filter [47] was employed to perform template matching and motion tracking (number of particles = 3000), in which the state variables of particle $i$ at frame $k-1$ involved the planar pose parameters of the template on the 2D image plane and their first derivatives ${\overrightarrow{S}}_{j,i}^{k-1}=\{{\overrightarrow{p}}_{j,i}^{k-1},{\dot{\overrightarrow{p}}}_{j,i}^{k-1}\}$. The particle state evolved from time steps of $I_{}^{k-1}$ to $I_{}^{k+1}$ via a linear kinetic model, which yielded a set of ${\overrightarrow{S}}_{j,i}^{k+1}=\{{\overrightarrow{p}}_{j,i}^{k+1},{\dot{\overrightarrow{p}}}_{j,i}^{k+1}\}$. Evolved particles were then resampled [47] in accordance with the normalized weight of the particles, which led to discarding particles with smaller weight. The weight of particle $i$, ${\omega }_i$, was defined by the probability density function value of a normal distribution function, evaluated using the mean absolute difference between ${\mathrm{M}}_j$ and $I_{}^{k+1}$:

$${\omega }_i=\mathrm{N}\left(e_i|\mu ,\sigma \right),$$

(1)

$$e_i=\frac{1}{\left|{\mathrm{M}}_j\right|}{\displaystyle \sum _{q\in {\mathrm{M}}_j}\left|D\left({\mathrm{M}}_j(q),{\overrightarrow{p}}_{j,i}^{k+1}\right)-I_{}^{k+1}(q)\right|},$$

(2)

where $| |$ denotes the number of pixels, and $D$ represents the function that displaces template ${\mathrm{M}}_j$ given a template pose ${\overrightarrow{p}}_{j,i}^{k+1}$. The final template 2D pose ${\overrightarrow{p}}_j^{k+1}=\left\{t_{x_j}^{k+1},t_{y_j}^{k+1},r_{z_j}^{k+1}\right\}$ was estimated by averaging the poses ${\overrightarrow{p}}_{j,i}^{k+1}$ of the resampled particles. After template registration, the bone template was updated by linearly interpolating the intensity values from $I_{}^{k+1}$ at the grid points of the template.

Increments in in-plane rotation $\Delta R_j^{2D}$ and translation $\Delta v_j^{2D}$ of the template (Figure 2C) between preceding and current frames $\left(I_{k-1}\to I_k\right)$ were approximated by halving the motion increments from preceding to succeeding frames $\left(I_{k-1}\to I_{k+1}\right)$ and were expressed in $G_{XYZ}$ following the equations:

$$\Delta v_j^{2D}=R_{G_{XYZ}}^{F_{xyz}^j}{\left[\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(t_{x_j}^{k+1}-t_{x_j}^{k-1}),\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(t_{y_j}^{k+1}-t_{y_j}^{k-1}),0\right]}^T,$$

(3)

$$\Delta R_j^{2D}=R_{G_{XYZ}}^{F_{xyz}^j}R(\varphi ),$$

(4)

$$\varphi =\frac{1}{2}\left(r_{z_j}^{k+1}-r_{z_j}^{k-1}\right),$$

(5)

where $R_{G_{XYZ}}^{F_{xyz}^j}$ describes the rotation between coordinate systems $G_{XYZ}$ and $F_{xyz}^j$. The above procedure was also carried out between $I^{k-2}$ and $I^k$, which approximated the in-plane rotation and translation increments from the other X-ray view. Given $\Delta R_j^{2D}$ and $\Delta V_j^{2D}$ from the two X-ray views and the registered pose in the preceding frame $T(R^{k-1},v^{k-1})$, $R_{treg}^k$ and $v_{treg}^k$ were estimated by:

$$R_{treg}^k=\left({\displaystyle \prod _j\Delta R_j^{2D}}\right)R^{k-1},$$

(6)

$$v_{treg}^k=v_{}^{k-1}+\Delta v_{}^{3\mathrm{D}},$$

(7)

where the 3D incremental translation $\Delta v_{}^{3D}$ is determined by minimizing the differences between its 2D projections and $\Delta v_j^{2D}$ on the two X-ray views and is defined as:

$$\Delta V_{}^{3D}=\arg \underset{\Delta V_{}^{3D}}{\min }{\displaystyle \sum _j{\left|\left(\left(\left(z_j\times \Delta v_{}^{3D}\right)\times z_j\right)\cdot \frac{\left|v_j^{k-1}-s_j\right|}{\left|v^{k-1}-s_j\right|}\right)-\Delta v_j^{2D}\right|}^2},$$

(8)

where $s_j$ is the position of the X-ray source of the $j\mathrm{th}$ X-ray unit, $v_j^{k-1}$ is the position of the perspective projection of $v_{}^{k-1}$ onto the $j\mathrm{th}$ image plane, and $z_j$ is the unit normal to the $j\mathrm{th}$ image plane (Figure 2C). The Levenberg–Marquardt nonlinear least-squares method was used to solve for (8) [48].

2.4.2. Frame Interpolation

To convert the interleaved biplane images into a synchronous image pair, the missing image at time point $I^k$ was replenished by the intermediate image between $I_{}^{k-1}$ and $I^{k+1}$, which was created using a motion compensated frame interpolation (MCFI) method [49] (Figure 1B). In the MCFI, for each 8 × 8 image block within the intermediate image, a bidirectional motion vector was estimated by minimizing the sum of differences between the enlarged 12 × 12 image blocks in $I_{}^{k-1}$ and $I^{k+1}$ in the direction of that motion vector (Figure 1B). After the estimation of the motion vectors, the pixel intensities of the intermediate image were computed using an overlapped block-motion compensation method [50]. The originally acquired fluoroscopic image $I_{}^k$ together with the interpolated image, $J^k$, from the other X-ray view thus formed a configuration of biplane fluoroscopic imaging (Figure 1C). It should be noted that the frame interpolation was only conducted on a rectangular region, defined by the union of successive templates, to minimize computations on image regions irrelevant to the target bone (Figure 1C).

2.4.3. Biplane 3D/2D Image Registration

A 3D/2D image registration procedure based on the forward projection model and biplane fluoroscopic images was implemented to precisely determine the pose parameters of the bones (Figure 3). A forward projection model $P$ [51] constructed with ray-tracing with trilinear interpolation [52] was utilized to quickly generate a digitally reconstructed radiography image (DRR, $D_j$) given the pose parameters $T(\overrightarrow{\mathrm{P}})$, the volumetric bone model $V$, and calibrated intrinsic parameters $Q_j$ (Equation (10)) (Figure 3). A set of pose increments $\Delta {\overrightarrow{\mathrm{P}}}_k={\left\{\Delta T_X^{},\Delta T_Y^{},\Delta T_Z^{},\Delta {\theta }_X^{},\Delta {\theta }_Y^{},\Delta {\theta }_Z^{}\right\}}^T$ from the template-registered pose $T_{treg}$ was determined by minimizing the weighted sum of criterion functions $\mathrm{CF}$ using a vectorized genetic algorithm (GA) [53], followed by the Nelder–Mead Simplex method [54]:

$$\Delta {\overrightarrow{\mathrm{P}}}_{}^k=\arg \underset{\Delta {\overrightarrow{\mathrm{P}}}_k}{\min }{\omega }_I\cdot \mathrm{CF}\left(D_{j=1\mathrm{or}2},I_{}^k\right)+{\omega }_J\cdot \mathrm{CF}\left(D_{j=2\mathrm{or}1},J_{}^k\right),$$

(9)

$$D_j=P\left(T({\overrightarrow{\mathrm{P}}}^k),V,Q_j\right),$$

(10)

$$T({\overrightarrow{\mathrm{P}}}^k)=T(\Delta {\overrightarrow{\mathrm{P}}}^k)T_{treg},$$

(11)

where $\mathrm{CF}$ is set to be the “gradient difference” similarity measure between the digitally reconstructed radiographs and the fluoroscopic images (Figure 3). The weighting factors ${\omega }_I$ and ${\omega }_J$ were set to be 1.0 and 0.8, respectively, as determined in a preliminary evaluation. The number of populations in GA was set to 50, as this was the largest number of DRRs that could be generated in parallel in our current implementation. The maximum number of generations in GA was empirically set to 60, because a larger number of generations did not improve the pose estimation accuracy in our preliminary evaluation. All procedures of the MIBFT, including template registration, frame interpolation, and 3D/2D image registration, were implemented using MATLAB (MATLAB 2018b, MathWorks Inc., Natick, MA, USA).

2.5. In Vitro Motion Experiment

This study was approved by China Medical University and Hospital Research Ethics Committee (No: DMR101-IRB1-139). The knee specimen used in the current in vitro experiment was obtained from a male donor who received a below-hip amputation for disorders other than the knee joint and had provided written informed consent. Access to information identifying the donor was not provided to the authors.

To evaluate the performance of the proposed MIBFT in measuring tibiofemoral kinematics with interleaved biplane fluoroscopy, we conducted a simulated knee-motion experiment on the knee specimen described as follows. Lengths of the femur and tibia of the fresh-frozen knee specimen were approximately 260 and 290 mm, respectively. A plastic connector was affixed onto the proximal end of the femur and distal end of the tibia with bone cement (Figure 4A). Four marker clusters, each attached with four infrared retroreflective markers with a 7 mm radius, were rigidly implanted into the femur and its connector. Three other clusters were implanted into the tibia and its connector. The knee specimen with the clusters and connectors was scanned using a clinical CT scanner (Optima CT660, GE Healthcare, Chicago, IL, USA) to acquire CT images of the whole construct (peak voltage = 120 kVp, tube current = 180 mA, pixel size = 0.709 mm × 0.709 mm, slice thickness = 0.625 mm, matrix size = 512 × 512 × 945).

A clinical biplane fluoroscopy system (Allura Xper FD20/20, Philips Medical Systems, Amsterdam, The Netherlands) was configured with source-to-detector distances of 1195 and 1050 mm for the two X-ray units and was operated with tube voltages ranging from 50 to 59 kVp, a tube current of 5 mA and a frame rate of 30 frames/s per X-ray unit, which allowed us to obtain 512 × 512 8-bit interleaved fluoroscopic images at an effective frame rate of 60 frames/s. Prior to fluoroscopic image acquisition, the specimen was mounted on a custom-made motion testing frame through the connector. While the femoral connector was rigidly mounted to the frame, the femur was kept immobile, and the tibia was allowed to manually rotate about the knee by the operator, and vice versa. Following the cadence from a metronome, the operator extended and flexed the knee by moving the end of the freely moving segment at five rotation speeds (0.25 cycle/s, 0.375 cycle/s, 0.5 cycle/s, 0.625 cycle/s, and 0.75 cycle/s), during which the interleaved biplane fluoroscopic images were acquired. These examinations were termed femoral and tibial isolated flexion/extension (FE). The connector was then detached from the frame, and the specimen was manually bent and extended by the operator to acquire fluoroscopic images of the tibiofemoral FE at four prescribed speeds (0.375 cycle/s, 0.5 cycle/s, 0.625 cycle/s, and 0.75 cycle/s). Overall, three trials per motion speed, each with at least one complete motion cycle, were collected for the femoral isolated FE, tibial isolated FE, and tibiofemoral FE scenarios. During the motion examinations, the 3D trajectories of the markers on the cluster frames were reconstructed using an optical motion capture system (MoCap) with 14 calibrated and synchronized infrared cameras (VICON 612, Oxford Metrics, Oxford, UK).

2.6. Evaluation of the MIBFT

2.6.1. Standard Reference Determination

The reconstructed marker coordinates from the MoCap system were expressed in $G_{XYZ}$ (Figure 4B) after a spatial transformation that was obtained by registering $G_{XYZ}$ to the coordinate system of the MoCap using a reported calibration procedure [55] with a custom-made object. The same set of markers presented on the CT images were also semimanually segmented and reconstructed into polygonal models.

The 3D coordinates of the markers were herein extracted by fitting spheres of corresponding radius to the polygonal models and were expressed in their corresponding anatomical frames. The standard references of the bone pose parameters were determined by coregistering the coordinates of the corresponding markers in the anatomical frame and in $G_{XYZ}$ in the least-squares sense [56]. The standard reference kinematics of the tibiofemoral joint were determined according to standard reference poses of the femur and tibia, and the tibiofemoral angles were expressed in a commonly used joint coordinate system [57]. The tibiofemoral translations were described as the linear displacement of the midpoint of the epicondyles with respect to the tibial anatomical frame. The rotation speeds of the tibiofemoral motions were also derived using the quantified flexion angles at each time instants.

2.6.2. Error Metrics

For each trial in the cadaveric experiment, manual adjustment of the bone pose followed by single-plane 3D/2D image registration [40] was only performed for the first image frame (k = 1). The MIBFT was then executed to determine the bone pose parameters in the following successive image frames (k = 2, 3, …, N). The pose parameters in the first and last image frames were discarded for the error analysis, because they were not included in MIBFT analysis due to the lack of both a preceding and succeeding frame. A 64-bit personal computer with an Intel CPU (Core i7-9700K, 3.6 GHz), an NVIDIA graphics card (GeForce RTX 2080), and 16 GB of RAM were used for all numerical computations.

The bone pose errors were quantified as the deviations of the MIBFT-determined pose parameters from the standard references. The absolute errors at each registration stage and for different tibiofemoral flexion speeds were analyzed. For each trial, the bias (i.e., mean of the errors), precision (i.e., standard deviation (SD) of the errors), and root-mean-squared errors (Rms errors) across all image frames were computed. Descriptive statistics of these error metrics over all trials were obtained. To express the overall positioning error of a bone, a commonly used metric, the mean target registration error (mTRE) [58], was also quantified for each image frame. The agreement between the MIBFT-determined tibiofemoral angles and translations and those obtained from the standard references were quantified by limits-of-agreement (LoA) plots and analysis [59].

3. Results

The Rms error of the MIBFT-determined femoral pose parameters ranged from 0.11 to 0.35 mm in translation and 0.18 to 0.47° in rotation (

Table 1. Means (standard deviation) of the root-mean-squared errors (Rms errors) of each motion component over all trials for the femur and tibia during isolated femoral and tibial flexion/extension (FE) and tibiofemoral FE. The averaged bias and precision over all trials are also shown.

Component	Error Metric	Femoral FE	Tibial FE	Tibiofemoral FE
		Femur	Tibia	Femur	Tibia
Tx (mm)	Rms error	0.28 (0.04)	0.23 (0.03)	0.15 (0.03)	0.18 (0.05)
	Bias ± precision	−0.07 ± 0.28	0.08 ± 0.22	−0.07 ± 0.13	0.08 ± 0.16
Ty (mm)	Rms error	0.35 (0.04)	0.22 (0.03)	0.11 (0.01)	0.20 (0.04)
	Bias ± precision	0.27 ± 0.23	0.12 ± 0.18	0.02 ± 0.10	0.13 ± 0.15
Tz (mm)	Rms error	0.29 (0.05)	0.26 (0.05)	0.18 (0.02)	0.22 (0.04)
	Bias ± precision	−0.13 ± 0.26	0.01 ± 0.26	−0.04 ± 0.17	0.00 ± 0.22
θx (°)	Rms error	0.20 (0.06)	0.30 (0.07)	0.19 (0.04)	0.31 (0.13)
	Bias ± precision	0.03 ± 0.20	0.12 ± 0.27	0.04 ± 0.19	0.14 ± 0.27
θy (°)	Rms error	0.28 (0.04)	0.43 (0.05)	0.47 (0.10)	0.49 (0.11)
	Bias ± precision	−0.03 ± 0.28	0.09 ± 0.42	−0.23 ± 0.41	0.07 ± 0.46
θz (°)	Rms error	0.20 (0.06)	0.29 (0.05)	0.18 (0.02)	0.33 (0.14)
	Bias ± precision	−0.03 ± 0.19	0.02 ± 0.29	0.03 ± 0.18	0.01 ± 0.33

). The corresponding ranges of Rms error for tibial poses were 0.18 to 0.26 mm in translations and 0.29 to 0.49° in rotation (Table 1). The absolute value of the average bias over all trials was less than 0.27 mm in translation and 0.23° rotation for both the femur and tibia. The corresponding values for the precision were 0.28 mm and 0.46° (Table 1). The errors in the pose parameters were successively reduced throughout the MIBFT process, in which the initial pose parameters taken from the registered pose in the preceding frame gave averaged absolute errors up to 2 mm in translations and 1.1° in rotations that were reduced to 0.6 mm and 0.5°, respectively, after 2D/2D template registration. The ultimate averaged absolute errors were down to 0.3 mm and 0.3° after the 3D/2D image registration (Figure 5).

The means (SDs) of the tibiofemoral FE speeds for trials with different rotation speeds were quantified as 42.6 (16.4)°/s, 65.7 (32.0)°/s, 86.7 (40.7)°/s, 121.0 (54.7)°/s, and 154.5 (68.6)°/s for the isolated femoral FE tasks; 51.6 (24.3)°/s, 74.2 (34.7)°/s, 96.9 (44.3)°/s, 122.5 (50.1)°/s, and 158.5 (68.6)°/s for the isolated tibial FE tasks, and 84.6 (41.2)°/s, 101.4 (34.3)°/s, 121.5 (62.8)°/s, and 145.0 (80.7)°/s for the tibiofemoral FE tasks. The influence of rotation speeds on the pose errors is shown in Figure 6, in which the maximum differences of the moving averaged absolute errors across all tested speeds were 0.16, 0.23, and 0.22 mm for the three translational components and 0.26, 0.22, and 0.26° for the rotational components. The moving averaged absolute errors were below 0.36 mm in translation and 0.44° in rotation (Figure 6).

When comparing the mTRE obtained from the proposed method with and without incorporating template registration, a substantial difference was found (Figure 7). With template registration, the MIBFT-determined pose parameters yielded final mTRE values of generally less than 2 mm (Figure 7A). When template registration was not utilized, a proportion of registrations unsuccessfully localized the bone in 3D space, as indicated by instances with final mTRE values larger than 6 mm (Figure 7A). It also led to gradually increased averaged final mTRE values with respect to initial mTRE values, as shown in Figure 7B.

The agreement of the MIBFT-determined tibiofemoral 6-dof kinematics with those of the standard reference values is shown in Figure 8. The MIBFT reproduced tibiofemoral kinematics that were in good agreement with the standard reference kinematics. The Bland–Altman plot gave mean differences for all translation components ranging from −0.17 to 0.22 mm and an LoA in the range of −0.67 to 0.88 mm. For flexion/extension, abduction/adduction, and internal/external rotation, the mean differences ranged from −0.04 to 0.38°, and the LoAs were in the range of −0.85 to 1.6°.

Based on the same experimental dataset, the MIBFT appeared to outperform typical model-based tracking using a single-plane image [24] and the previously reported MCPI method [40] in terms of the Rms errors of skeletal poses (Figure 9). As indicated by the boxplot, the median, 25 and 75% percentile of the Rms error distribution for each motion component were lowest for the MIBFT method (Figure 9). The present study demonstrated bias, precision, and Rms errors that were smaller or within the ranges of accuracy of those from model-based tracking techniques with synchronous or clinically asynchronous biplane X-ray imaging previously reported in the literature (

Table 2. Comparison of reported accuracies in the measurement of tibiofemoral kinematics obtained from biplane model-based tracking in the literature and the present study. The accuracy values are represented in the form of bias (precision) [Rms error].

Study	Equipment e	Activity	Add/Abd (°) (X-axis)	I/E Rot (°) (Y-axis)	Flex/Ext (°) (Z-axis)	AP Tran. (mm) (X-axis)	PD Tran. (mm) (Y-axis)	LM Tran. (mm) (Z-axis)
Li et al. [35] a	Synchronous BXI	Dynamic	0.31 (0.72)	−0.16 (0.61)	0.37 (0.91)	0.24 (0.16)	−0.11 (0.18)	−0.13 (0.18)
Anderst et al. [30]	Synchronous BXI	Dynamic	−0.11 (0.30) [0.54]	1.01 (0.62) [1.44]	−0.3 (0.94) [1.75]	−0.68 (0.74) [1.54]	−0.16 (0.37) [0.69]	−0.49 (0.31) [0.69]
Ohnishi et al. [60] b	Clinical BXI	Static	[0.34]	[0.55]	[0.62]	[0.51]	[0.49]	[0.53]
Giphart et al. [31] c	Synchronous BXI	Dynamic	−0.03 (0.61)	−0.05 (0.69)	0.03 (0.43)	−0.02 (0.49)	0.08 (0.60)	0.27 (0.71)
Stentz-Olesen et al. [32] d	Synchronous BXI	Dynamic	0.05 (0.71)	−0.17 (0.77)	−0.11 (0.45)	−0.05 (0.56)	−0.23 (0.38)	0.16 (0.68)
Present study	Clinical BXI	Dynamic	−0.04 (0.26) [0.26]	0.38 (0.63) [0.73]	0.02 (0.32) [0.32]	0.14 (0.24) [0.28]	−0.17 (0.26) [0.31]	0.22 (0.34) [0.40]

^a Data were extracted from Knee 1. ^b Details regarding the definition of the joint axes were not provided in the study. ^c Data were extracted from an isolated extension trial. ^d Precision values were taken as the standard deviation of errors estimated from the Bland–Altman plots. ^e Abbreviation “BXI” indicates the biplane X-ray imaging system.

).

4. Discussion

A new model-based tracking scheme was developed to accommodate interleaved biplane fluoroscopy for measuring tibiofemoral kinematics. It was primarily accomplished by synthesizing an intermediate image with a frame interpolation approach to integrate asynchronous information provided by the clinical biplane fluoroscopy to yield a pseudosynchronous pair of biplane images. In previously reported techniques with the same objective, the approaches have normally relied on the data generated from the interpolation of 3D skeletal kinematics [39,40], which must be obtained from an additional image registration procedure based on a single-plane fluoroscopic image. This therefore led to two image analysis steps for an entire kinematics reconstruction process. Successful single-plane image registration must be ensured to generate precise kinematic interpolation, which subsequently affects the accuracy of the measured bone poses. In contrast, the present approach requires only a single image analysis step. In combination with template registration for initial pose estimation, fully automated image registration for a complete image set was achieved.

Directly comparing the present accuracy values to those reported with different biplane fluoroscopy-based model-based tracking approaches was difficult, as the imaging system used, experimental conditions of the motion tasks, and the reported error metrics were different (Table 2). Nonetheless, a comparison of the primary error metrics calculated for the determination of tibiofemoral kinematics, such as bias, precision, and Rms error (Table 2), was made. It appeared that the present MIBFT with a clinical biplane X-ray imaging system was capable of producing submillimeter and subdegree accuracy in terms of the compared error metrics and was well within the values generated by techniques based on custom-made imaging systems, which supports the applicability of the clinical asynchronous X-ray imaging system for precise functional examinations of the knee. On the same dataset, we further implemented and evaluated a single-plane MBT [24] and the relevant MCPI algorithm specifically devised for clinical asynchronous biplane fluoroscopy [40]. The present MIBFT was shown to not only be capable of addressing the out-of-plane translation problem in single-plane MBT and reducing errors in out-of-plane rotation, as shown in MCPI [40], but further reduce errors in other relatively accurate motion components (Figure 9).

Among the 6-dof reconstructed tibiofemoral kinematics, the internal/external rotation was less accurate than the other motion components (Table 2) and was the only one that showed an LoA beyond 1° (Figure 8). Similar results were also found in previous studies [31,32] (Table 2) conducted on the synchronous biplane imaging system. This is reasonable, as the nature of the cylindrical shape of both the femur and tibia leads to a less sensitive, and thus less accurate, model rotation around the long axis during image registration [32]. On the other hand, the specimen leg was hung on a frame in the experimental setup used in the present study, which caused the long axes of the bones to be most near-parallel to the y-axis of the fluoroscopic coordinate system (Figure 1C), even during motion tasks. The less accurate y-axis rotation of the bone (Table 1), as a result of being an out-of-plane motion component from the perspective of both detectors [40], would also propagate the errors to the long-axis rotation of the bones and the description of the tibiofemoral internal/external rotation.

Template registration on video fluoroscopy has been applied for the analysis of the planar displacement of vertebrae during motion in three-parameter space [13,61,62,63]. While the pose parameters may be accurately reproduced through a semirigid template [61] or with an exhaustive search strategy to ensure a global optimum [13,63], increased computational effort and run time are inevitable. Since accurate estimation of template poses was not the major aim of the study, a quick and simple approximation of template displacement between successive frames using a particle filter was shown to be a valid alternative. Integrating the planar displacements from two perspectives, the estimated 3D bone displacement appeared to greatly compensate for the deviations in bone poses between successive image frames (Figure 5). This diminishes the search space of pose parameters required in the subsequent 3D/2D image registration and therefore increases the probability of finding the optimal solution, allowing for measurement of faster motions involving a greater between-frame bone displacement (Figure 7).

The influences of the motion speeds on the determination of bone poses were analyzed in the present study, as the motion speed was found to affect the quantification of joint motion via fluoroscopy [26]. The average and instantaneous maximum tibiofemoral FE speeds were approximately 158 and 250°/s, respectively (Figure 6), which were higher than the values reported in [31] and [64] (i.e., 45–60°/s) and comparable to those reported in [65] and [26] (i.e., 202–225°/s). In addition to the factors related to the equipment used (e.g., fluoroscopy), the quality of the synthetic images obtained from the frame interpolation should be considered, especially when the motion activity is fast. A higher motion speed could result in a more pronounced object displacement and nonrigid deformation between successive images, which would directly increase the required search space for the motion vectors and probably deteriorate blocking artifacts caused by inaccurate optical flow. Recently reported deep-learning-based approaches (e.g., [66]) showed great potential in reproducing high-quality interpolation images with minimal artifacts; such approaches can be feasibly integrated in the present MIBFT and may be taken into account in future applications. Nonetheless, the current implementation of a block-based frame interpolation algorithm yielded acceptable interpolation images for biplane image registration, and their influence on the quantification of high-speed kinematics was limited to approximately 0.3 mm and 0.3° (Figure 6).

5. Conclusions

A new fully automated model-based tracking scheme in conjunction with clinical asynchronous biplane fluoroscopy has been developed and evaluated for the accuracy on the measurement of dynamic tibiofemoral motions. Less than 0.4 mm Rms error in measuring tibiofemoral translation and 0.8° Rms error in measuring rotation were achieved, comparable to those of a variety of model-based tracking techniques using custom-made synchronous biplane imaging systems. This may encourage the further use of clinical imaging systems for the noninvasive and precise examination of dynamic joint functions in clinical practice and extend the application of MBT for kinematic measurement investigations in a budge-efficient way. The presented approach may be also integrated into image-guided intervention procedures to spatially align the pre-interventional image data with the intra-intervention data.