Integration of statistical shape modeling and alternating interpolation-based model tracking technique for measuring knee kinematics in vivo using clinical interleaved bi-plane fluoroscopy

Subjects

Twelve healthy male volunteers (age: 30.3 ± 17.1 years; height: 172.4 ± 7.0 cm; body mass: 64.5 ± 9.5 kg) without any neuromusculoskeletal disease or surgical history of the lower limbs participated in the current study. The subjects were fully informed of the experimental protocol approved by the Institutional Review Board and gave their written consent (China Medical University Hospital Research Ethics Committee, No: CMUH107-REC2-078). All subjects were scanned by CT at a voxel size of 0.709 mm ×0. 709 mm ×0.625 mm (Optima CT660; GE Healthcare, Chicago, IL, USA) for the reconstruction of subject-specific volumetric models of the knee, which were not involved in building the knee SSM reported by (Lu et al. (2021)).

Experimental protocol

The subjects performed active flexion and extension of the left knee with the right leg supporting the body weight. In addition to the tested tasks, all subjects performed static standing for subject calibration, from which the knee pose was obtained as the baseline for subsequent analysis. A clinical bi-plane fluoroscopy system (Allura XPER FD 20/20; Philips Medical Systems, Best, Netherlands) was used to acquire the interleaved bi-plane X-ray images of the knee during tested tasks at a resolution of 512 × 512, a grayscale of 8-bit, and an effective frame rate of 60 fps with a constant time interval of 1/60 s between two X-ray units. The two fluoroscopic detector units were positioned orthogonally with each other. Each fluoroscopic unit was modeled as an ideal perspective projection of a point source of X-ray. Before data acquisition, a well-established experimental calibration procedure was performed to obtain the intrinsic and extrinsic parameters for each projection model of the fluoroscopy units (Lin et al., 2014b; Lin et al., 2018). Four lead skin markers were attached to the distal thigh and proximal shank, and one lead marker was attached to the patella to determine the spatial transformation between the X-ray image pairs. For each pair of fluoroscopy images, the 3D coordinates of the lead markers were first determined using radio-stereometric analysis (Karrholm, 1989). The obtained coordinates of the lead markers from any two image pairs were then co-registered to determine the transformations between the two image pairs. Each subject stood with the tested foot on a rotating plate and the left knee located at the isocenter of the bi-plane imaging system during subject calibration. Multiple views of fluoroscopic images were acquired by rotating the tested lower limb vertically around the isocenter of the bi-plane fluoroscope for 0°, 30°, 45°, and 60°. From these positions, three combinations of asynchronous X-ray image pairs were used for the subsequent reconstruction of personalized knee model: (1) one image pair (two orthogonal images from the 0° position); (2) two image pairs (four images from the 0° and 45° positions); (3) three image pairs (six images from the 0°, 30°, and 60° positions).

Statistical shape modeling of the knee

The general procedure of the knee SSM included (1) obtaining a set of CT-derived training shape models, (2) choosing a reference model with a predefined surface mesh; (3) establishing shape (mesh) correspondence between individual training models by transforming the reference model to individual training ones; and (4) determining the mean model and primary modes of shape variations using Principal Component Analysis (PCA). The training shape models for the SSM were reconstructed from the CT data of the distal femur and the proximal tibia from 60 healthy Chinese males (Lu et al., 2021). The mesh of the femoral reference model had 3524 vertices, and that of the tibial reference model had 1901 vertices. Individual training models were obtained by co-registering the reference and CT-based models via shape correspondence using the iterative closest point (ICP) (Besl & McKay, 1992) and coherent point drift (CPD) (Myronenko & Song, 2010) methods, and shape alignment using Generalized Procrustes analysis (GPA) (Cootes & Taylor, 2004). The variations of the individual shape models were then decomposed to a set of principal components (eigenvectors) using principal component analysis (PCA) (Wold, Esbensen & Geladi, 1987). Thus, each training shape model could be described as the mean model superimposed with a linear combination of the principal components.

For each of the current 12 knee joints, the subject-specific bone shape model was generated using the trained SSM with multiple asynchronous 2D images of the bone via a two-phase optimization scheme, referred to as SSM-reconstructed model (Lu et al., 2021). In the two-phase optimization scheme, the first phase involved searching for the optimum pose (i.e., six degrees of freedom of the bone) and shape (i.e., parameters for the first 10 principal components) that maximizes the similarity between the pseudo digitally reconstructed radiographs (DRRs) of the SSM-reconstructed bone model and the multiple asynchronous 2D fluoroscopy images (Lu et al., 2021). Generally, a DRR was generated by casting rays from the X-ray point source towards the image plane through the CT-reconstructed model, from which the interior voxel values (i.e., the ideal attenuation coefficients) of the bone encountered by each of these rays were accumulated to derive the grayscale values of the DRR (Engel et al., 2006). For generating the DRRs of an SSM-reconstructed bone model, the model was transformed to a pseudo volumetric one by a voxelization process, which was achieved by first finding intersections between the shape model and a set of parallel virtual 2D transverse slices (pixel size of 0.5 mm ×0.5 mm and inter-slice distance of 0.5 mm) (Patil & Ravi, 2005). The resulting virtual voxels interior to the shape model were assigned a constant of 700 to simulate the Hounsfield unit (HU) value of bone, while voxels outside the contours were assigned −1.000 to simulate air. As a result, the DRRs of the pseudo volumetric SSM-reconstructed models in space were generated for a given bone pose and shape. In the second phase, the shape model and its pose obtained in the first phase were further refined. The shape model was further refined by including ten additional shape parameters (i.e., 11th to 20th principal components) to better match the fluoroscopic images taking the pose and shape parameters obtained in the first phase as fixed parameters, giving the final shape of the SSM-reconstructed bone model (Lu et al., 2021). The genetic algorithm (Goldberg & Holland, 1988) was employed to search for the optimal poses and shape coefficients parameters. By partitioning the shape and pose parameters, the two-phase optimization approach enabled the consideration of more principal components with increased accuracy but reduced computational effort compared to a single-phase optimization with 10 principal components (Lu et al., 2021).

3D/2D registration using AIMT with SSM-reconstructed and CT-based models

We utilized a validated AIMT approach (Lin et al., 2020) to accomplish the fully automated kinematics tracking of the knee, from which the SSM-reconstructed and CT-reconstructed subject-specific bone models were separately registered to the interleaved bi-plane fluoroscopic images to obtain the six-degree-of-freedom pose parameters (Fig. 2). The AIMT technique consists of three stages, namely (1) 2D/2D template registration, (2) motion-compensated frame interpolation, and (3) bi-plane 3D/2D image registration. The first stage analysis aimed to estimate the bone model’s pose parameters for the current frame from the registered pose in the preceding frame. This was accomplished by estimating the 2D pose increment between the preceding and current image frames using template registration with the particle filter (Arulampalam et al., 2002). The 2D pose increments from the two fluoroscopic views were then summed up to obtain the 3D pose increment used to estimate the pose parameters of the bone model in the current frame (Lin et al., 2020). In the second stage, the interleaved bi-plane images were converted into pseudo synchronous bi-plane image pairs. To this end, an intermediate image frame (Iⁿ) between any two consecutive image frames I ⁿ⁻¹ and I ⁿ⁺¹ in the same fluoroscopic view was synthesized using a motion-compensated frame interpolation method (Zhai et al., 2005) (Figs. 2 and 3. In the final stage, a bi-plane model-based 3D/2D image registration scheme based on the forward projection model was implemented to precisely determine the pose parameters of the bones (Fig. 4). The 3D pose parameters of the bone were determined by maximizing the normalized cross-correlation between gradients (i.e., so-called gradient correlation) (Penney et al., 1998). Optimization procedure for 3D/2D image registration was performed to minimize a metric of gradient correlation (GC, f_GC) that quantified the similarity of pixel intensity between the model-projection DRRs (I_DRR) and the pseudo-bi-plane fluoroscopic images (I_Fluoro). DRR image (I_DRR) and fluoroscopic image (I_Fluoro) are transformed by applying horizontal (i) and vertical Sobel templates (j) to obtain four gradient images, namely dI_DRR/di, dI_DRR/dj, dI_fl/di, dI_fl/dj. GC is defined with normalized cross correlation between I_Fluoro and I_DRR as the follows Eq. (1). The sum of −f_GC was taken as the cost function to be minimized by using N number of alternating images (N = 2) with design variables of six degrees of freedom (x, y, z, α, β, γ) Eq. (2). The genetic algorithm was employed to search for the optimal poses (Goldberg & Holland, 1988). Noted that the pose parameters reproduced using the CT-derived bone models were taken as the standard reference for evaluating the proposed method. ${\displaystyle f_{GC}=\frac{\sum _i\left[d{I}_{fl}/di-{\overline{d{I}_{fl}/di}}_{Fluoro}\right]\sum _i\left[d{I}_{DRR}/di-{\overline{d{I}_{DRR}/di}}_{DRR}\right]}{\sqrt{\sum _i{\left[d{I}_{fl}/di-{\overline{d{I}_{fl}/di}}_{Fluoro}\right]}^2}\sqrt{\sum _i{\left[d{I}_{DRR}/di-{\overline{d{I}_{DRR}/di}}_{DRR}\right]}^2}}}$ (1)${\displaystyle +\frac{\sum _j\left[d{I}_{fl}/dj-{\overline{d{I}_{fl}/dj}}_{Fluoro}\right]\sum _j\left[d{I}_{DRR}/dj-{\overline{d{I}_{DRR}/dj}}_{DRR}\right]}{\sqrt{\sum _j{\left[d{I}_{fl}/dj-{\overline{d{I}_{fl}/dj}}_{Fluoro}\right]}^2}\sqrt{\sum _j{\left[d{I}_{DRR}/dj-{\overline{d{I}_{DRR}/dj}}_{DRR}\right]}^2}}}$ where ${\overline{d{I}_{fl}/di}}_{Fluoro}$, ${\overline{d{I}_{DRR}/di}}_{DRR}$, ${\overline{d{I}_{fl}/dj}}_{Fluoro}$, ${\overline{d{I}_{DRR}/dj}}_{DRR}$ are the mean values of the transformed images. (2)${\displaystyle f_{\mathit{GC},\mathit{all}}=arg\underset{f_{\mathit{GC}}}{min}\sum _{p=1}^N{\left(x,y,z,\alpha ,\beta ,\gamma \right)}_p}$

From the registered poses of the bones, the knee kinematics were then calculated. The anatomical coordinate systems (ACS) of the CT-based models of the femur and tibia were determined by the 3D geometry features of the bone models (Miranda et al., 2010), with the positive x-axis directed anteriorly, positive y-axis directed superiorly and positive z-axis directed to the right. The CT-based model was co-registered with the corresponding SSM-reconstructed bone model via the ICP method, and the ACS of the CT-based model was then assigned to the registered SSM-reconstructed model. This approach minimized the possible discrepancies in the coordinate systems between SSM- and CT-based models. The knee kinematics was described as the tibial pose in the femoral ACS. The joint angles were computed using a z − x − y Cardanic sequence (Grood & Suntay, 1983), giving flexion/extension, abduction/adduction and external/internal rotations.

Evaluation of the AIMT technique with SSM-reconstructed bone model

For the quantification of the shape differences between the SSM-reconstructed model and the corresponding CT-reconstructed model, the models were first spatially aligned to each other via the ICP method (Besl & McKay, 1992), and the point-to-surface distance 𝔢 for each point p on the surface of the SSM-reconstructed model was then calculated as its shortest Euclidean distance to a corresponding point p′ on the surface S of the CT-reconstructed bone model as follows (Cignoni, Rocchini & Scopigno, 1998). (3)${\displaystyle e\left(\mathrm{p},\mathrm{S}\right)=\left|\underset{p\prime ϵS}{min}d\left(p,p\prime \right)\right|}$

The root-mean-squared values of 𝔢 over the entire surface points (RMSe) were obtained as a measure of the shape error of the SSM-reconstructed model.

For the quantification of the differences between the registered bone kinematics during the tested dynamic activity using SSM-reconstructed and CT-reconstructed models, the mean target registration error (mTRE) (Van de Kraats et al., 2005) was calculated as the mean of 𝔢 values over the entire SSM-reconstructed model surface for each image frame. The mean of mTRE (mmTRE) was then obtained by averaging the mTRE of each frame over the dynamic motion cycle. The peak mTRE (pmTRE) over the dynamic motion cycle was also obtained. In this study, a successful registration was defined as mTRE less than 1.5 mm.

The differences between the registered joint poses using the SSM-reconstructed and CT-reconstructed models were quantified by the mean absolute differences (MAD) in each of the six kinematic components (i.e., three translations and three rotations) over the dynamic motion cycle. The peak absolute differences (PAD) over the dynamic motion cycle was also obtained for each kinematic component. The mean and standard deviation (SD) of the MAD across all the subjects gave the bias and precision of the model used, respectively.

Computational costs (time) in reconstructing SSM-based models from 1, 2 and 3 static image pairs and registering these SSM-based models to interleaved biplane fluoroscopic images over the tested dynamic activity were also obtained on a 3.7 GHz, 32 GB computer running a MATLAB implementation.

Statistical analysis

One-way ANOVA was used to test the differences in mmTRE, MAD in six-degree-of-freedom components, computational costs in reconstructing SSM-based models from 1, 2 and 3 image pairs and registering SSM-based models to fluoroscopic images with paired t-tests by using Bonferroni correction for post hoc pairwise comparisons. All comparisons were performed using SPSS 20.0 (SPSS, IBM, Armonk, New York, NY, USA). The significance level was set at α = 0.05.