A PCA-based bio-motion generator to synthesize new patterns of human running

Data gathering

An experimental phase was carried out with the aim of gathering a database of the running movements that we needed to include in the biomechanical model. The data consisted of the 3D joint angles of the main body joints.

The articulated human body model used in our study comprised 21 segments and 20 joints distributed throughout the body. Lower limb: hip, knee, ankle and metatarsophalangeal joint; upper limb: shoulder, elbow and wrist; trunk: pelvis, L5, L3, T12, T8 and Neck.

The positions of the joints of the human body model were defined in a recursive mode with respect to the origin joint (father) of the related segment. This methodology is based on the BioVision Hierarchical data (BVH) format (Meredith & Maddock, 2001).

Each subject was kinematically characterised by means of 64 variables, defined as follows:

• 1.

  Vertical (Z) position of the root segment, in our case the hips.

• 2.

  Tri-dimensional orientation (X, Y, Z) of the total amount of segments with respect the root. 21 × 3 = 63 variables.

Study sample

Eighteen people composed the study sample, with the same number of male and female. Their ages ranged from 21 to 44 years (average age: 31 years). They were selected according to some specific parameters, trying to cover a wide range in the anthropometric characteristics of height and body mass index (Table 1). Ethical approval was obtained from the ethics committee of the Universitat Politècnica València. All participants gave written informed consent.

Measurement and protocol

The measurements were performed using commercial equipment based on 17 inertial sensors: MOVEN studio. The commercial system has been validated by previous studies (Zhang et al., 2013; Thies et al., 2007). A sampling frequency of 120 Hz was used. This system showed a very high sensitivity to electromagnetic fields. For this reason, the measurements of running trials were done outdoors in a location free of electromagnetic pollution.

Experimental procedure

For the purpose of controlling the pace of running, a 20-metre-long corridor, delimited with cones every 5 m was set up. Thus, we obtained four areas, one area of acceleration, two of constant speed and a final deceleration area. Running at constant speed presents a periodic timing in which the period depends on the velocity (Novacheck, 1998). Nevertheless, acceleration and deceleration periods are out of phase and the duration of cycles is variable. Therefore, the running cycles used to create our model were selected within the area of constant speed.

In the case of running, the pattern of the movement changes with velocity (e.g., stride length, maximum joint angles, etc.). For this reason, each runner completed six running trials at different speeds. Initially, subjects started running at normal speed. In the second measurement subjects ran at their maximum speed. The third and fourth trial were performed at a pace between normal and maximum speed. The fifth trial was performed at the minimum speed at which each runner was able to run, on the edge between walking and running. Running defined in this case as when there is no phase of bipodal support (Biewener et al., 2004). The last trial was performed at an average speed between the lower and the normal speed.

This procedure allowed us to obtain six observations representing the whole range of speeds that each subject could execute.

Mathematical procedure

The methodology used in our study comprised five steps:

1-Reduction of intra-personal variability: joint angles are periodic by nature. We took the most representative single stride for the purpose of reducing variability and dismiss the phases with no consistent speed, such as acceleration and deceleration steps. The selected stride was picked in the middle of the running sequence, guaranteeing constant speed.

2-Time normalization: time normalization is usually employed for the temporal alignment of cyclic data obtained from different trials with different duration. In our approach, the number of samples for each stride depended on the velocity of the running. At this point we normalised the variable “time” applying an interpolation technique based on cubic splines through the measured values of the whole sequence of samples. This technique enables the normalization of all the measurements to the percentage of the running cycle. The cubic spline was applied to normalize at 50 equispaced time intervals per each variable (see the example in Fig. 1). The application of the cubic spline to the 64 kinematical variables makes a total amount of 50 × 64 (3,200) data per each subject.

3-Data cleaning: at this point a detailed checking and cleaning of inaccuracies of the kinematical data was conducted. These type of inaccuracies were caused mostly by the measurement system. The prevention of errors at this point is preferable to their later correction once the model has been created. All the measurements have been manually analysed thoroughly by an experienced examiner. The identified inaccuracies were treated as follows:

• (a)

  Angular offsets: this common error usually appears during the standing posture and can affect the later joint angles registers during the trial (running) (Mills et al., 2007). Offsets have been corrected manually eliminating (adding or subtracting) the difference between the initial angle observed and the expected angle of the body segment at this position.

• (b)

  Positioning error: due to the fact that our measurement system, based on inertial sensors, uses the earth’s magnetic field to determine the reference position of each subject, it is quite common to find subjects with slight differences in their initial reference positions. In this case, we have proceeded by correcting the reference system aligning it with the direction of running forward. Thus, we can guarantee that all measurements are equally oriented.

• (c)

  Non-physiological angles: some errors in the registration of joint angles were detected in the database. These inaccuracies came from errors of the inertial sensors. In this case, it was not possible to correct the error effectively, thus we proceed by eliminating these observations from the database and repeating the measurement.

4-Dimensionality reduction: the database of all measurements was combined in a single matrix W. The initial number of observations is 108 (18 subjects × 6 velocities = 108). But three measurements fail, therefore W has 105 rows (observations) and 3,200 columns (50 equispaced time intervals × 64 kinematical variables). Motion data of each observation is enclosed in the rows of the matrix $W=\left(w_i\right),i=1,\dots ,105$.

Before the creation of the bio-motion generator, by means of a regression model, it was needed to reduce the dimensionality of the motion data. Computing a PCA on the running data (contained in matrix W), resulted in a decomposition of the data matrix W into an average running vector w₀ and 3,200 weighed components, arranged in a 3,200 × 3,200 matrix V: (1)${\displaystyle W=W_0+\alpha \cdot V}$ where W₀ is a 105 × 3,200 matrix with all rows equal to w₀ and $\alpha =\left({\alpha }_i\right)$ with i = 1, …, 105 is a 105 × 3,200 matrix of PCA scores. Each observation w_i was thus expressed by a linear combination of scores α_i and PCA components (columns of matrix V). Components represented factors related to gender, anthropometrical traits and running speed. And scores represented individual characteristics of each runner and performance of the running trial related to the previous factors. PCA components are arranged in descendent order of explained variance of the original matrix data. Thus the first columns in matrix V retained most of the information in the data sample and it was possible to select a reduced number of components c. (2)${\displaystyle w_{ic}=w_0+{\alpha }_{ic}\cdot V_c.}$

Above w₀ denoted the average of all the running samples. The matrix V_c contained the first components. α_ic represents the c scores of each observation of the database in the reduced dimension space formed by the selected components. As score values change from negative to positive values, the movement of the runner changes from men to women; high BMI to low BMI; high speed to low speed, etc.

The decision of how many components to retain was a critical issue in the exploratory factor analysis. To perform this decision we used the methodology of Parallel Analysis (PA) (Hayton, Allen & Scarpello, 2004). PA is a Monte-Carlo based simulation method that compares the observed eigenvalues (components) with those obtained from randomized normal variables. A component is retained if its explained variance or information is higher than the information provided by the eigenvectors derived from the random data.

5-Regression model: one of the objectives of our work was to generate a statistical model capable of synthetizing new realistic running motion from a set of desired data: age, gender, height, weight and velocity, also called 1D data. Accordingly, to devise the bio-motion generator, we needed to establish the correlation between the 1D data and PCA scores, which provide the signature of each motion.

The correlation was obtained as a regression model, combining a partial least squares (PLS) regression model as a first step and a linear regression model (LRM) as a second step.

PLS methodology is explained in Wold (2006) and Geladi & Kowalski (1986). This type of regression model is suitable for the kind of data involved in the bio-motion generator since the input data of the model is strongly correlated (anthropometrical information).

The PLS regression model takes the 1D data—age, height, weight and velocity—as input information and produces a set of PCA scores as output. The LRM model was applied to these output PCA scores to reflect the influence of gender in the PCA scores.

In the first step, we estimated a PLS model considering anthropometrical data and velocity of the movement as independent variables and the PCA scores as dependent variables. The general formula of a PLS model is: (3)${\displaystyle Y-Y_0=B\cdot \left(X-X_0\right)+E}$ where Y is the matrix of dependent variables, X is the matrix of independent variables, X₀ and Y₀ are the matrices of mean data, B is the coefficient matrix of the PLS model and E is the prediction error matrix.

Having four different 1D measurements for each subject x_i = [age, height, weight, velocity], we built a 105 × 4 matrix X formed by the concatenation of the vectors x_i. Notice that since matrix =α , it is constituted by PCA scores and the mean matrix Y₀ is the zero matrix. Finally, we estimated the coefficient matrix B from sample data X and α with the PLS algorithm and substituted in the equation above, obtaining the following model: (4)${\displaystyle \alpha =B\cdot \left(X-X_0\right)+E.}$ Where X₀ is the matrix of mean anthropometrical data and velocity, and E is the prediction error matrix. PLS decomposes the independent and dependent variables in component spaces in order to obtain their correlation. The number of significant PLS components in the model was selected in a leave-one-out procedure and according to the explained variance (R²) criteria.

Secondly, the influence of gender was modelled with a LRM of the prediction error matrix E with coefficients a =(a₁, …, a_c) and b =(b₁, …, b_c), where c is the number of retained PCA components: (5)${\displaystyle \hat{E}=\mathit{a}+\mathit{b}\cdot gender}$ where gender = 0 for men and gender = 1 for women.

This way, the motion information related to gender which is part of the PLS error matrix E, and uncorrelated with the prediction derived from the PLS regression, was modelled. Notice that a_j and b_j were considered zero whenever their F-value was below a desired level of statistical significance of 95%.

Once we have obtained B, a and b, whenever we want to synthesize a running motion from new anthropometrical data and velocity, we obtain the corresponding scores $\hat{\alpha }$ of the new realistic running motion by the following formula: (6)${\displaystyle \hat{\alpha }=\left[\mathit{a}-B\cdot X_0\right]+\left[\begin{array}{c}B\hfill \\ \mathit{b}\hfill \end{array}\right]\cdot \left[\begin{array}{cc}X\hfill & gender\hfill \end{array}\right]}$ where X₀ is the matrix of mean anthropometrical data and velocity.

Validation methodology of the bio-motion generator

To validate the five-step methodology described to develop the bio-motion generator we propose a comparison between each recorded observation and the prediction of running motion generated by the model by means of the ‘leave-one-out’ procedure. The recorded observation is considered the true angle curve of the running motion. The predicted motion is estimated using the ‘leave-one-out’ validation technique; that is, not including that observation in the bio-motion generator. We wish to determine if both curves are reproducible and sufficiently similar to consider that they represent the same motion. For this purpose, we use the Intraclass Correlation Coefficient (ICC) as a measurement of the reliability, and the Standard Error of Measurement (SEM) as a direct measurement of the global error between true and predicted angles. Theoretically, the ICC is defined as the ratio between the true variance and the predicted variance. The ICC varies between 0 and 1 and can be interpreted as the proportion of variance due to the methodology (true versus predicted data) in the total variance. An ICC greater than 0.8 is generally considered to be good (Fleiss, 1999). The ICC is determined between the measured or true curve (T_c) and the estimated curve (E_c) provided by the bio-motion generator. The ICC is determined from the variance of both curves (T_c) and (E_c) following the next equation: (7)${\displaystyle \mathrm{ICC}=\frac{{\sigma }_{E_c}^2}{\left({\sigma }_{T_c}^2+{\sigma }_{E_c}^2\right)};}$

On the other hand, SEM represents the existent difference between observed (T_c) and estimated curves (E_c) determined with the bio-motion generator, and provide an indication of the real magnitude of the error. (8)${\displaystyle \mathrm{SEM}={\sigma }_c\cdot \sqrt{\left(1-\mathrm{ICC}\right)};}$ Where σ_c is the combined standard deviation of the true scores (T_c) and observed scores (E_c). And S_E is the combined standard deviation of the true scores and observed scores.

We have obtained the SEM for each pair of true and predicted angles for the three spatial directions in all the joints that form the human model. For that reason, we have represented the SEM by its descriptive statistics (mean, std., 5-percentile and 95-percentile).

Parallel analysis

The results of the PA (Fig. 2) have been obtained with the explained variance of the main components extracted from the original data and the same obtained from randomized data. The intersection point of both curves indicates the optimal number of components to extract from the PCA. The original number of dimensions was 72 (3 related to the pelvis translation + 69 related to the body segments orientation). The results of the PA recommend retaining the first 12 eigenvalues, which explain the 88.16% of the total variance. Thus, the PCA allowed a percentage of data reduction of 83%, from 72 variables to 12 weighed components.

Regression model

As it has been explained in the methodology, the regression model consists of two parts, the first including the anthropometrical data (PLS) and the second the gender (LRM). The dependent variables of the PLS are the scores of the first 12 principal components (PC) of the kinematical running motion. Therefore, they are uncorrelated and the optimal number of PLS components are separately determined for each PC score (PC 1... PC 12) according to its adjusted R² plot (Fig. 3). PLS components are retained until their R² curve exhibits a decrease or a non-significant increase. Thus, for instance, two PLS components are retained for PC 1, whereas no components are considered for PC 7 and PC 9. Notice that for those PC with 0 retained components, the PLS model provides their mean value as output. This way, the motion information associated to those PC which is provided by the PLS model is the average motion.

With regard to the LRM, which analyses the influence of gender in the kinematics of running, the PCA scores which are significantly affected by gender are PC3, PC7, PC8, PC9 and PC11 (Table 2). The prediction obtained in the first step of the model is improved by the influence of gender on these PC. PC 7 and PC 9 are only affected by gender, since their number of retained PLS components was 0.

Validation of the bio-motion generator

The results of the reliability study, computed from the 90 observations and the same calculated by means of the leave-one-out technique, showed that the mean and standard deviation of ICC, was 0.91(0.04) with a 5 percentile of 0.829 and 95 percentile of 0.971. Only one subject exhibit an ICC lower than 0.8 in two observations (Fig. 4).

The SEM between the real and the predicted angles determined with the leave-one-out model showed a mean (std.) of 4.16°(6.80°), with a 5 percentile of 0.41°and 95 percentile of 14.23°(Fig. 5).