A simulation study to investigate an extension to the point cluster technique

Joint kinematics are an important and widely utilized metric in quantitative human movement analysis. Typically, trajectory data for skin-mounted markers are collected using stereophotogrammetry, sometimes referred to as optical motion capture, and processed using various mathematical models to estimate joint kinematics (e.g., angles). Among the various sources of noise in optical motion capture data, soft tissue artifacts (STAs) remain a critical source of error. This study investigates the performance of the point cluster technique (PCT), an extension of the PCT using perturbation theory (PCT-PT), and singular value decomposition least squares (SVD-LS) method (as a reference) for 100 different marker configurations on the thigh and shank during treadmill walking. This study provides additional evidence that the PCT method is significantly limited by the underlying mathematical constraints governing its optimization process. Furthermore, the results suggest the PCT-PT method outperforms the PCT method across all performance metrics for both body segments during the entire gait cycle. For position-based metrics, the PCT-PT method provides better estimates than the SVD-LS method for the thigh during majority of the stance phase and provides comparable estimates for the shank during the entire gait cycle. For knee angle estimates, the PCT-PT method provides equivalent results as the SVD-LS method.

family of methods 22 , some of which have been found to be unreliable solutions to mitigating STAs 23 .This method incorporates joint constraints in a cost function via a mathematical model describing a ball-and-socket joint, which means this method is somewhat limited by treating every joint as a ball-and-socket 9 .For example, Lu et al. 21utilized a sinusoidal function proposed by Cheze et al. 11 to artificially represent STAs in the evaluation of their multi-link musculoskeletal model that exploits this ball-and-socket assumption.While 21 showed minimal errors in axial rotation as well as abduction/adduction rotations at the joints, a subsequent validation study that used dual-plane fluoroscopy as ground truth revealed significant errors (i.e., mean errors of 10 degrees for rotations and 10-15 mm for translations) 24 .
Unlike the SVD-LS and global optimization methods, the PCT is grounded in physics and is free of potentially controversial joint constraint assumptions.The PCT assumes a body segment is a rigid body, which means the body segment's moment of inertia in a body-fixed frame of reference must be constant.Assuming each marker to have unit mass, the eigenvectors, eigenvalues, and center of mass (CM) of the subsequent mass distribution are calculated at each time step.Then, the mass of each marker is systematically adjusted such that the difference between norm of current time step's eigenvalues and those of a reference (i.e., those from a static calibration) is minimized.The segment pose is then given by the eigenvectors and CM of the newly adjusted mass redistribution.Note that this method is specifically designed to address non-rigid body STA components.
Importantly, details regarding how to optimize the algorithm's performance are not well understood.For example, an important assumption underlying the PCT is a uniform marker distribution, but the characteristics governing the nature of the uniform marker distribution (e.g., minimum marker separation) are presently unknown 15 .While this method is theoretically compatible with N> 3 markers, the proof-of-concept study demonstrated the validity for N = 8 markers only 15 .To address the fact that the PCT tends to fail in cases where marker clusters possess axes of rotational symmetry 15 , Carman et al. 25 carried out a more extensive investigation into marker configurations that caused the PCT to fail during movement trials.The PCT results suggested that anatomical frame axis alignment agreed very well with the SVD-LS results for the same settings.However, this study did not have ground truth data for comparison, used a similar N = 8 marker configuration 15 , and only con- sidered the shank -a body segment known to have minimal STAs.As aforementioned, Cereatti et al. 19 conducted a comparative assessment of the performance of various bone-pose estimators for a typical flexion-extension pattern during a gait cycle, which was corrupted with instrumental noise and STAs.Two marker configurations ( N = 4 and N = 12 markers) were used for both the thigh and shank.Under their experimental conditions, the PCT failed to produce any results for N = 4 markers and produced inconsistent results for N = 12 markers.The authors concluded that these findings were a result of the PCT algorithm optimizing the eigenvalue norm as opposed to individual eigenvalues themselves 19 .
While the PCT holds promise, past research also provides clear evidence of the dependence of the PCT's performance on the number and configuration of markers.Thus, the purpose of this study is twofold.First, this study describes and evaluates a novel bone pose estimation method that extends PCT by utilizing perturbation theory, optimizing individual eigenvalues, and using a mass redistribution approach that is independent of analytical constraints.Next, a comparative assessment of the traditional PCT, the SVD-LS, and the novel methods for various marker configurations is conducted.

PCT with perturbation theory
Assuming the markers have unit mass, the inertia tensor, I I I(t) , is given by Eqn. (1): where x j (t) , y j (t) , and z j (t) is the three-dimensional global position of the jth marker.Let (t) = [ 1 (t), 2 (t), 3 (t)] ⊺ denote the eigenvalues and E E E(t) = [e e e 1 (t), e e e 2 (t), e e e 3 (t)] denote the corresponding eigenvectors associated with I I I(t) .In the presence of STAs, I I I(t) and its corresponding eigendecomposition, (t) and E E E(t) , are time varying.Eigenvalue perturbation theory 26,27 can effectively "denoise" a matrix like the inertia tensor when implemented iteratively.Let I I I 0 represent the inertia tensor and 0 represent its eigenvalues at the initial time step, t 0 , which are equivalent to those calculated from a static calibration.Since the inertia tensor is necessarily symmetric, a perturbation matrix, δI I I C , can be defined using Eqn.(2): The change in the current jth eigenvalue, δ C,j , can be computed from perturbation theory 26,27 via Eqn. (3): (1) where e e e C,j is the current jth eigenvector.As detailed in "Appendix A" in the Supplementary Information, expand- ing the right hand side of Eqn.(3), multiplying the terms out, applying the result to each eigenvalue, and rearranging yields Eqn.(4): The current perturbation vector, δK K K C , can be obtained by solving Eqn.(4) by minimizing the L 2 norm least squares solution of Ax = b (i.e., lsqminnorm function in MATLAB 28 ), after which δI I I C can be calculated.The algorithm then updates the current eigenvalues, C , by minimizing the norm of the difference, � C − 0 � , subject to constraints.The current inertia tensor, I I I C , along with its eigendecomposition, C and E E E C , are initialized using the raw marker data polluted with STAs.A step size, s, is initialized to be 1, and the current eigenvalue offset, δ C , is initialized to be s • ( 0 − C ) .After determining δI I I C from the solution to Eqn. (4), the current eigenvector offset, δE E E C , is obtained utilizing Eqn. ( 5): Next, the conditions described by Eqns.( 6) and ( 7) are evaluated for each eigenvector: If these conditions are met, the estimates are updated.Otherwise, s and δ C are halved and the process is repeated until s < 2 −13 .If � 0 − C � < 10 −7 and the e e e C,j are mutually orthogonal, the final estimates are used to represent the current inertia tensor that has been appropriately denoised via perturbation theory.
Next, the task is to recover the mass distribution representing the "true" global marker positions, . Let I I I F (t) represent the denoised inertia tensor and m j represent the constant non- unit mass of the jth marker.At some time t k , the relationship between the inertia tensor, marker positions, and mass distribution, m m m = [m 1 , . . ., m N ] , is given by Eqn.(8): Eqn. (8) effectively contains 6 nonlinear equations at some time, t k , as shown in Eqn. ( 9): Stacking equations for every time step yields 6N T equations, where N T is the total number of time steps.The Levenberg-Marquardt algorithm 29 can solve this system of nonlinear equations to obtain the non-unit mass distribution, m m m , which is then used to obtain the unconstrained center of mass (CM), T T T UC (t) .The CM is then constrained with a scaling factor α such that the distance between the CM and the centroid (CM for unit mass distribution) is bound for each time instant (see "Appendix B" in the Supplementary Information).The constrained CM is given by Eqn. ( 10): where c c c(t) represents the centroid.Let �T T T(t) represent the offset of the constrained CM with respect to the centroid projected onto the first denoised eigenvector, e e e F,1 , given by Eqn. ( 11): If �T T T(t) < 0 for the entire gait cycle, the constrained CM is treated as the final CM estimate, T T T(t) .However, if �T T T(t) > 0 for the entire gait cycle, then the constrained CM is reflected about the centroid to obtain the final CM estimate, T T T(t) .If �T T T(t) changes sign, the roots of �T T T(t) must be examined (see "Appendix C" in the (4)    . . .

Simulation study
The data utilized in this study was from an extensive open-source repository focused on providing researchers with the means to characterize STAs for a variety of activities and participants 31 .A conical frustum was created using a previously described anthropometric model 32 (Fig. 1).Virtual marker configurations were generated based on locations for which STA data are available for a particular body segment 31,33 .Consistent with common marker sets 34,35 , the marker areas were divided into quadrants for which one virtual marker was placed randomly within each quadrant.A marker configuration was rejected if the markers were not adequately spaced.The result is a uniform distribution for N = 4 markers that is sufficient for estimating a body segment's pose 36 .Following the virtual marker placement, anatomical landmarks are also generated for the thigh (medial epicondyle (ME), lateral epicondyle (LE), femoral head (FH), and greater trochanter (GT)) and shank (medial malleolus (MM), lateral malleolus (LM), fibular head (HF), and tibial tuberosity (TT)).A total of N c = 100 configurations were generated.
Fortunately, STA and corresponding anatomical frame pose data has been made available for a variety of activities and subjects 31 .Here, treadmill walking data for 3 full gait cycles of the left leg for a 75 year old male participant was utilized 33 .Ground truth data are provided by two fluoroscopes configured to acquire X-ray images at a maximal frame rate 33 .First, the markers used in that data collection were assigned to a quadrant consistent with Fig. 1.For each quadrant, the STA data was obtained by averaging the raw STA data over all markers.Gait cycles were extracted from heel strikes, which was then used to temporally normalize both the ground truth and optical motion capture data.The ground truth data (collected at 30 Hz) was then interpolated such that the time stamps match that of the optical motion capture data (collected at 240 Hz).Finally, the STA data for the markers in a given quadrant on a given body segment were averaged for a characteristic STA profile for any virtual marker placed in that quadrant.The raw anatomical frame data provided by 33 were also pre-processed to obtain a kinematic input to drive the model.

Figure 1.
Conical frustums for the thigh (left) and the shank (right) with their respective reference marker configurations.Black spheres denote the anatomical landmarks (ME, LE, FH, GT, MM, LM, HF, TT).Red spheres denote the virtual skin-mounted markers, which can only be placed in the green shaded areas.The X (red), Y (green), Z (blue) vectors correspond to the anterior-posterior, medial-lateral, and inferior-superior directions, respectively.The inter-marker angular distance illustrated here is 60 • and height difference is greater than 15 cm.The parameters ( r 1 , r 2 , h, θ ) for the thigh are (5.5 cm, 9.5 cm, 25 cm, 60 • ) and for the shank are (4 cm, 5.5 cm, 17.5 cm, 60 • ).www.nature.com/scientificreports/

Evaluation metrics
The quadrant specific STA noise models were applied to each marker configuration (denoted by R) to generate the noisy marker configuration data.The PCT, SVD-LS, and novel PCT-PT methods were used to obtain the reconstructed configurations (denoted by r).The reconstruction error for the k th configuration using the method, M, for the segment, S, was evaluated by the center of mass reconstruction offset (TRO), eigenvector reconstruction offset (eRO), and anatomical landmark reconstruction offset (ALRO), which are computed via Eqns.( 12)-( 14), respectively: where a p (t) represents the global position vector for the p th anatomical landmark.For eRO, results will be pre- sented for each eigenvector.Descriptive statistics (i.e., the mean and standard deviation) were computed for all proposed performance metrics.For the sake of completeness, the reconstructed anatomical frames 37 were further evaluated with two additional performance metrics -the anatomical frame origin offset (AFOO) and the knee angle 38 offset (KAO).The AFOO is a nonlinear combination of the errors of TRO and eRO for an individual body segment whereas the KAO is a nonlinear combination of the errors of TRO and eRO for both body segments.

Results
Consistent with the findings of previous studies 19,25 , the PCT method produced discontinuous results for many of the generated virtual marker configurations.For the interested reader, "Appendix D" in the Supplementary Information describes how these discontinuities in the PCT results are a direct consequence of how the method was developed mathematically.Thus, only the subset of marker configurations that produced continuous time series were used to calculate descriptive statistics for the PCT method to capture what the method is capable of when physiologically consistent results are produced.Specifically, the number of continuous configurations for the PCT method were 23 and 17 for the thigh and shank, respectively.For the results to follow, the average was taken over all marker configurations for the SVD-LS and PCT-PT methods while the average was taken over the marker configurations that produced valid results for the PCT method.
The SVD-LS center of mass reconstruction offset (TRO) results did not have an envelope for either body segment.This finding suggests that the TRO for the SVD-LS method is independent of which marker configuration was used.As is illustrated in Fig. 2, the envelope for the PCT-PT method fell below the SVD-LS curve for the thigh and was comparable to SVD-LS curve for the shank during majority of the stance phase.For the majority of the gait cycle, the PCT-PT method also fell below the average curve for the PCT method for both body segments.Unsurprisingly given the inconsistency in the results, the PCT results had the largest and least smooth envelope for both the thigh and shank.
Next, the SVD-LS and PCT-PT eigenvector reconstruction offset (eRO) results were identical.The eRO results for these two methods were the smallest for the first eigenvector for the thigh and the smallest for the third eigenvector for the shank as is illustrated in Fig. 3.Note that the offsets for the other two eigenvectors for both body segments were also small (i.e., less than 5 degrees).The average eigenvector offset curve for the PCT-PT method is below the average eigenvector offset curve for the PCT for both body segments.The envelopes for the first and second eigenvectors also have minimal overlap during the majority of the stance phase for the thigh and during the majority of the gait cycle for the shank.Finally, the anatomical landmark reconstruction offset (ALRO) for the SVD-LS method had an envelope that differed from the TRO results for the thigh (Fig. 4) and shank (Fig. 5).The average curve for the PCT-PT method fell below the envelope of the SVD-LS method during the majority of the stance phase for the thigh and was comparable to the envelope of the SVD-LS method for the shank.For both body segments, the envelope for the PCT-PT method was below the average curve for the PCT during majority of the gait cycle.For the thigh (shank), the envelopes for all three methods overlap less (more) frequently for the epicondyles (malleoli) than the other two anatomical landmarks.The results for all performance metrics are reported in Tables 1 and 2 for the thigh and shank, respectively.
For the anatomical frame origin offset (AFOO), the envelope for the PCT-PT method was below the envelope for the SVD-LS method during majority of the stance phase for the thigh and during the majority of the entire gait cycle for the shank (Fig. 6).The envelope for the PCT-PT method was also below the envelope for the PCT  method for the majority of the gait cycle.The PCT method had the largest envelope while the average curve had the largest offset for both body segments.
For all knee joint angles, the reconstructed values for the PCT-PT and SVD-LS methods applied to virtual marker trajectories with representative STA data agree well with the ground truth curve for majority of the gait cycle (Fig. 7).For flexion-extension, both methods have excellent agreement for the entire gait cycle.By contrast, the PCT method overestimates flexion-extension, deviates more from the ground truth curve for internalexternal rotation, and exhibits larger envelopes across all three angles.

Discussion
This study provides additional evidence that the PCT method is significantly limited by the underlying mathematical constraints governing its optimization process 15,19,25 .As a direct result of these constraints 15 , this method frequently yields physiologically meaningless results (i.e., discontinuities), specifically following the mass redistribution optimization.As such, the performance of the PCT method was shown to vary significantly depending on the marker configuration used, which is consistent with previous studies 19,25 .Recall the PCT method yielded continuous results for only 23 marker configurations for the thigh and 17 marker configurations for the shank.www.nature.com/scientificreports/ The novel PCT-PT method significantly improves upon the PCT method by producing physiologically meaningful results for all marker configurations that are consistently superior on average with smaller envelopes.Surprisingly, the SVD-LS center of mass offset (TRO) results did not have an envelope (i.e., standard deviations equal to zero in Tables 1 and 2), implying that the TRO for this method is independent of the marker configuration.It is likely that this finding is a consequence for using the same quadrant-specific representative STA data for the virtual markers, which is also likely artificially inflating the performance of the algorithms presented here as compared to results reported in the literature 33,39,40 .However, this lack of variance was not true for the SVD-LS eigenvector reconstruction offset (eRO) results.Interestingly, the eigenvector estimates for the PCT-PT and SVD-LS methods were identical.The eigenvectors for the PCT-PT method are obtained by solving a least-squares fitting problem after translating the global frame to the center of mass 30 .The SVD-LS method 16 essentially solves a least-squares problem to obtain the centroid (center of mass with unit mass distribution) and the eigenvectors simultaneously.The fact that the PCT-PT and SVD-LS methods yield the same rotation matrix estimates implies that the least-squares solution for the rotation matrix is also independent of the center of mass estimates.www.nature.com/scientificreports/frustrum models utilized by the PCT-PT method.With respect to marker configurations, other future work could investigate how the PCT-PT method performs on widely used marker configurations like the Helen-Hayes marker set 34 and Vicon plug-in gait model 35 .It should also be noted that SVD-LS, PCT, and PCT-PT are all methods designed to eliminate the non-rigid body components of STA.Previous work 10 has shown that these methods do not improve upon each other (i.e., applying multiple non-rigid body methods does not mitigate any rigid body STA).As such, future work should also consider how to mitigate rigid body components of STA, which have been shown to be related to joint kinematics 17 .

Conclusion
The novel PCT-PT method has been shown to have superior performance compared to the PCT method from which it builds.This study also provides a detailed description regarding why the PCT method so frequently produces discontinuous (i.e., physiologically impossible) results.For position-based metrics, it has also been shown to produce better reconstruction estimates than the SVD-LS method for the thigh during majority of the stance phase and produce comparable estimates for the shank during the entire gait cycle during treadmill walking.It has also been shown to produce the same eigenvectors and consequently, knee angle estimates as the SVD-LS method.Future work should consider application of this method for non-uniform marker placements, collection of varied STA data and application of the method developed here to that data.

( 3 )
δ C,j = e e e C,j T δI I I C e e e C,j

Figure 2 .
Figure 2. Center of Mass reconstruction offset (TRO) for the (a) thigh and (b) shank over the entire gait cycle.The average TRO line (solid) and standard deviation envelope (shaded) is shown for the PCT (red), PCT-PT (blue), and SVD-LS (green).The shaded window represents the stance phase where the vertical black line bounding the shaded window on the right corresponds to toe-off.

Figure 3 .
Figure 3. Eigenvector Reconstruction Offset (eRO) for the (a) thigh and (b) shank over the entire gait cycle.The average eRO line (solid) and standard deviation envelope (shaded) is shown for the PCT (red), PCT-PT (blue), and SVD-LS (green) for the first eigenvector (top), second eigenvector (center), and third eigenvector (bottom).The shaded window represents the stance phase where the vertical black line bounding the shaded window on the right corresponds to toe-off.

Figure 4 .
Figure 4. Anatomical Landmark Reconstruction Offset (ALRO) for the thigh over the entire gait cycle.The subplots illustrate the results of the four anatomical markers corresponding to the thigh.The average CPRO line (solid) and standard deviation envelope (shaded) is shown for the PCT (red), PCT-PT (blue), and SVD-LS (green).The shaded window represents the stance phase where the vertical black line bounding the shaded window on the right corresponds to toe-off.

Figure 5 .
Figure 5. Anatomical Landmark Reconstruction Offset (ALRO) for the shank over the entire gait cycle.The subplots illustrate the results of the four anatomical markers corresponding to the shank.The average line (solid) and standard deviation envelope (shaded) is shown for the PCT (red), PCT-PT (blue), and SVD-LS (green).The shaded window represents the stance phase where the vertical black line bounding the shaded window on the right corresponds to toe-off.

Figure 6 .
Figure 6.Anatomical Frame origin offset (AFOO) over the entire gait cycle.The average line (solid) and standard deviation envelope (shaded) is shown for the PCT (red), PCT-PT (blue), and SVD-LS (green) for the thigh (top) and shank (bottom).The shaded window represents the stance phase where the vertical black line bounding the shaded window on the right corresponds to toe-off.

Table 1 .
Performance metrics for each algorithm for the thigh.The mean and standard deviations are calculated across all marker configurations for the full gait cycle (GC), stance phase (StP), and swing phase (SwP).

Table 2 .
Performance metrics for each algorithm for the shank.The mean and standard deviations are calculated across all marker configurations for the full gait cycle (GC), stance phase (StP), and swing phase (SwP).