Comparison between Helical Axis and SARA Approaches for the Estimation of Functional Joint Axes on Multi-Body Modeling Data

: Functional methods usually allow for a ﬂexible and accurate representation of joint kinematics and are increasingly implemented both for clinical and biomechanics research purposes. This paper presents a quantitative comparison between two widely adopted methods for functional axis estimation, that is, the helical axis theory and the symmetrical axis of rotation approach (SARA). To this purpose, a multi-body model was developed to simulate the lower limb of a subject. This model was designed to reproduce different motion patterns, that is, by selecting the active degrees of freedom of the simulated ankle joint. Thanks to virtual markers attached to each segment, the multi-body model was used to generate simulated motion capture data that were then analyzed by instantaneous helical axes and SARA algorithms. To achieve a synthetic representation of joint kinematics, a mean helical axis and an average SARA functional axis were estimated, along with dispersion parameters and rms distance data that were used to quantitatively assess the performance of each method. The sensitivity of each algorithm to different combinations of range and speed of motion, scattering of marker clusters, sampling rate, and additive noise on markers’ trajectories, was ﬁnally evaluated.


Introduction
Human motion analysis is particularly relevant nowadays both for clinical purposes and for biomechanics research. For instance, loading conditions in human joints can be investigated through musculoskeletal analyses, to evaluate the effects of a therapeutic or surgical treatment, or to guide the clinician in the diagnosis of postural irregularities. Although different systems are available today to perform the acquisition of the human movement, such as inertial measurement units whose applications can be increasingly observed in the recent literature [1], most of the analyses are still conducted, especially at clinical level, through motion capture (stereo-photogrammetry) systems [2]. These systems typically involve passive reflective markers that are fixed on the surface of body segments (such as the trunk or a limb) and that are tracked by multiple cameras in a controlled environment. The number and the positioning of the markers selected for each body segment have an impact on the quality of the measurement and on the following analyses [3,4]. In particular, the marker reference positions are often used to implement biomechanical models that enable the estimation of the kinematics and dynamics related to a single segment or, most likely, to a group of segments up to the whole body. Although commercial motion capture systems can already implement the most accepted and widely recognized marker positioning protocols and related biomechanical modeling [5][6][7][8], there is still interest in the development of new techniques and models for the investigation of human motion, especially when subject-specific analyses are required.

The Multi-Body Model
The core of this work is represented by the definition of a multi-body model, used as a generator for noisy marker trajectories resembling an actual motion capture acquisition of a gait analysis. The physical characteristics implemented in the model, as well as the range of motion and many other features that will be discussed, were selected to be comparable with the conditions of real biomechanical analyses to avoid misinterpretation of the results due to non-realistic joint configurations. Moreover, simulation parameters such as the sampling rate and measurement noise were defined to match the conditions of real experimental trials and to achieve a realistic representation of human motion but in a controlled environment. The model was developed in Simscape ™ (The MathWorks, Natick, MA, USA), version 5.2 (R2021b).

Definition of Shank and Foot Segments
The model consists of two rigid links, that is, the shank and the foot, whose dimensions were selected from average data. Although the geometry does not actually contribute to the calculations in the model, due to the rigid-body approximation, a prismatic solid body was associated with each segment, with the dimensions presented in Table 1. The center of mass of each link was coincident with the centroid of each prismatic solid body. The geometry was necessary to help in the selection of four points of interest for each segment that corresponds to four reflective markers located on the skin of a subject. The location of the markers was defined arbitrarily on the frontal surface of the shank and on the top surface of the foot, as shown in Figure 1a. Those positions defined a default configuration of each marker set (or cluster). In this configuration, the longitudinal distance between proximal and distal markers was set equal to half the length of the corresponding segment, whereas the lateral distance between each couple of lateral and medial markers was set equal to the total width of the corresponding segment. Medial and lateral positions refer to a right lower limb. The markers' positioning followed a similar configuration used by the authors in a previous study, in which rigid plastic shells were used to limit rigid and non-rigid cluster deformation during in vivo analyses [19].  Among each marker cluster, one of the markers was placed out of the plane includin the other three markers to account for any variability in marker positioning that normal occurs during experimentation. Then, the scattering of the four markers belonging to eac cluster was governed by a multiplicative factor (kS for the shank, kF for the foot) that w used to change the distance between each marker and the corresponding cluster centro as compared to the default configuration (kS = kF = 1). These factors (or gains) respective affected the location of each marker in the frontal plane of the shank or on the top surfa of the foot. Among each marker cluster, one of the markers was placed out of the plane including the other three markers to account for any variability in marker positioning that normally occurs during experimentation. Then, the scattering of the four markers belonging to each cluster was governed by a multiplicative factor (k S for the shank, k F for the foot) that was used to change the distance between each marker and the corresponding cluster centroid as compared to the default configuration (k S = k F = 1). These factors (or gains) respectively affected the location of each marker in the frontal plane of the shank or on the top surface of the foot.

Joint Kinematics Setup
Having defined the two segments, two joints were implemented to simulate the knee and the ankle. Although the paper focuses on ankle joint motion, the knee was introduced only to achieve a more realistic representation of markers' trajectories in the global coordinate system (GCS). Due to its minor role in the methodology here discussed, the knee joint is approximated as a simple revolute, with a single rotational DoF acting between two frames: the first one translates horizontally along the X-axis of the GCS at constant speed, set equal to 1.4 m/s (average walking velocity); the second frame is integral with the shank and located at the knee center of rotation. By model definition, the knee AoR is parallel to the Z-axis of the GCS. As stated above, the translational DoF assigned to the knee joint allowed to generate marker coordinates whose magnitude (in mm) could result to be compatible with the actual trajectories recorded during a gait analysis in a movement analysis laboratory.
As far as the ankle is concerned, the model allowed relative motion between a first (base) frame integral with the shank and located at the distal end of the segment and a second (follower) frame integral with the foot and located at the proximal end of the corresponding segment. In this case, in order to achieve a more complex and detailed description of the natural joint kinematics, a bushing joint (6 DoFs) was selected. Based on this approach, the simulated ankle joint enables three independent translations and three independent rotations along the base frame axes. Therefore, the model enables different motion configurations, for instance, by selecting only a single active DoF or multiple DoFs (up to 6).
In order to run the model and generate the markers' trajectories, the complete joint kinematics must be provided as an input; hence, the rotation angle, angular velocity, and angular accelerations of each joint had to be defined in the first place. To this purpose, the three ankle rotation angles were defined as piecewise-defined functions of the gait cycle, by considering ranges of motion and curve trends respectful of natural joint kinematics. The translational DoFs along the antero-posterior, longitudinal, and lateral directions were defined with a similar approach; however, in this case, the ranges of motion were selected arbitrarily due to the lack of specific information from the literature. However, to avoid unrealistic joint configurations, the ranges of the translational DoFs were selected to achieve AoR translations coherent with the linear dispersion of helical axes observed on cadaver specimens [27]. Each ankle DoF signal could be easily varied in terms of magnitude and duration, hence enabling the possibility to test the model for different ranges of motion and speed. Figure 1b,c show the patterns of ankle and knee motion for a general tri-planar motion configuration, which considers 6 active DoFs at the ankle joint. Unlike the ankle joint, the knee, modeled as a simple revolute joint, only requires the definition of the flexion/extension DoF to run the model. The same pattern of knee motion shown in Figure 1b was considered for all the simulations, which always included three full gait cycles. As highlighted in Figure 1b,c, all the DoFs were smoothed by a low-pass filter (Butterworth, 4th order) before they were provided as input to the model to remove the effect of undesired spikes in the velocity and acceleration signals.
To summarize, the setup of joint kinematics in the model allows for customizing: 1.
The number of active DoFs allowed by the ankle joint (modeled as a bushing); 2.
The range of motion of each active DoF; 3.
The duration of each gait cycle, hence the speed of the motion patterns.

Generation of the Markers' Coordinates
In order to collect the coordinates in the GCS, a transform sensor was included in the model for each marker. This component measures time-dependent relationships, as the position, between two frames. For each marker, the respective sensor measures the position of the marker frame relative to the GCS frame, which is fixed. In this way, the global coordinates of each marker could be assessed, similarly to what is performed during a common motion capture acquisition. A fixed step size equal to 1 ms was selected to run the model, corresponding to a native sampling frequency of marker's data equal to 1000 Hz. All the input for the model (DoFs signals) was provided at the same sampling frequency to ensure consistency. Commercial motion capture systems typically work at much lower sampling rate due to limitations in cameras' refresh rate or simply to limit the amount of data collected for long-lasting observations. Therefore, to achieve a realistic representation of the dataset, the original marker's coordinates produced by the model were downsampled. The default target frequency f s was set to 100 Hz, even though the model enables any arbitrary selection of this parameter.
Finally, to simulate the effect of measurement noise from the instrumentation, as well as the artifacts produced by skin sliding and soft tissue, the downsampled coordinates were altered by additive white Gaussian noise. The simplified approach considered was the same adopted by Ehrig et al. [26], which consists of two different sources for the noise signals that can be added to each marker trajectory. In particular, white Gaussian noise is primarily added to each coordinate to simulate random accuracy errors related to the instrumentation and to marker cluster non-rigid deformation. Then, a single white Gaussian noise signal is generated and added to all the markers included in the same cluster. This second source of (common) noise simulates those accuracy errors related to skin sliding that should affect the positioning of each marker within the cluster in the same way (rigid transformations, [15]). This approach is rather a simplified one, since soft tissue artifacts, in particular, have quite different dynamics and depend on the motor task considered [14,16]. However, this method can be considered sufficient to achieve a more realistic representation of the dataset. To define the magnitude of marker noise, the signal-to-noise ratios (SNRs) for the common (SNR cmn ) and the random (SNR rnd ) sources were included as configurable parameters in the model and can be selected arbitrarily. By default, they are both set at 100 dB. An example of the corresponding noise signals is presented in Figure 2. Figure 1d shows an example of noisy trajectories of shank and foot markers in the GCS, at the end of the generation process. Appl

AoR Estimation
In this section, the methodologies considered for the estimation of the AoR are presented. For each method, a brief description of the algorithm is reported, and references are provided. Moreover, a discussion about the averaging and the assessment of the dis-

AoR Estimation
In this section, the methodologies considered for the estimation of the AoR are presented. For each method, a brief description of the algorithm is reported, and references are provided. Moreover, a discussion about the averaging and the assessment of the dispersion and accuracy of the results is given.

Instantaneous and Mean Helical Axis
The algorithm for IHA estimation used is based on the original work of Woltring et al. [23], and it has been already outlined in previous works from the authors [10,19]. It considers a screw-like helical motion assumption. Firstly, this procedure requires the transformation of shank and foot markers' global trajectories in the local coordinate system (LCS) of one of the segments considered, that is, the shank. A representation of the shank LCS used to evaluate IHAs is presented in Figure 3a, whose elements are as follows:

4.
O S , the origin located at the mid-point of the SDistMed-SDistLat segment; 5.
w, axis orthogonal to the plane (u,v).
Then, a rigid body fitting algorithm based on singular value decomposition (SVD) is run to calculate an equivalent rigid body model of the foot. The latter is needed to estimate the six DoFs (three translational, three rotational) of the foot relative to the shank LCS that are then filtered through a self-tuning, 0-phase-shift, 2nd-order Butterworth digital filter [37]. The filtered DoFs and their time derivatives are finally used to calculate the axis parameters n i and s i that are the direction vector and the position vector representative of the IHA at the i-th time frame. It is well known that the calculation of the axis parameters is affected by the magnitude of the angular velocity ω i . It can be evaluated starting from an estimate of the attitude vector ϑ i , used to define the rotational DoFs, and its time derivative. Since low ω i values compromise the calculation of the IHA, a threshold of 0.25 rad/s is often considered to filter out the time frames that may lead to ill-posed estimation of the helical axis.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 8 of 26 [37]. The filtered DoFs and their time derivatives are finally used to calculate the axis parameters ni and si that are the direction vector and the position vector representative of the IHA at the i-th time frame. It is well known that the calculation of the axis parameters is affected by the magnitude of the angular velocity ωi. It can be evaluated starting from an estimate of the attitude vector ϑi, used to define the rotational DoFs, and its time derivative. Since low ωi values compromise the calculation of the IHA, a threshold of 0.25 rad/s is often considered to filter out the time frames that may lead to ill-posed estimation of the helical axis. Given the full representation of the joint kinematics provided by a set of IHAs, a mean helical axis (MHA) can be calculated by a weighted least squares algorithm as detailed in the literature [10,38,39]. According to Ehrig et al. [34], the most appropriate choice of the weight is: with ϑi as the angle of rotation corresponding to the rotational DoFs for the i-th time frame.
Other formulations consider different weighting factors [23,32,34], as follows: Given the full representation of the joint kinematics provided by a set of IHAs, a mean helical axis (MHA) can be calculated by a weighted least squares algorithm as detailed in the literature [10,38,39]. According to Ehrig et al. [34], the most appropriate choice of the weight is: with ϑ i as the angle of rotation corresponding to the rotational DoFs for the i-th time frame. Other formulations consider different weighting factors [23,32,34], as follows: To evaluate the effect of the different weights, the MHA estimation was performed considering all the three solutions presented in Equations (1)-(3).

Symmetrical Axis of Rotation Approach (SARA)
The SARA algorithm used in this work is based on the original work of Ehrig et al. [26] and considers two marker sets that are simultaneously moving relative to the GCS. To run the algorithm, an LCS is assigned to each marker cluster, that is, to each segment. The same shank LCS used to run the IHA algorithm was used, whereas the foot LCS was defined as shown in Figure 3b and includes:

1.
O F , the origin located at the centroid of the foot marker cluster; 2.
c, axis orthogonal to the plane (a,b).
By mathematical formulation, SARA considers that the estimated AoR, used to represent the relative motion between each LCS and the GCS, remains stationary with respect to the corresponding segment. Rigid body transformations are performed to evaluate the rotation matrices and translation vectors associated with each segment, which relate the "fixed" reference marker set represented in each LCS to the corresponding time-varying configuration in the GCS. This means that, for the shank segment, a rotation matrix R i and a translation vector T i must be estimated, whereas, for the foot segment, separate matrix S i and vector D i are considered. Each matrix and vector are evaluated for any time, similarly to the IHA parameters. As detailed in the literature [26,40], this calculation is performed by minimization of an objective function. This problem, which is equivalent to a linear least squares one, is solved via SVD and yields the functional axis parameters for each LCS, that is, a position vector c, which defines a reference point on the axis, and an orientation vector v. These parameters, due to the nature of the algorithm, are time-independent; therefore, two couples of constant vectors (c shank , v shank ) and (c foot , v foot ) are obtained as a consequence of the SVD. Each couple identifies the position and the orientation of the functional axis in the corresponding LCS. Due to measurement noise, when transformed in the GCS by means of the (R i , S i ) matrices and the (T i , D i ) vectors, the axes identified by each couple do not coincide, and therefore an average is performed for each time frame.
To ease the comparison between SARA and IHA approaches, the axes calculated by both methods were represented in the shank LCS. Unlike the IHA approach, which is subjected to an angular velocity threshold to avoid ill-posed estimation, SARA is able to perform the AoR evaluation for all the time frames available.

Assessment of Functional Axes' Dispersion and Accuracy
Compared to the SARA approach, given the different mathematical formulation, the IHA methodology yields the estimation of generally quite scattered set of axes that require some averaging, as already discussed regarding the MHA definition. However, it is interesting to assess how much the axes are dispersed around an average position and orientation. For this reason, linear (d eff ) and angular (χ eff ) dispersion parameters were calculated as a measurement of the rms distance and orientation between the MHA and the IHAs of the corresponding dataset [38].
Although it does not retain the same physical meaning, due to the different formulation, the same methodology was applied to the SARA dataset so as to evaluate an average AoR estimated by SARA algorithm. Indeed, due to the matrices (R i , S i ) and vectors (T i , D i ) being time-dependent, variability in the positioning and direction of the AoR must be expected for the SARA method too.
While dispersion parameters allow one to assess the variability in joint kinematics, they do not provide a measurement of the AoR estimation accuracy by themselves. Indeed, the axes might result to be quite low scattered but with a significant bias error in the orientation or position relative to the real AoR. For this reason, a measurement of the distance between the exact AoR, imposed by joint kinematics setup in the multi-body model, and the estimated AoR (both with helical axis and SARA algorithms) was performed, similarly to previous works in the literature [26]. To do so, the intersections between the estimated AoR and two quasi-sagittal planes perpendicular to the v axis of the LCS (see Figure 3a) and located at a distance of ±50 mm from the origin were evaluated. Similarly, the intersections of the exact AoR with the two planes were also calculated. Then, the distance between each couple of intersection points was computed in both quasi-sagittal planes (d Lat for the lateral and d Med for the medial side), yielding the final value of: where ∆ is the rms distance between exact and estimated AoRs, representing an overall descriptor of the AoR estimation accuracy. The large distance (100 mm) selected between the two quasi-sagittal planes was chosen to stress any misalignment between the exact AoR and the estimated one. Moreover, the two planes have also an anatomical interpretation since they could approximately represent two lateral and medial sagittal planes enclosing the ankle complex and, in particular, the malleoli, which are relevant anatomical references in the determining of the ankle joint rotation axis.

General Concepts
In this section, the design of the simulations performed on the multi-body model is presented. While each of the following subsections refers to a certain set of simulations (S1, S2, S3, and S4), based on a specific combination of parameters that will be discussed, four common motion configurations were considered throughout the whole analysis. As stated in Section 2.1.2, the same knee flexion/extension pattern was chosen for all the simulations, and therefore the expression "motion configuration" refers only to the ankle joint. In particular, the configurations tested are the following:

1.
A: 1 DoF, flexion/extension about a fixed axis, parallel to the v axis; 2.
B: 2 DoFs, flexion/extension about a moving axis, parallel to the v axis and translating horizontally in the antero-posterior direction (directed backward); 3.
C: 2 DoFs, screw-like helical motion, realized by means of flexion/extension occurring about an axis parallel to the v axis, combined with a translation (directed laterally) along the same axis; 4. D: 6 DoFs, tri-planar motion, flexion/extension occurs about an axis that has 3 translational DoFs, same trends as the ones selected for B and C configurations, and has an average inclination in the frontal plane of 10 deg (the medial tip lower than the lateral one). To simulate the limited contributions of the two remaining rotational DoFs, a small oscillation of ±1 deg about the average axis position was introduced, with the same curve trend of the flexion/extension DoF.
Not all the simulations presented in the following were repeated for each motion configuration for the sake of clarity and conciseness. In particular, A e D configurations were always considered since they represent, in our belief, the two most significant motion conditions among the ones tested. The performance of the algorithms for the other configurations (B and C) were discussed in the first set of simulations (S1), which is still the most extended among the ones presented in this paper, to allow some discussion about those specific patterns of motion. On the other hand, configuration A is the most simplistic one, and it allows for proving the performance of each algorithm in the most restrained environment. Configuration D, although with some approximations, allows for testing the algorithms in a more realistic environment, in which the joint under investigation has full tri-planar motion. Such configuration is the one depicted in Figure 1. Translational and rotational motion patterns for the other configurations (A, B, and C) were derived from the 6-DoFs tri-planar motion by "turning off" the unnecessary DoFs. In all the motion configurations tested, three full gait cycles were considered.
The noisy marker coordinates generated by the model, for each motion configuration and simulation parameter selected, were then elaborated through the AoR estimation algorithms shown in Sections 2.2.1 and 2.2.2. In particular, the helical axis estimation was performed for the three different weighting factors shown in Equations (1)-(3) and summarized in Table 2. Therefore, a total of four estimation algorithms were run for each multi-body simulation considered. For each repetition, dispersion parameters (d eff , χ eff ) and rms distance ∆ were evaluated, as detailed in Section 2.2.3. Table 2. Helical axis weighting and corresponding nomenclature used in the simulations.

Nomenclature
Weight w i References

First Set, S1: Range and Speed of Motion Patterns
In the S1 set, different ranges of motion (ROMs) were defined for each configuration, as follows: In particular, the range of motion for configuration A only refers to the flexion/extension angle, being the only DoF available. For configurations B and C (2 DoFs), the default range (+20, −30) deg was selected for the flexion/extension, whereas different ROMs were chosen for the translational DoF (horizontal displacement for B, lateral displacement for C configuration). The range of motion of flexion/extension for configuration D was varied similarly to configuration A, whereas the ROMs of the three translational DoFs were set as constants and equal to 2 mm. For all the configurations, positive flexion/extension angles refer to dorsiflexion, negative angles to plantarflexion.
The speed of each motion pattern was defined indirectly by tuning of the gait cycle duration t cyc . Three values of cycle duration were considered, equal to 1.6 s, 1.2 s, and 0.8 s, that correspond to increasing speed values labeled as V1, V2, and V3.

Second Set, S2: Scattering of the Markers within Each Cluster
In the S2 set, the scattering of the markers in each cluster was changed by tuning of the gains k S and k F , as detailed in Section 2.1.1. Three values of each gain were selected, equal to 0.75, 1, and 1.25. The 9 combinations of the gains were tested for configurations A and D, with the default range of motion of flexion/extension equal to (+20, −30) deg. A default cycle duration equal to 1.2 s was considered.

Third Set, S3: Speed of Motion Patterns and Sampling Frequency
Since the number of time frames obtained for each cycle can affect the performance of the AoR estimation algorithms, due to different time discretization of a kinematic pattern, several sampling frequencies f s were tested in the S3 set. These frequencies correspond to the target frequency of the downsampling applied to the original trajectory data produced by the model, which ran at 1000 Hz. The four values selected were equal to 100, 250, 500, and 1000 Hz. For a fixed sampling frequency, the discretization of joint motion also depends on its speed. For this reason, as considered in the S1 set, the speed of each motion pattern was varied by tuning of the gait cycle duration t cyc , yielding the same three increasing values of speed labeled as V1, V2, and V3. The simulations were performed for A and D configurations, with the default range of motion of flexion/extension equal to (+20, −30) deg.

Fourth Set, S4: Common and Random Sources of Additive Noise on Markers' Trajectories
In the S4 set, the effect of the different sources of accuracy errors was observed by tuning of the values SNR cmn and SNR rnd to evaluate the sensitivity of the AoR estimation algorithms to common (rigid) and random (non-rigid) additive noise on marker's trajectories. Four values were tested for each gain, that is, 100, 90, 80, and 70 dB. The simulations were performed for A and D configurations, with the default range of motion of flexion/extension equal to (+20, −30) deg. A default cycle duration equal to 1.2 s was considered.

Results
The results of the simulations detailed in Sections 2.3.2-2.3.5 are here presented. Each subsection refers to a specific set of simulations. For each set of simulations, the angular and linear dispersion parameters, as well as the rms distance values between the exact and estimated AoRs, were calculated. This process was repeated for each algorithm tested, that is, the IHA, with three different weighting factors (Table 2), and SARA. In the following, the main results of these analyses are presented, while the complete sets of results are reported in Appendix A for the sake of conciseness. Figure 4 shows an example of IHA sets resulting from the analyses, along with the resulting dispersion parameters and MHA, when a constant weighting factor is selected. The same calculation was performed for all the ROM/speed combinations and weighting factors, and then the results were compared to the outcomes of the SARA approach. Figure 5 shows a direct comparison between the MHAs obtained with the different weighting factors and the functional axes calculated with the SARA algorithm for motion configuration A. The exact AoR is also shown to ease the interpretation of the results. The same comparison, performed for configuration D, is shown in Figure 6, in two different views. All the plots refer to the shank LCS. The full set of the dispersion parameters and the rms distance resulting from the analysis is reported in Appendix A, Table A1. The number of frames (N f ) considered by each algorithm is also presented. Appl

Results
The results of the simulations detailed in Sections 2.3.2-2.3.5 are here presented. Each subsection refers to a specific set of simulations. For each set of simulations, the angular and linear dispersion parameters, as well as the rms distance values between the exact and estimated AoRs, were calculated. This process was repeated for each algorithm tested, that is, the IHA, with three different weighting factors (Table 2), and SARA. In the following, the main results of these analyses are presented, while the complete sets of results are reported in Appendix A for the sake of conciseness. Figure 4 shows an example of IHA sets resulting from the analyses, along with the resulting dispersion parameters and MHA, when a constant weighting factor is selected. The same calculation was performed for all the ROM/speed combinations and weighting factors, and then the results were compared to the outcomes of the SARA approach. Figure 5 shows a direct comparison between the MHAs obtained with the different weighting factors and the functional axes calculated with the SARA algorithm for motion configuration A. The exact AoR is also shown to ease the interpretation of the results. The same comparison, performed for configuration D, is shown in Figure 6, in two different views. All the plots refer to the shank LCS. The full set of the dispersion parameters and the rms distance resulting from the analysis is reported in Appendix A, Table A1. The number of frames (Nf) considered by each algorithm is also presented. Given the constant sampling rate equal to 100 Hz, an increase in motion speed corresponds to a reduction in the number of frames available to represent each gait cycle. When performing on the lowest number of samples, the IHA algorithm exhibits a relevant decrease in the general performance, even in configuration A, as signaled by the increase in   Given the constant sampling rate equal to 100 Hz, an increase in motion speed corresponds to a reduction in the number of frames available to represent each gait cycle. When performing on the lowest number of samples, the IHA algorithm exhibits a relevant decrease in the general performance, even in configuration A, as signaled by the increase in the linear dispersion. Especially for configurations A and B, this effect seems even clearer when the largest ROMs are considered, with the linear dispersion increasing up to 16.7 mm (conf. A) and 27.8 mm (conf. B). Curiously, the linear dispersion of the IHA set observed for configuration B is orthogonal to the horizontal displacement of the exact AoR imposed by the model, as shown in Figure 4b. Regarding the helical motion configuration (C), the IHA algorithm performs similarly for any of the ROMs considered; however, the linear dispersion is still increased and exceeds 10 mm at the highest speed. For the tri-planar motion configuration (D), the IHA algorithm shows similar results in terms of linear dispersion when the narrower ROMs are selected, regardless of the speed selected. However, for the widest ROM, a similar increase in linear dispersion is obtained at highest speed. As far as the rms distance ∆ is concerned, the same dependency with the ROM/speed combinations shown for the linear dispersion parameters is observed.   As regards the angular dispersion of the IHA sets, the results do not show a significant correlation with the ROM/speed combination selected. In particular, the angular dispersion is always lower than 1.4 deg for configuration A, 0.8 deg for configuration B, 1.1 deg for configuration C, and 7 deg for configuration D. The increase in angular dispersion in tri-planar motion is expected due to the increased variability of the exact AoR in that configuration. Those values refer to the worst scenario, whereas, in the average ROM/speed combination, the angular dispersion is lower than 0.6 deg for configurations A and B, 0.8 deg for configuration C, and 3.5 deg for configuration D. The largest dispersion values are always obtained when a constant weighting factor is selected. The use of the other different weighting factors definitely reduces the angular dispersion of the IHA sets; however, it also has a negative impact on the estimation of the MHA location, as signaled by the increased rms distance. This effect is observed in most of the ROM/speed combinations selected; however, it is far more evident in A and C configurations.

Results of S1 Set: Range and Speed of Motion Patterns
Compared to the helical axis method, the SARA algorithm shows far more stable results for most of the ROM/speed combinations selected, as expected. This result is particularly evident in A, B, and C configurations, both by considering the linear and angular dispersion values. The linear dispersion reaches 0.04 mm for the A and C configurations, both of which consider a stationary AoR along the u and w axes of the shank LCS, whereas it increases up to 0.52 mm and 2.54 mm for B and D configurations, respectively. Similar results are obtained as far as the angular dispersion is concerned, achieving values lower than 0.3 deg for all the ROM/speed combinations and motion configurations tested. However, the most significant improvement shown by the SARA methodology is related to the estimation accuracy of the AoR, as signaled by the low rms distance. For most of the ROM/speed combinations, the rms distance ∆ measured between the SARA-estimated and the exact AoRs is lower than the one calculated by considering the MHA. This result is particularly evident for A (see Figure 5) and C configurations (∆ lower than 0.04 mm); however, it is confirmed for most of the trials performed by choosing the B configuration. Interestingly, the results obtained for the tri-planar motion configuration (D, see Figure 6) are comparable with the ones yielded by the helical axis methodology, with the latter performing even better in a few ROM/speed combinations. In particular, the IHA algorithm shows a limited improvement in the AoR estimation accuracy for the narrowest and average ROMs, regardless of the speed selected.

Results of S2 Set: Scattering of the Markers within Each Cluster
The dispersion parameters and the rms distance resulting from the analysis, for the different gains k S and k F tested, are reported in Appendix A, Table A2.
As far as the linear and angular dispersion parameters are concerned, no significant dependency with the different scattering is observed. In any of the combinations selected, the linear dispersion is lower than 2.5 mm when a simple 1-DoF motion configuration (A) is considered, and it increases up to about 11.3 mm when tri-planar motion (D) is imposed. Similarly, angular dispersion is always lower than 0.75 deg in configuration A and lower than 3.5 deg in configuration D.
As far as the rms distance is concerned, the SARA outperforms the helical axis for all the combinations tested in configuration A, whereas the helical axis algorithm shows better estimation accuracy for four out of nine trials performed in configuration D. Unlike the dispersion parameters, an increased variability in the estimation accuracy depending on the different scattering combinations is observed for the 1-DoF configuration. For k S = 0.75, the least-accurate results are obtained when the largest scattering of the foot cluster is selected. In the other conditions (k S = 1 or k S = 1.25), the least-accurate estimation is obtained when the lowest and average scattering of the markers in the foot cluster is considered. However, the absolute differences between the rms distance values resulting from each trial are quite limited and do not signal a particular trend. In configuration D, the ∆ values obtained are almost constant and do not seem to be correlated to the several scattering combinations tested.

Results of S3 Set: Speed of Motion Patterns and Sampling Frequency
The dispersion parameters and the rms distance values resulting from the analysis, for the different speed and sampling frequency values selected, are presented in Appendix A, Table A3. The number of frames (N f ) considered by each algorithm is also presented.
When the lowest (V1) and average (V2) speed are considered, a limited improvement in the estimation of AoR positioning is observed by rising the sampling frequency for configuration A, whereas no significant improvement is shown for configuration D. In particular, for the 1-DoF motion configuration, the linear dispersion and rms distance are halved at 1000 Hz (about 1 mm), while the angular dispersion falls below 0.3 deg (starting from 0.9-1.64 deg at 100 Hz).
At highest speed (V3), the improvement in linear dispersion and rms distance obtained at a high sampling rate is evident for both motion configurations, although it is very clear especially for configuration A. The linear dispersion obtained at 100 Hz (10.6 mm) is reduced to 2.46 mm at 250 Hz, and it falls below 2 mm at higher frequencies. Similarly, the angular dispersion is generally reduced and achieves a very low magnitude (below 0.1 deg) at maximum frequency. The same behavior is observed for configuration D; however, the improvement in estimation accuracy is less significant.
The AoR estimation based on the SARA method shows an almost constant level of accuracy for most of the trials analyzed, which is generally higher than the one observed for the IHA approach. Figure 7 shows the estimation of the angular velocity performed by the IHA algorithm and its comparison with the angular velocity set as input of the multibody model, by choosing motion configuration A. As shown in Figure 7a, the evaluation accuracy of the angular velocity at low speed (V1) is very high, as the AoR estimation. Figure 7b shows the same signals when a higher speed (V3) is selected. In this condition, by assuming a constant sampling rate, the angular velocity is predicted with less accuracy as compared to the result shown in Figure 7a. This outcome affects the accuracy of AoR estimation performed by means of the helical axis approach, whereas it has no impact on SARA methodology.

Results of S4 Set: Common and Random Sources of Additive Noise on Markers' Trajectories
In Appendix A, Table A4, the dispersion parameters and the rms distance resulting from the analysis are collected for different combinations of common and random noise signals added to the markers' trajectories provided by the multi-body model. Figures 8  and 9 show a direct comparison between the MHAs obtained with the different weighting factors and the functional axes calculated with the SARA algorithm, respectively, for motion configurations A and D. The exact AoR is also shown to ease the interpretation of the results. All the plots refer to the shank LCS.
The results highlight the significant impact of marker noise on the accuracy of AoR estimated by helical axis methodology. This effect seems to be particularly highlighted in the transition between SNRcmn = 100 dB and SNRcmn = 90 dB, and it increases at low SNRrnd values. For SNRcmn values lower than 90 dB, the addition of a random noise signal does not seem to increase the linear dispersion of the IHA sets in a significant way. The same behavior is highlighted both in A and D motion configurations. On the other hand, the angular dispersion is generally increased by lowering the SNRs of both source signals, and this effect is even more relevant for configuration A. Surprisingly, although linear and angular dispersion parameters show a significant increase even at relatively low levels of noise, the accuracy of AoR estimation by means of the MHA is generally not compromised for SNRcmn higher than 80 dB, regardless of the SNRrnd values selected. This effect is more evident for configuration A as compared to the tri-planar motion configuration (D).
By comparing the results presented so far with the outcomes of the SARA approach, the latter clearly outperforms the helical axis methodology for all the noise conditions tested (see Figures 8 and 9). In particular, the rms distance between the exact and esti-

Results of S4 Set: Common and Random Sources of Additive Noise on Markers' Trajectories
In Appendix A, Table A4, the dispersion parameters and the rms distance resulting from the analysis are collected for different combinations of common and random noise signals added to the markers' trajectories provided by the multi-body model. Figures 8  and 9 show a direct comparison between the MHAs obtained with the different weighting factors and the functional axes calculated with the SARA algorithm, respectively, for motion configurations A and D. The exact AoR is also shown to ease the interpretation of the results. All the plots refer to the shank LCS.
The results highlight the significant impact of marker noise on the accuracy of AoR estimated by helical axis methodology. This effect seems to be particularly highlighted in the transition between SNR cmn = 100 dB and SNR cmn = 90 dB, and it increases at low SNR rnd values. For SNR cmn values lower than 90 dB, the addition of a random noise signal does not seem to increase the linear dispersion of the IHA sets in a significant way. The same behavior is highlighted both in A and D motion configurations. On the other hand, the angular dispersion is generally increased by lowering the SNRs of both source signals, and this effect is even more relevant for configuration A. Surprisingly, although linear and angular dispersion parameters show a significant increase even at relatively low levels of noise, the accuracy of AoR estimation by means of the MHA is generally not compromised for SNR cmn higher than 80 dB, regardless of the SNR rnd values selected. This effect is more evident for configuration A as compared to the tri-planar motion configuration (D).     By comparing the results presented so far with the outcomes of the SARA approach, the latter clearly outperforms the helical axis methodology for all the noise conditions tested (see Figures 8 and 9). In particular, the rms distance between the exact and estimated AoRs remains at very low values (below 5 mm) for almost all the trials, whereas it reaches about 20 mm (in the worst scenario) when the helical axis method is used. Figure 10 shows the estimation of the angular velocity performed by the IHA algorithm and its comparison with the angular velocity set as input of the multi-body model by choosing motion configuration A and different SNR cmn values. A progressively worse evaluation of the angular velocity is obtained by increasing the SNR of the common noise signal, affecting the accuracy of the estimation algorithm.

Discussion
This paper proposed a comparison between the helical axis and SARA approaches for the estimation of functional joint axes, starting from simulated motion capture data provided by the multi-body modeling approach. Compared to similar studies from the literature that considered simulated motion data [26], the multi-body model provided additional flexibility to the configuration of the motion patterns, thus allowing an extended comparison between the performance of the two methodologies considered. In particular, this work provided a quantitative analysis in terms of dispersion and estimation accuracy of the functional joint axes estimated by the two algorithms that can improve the knowledge of the algorithms themselves, in support to other studies that compared the methods by a mathematical perspective [34]. Moreover, the effect of the different weighting factors in MHA estimation has not really been discussed by objective evaluation of the results, as presented in the current work. Reichl et al. [40] presented a comparison between finite helical axis and SARA on simulated motion data; however, they focused on a simpler knee joint model based on a compound hinge joint with 2 DoFs and did not provide any discussion about marker positioning, speed of motion patterns, and sampling rate. A similar analysis was conducted by the authors on experimental data [35]; however, this approach lacks the possibility to change the conditions of the trials with the same flexibility shown in the current paper.

AoR Estimation by Means of Helical Axis is Significantly Affected by Time Discretization
As reported in Section 3.1, by varying the range and speed of motion patterns, a significant dependency between estimation accuracy and the number of time frames available was observed for the IHA algorithm. In particular, this phenomenon was highlighted when high speed and range of motion are considered. Indeed, given the constant sampling rate of the simulated motion capture data, those conditions correspond to the worst time discretization of the trajectories used to reconstruct joint motion. As confirmed by the simulations presented in Section 3.3, the improvement in IHA estimation accuracy, when higher than 100 Hz sampling frequencies are selected, is obvious, as signaled by the low values of linear dispersion deff and rms distance Δ achieved (Table A3). Of course, high sampling frequency values such as 1000 Hz are not desirable during real experimentation due to the limitations of the instrumentation and data acquisition systems, though

Discussion
This paper proposed a comparison between the helical axis and SARA approaches for the estimation of functional joint axes, starting from simulated motion capture data provided by the multi-body modeling approach. Compared to similar studies from the literature that considered simulated motion data [26], the multi-body model provided additional flexibility to the configuration of the motion patterns, thus allowing an extended comparison between the performance of the two methodologies considered. In particular, this work provided a quantitative analysis in terms of dispersion and estimation accuracy of the functional joint axes estimated by the two algorithms that can improve the knowledge of the algorithms themselves, in support to other studies that compared the methods by a mathematical perspective [34]. Moreover, the effect of the different weighting factors in MHA estimation has not really been discussed by objective evaluation of the results, as presented in the current work. Reichl et al. [40] presented a comparison between finite helical axis and SARA on simulated motion data; however, they focused on a simpler knee joint model based on a compound hinge joint with 2 DoFs and did not provide any discussion about marker positioning, speed of motion patterns, and sampling rate. A similar analysis was conducted by the authors on experimental data [35]; however, this approach lacks the possibility to change the conditions of the trials with the same flexibility shown in the current paper.

AoR Estimation by Means of Helical Axis Is Significantly Affected by Time Discretization
As reported in Section 3.1, by varying the range and speed of motion patterns, a significant dependency between estimation accuracy and the number of time frames available was observed for the IHA algorithm. In particular, this phenomenon was highlighted when high speed and range of motion are considered. Indeed, given the constant sampling rate of the simulated motion capture data, those conditions correspond to the worst time discretization of the trajectories used to reconstruct joint motion. As confirmed by the simulations presented in Section 3.3, the improvement in IHA estimation accuracy, when higher than 100 Hz sampling frequencies are selected, is obvious, as signaled by the low values of linear dispersion d eff and rms distance ∆ achieved (Table A3). Of course, high sampling frequency values such as 1000 Hz are not desirable during real experimentation due to the limitations of the instrumentation and data acquisition systems, though significant improvements in estimation accuracy are shown even for more feasible solutions, such as 250 and 500 Hz. Moreover, high discretization would result in a large number of time frames that can affect the efficiency of the algorithms in applications where the computational time is relevant. As shown in Section 3.3, the SARA approach does not seem to be significantly influenced by the number of time frames, regardless of the speed selected. This result is quite interesting and widens the field of application of this methodology to fast and wide motion patterns, even for relatively low sampling frequency. The dependency between AoR estimation accuracy, in terms of linear dispersion, and time discretization of the motion patterns is highlighted by the results shown in Figures 11a and 12a  selected. This result is quite interesting and widens the field of application of this methodology to fast and wide motion patterns, even for relatively low sampling frequency. The dependency between AoR estimation accuracy, in terms of linear dispersion, and time discretization of the motion patterns is highlighted by the results shown in Figures 11a  and 12a for representative conditions.  The results presented in Section 3.1 show that, when the linear dispersion of the IHA set is high, mainly due to large ROMs and speed, the rms distance is also increased (see Figures 11b and 12b). This result signals that, in general, an increase in the linear dispersion has an impact on the estimation accuracy of the MHA position, regardless of the weighting factor chosen. As far as the angular dispersion χeff is concerned, it seems to be independent of the several ROM/speed combinations tested, but it is affected far more significantly by the weight selected for MHA estimation. In particular, the correction proposed by Ehrig et al. [34] shows the most limited angular dispersion among the several solutions, but it has an impact on the accuracy of AoR estimation in terms of its positioning (see Table A1, IHA (w.3) columns). The corresponding increase in rms distance was  selected. This result is quite interesting and widens the field of application of this methodology to fast and wide motion patterns, even for relatively low sampling frequency. The dependency between AoR estimation accuracy, in terms of linear dispersion, and time discretization of the motion patterns is highlighted by the results shown in Figures 11a  and 12a for representative conditions.  The results presented in Section 3.1 show that, when the linear dispersion of the IHA set is high, mainly due to large ROMs and speed, the rms distance is also increased (see Figures 11b and 12b). This result signals that, in general, an increase in the linear dispersion has an impact on the estimation accuracy of the MHA position, regardless of the weighting factor chosen. As far as the angular dispersion χeff is concerned, it seems to be independent of the several ROM/speed combinations tested, but it is affected far more significantly by the weight selected for MHA estimation. In particular, the correction proposed by Ehrig et al. [34] shows the most limited angular dispersion among the several solutions, but it has an impact on the accuracy of AoR estimation in terms of its positioning (see Table A1, IHA (w.3) columns). The corresponding increase in rms distance was highlighted especially for motion configurations A (1 DoF, flexion/extension) and C (2 The results presented in Section 3.1 show that, when the linear dispersion of the IHA set is high, mainly due to large ROMs and speed, the rms distance is also increased (see Figures 11b and 12b). This result signals that, in general, an increase in the linear dispersion has an impact on the estimation accuracy of the MHA position, regardless of the weighting factor chosen. As far as the angular dispersion χ eff is concerned, it seems to be independent of the several ROM/speed combinations tested, but it is affected far more significantly by the weight selected for MHA estimation. In particular, the correction proposed by Ehrig et al. [34] shows the most limited angular dispersion among the several solutions, but it has an impact on the accuracy of AoR estimation in terms of its positioning (see Table A1, IHA (w.3) columns). The corresponding increase in rms distance was highlighted especially for motion configurations A (1 DoF, flexion/extension) and C (2 DoFs, helical motion), whereas it was quite negligible for the other configurations that share linear translation DoF of the AoR along the antero-posterior direction. The presence of translational motion occurring in the sagittal plane might mix up the analysis and should be investigated in the future. However, the results obtained by choosing the different weighting factors do not clearly support the use of one weight rather than the other. In general, the corrections help in achieving a more confined set of helical axes; however, they may compromise the identification of the correct location of the AoR.

SARA Generally Provides Better Estimation Accuracy Due to Its Intrinsic Stability
Compared to the helical axis approach, SARA highlighted the best performance in terms of rms distance between the estimated AoR and the exact one, which is the parameter used to assess the accuracy of the AoR estimation. In particular, this result was confirmed for A, B, and C configurations for all the conditions tested, whereas both algorithms performed in a similar way for the tri-planar motion configuration (D, as shown in Figures 11b and 12b). As far as the linear and angular dispersion parameters are concerned, SARA showed the lowest values in most of the conditions observed; however, this result was expected since it is directly connected to its formulation being quite different than the IHA's. As presented in Section 2.2.2, SARA initially considers two separate LCSs, respectively integral with the shank and foot, to evaluate the relative motion between each segment and the fixed GCS. In this calculation, each AoR is considered stationary with respect to the corresponding segment. Then, the functional axis obtained at the end of the procedure is estimated as the average of the AoRs calculated for each segment. Compared to the SARA approach, the helical axis does not consider any of those hypotheses for determining the relative motion between two rigid segments, yet the additional freedom of this methodology plays a key role in the "stability" of the results. For this reason, the variability in the IHA sets is expected to be more pronounced than the one observed as the outcome of the SARA algorithm. This aspect must be taken into account when the two methods are compared. For tri-planar motion (see Table A1), the angular dispersion obtained for the IHA approach is quite a good approximation of the actual orientation variability of the joint AoR imposed by the model; however, this is not confirmed for all the ROM/speed tested. On the contrary, linear dispersion is always an overestimation of the actual linear variability of the joint AoR. This confirms the weakness of the helical axis approach in determining the correct positioning of the AoR, already highlighted in the literature [34].

SARA Is More Suitable Than IHA for Kinematic Analyses Based on Noisy Markers' Trajectories
As pointed out in the literature [40], the intrinsic variability of the helical axis approach also affects its sensitivity to noise. As confirmed by the results presented in Section 3.4, in which different levels of additive noise were considered to alter the markers' trajectories as provided by the multi-body model, the performance of the IHA algorithm considerably worsen when SNR cmn values lower than 80 dB are selected. In particular, the AoR estimation accuracy seems to be mainly affected by the source of the common noise signal, which is used to simulate non-rigid deformations of the marker clusters such as the one related to the soft tissue artifact. Oddly, the rms distance measured at the maximum level of the common noise signal decreases while lowering SNR rnd . This result might be due to the high level of random noise, which could contribute to the overall noise signal by subtracting to the common one. The SARA algorithm, on the other hand, can benefit from the lack of numerical differentiation and filtering that IHA methodology requires [34,37]; therefore, it might be more suitable for any application based on position-based rather than velocity-based datasets, such as those provided by motion capture analysis. Numerical differentiation can indeed intensify the effect of noise on marker coordinates and compromise the reconstruction of joint kinematics. In addition, unlike the IHA methodology, the SARA algorithm does not require the calculation of the foot DoFs relative to the shank LCS, thus removing the impact of the errors associated with coordinate transformations. In support of these considerations, the AoR estimation performed by the SARA algorithm was always more accurate than the one obtained by application of the IHA methodology, even for high levels of noise.

Different Scattering of the Marker Clusters Has a Limited Impact on AoR Estimation Accuracy
Finally, the simulations presented in Section 3.2, which considered different scattering of the marker clusters on the surface of each segment, showed that no significant dependency between dispersion parameters and marker disposition was observed for both algorithms. The larger variability in estimation accuracy was obtained for the 1-DoF motion configuration (A); however, it was not possible to derive a correlation with the marker scattering in each cluster. This aspect should be taken into account during experimentation since marker positioning has an impact on bone pose estimation, which could not be investigated with the modeling approach considered in this study.

Limitations
The main limitations of this work are represented by the approximations in the ankle joint motion patterns that were provided as input for the multi-body model. In particular, the three rotational DoFs of the bushing joint do not accurately match the actual physiological ones. Unlike the anatomical interpretation of joint axes (antero-posterior, longitudinal, medio-lateral), the modeling of the joint in the multi-body approach considers three mutually orthogonal axes. Moreover, the DoFs associated with any joint axes are independent by mathematical formulation of the bushing. This feature leads to a simplified representation of the joint DoFs, which result to be decoupled in the model, whereas they cannot be isolated in the real scenario. For these reasons, the simulated flexion/extension, adduction/abduction, and inversion/eversion motion patterns can only approximate the real joint DoFs and have limited meaning from the anatomical/physiological point of view. However, this approximation was considered reasonable for the aim of the current work, which focused on the objective evaluation of the helical axis and SARA techniques in a realistic scenario, rather than on the accurate representation of ankle joint motion through multi-body modeling. Moreover, this simplification greatly facilitates the formulation of the model.

Conclusions
This paper showed a comparison between different approaches for the estimation of functional joint axes. The work was conducted at a simulation level by developing a multibody model resembling the lower limb of a subject, which was used to generate simulated motion capture data. The helical axis and SARA approaches were used to estimate the AoR of the simulated ankle joint for several motion patterns and testing conditions. The main conclusions are reported in the following:

1.
For the 1-DoF and 2-DoFs configurations, the SARA proves to be more accurate than the helical axis methodology, regardless of the different conditions considered; 2.
For the 6-DoFs configuration, the SARA outperforms the helical axis approach when a high level of noise is considered or when the number of samples used to discretize the gait cycle is limited; 3.
The IHA methodology allows a description of joint kinematics that is intrinsically less constrained than the one enabled by the SARA. This could be useful when a more detailed, time-dependent representation of joint kinematics is necessary.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement: Data available on request.
Acknowledgments: The author wishes to thank Laura Bellaera for her valuable contribution in the code implementation of the SARA algorithm during her M.Sc. thesis.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A
This section includes the dispersion parameters and rms distance values obtained for each set of simulations (S1 . . . S4) and for the several algorithms tested.  Table A4. Dispersion parameters (d eff , in mm, and χ eff , in deg) and rms distance (∆, in mm) between exact and estimated AoRs for different common (SNR cmn ) and random (SNR rnd ) signal-to-noise ratios (in dB) of the corresponding sources of noise signals. Results are given for A and D motion configurations.