A marker refers to a small spherical tool used in photogrammetry to acquire spatial information of a specific point in the human body. When a motion is performed after attaching numerous markers to the body, two or more high-speed cameras can be used to capture such motion. An analysis of body motion using the three-dimensional (3-D) positions of the markers acquired from a series of image frames is the most commonly used method in biomechanics.

The mechanical principle of assuming an object as a particle or a rigid body has been applied, as is, on human body motion. A multi- body system model that connects various segments with joints has also been applied to the human body. In addition, the position vectors of markers acquired from images have been used to construct mathe- matical models of variables of interest. The quality of mechanical analysis of human body motion depends on appropriate selection and com- bination of physical and mathematical models.

The observed movements of markers can be viewed as the effect of joint motion by muscle contraction being the primary cause. Here, joint motion is systematic or fixed effects. Moreover, the uncertainty and unavoidable error that occur during the measurement process (Taylor, 1997) can also be major causes that affect the measured values of the markers, which is referred to as random effects. A statistical model is defined as a model that contains mean and variance of resulting values from systematic and random effects, and the associated assumptions. Joint motions in the human body consist of rotational motion, and its elements such as angle, angular velocity, and angular acceleration are described by a transformation matrix. Rotational motions are centered on an axis of rotation (AoR) or center of rotation (CoR). Joint centers such as AoR and CoR are used in various ways in bio- mechanics, including calculation of momentum and defining a coordinate system, while they also represent a major point of interest in orthopedics for diagnosis of joint functions (Kim, 2013b). Therefore, a model for observed markers is more valuable when it is expressed as a functional relationship between joint center and transformation matrix. Touching or direct measurement of a joint center is not possible; thus, it can only be estimated by using the appropriate joint motion.

In biomechanics, some studies that estimate the joint center (Ehrig, Taylor, Duda, & Heller, 2006; Gamage & Lasenby, 2002; Piazza, Erdemir, Okita, & Cavanagh, 2004; Siston & Delp, 2006) can be considered as linear models. However, systematic modeling were not carried out because procedures for determining the uncertainty of estimation para- meters and goodness of fit of the model were not performed. A statistical model is valuable for concisely summarizing data as major systematic effects, and the size and nature of unexplained variations (McCullagh & Nelder, 1989). Therefore, an appropriate model for marker positions should be established and the model should be used to estimate the parameters of interest while also testing the uncertainty of the parameters and validity of the model.

The present study examined the statistical model of segment marker positions observed from images of human body joint motion. The objective of the present study was to examine linear models related to CoR estimation and to identify the strengths and weaknesses of the models for applying an appropriate general linear model that shows marker vectors in response to random and fixed effects that include CoR. The study also aimed to determine the goodness of fit of joint motion model for photogrammetry data by examining the validity of the model and the uncertainty of parameters estimated from assump- tions in the study model.

1. A linear model for markers

The kinematical model of observed markers attached to a distal segments with respect to a proximal segment coordinate system, including error e, can be expressed as shown below (Figure 1):

represents the 3-D position vector obtained from the image frame (j=1,⋯,n) of the marker (i=1,⋯,m). The term (omicron) denotes the CoR position vector expressed with respect to the proximal segment coordinate system, while is the position vector of the marker expressed in the distal segment coordinate system. Meanwhile, represents the transformation matrix that transforms the distal segment to the proximal segment components. Both proximal and distal coordinate systems are rectangular coordinate systems and are fixed to the segments and move together with those segments. When looking from the proximal segment toward the distal segment, rotational motion occurs at the center . If the human body segment is assumed to be a rigid body, and are vectors that are fixed in all phases of motion when looking from each coordinate system. In biomechanics, CoR is of greater interest than . and must be estimated from joint motion. Generally, 3 or more markers are attached to each segment to acquire the 3-D motion information of the segment.

Figure 1. A marker position on the distal segment with respect to the proximal segment coordinate system

Modify Delete

The origin of the statistical model for the segment markers can be found in the study of Chan (1965). If there are the measured coordinate data for a circle on a plane (j=1,⋯,n), then

where (ξ,  η) represents the center point of the circle, τ is the radius, and is the angle formed with the x-axis of the coordinate system. The terms and  are errors, which were assumed to be normally distributed with the mean and variance of 0 and , respectively; the error terms are independent of each other. Here, the parameters of interest are the center point (ξ, η) and radius τ.

It should be noted that all components in Equation 2 were expressed by using the same coordinate system components. If we represent the realized values of the random variables and as , the center of the circle as , and the distance from the center of the circle to the measured values as vector , there are two basic forms of the least squares methods that fit the circle.

‖⋅‖ is the norm of the vector. The first equation was proposed by Robinson (1961), and the second equation was by Delogne (1972) and Kåsa (1976). The first equation minimizes the normal distance between the data and the circle. This distance is a geometric distance; thus, the method is called the geometric fit method (Gander, Golub, & Strebel, 1994). Other methods are called algebraic fit method, which usually minimize the square of the algebraic expression representing with parameters (Kim, 2011). The second equation shown above falls into this category. The least-squares estimate of a geometric fit function is the same as the maximum likelihood estimate (Chan, 1965), and the solution must be derived by iteration method with nonlinear equations. On the other hand, the algebraic fit method can derive the solution directly from linear equations (closed form). The 2-D models of Equations 2, 3, and 4 can be expanded into a 3-D spherical model.

One should note two points about in Equation 1. First, in Equations 3 and 4 is a scalar, whereas is a vector. Second, unlike ο, , and , was expressed with respect to the distal coordinate system. For an analysis that assumes the human body to be a rigid body, may not be important as ο. However, if the human body is assumed to be a deformable body, not a rigid body, then takes on greater importance, and it must be expressed with components of a distal coordinate system, which is a local coordinate system (Figure 1). What is more important is that even in the mechanics of rigid body, expressing with components of a local coordinate system is more advantageous. Moreover, estimating as a vector value is valuable for analysis of 3-D human body motion.

When in Equation 1 is given, vectors ο and become unidentifiable. A method for resolving this is to transform to a distance value, that is, a scalar value.

Equation 5 comes from triangular inequalities, and A is an orthogonal matrix ; therefore, . If error is assumed to be small, then Equation 5 approaches

Meanwhile, minimizing the square distance, not the displacement of error vector, corresponds to Equations 3 and 4. Here, is independent to the coordinate system (). Vector cannot be estimated with Equations 3 and 4. Even if Equation 1 is approached by using a covariance structure analysis of the measurement model, it cannot be estimated. Since this is linear transformation, infinitely many solutions exist due to no appropriate constraints. Besides, Späth (1998) proposed a model that uses components of a spherical coordinate system, not a 3-D Cartesian coordinate system. However, even this model cannot describe rotational motion with three full degrees of freedom.

To estimate τ_i, which is expressed by a proximal segment coordinate system, a method called transformation method (Ehrig et al., 2006; Piazza et al., 2004; Siston & Delp, 2006) should be used. This is the segment coordinate system being given as transformation matrix in Equation 1. When and in Equation 1 are given, the equation can be viewed as a general linear model (Monahan, 2008; Ravishanker & Dey, 2002), and both ο and may be estimated. The transformation matrix may be comprised of Eulerian angles, directional cosines, and quaternions, in addition to the coordinate system of the segments (Ginsberg, 1998; Goldstein, 1980). The estimator of transformation method is an unbiased estimator with smaller variance than that of fitting method (Kim, 2013a). Thus, transformation method is recommended.

The following shows Equation 1 expressed as a matrix equation conventionally used in general linear models.

With y as the response variable,

where c refers to arranging by columns (McCulloch, Searle, & Neuhaus, 2008) and m number of markers consisting of n number of vectors were arranged to form a single column (3N × 1). In the same manner, if matrix and vectors are expressed as

then, X becomes a model matrix with a size of 3N × 3(m+1). If Kronecker product ⊗ is used, then

where is the vector , N number of components are all '1', and and are identity matrix, size of 3 and m, respectively. The term X can be called explanatory variables. β is a vector of para- meters related to the mean of y, where

Equation 7 represents a general linear model with joint center ο, and ,⋯, are the effects of m number of markers from rotational motion, where the expected value and variance of y is E(y)=Xβ and Var(y)=Σ, respectively. Therefore, the expected value and variance of error e becomes E(e)=0 and Var(e)=Σ. Σ is a 3(N×N) covariance matrix that has 3N(N+1)/2 of parameters that need to be estimated. Therefore, general assumptions are made on error e in a general linear model. These are an independence assumption that all N number of error vectors are independent to each other and a homoscedasticity assumption that all error vectors have the same variance of .

Transformation method estimates β under these assumptions. Equation 13 also contains the assumption that error vector components x,y,and z are independent each other while all have the same variance of . We revised such assumptions. The components x,y,and z are independent from each other, but they have different variances. Here, the error covariance matrix is

where V is

Σ is a diagonal matrix consisting of 3 parameters. If Σ is known, a generalized least-squares estimator of β can be used. However, it is unknown and we should estimate it. An ordinary least-squares estimator

is an unbiased estimator even under heteroscedasticity conditions (E(β)=β). However, variance is biased toward

and does not satisfy the efficiency or the consistency of the covariance matrix, which makes it inappropriate for estimation (Diggle, Heagerty, Liang, Zeger, 2002). Accordingly, the present study used two different methods to derive the variance of estimator . The first method was to estimate , , and separately for each component x,y,and z, as shown in Equation 15, while the second method was to use the hetero- scedasticity-consistent covariance matrix estimator by MacKinnon and White (1985). This is a method that derives the residual by using a generalized least-squares estimator and then taking it into consid- eration for the variance of . The estimation vector obtained from substituting for β in Equation 7 is

Here, is a projection matrix that transforms into an estimation space. If residual between y and is to be

then the variance of , by using covariance matrix obtained by substituting least-squares estimators ,,and of components , , and in Equation 15 into Equations 14 and 15, becomes

The variance of using HC3 (heteroscedasticity-consistency) by MacKinnon and White (1985) is as follows:

where, is

The term is a component of residual , while (k=1,⋯,3N) is a diagonal component of . HC3 is a robust method for deriving the standard error of estimation vectors, showing excellent results under various conditions (Long & Ervin, 2000). Residual variance is

The estimated covariance matrix , residual , and estimated residual variance were used to derive various types of residuals (Draper & Smith, 1998) for an analysis of model validation.

2. Empirical data analysis

For applying a general linear model to measured marker position vectors in human body joint motions, the present study chose the right hip joint. The hip joint has 3 degrees of freedom, consisting of flexion/ extension, adduction/abduction, and internal/external rotation. Joint center was estimated based on hip joint motions. A total of 5 male participants were selected for the experiment and their physical char- acteristics consisted of height of 179±2 cm, weight of 73±3 kg, and age of 20.6±.89 years. The participants did not have any history of musculoskeletal disease and voluntarily participated in the experiment. The experiment was conducted with the approval from the institutional review board at K University. A total of 7 markers were attached to anatomical landmarks in the participant's pelvis and femur as follows: two on the anterior left and right anterior superior iliac spine (ASIA), two on the posterior left and right posterior superior iliac spine (PSIS) on the pelvis, and one each on the greater trochanter (GT), lateral epicondyle (LE), and medial epicondyle (ME) on the femur. Moreover, 9 additional markers were attached, 3 each on the medial, anterior, and lateral sides of the femur. The total number of markers used in the experiment was 16, with 7 markers attached on landmarks used to define the coord- inate system for each segment and the remaining 9 markers used for analyzing the model used in the present study (Figure 2).

Figure 2. Movements of a hip joint for estimating the CoR and marker positions

Modify Delete

The participants performed 3 types of motions used in precedent studies (De Rosario, Page, Besa, & Valera, 2013; Ojeda, Martínez-Reina, & Mayo, 2014). The starting point of the motion was a neutral position. The motions consisted of 30° flexion - neutral position - 30° extension - neutral position as the first type of motion; 30° abduction - neutral position as the second type of motion; and 30° flexion - circum- duction to 30° extension point - neutral position as the third type of motion. The 3-D position vectors of the markers were acquired in centimeters by using the Optitrack Motion Capture System (Nationalpoint, USA). The camera used was Prime13, and 7 cameras were used at a sampling rate of 120 frames per second (FPS). The raw data for 3-D position vectors did not undergo any filtering process. The coordinate systems for the pelvis and femur were defined by using a method described by Wu et al. (2002), and GT was used instead of the hip joint center for the femur. The pelvis and femur coordination systems are represented in Figure 1 as proximal and distal segment coordinate systems, respectively. In the coordinate systems, the +x, +y, and +z axes represented anterior, superior, and right directions. The 9 markers (m=9) used in the analysis were transformed into components of the pelvis coordinate system, same as r_ijin Equation 1. MATLAB ver 6.1 (The Math Works) was used for data analysis and model-based estimation.

Estimated ο ̂ was the position vector of the hip joint center in the pelvis coordinate system, with the center point of the left ASIS (LASIS) and right ASIS (RASIS) as the origin, while was the position vector of the marker in the femur coordinate system located on the hip joint center (Figure 2). In particular, do not represent the positions of the markers attached to landmarks but rather the markers attached to the medial/lateral and anterior sides on the point where the femur was divided roughly into 3 parts; thus, the estimated positions varied by participants. Consequently, it did not appear in the results. Estimates of parameters and were as shown in Table 1.

The hip joint center estimate vector was slightly different for each participant due to difference in size. The ranges for x,y, and z axes were -2.423 ~ 2.264, -15.464 ~ -12.187, and 11.066 ~ 12.153, respectively, meaning the x axis showed the highest deviation. The ordinary least-squares variance estimate of was .983~3.063, showing differences between the participants. Among most of the variance estimate for each component, the x axis component was found to be the biggest. In participants 1 and 5, the size appeared in the order of x > z > y, whereas the order was x>y>z in participants 2, 3, and 4. The standard error of from three different methods (,,) showed similar results.

With respect to the residuals for diagnosis of the model, the five participants showed mostly similar patterns, and the analysis used the residual of participant 5. In Figure 3, the graphs showed the residuals for the x,y,and z components, where the left side showed the residuals with fitted values and the right side the residuals with time. On the right side, non-continuous forms appeared near frames 360 and 550, which was due to all three motions being used in the analysis, these points were where these motions were distinguished.

Figure 3. Plots of standardized residuals with fitted values and frames

Modify Delete

The figure showed that the residuals were not scattered, but rather, they showed a continuous form, and such pattern was more prominent in the images on the right. This can be definitely viewed as a data structure with autocorrelation between the measured values. With respect to homoscedasticity, as shown in the images on the right, x,y,and z components are determined to show homoscedasticity, but the size of the residuals varied according to the motion performed. It can be surmised that the size may vary according to the type of motion and marker position.

Figure 4 showed a box-and-whisker and Q-Q plots of the residuals. On the left, x showed the largest deviation, while outliers appeared on the +x and –y directions. On the right for normality testing, x,y,and z all showed a symmetrical distribution relative to 0 and showed a shape with a slightly longer tail than the normal distribution. The lengths of the tails seemed to be short, and the distribution was balanced on both sides; thus, we believed that the influence on the estimates may be ignored. For further evaluating the variance of the residuals, graphs were plotted, as shown in Figure 5, with the observed and estimated vectors divided into x,y,and z components. The x axis was an anterior/posterior axis, which showed the largest displacement in flexion/extension motion. From the graphs, large differences between the observed and estimated vectors were found at the phases of the directional changes (slope = 0), which was believed to be the primary cause of the variance estimates appearing large.

Figure 4. Box and whisker and normal plots of the residuals of each of the x, y, z components

Modify Delete

Figure 5. Observed and fitted positions for 9 markers of each of the x,y,zcomponents. Thin lines (━) are observed positions, and thick lines (━) are fitted positions. L, A, and M indicate lateral, anterior, and medial, respectively (Figure 2).

Modify Delete

A statistical tool is needed to solve the problem of estimating para- meters of interest with measured values that include errors. Solving a given problem by using statistics is a difficult task with unavoidable errors, but as long as it is expressed in statistical terms, the solution is predetermined (Faraway, 2005). The synovial joint, which is responsible for most human body movements, is mostly assumed to be a structure where rotational motion is present. A rotational motion may be described by a rotational matrix, while CoR, where instantaneous AoR that generate rotational motion pass through, is significant in biomechanics. Moreover, human body segments are generally assumed to be rigid bodies; thus, the relative distance between markers attached to these segments is considered a constant value that does not change during motion. The significance of the present study can be found in the fact that it combined the measured values that include such errors and characteristics of the human body with statistical tools in expressing the relationship between the mean marker positions observed during joint motion and CoR as a linear model for the analysis.

The model used in the present study differed slightly from general linear models in two aspects. First, the explanatory variable was the same for all markers (Equations 7 and 11). Based on the assumption of rigid body, the position of the marker attached to a distal segment is fixed, and all markers perform the same rotational motion due to the nature of a rigid body; thus, the same transformation matrix is applied. Second, the number of parameters being estimated increases according to the number of markers used (Equation 12). It is possible to estimate ο and τ with a single marker. However, 3 or more markers are generally used in experiments to acquire information on translation and rotational motions of a segment. Here, additional markers are additional parameters that need to be estimated, and the structure has a different value of , depending on where the markers are attached. If too many markers are used, then numerical analysis involving estimation becomes inefficient.

In biomechanics, ordinary least squares are used often (Kim, 2010). When no association between the measured vectors and the same variance are assumed, according to the Gauss-Markov theorem, the ordinary least-squares estimator is the best linear unbiased estimator (BLUE). The present study assumed that including vector components x,y,and z, m number of markers and n number of measures vectors within-markers were independent from each other. Moreover, com- ponents x,y,and z of an orthogonal coordinate system were assumed to have different variances from each other. This was based on the belief that the variance of x,y,and z components would vary according to the motion being performed. Consequently, the sizes of the variances differed. However, this did not have an effect on the standard error of estimated vectors of the joint center. As shown in Figure 3, the number of frames was approximately 900, and the number of frames that can be considered as sample is believed to be large. In the present study, images were acquired with a sampling rate of 120 FPS. If the sampling rate is increased, then the amount of data also increases and it will take more time to calculate. Therefore, deciding on an appropriate sampling rate can also have an impact on the validity of a study.

As seen in residual plots, data were continuous and showed a form that changed according to the motion performed. This was because the values of the image data measured continuously over time were longitudinal data that were not independent from each other. We believe that estimations should be performed with taking into account such data characteristics.

The present study applied a general linear model for estimating of CoR during joint motion to investigate the uncertainty of estimation vectors and validity of the model, based on which the following con- clusions were reached. The standard errors of the estimated vectors did not appear to be large among the models that assumed homosced- asticity and heteroscedasticity of errors. Residuals showed a pattern of changing according to motion, appearing especially large when the direction changed. Therefore, we believe that the current model is limited and should have the specific structure for variance-covariance matrix. Hence, consideration should be given to revising the model to a mixed-effects model that assumes markers with random effects or a model that explains the effects of soft tissue artifacts, which represent movement of soft tissues such as the skin, generated when motion is created.