A Simple Algorithm for Assimilating Marker-Based Motion Capture Data During Periodic Human Movement Into Models of Multi-Rigid-Body Systems

1. Introduction

Human movement analysis is performed in various fields, including biomechanics, physiology, orthopedics, neurology, and sports science, playing an important role for understanding physical functions, motor control, and motor dysfunctions (e.g., see Harris and Smith, 1996; Winter, 2009; Lu and Chang, 2012 for a review). Optoelectronic stereo-photogrammetry, referred simply to as the motion capture system, is the most popular method used in the human movement analysis, in which spatio-temporal changes in positions of reflective markers placed on anatomical landmarks of the human body, allow us to describe the body movement quantitatively in the computer. In this process, the human body is often modeled by a multi-link-rigid-body system, and reflective-marker-based motion capture data are assimilated into the model for characterizing kinematics and kinetics of the movements quantitatively as well as mechano-dynamic analysis of the movement such as the inverse dynamics analysis to estimate joint torque profiles along the movement (Sibella et al., 2003; Ren et al., 2008; Alexander and Schwameder, 2016).

A goal of the data assimilation process is to make the motions captured from an experimental subject physically consistent with those of the model, while they are as close as to the actual motions of the subject (Cappozzo et al., 2005). This is because, if no assimilation processing is performed, the original, non-assimilated motion data exhibit several inconsistencies with the model, such as temporal variations in the length of each link, which should not happen under the rigid-body-model assumption. The underlying cause lowering an accuracy of the assimilation is the fact that the human body is not a simple multi-rigid-body system: skeletal limbs are covered by deformable muscular-tissues, the trunk, foot, etc., are composed of a number of small bones despite a common modeling simplification with single or a few links, and no joints are rotary nor spherical. Nevertheless, the multi-rigid-body modeling is still useful for the quantitative human movement analysis (Dicharry, 2010; Rusaw and Ramstrand, 2011). Moreover, for practical applications, it is preferable that the accurate data assimilation can be performed with a small number of markers used for the motion capture (Simon, 2004).

In this study, we consider the minimum numbers of markers necessary for identifying a position and a posture of each link, and movements of the markers on the skin relative to the underlying skeletal system, referred to as the soft tissue artifact (STA), as one of the major causes that lower accuracy of the assimilation (Camomilla et al., 2017). STA has its origin in wobbling of the soft tissue, stretching of the skin, and activation of the muscles. STA-characteristics have been investigated by comparing movements of markers attached on the skin with those of markers attached to intra-cortical pins fixed directly into the bones (Lafortune et al., 1992; Benoit et al., 2006; Andersen et al., 2012; Blache et al., 2017), to a splint fixed to the bone (Cappozzo et al., 1996), and to the Percutaneous Skeletal Tracker (Holden et al., 1997). Trials for measuring movements of the bone using X-ray (Sati et al., 1996; Li et al., 2017) and MRI have been also carried out (Ryu et al., 2006). Moreover, a number of computational methodologies have been proposed to minimize and/or compensate the effects of STA, which include the multiple anatomical landmark calibration (Cappello et al., 1997, 2005), the dynamic calibration (Lucchetti et al., 1998), the point cluster technique (Andriacchi et al., 1998; Alexander and Andriacchi, 2001), the multi body optimization (Lu and O'Connor, 1999; Yoshikawa et al., 2012, 2013; Clément et al., 2015; Richard et al., 2017), approaches with Kalman filters (Cerveri et al., 2003, 2005; Halvorsen et al., 2008; Bonnet et al., 2017b), an optimization using the ground reaction forces (Riemer et al., 2008), the linear 3D interpolation and approximation methods (Dumas and Cheze, 2009), among others. However, reliable yet simple assimilations of skeletal motions have not been achieved satisfactorily (Leardini et al., 2005; Peters et al., 2010; Camomilla et al., 2017).

Primary causes of the difficulty to minimize and/or compensate the effects of STA are associated with the fact that the patterns of the artifacts are task-dependent, although the STA is reproducible within a given specific task, but not among subjects (Leardini et al., 2005). Since most of the methods proposed in literatures so far need to be task-adjusted, they require a precise calibration and tuning in numerical compensations specifically for each subject and each task. Thus, it would be very useful if an algorithm focuses only on the relationship between STA and skeletal motions in common across movements, independent of tasks.

Previously, we studied motions of markers attached on the anatomical landmarks during human walking, and showed that STA profiles were also periodic (Inoue et al., 2016). In this paper, we propose a novel and simple algorithm to assimilate motion-captured periodic human movements, represented by positions of markers, that are affected by STA into models of multiple rigid-link systems. In section 2, we describe the proposed algorithm using a simple example of a two-link model that moves in the two-dimensional space. We then performed a numerical experiment using a model of human walking to assess the efficiency of the algorithm, for which methods for the simulation are summarized in section 3. Section 4 reports the results of the numerical experiment. This is a preliminary study that aims to develop a theoretical foundation of the algorithm in three-dimensional system. We discuss whether the algorithm can be applied to posture estimation in three-dimensional space in section 5.

2. The Proposed Algorithm

2.1. Problem Setting

In this section, we consider temporal changes in a posture of a planar two-link system that moves periodically. This is just for simplicity, and the algorithm proposed here can be applied to systems with more links as shown later in this paper, and could extended for three-dimensional systems. The two-link model consists of two rigid links, referred to as the link-A and the link-B as shown in Figure 1A, which are connected to each other by a pin-joint. We consider a situation such that a soft tissue (the gray area in Figure 1A) surrounds the rigid-link system. The minimum number of markers necessary for specifying the position and posture of each link is two for the planar system. Thus two markers, referred to as the marker-i1 and the marker-i2 for i={A,B}, are attached on the surface of the soft tissue. We assume that those markers are attached accurately to the landmarks of each link (the landmarks-A1 and A2 for the link-A and the landmarks-B1 and B2 for the link-B) when the system is at rest without any movement. However, each marker may excurse from the landmark, causing STA during periodic movements of the system.

FIGURE 1

Figure 1. Two-link model, which is used to explain posture estimation algorithm. The two-link model consists of two rigid links, and surrounded by soft tissue. In each panel, blue circles represent markers located exactly on the landmarks of the two-link model, and red circles represent markers captured by motion capture system, which excurse due to STA. (A) Conceptual diagram of two-link model. Brown rectangular areas are rigid links, and gray shaded region is soft tissue. (B) Position vectors of landmark markers [${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}$, i = {A, B}, j = {1, 2}], and captured markers [${\text{m}}_{ij}^{(\mathrm{\text{G}})}$]. (C) Position vectors of local coordinate system [the link coordinate (${\bar{\text{o}}}_i^{(\mathrm{\text{G}})}$) and the landmark coordinate (${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}$)] and coordinate axes of those in global coordinate system [$x_i^{(\mathrm{\text{L}})}$-$y_i^{(\mathrm{\text{L}})}$ and $x_i^{(\widetilde{\mathrm{\text{M}}})}$-$y_i^{(\widetilde{\mathrm{\text{M}}})}$]. (D) Position vectors of the marker coordinate system in global coordinate system [${\text{o}}_i^{(\mathrm{\text{G}})}$] and in the landmark coordinate system [${\text{o}}_i^{(\widetilde{\mathrm{\text{M}}}-i)}$], and coordinate axes of that [$x_i^{(\mathrm{\text{M}})}$-$y_i^{(\mathrm{\text{M}})}$].

We denote by ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ the spatio-temporal position of the marker-ij in the global coordinate system, where the superscript (G) represents the global coordinate X-Y as in Figure 1B for i={A,B} and j={1,2}. The integers n = 1, 2, ⋯  indicate the data numbers or the sampling times. ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ corresponds to an experimentally motion-captured positions of the marker-ij contaminated by unknown STA. Now, the positions of the markers located exactly on the landmarks are denoted as ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ (Figure 1B), which cannot be measured experimentally, but can be estimated if we can estimate the unknown STA time-profiles as ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]={\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]-\mathrm{\text{STA}}$. The problem here can be stated as follows: the motion capture experiments provide time-series data of ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$, from which we would like to estimate the positions and postures of the link-A and the link-B at each time instant so that they are consistent with the assumption for the two-rigid-link model. This may be possible if we can estimate the STA profiles and then ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$, as described in this sequel.

2.2. Definitions of Marker-Positions and the Multi-Link System

We define a local coordinate system, referred to as the link-coordinate system, which is fixed on the link-i as in Figure 1C to describe the position and posture of the link-i. The link-coordinate system of the link-i is denoted by $x_i^{(\mathrm{\text{L}})}$-$y_i^{(\mathrm{\text{L}})}$, where the superscript (L) represents the link-coordinate system. Position of the origin and the posture of the link-coordinate changes in the global coordinate system together with translation and rotation of the link-i. The origin of the link-coordinate system ${\bar{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ represents the global position of the link-i, and the coordinate transformation matrix ${\bar{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$, consisting of normalized orthogonal basis vectors in the directions of $x_i^{(\mathrm{\text{L}})}$- and $y_i^{(\mathrm{\text{L}})}$-axes, represents the posture of the link-i. We set the link-coordinate system of the link-i such that the origin is located at the centroid of the triangle having the vertices at the joint and two landmarks (${\tilde{\text{m}}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{m}}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]$). As in Figure 1C, the $y_i^{(\mathrm{\text{L}})}$-axis of the link-coordinate system of the link-i is defined so that it is parallel to the vector directing from ${\tilde{\text{m}}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]$ to ${\tilde{\text{m}}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]$, and then the $x_i^{(\mathrm{\text{L}})}$-axis is defined so that it becomes orthogonal to the $y_i^{(\mathrm{\text{L}})}$-axis¹. It is essential for the proposed algorithm to notice that ${\bar{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\bar{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ of the link-coordinate system defined here cannot be determined easily from the motion-captured marker data ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$, since three vertices at the joint, ${\tilde{\text{m}}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{m}}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]$ are all unknown. That is, the maker positions ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ may show temporal deviations from the landmark positions ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ due to unknown STA, and thus the landmark positions ${\tilde{\text{m}}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{m}}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]$ are not available directly from the marker positions. Obviously, the joint position needs to be estimated somehow from the STA-affected marker positions.

If all markers were strictly fixed on the landmarks of each link, i.e., if the motion-captured marker positions were not affected by STA, the link-coordinate systems could be estimated easily, since the landmark-marker positions ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ are exactly the same as the motion-captured marker positions ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$, and the position of the joint could be calculated as described below.

2.2.1. Estimation of the Joint Position From Non-STA-Affected Markers

For estimating the joint position in the global coordinate system, we define another local coordinate system based on the landmark positions ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ (Figure 1C), referred to as the landmark-coordinate system and denoted by $x_i^{(\widetilde{\mathrm{\text{M}}})}$-$y_i^{(\widetilde{\mathrm{\text{M}}})}$. The origin of the landmark-coordinate system in the global coordinate system, ${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ for the link-i, is set as the midpoint of the landmark-i1 and the landmark-i2, and the $x_i^{(\widetilde{\mathrm{\text{M}}})}$-axis of the landmark-coordinate system relative to the global coordinate system, which determines the posture matrix ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$, is set as the normalized vector directing from the landmark-i1 to the landmark-i2:

$\begin{array}{cc}\begin{array}{rcl}{\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]& =& \frac{{\tilde{\text{m}}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]+{\tilde{\text{m}}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]}{2},\end{array}& (1)\end{array}$

$\begin{array}{cc}{\tilde{A}}_i^{(\text{G})}[n]=\frac{1}{\left|{\tilde{m}}_{i\text{-diff}}^{(\text{G})}[n]\right|}\left[\begin{array}{cc}{\tilde{m}}_{i\text{-diff}}^{(\text{G})}[n]& \underset{\left|}{\overset{\right|}{{\displaystyle |}}}{R}_{\frac{\pi }{2}}{\tilde{m}}_{i\text{-diff}}^{(\text{G})}[n]\end{array}\right],& (2)\end{array}$

where

$\begin{array}{l}{\tilde{m}}_{i\text{-diff}}^{(\text{G})}[n]={\tilde{m}}_{i2}^{(\text{G})}[n]-{\tilde{m}}_{i1}^{(\text{G})}[n],\\ R_{\frac{\pi }{2}}=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right].\end{array}$

The joint position in the landmark-coordinate system of the link-i is constant in time during any movement of the model, because the skeletal part of the link-i (each of the brown rectangular areas in Figure 1A) is assumed to be rigid body and the joint is pinning such two rigid bodies. Thus, we denote it by a time-independent constant vector ${\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, where the superscript ($\widetilde{\mathrm{\text{M}}}$-i) means that the vector is represented in the landmark-coordinate system for the link-i. The joint position in the global coordinate system, denoted as the time-varying vector ${\tilde{\text{j}}}_i^{(\mathrm{\text{G}})}\left[n\right]$, can be described by the position and the posture of the landmark-coordinate system of the link-i as follows:

$\begin{array}{cc}\begin{array}{r}{\tilde{\text{j}}}_i^{(\mathrm{\text{G}})}\left[n\right]={\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]+{\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]{\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}.\end{array}& (3)\end{array}$

Since the link-A and the link-B are connected by the joint at a single point in the global coordinate system, referred to as the ${\tilde{\text{j}}}^{(\mathrm{\text{G}})}\left[n\right]$, the joint position described in the link-A and that in the link-B by Equation (3) must be equal to each other for any time instant, i.e.,

$\begin{array}{cc}\begin{array}{r}{\tilde{\text{j}}}^{(\mathrm{\text{G}})}\left[n\right]={\tilde{\text{o}}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]+{\tilde{\text{A}}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]{\tilde{\text{j}}}_{\mathrm{\text{A}}}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-A}})}={\tilde{\text{o}}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]+{\tilde{\text{A}}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]{\tilde{\text{j}}}_{\mathrm{\text{B}}}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-B}})}.\end{array}& (4)\end{array}$

This is the joint constraint that assures the two adjacent links are connected at a single joint. The two unknown constant vectors ${\tilde{\text{j}}}_{\mathrm{\text{A}}}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-A}})}$ and ${\tilde{\text{j}}}_{\mathrm{\text{B}}}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-B}})}$, representing the joint position in the landmark-coordinate system of the link-A and the link-B, can be obtained by solving Equation (4) at any two different time instances n1 and n2 as

$\begin{array}{cc}\left[\begin{array}{c}{\tilde{j}}_{\text{A}}^{(\tilde{\text{M}}\text{-A})}\\ {\tilde{j}}_{\text{B}}^{(\tilde{\text{M}}\text{-B})}\end{array}\right]={\left[\begin{array}{cc}{\tilde{A}}_{\text{A}}^{(\text{G})}[n1]& -{\tilde{A}}_{\text{B}}^{(\text{G})}[n1]\\ {\tilde{A}}_{\text{A}}^{(\text{G})}[n2]& -{\tilde{A}}_{\text{B}}^{(\text{G})}[n2]\end{array}\right]}^{-1}\left[\begin{array}{c}{\tilde{o}}_{\text{B}}^{(\text{G})}[n1]-{\tilde{o}}_{\text{A}}^{(\text{G})}[n1]\\ {\tilde{o}}_{\text{B}}^{(\text{G})}[n2]-{\tilde{o}}_{\text{A}}^{(\text{G})}[n2]\end{array}\right].& (5)\end{array}$

Then the joint position in the global coordinate system, ${\tilde{\text{j}}}^{(\mathrm{\text{G}})}\left[n\right]$, can be obtained by Equation (3). In this way, using the joint position ${\tilde{\text{j}}}^{(\mathrm{\text{G}})}\left[n\right]$ in the global coordinate system and the positions of two landmarks for each link ${\tilde{\text{m}}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{m}}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]$ that are experimentally available in the special case with ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]={\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$, we can determine the position ${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and the posture ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ of the link-i, completing the assimilation almost directly.

However, in practice, unlike the case without STA considered here, the motion-captured marker positions are not the same as those of the landmark positions.

2.3. A Naive Algorithm for Simply Minimizing the Joint-Constraint-Errors

Previously, we studied a preliminary algorithm to estimate skeletal movement of human body from captured positions of markers, which are attached to anatomical landmarks (Yoshikawa et al., 2012, 2013), referred to as the naive algorithm. Performances of the naive algorithm would be compared with those of the new algorithm proposed in this paper. The naive algorithm could reduce the influences of STA effectively through a multi-body optimization process. However, the obtained solutions (the estimated motions) do not necessarily coincide with the actual motions of the multi-link system. In this paper, we propose a novel but still simple algorithm to better assimilate the motion-captured data into models of multi-rigid-links during human periodic motion by estimating the unknown STA profiles.

The naive algorithm does not utilize a notion of the landmarks and the landmark-coordinate systems, and the local coordinate systems are defined based on the positions of the marker-i1 and the marker-i2 in a similar way to the landmark-coordinate system by removing the tilde-signs from the equations, because the tilde-signs are reserved for the landmark-markers. We refer to this “pseudo” local coordinate system as the marker-coordinate system, where the origin and the posture of the coordinate system (${\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]$) are defined as follows.

$\begin{array}{cc}\begin{array}{rcl}{\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]& =& \frac{{\text{m}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]+{\text{m}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]}{2},\end{array}& (6)\end{array}$

$\begin{array}{cc}{A}_i^{(\text{G})}[n]=\frac{1}{\left|m_{i\text{-diff}}^{(\text{G})}[n]\right|}\left[\begin{array}{cc}{m}_{i\text{-diff}}^{(\text{G})}[n]& \underset{\left|}{\overset{\right|}{{\displaystyle |}}}{R}_{\frac{\pi }{2}}{m}_{i\text{-diff}}^{(\text{G})}[n]\end{array}\right],& (7)\end{array}$

where

$\begin{array}{l}{\text{m}}_{i\mathrm{\text{-diff}}}^{(\mathrm{\text{G}})}\left[n\right]={\text{m}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]-{\text{m}}_{i1}^{(\mathrm{\text{G}})}\left[n\right].\end{array}$

The marker-coordinate system of the link-i cannot be fixed in the link, due to STA, which is why we call it “pseudo” local coordinate system. Thus, the marker-coordinate system of the link-i may exhibit temporal deviations from the landmark-coordinate system, and the relative position between the marker-coordinate system of link-i (${\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]$) and the actual joint position (${\tilde{\text{j}}}^{(\mathrm{\text{G}})}$) changes in time. Therefore, unlike in the landmark-coordinate system, the joint position vector in the marker-coordinate system ${\text{j}}_i^{(\mathrm{\text{M-}}i)}$ inevitably changes in time. Nevertheless, in the naive algorithm, the joint position in the marker-coordinate system of the link-i is estimated as the time-independent constant vectors ${\text{j}}_i^{(\mathrm{\text{M-}}i)}$ (i=A and B) that minimizes the cost function as

$\begin{array}{cc}\begin{array}{r}{\displaystyle \underset{{\text{j}}_{\mathrm{\text{A}}}^{(\mathrm{\text{M-A}})},{\text{j}}_{\mathrm{\text{B}}}^{(\mathrm{\text{M-B}})}}{argmin}}{\displaystyle \sum _n}|({\text{o}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]+{\text{A}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]{\text{j}}_{\mathrm{\text{A}}}^{(\mathrm{\text{M-A}})})-({\text{o}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]+{\text{A}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]{\text{j}}_{\mathrm{\text{B}}}^{(\mathrm{\text{M-B}})})|.\end{array}& (8)\end{array}$

In many cases, there was no pair of constant vectors for the joint positions (${\text{j}}_{\mathrm{\text{A}}}^{(\mathrm{\text{M-A}})}$ and ${\text{j}}_{\mathrm{\text{B}}}^{(\mathrm{\text{M-B}})}$) that make the cost function zero. In other words, the optimal joint position obtained for the time-varying marker-coordinate system for the link-A could not coincide to that for the link-B for all time instances. Thus, in the naive algorithm, the joint position j^(G)[n] in the global coordinate system is estimated as the midpoint of the optimal joint positions obtained for each of the two marker-coordinate systems as

$\begin{array}{cc}\begin{array}{r}{\text{j}}^{(\mathrm{\text{G}})}\left[n\right]=\frac{\left({\text{o}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]+{\text{A}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]{\text{j}}_{\mathrm{\text{A}}}^{(\mathrm{\text{M-A}})}\right)+\left({\text{o}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]+{\text{A}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]{\text{j}}_{\mathrm{\text{B}}}^{(\mathrm{\text{M-B}})}\right)}{2}.\end{array}& (9)\end{array}$

The position and the posture of each link are then estimated by using j^(G)[n], ${\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]$, and ${\text{A}}_{\mathrm{\text{i}}}^{(\mathrm{\text{G}})}\left[n\right]$. In this way, the marker positions and the joint position are both affected by STA in the naive algorithm even after the assimilation process, although the assimilated motions of the two-rigid-link model maximally satisfy the joint constraint.

2.4. The Proposed Algorithm for Periodic Movements

In the algorithm proposed here, unknown positions of the landmark markers and the joint in the global coordinate system are estimated based on the evidence that STA varies periodically during periodic movement, and then the position and the posture of each link are estimated. That is, according to the knowledge that STA profiles during a given periodic motion are also periodic (Inoue et al., 2016), the unknown STA components of the motion-captured marker positions in the unknown landmark-coordinate system are formally represented by a Fourier series. Interestingly, we show that the corresponding Fourier coefficients for the unknown STA can be determined based only on the joint constraint and the periodicity in the distance between two STA-affected markers.

The Fourier coefficients of STA in the landmark-coordinate system is defined first. The STA-affected marker positions in the landmark-coordinate system is denoted as ${\text{m}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$. Note that ${\text{m}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ is not the same as ${\tilde{\text{m}}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ that is the position of non-STA-affected landmark marker in the landmark coordinate system. ${\text{m}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ can be described formally as

$\begin{array}{cc}\begin{array}{r}{\text{m}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]={(-1)}^{g_j}\left[\begin{array}{c}\frac{C_i}{2}\\ 0\end{array}\right]+{\text{e}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right],\end{array}& (10)\end{array}$

with

$g_j=\{\begin{array}{ll}1,\hfill & \text{for }j=1\hfill \\ 2,\hfill & \text{for }j=2\hfill \end{array},$

where C_i is the distance between two landmark markers on the link-i, and ${\text{e}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ are the STA profiles represented in the landmark-coordinate system, although we do not know the landmark positions and STA for any instant of time n. The first term on the right-hand-side of Equation (10) represents the fixed position of the landmark-ij in the landmark-coordinate system of the link-i. Since the origin of the landmark-coordinate system is located at the midpoint of two landmarks (Equation 1) and the $x_i^{(\widetilde{\mathrm{\text{M}}})}$-axis directs from the landmark marker-i1 to the landmark marker-i2 (Equation 2), the absolute value of the x-component of the first term is one-half of the distance between two landmarks. The y-component of the first term is zero by definition of the $y_i^{(\widetilde{\mathrm{\text{M}}})}$-axis of the coordinate system.

We expand the STA for the marker-ij in the landmark-coordinate system into the following Fourier series:

$\begin{array}{cc}\begin{array}{l}{e}_{ij}^{(\tilde{\text{M}}\text{-}i)}[n]=\left[\begin{array}{c}{\displaystyle {\sum }_{k=1}^K\left\{a_{ij,x,k}\cos \frac{2\pi kn}{N}+b_{ij,x,k}\sin \frac{2\pi kn}{N}\right\}}\\ {\displaystyle {\sum }_{k=1}^K\left\{a_{ij,y,k}\cos \frac{2\pi kn}{N}+b_{ij,y,k}\sin \frac{2\pi kn}{N}\right\}}\end{array}\right]\\ \equiv \left[\begin{array}{cc}P[n]& 0\\ 0& P[n]\end{array}\right]\left[\begin{array}{c}{q}_{ij,x}^{(\tilde{\text{M}}\text{-}i)}\\ q_{ij,y}^{(\tilde{\text{M}}\text{-}i)}\end{array}\right],\end{array}& (11)\end{array}$

with

$\begin{array}{lll}\text{P}\left[n\right]& =& \left[\begin{array}{cccccc}cos\frac{2\pi 1n}{N}& \dots & cos\frac{2\pi Kn}{N}& sin\frac{2\pi 1n}{N}& \dots & sin\frac{2\pi Kn}{N}\end{array}\right],\\ {\text{q}}_{ij,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}& =& {\left[\begin{array}{cccccc}{a}_{ij,x,1}& \dots & a_{ij,x,K}& b_{ij,x,1}& \dots & b_{ij,x,K}\end{array}\right]}^{\mathrm{\text{T}}},\\ {\text{q}}_{ij,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}& =& {\left[\begin{array}{cccccc}{a}_{ij,y,1}& \dots & a_{ij,y,K}& b_{ij,y,1}& \dots & b_{ij,y,K}\end{array}\right]}^{\mathrm{\text{T}}},\end{array}$

where N represents the one cycle data length of the periodic motion, K is the order of Fourier series expansion, P[n] is the row vector composed of the cosine and the sine basis functions, and {a_{ij, x, k}, b_{ij, x, k}} and {a_{ij, y, k}, b_{ij, y, k}} are the Fourier coefficients that define the column vectors ${\text{q}}_{ij,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{ij,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$.

Let ${\text{o}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ and ${\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ be the position and the posture of the marker-coordinate system for the link-i, respectively (Figure 1D), relative to the landmark-coordinate system (${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$). Using the definition of the marker-coordinate system (Equations 6 , 7) and Equations (10, 11), the excursion of the origin and the tilt of axes of the marker-coordinate system in the landmark-coordinate system, i.e., ${\text{o}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ and ${\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$, can be expressed by using the Fourier-expanded STA as follows.

$\begin{array}{cc}{o}_i^{(\tilde{\text{M}}\text{-}i)}[n]=\frac{1}{2}\left[\begin{array}{cc}P[n]& 0\\ 0& P[n]\end{array}\right]\left[\begin{array}{c}{q}_{i1,x}^{(\tilde{\text{M}}\text{-}i)}+q_{i2,x}^{(\tilde{\text{M}}\text{-}i)}\\ q_{i1,y}^{(\tilde{\text{M}}\text{-}i)}+q_{i2,y}^{(\tilde{\text{M}}\text{-}i)}\end{array}\right],& (12)\end{array}$

$\begin{array}{cc}{A}_i^{(\tilde{\text{M}}\text{-}i)}[n]=\left[\begin{array}{rr}\hfill \cos \theta [n]& \hfill -\sin \theta [n]\\ \hfill \sin \theta [n]& \hfill \cos \theta [n]\end{array}\right],& (13)\end{array}$

where

$\begin{array}{cc}\begin{array}{r}\theta \left[n\right]={tan}^{-1}\frac{\text{P}\left[n\right]\left({\text{q}}_{i2,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}-{\text{q}}_{i1,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\right)}{C_i+\text{P}\left[n\right]\left({\text{q}}_{i2,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}-{\text{q}}_{i1,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\right)}.\end{array}& (14)\end{array}$

In this sequel, we show that C_i and the Fourier coefficients of ${\text{e}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ (i.e., STA in the landmark-coordinate system), ${\text{q}}_{ij,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{ij,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, can be determined based only on the joint constraint and the periodicity-assumption of STA. This means that we can obtain STA-free, corrected time-courses of the marker positions, as ${\text{m}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]-{\text{e}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ in the landmark-coordinate system, which can be transformed into the landmark-marker positions ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ in the global coordinate system. Once we determine ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ successfully, the position and the posture of each link (${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$), the joint positions ${\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ in the landmark-coordinate system, and the joint position ${\tilde{\text{j}}}^{(\mathrm{\text{G}})}\left[n\right]$ in the global coordinate system can be obtained using Equations (1, 2, 5), and Equation (3), respectively.

2.4.1. The Marker-Coordinate System Relative to the Landmark-Coordinate System

Any local vector u^(M-i)[n] in the marker-coordinate system of the link-i can be represented as

$\begin{array}{cc}\begin{array}{r}{\text{u}}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]={\text{o}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]+{\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]{\text{u}}^{(\mathrm{\text{M-}}i)}\left[n\right],\end{array}& (15)\end{array}$

in the landmark-coordinate system, which is also represented in the global coordinate as

$\begin{array}{cc}{u}^{(\text{G)}}[n]={\tilde{o}}_i^{(\text{G})}[n]+{\tilde{A}}_i^{(\text{G})}[n]o_i^{(\tilde{\text{M}}\text{-}i)}[n]+{\tilde{A}}_i^{(\text{G})}[n]A_i^{(\tilde{\text{M}}\text{-}i)}[n]u^{(\text{M-}i)}[n]& (16)\end{array}$

or, alternatively

$\begin{array}{cc}{u}^{(\text{G)}}[n]=o_i^{(\text{G})}[n]+A_i^{(\text{G})}[n]u^{(\text{M-}i)}[n],& (17)\end{array}$

by the coordinate transformation directly from the marker-coordinate system (with its origin at ${\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and the posture ${\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]$ in the global coordinate system) to the global coordinate system.

Comparing the first and the second terms of the right-hand side of Equation (16) with the first term of the right-hand side in Equation (17), we have

$\begin{array}{cc}\begin{array}{r}{\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]={\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]+{\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]{\text{o}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right],\end{array}& (18)\end{array}$

and comparing the third term of the right-hand side in Equation (16) with the second term of the right-hand side in Equation (17), we have

$\begin{array}{cc}\begin{array}{r}{\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]={\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]{\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right].\end{array}& (19)\end{array}$

By rearranging Equations (18) and (19), the landmark-coordinate system in the global coordinate system can be related to the marker-coordinate system as

$\begin{array}{cc}\begin{array}{rcl}{\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]& =& {\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]-{\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]{\left({\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]\right)}^{-1}{\text{o}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right],\end{array}& (20)\end{array}$

$\begin{array}{cc}\begin{array}{rcl}{\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]& =& {\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]{\left({\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]\right)}^{-1}.\end{array}& (21)\end{array}$

Equations (20) and (21) imply that we can obtain the landmark-coordinate system in the global coordinate system, if, somehow, we can obtain the position and the posture of the marker-coordinate system relative to the landmark-coordinate system, since ${\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]$ for the marker-coordinate system are available from the motion-captured marker positions using Equations (6) and (7).

2.4.2. The Cost Function for Ensuring the Joint Constraint

To estimate the Fourier coefficient of STA, ${\text{q}}_{ij,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{ij,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, a cost function to ensure the joint constraint, such as in Equation (8) for the naive algorithm, is considered also in the proposed algorithm. We show that, through the optimization process for the cost function, a time-invariant constant vector ${\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ representing the joint position in the landmark-coordinate system for each link can also be obtained simultaneously together with ${\text{q}}_{ij,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{ij,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$.

Using Equations (20) and (21), the joint position for the link-i in the global coordinate system defined by Equation (3) can be rewritten with the constant joint position vector in the landmark-coordinate system ${\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ as

$\begin{array}{cc}\begin{array}{l}{\tilde{j}}_i^{(\text{G})}[n]\\ ={\tilde{o}}_i^{(\text{G})}[n]+{\tilde{A}}_i^{(\text{G})}[n]{\tilde{j}}_i^{(\tilde{\text{M}}\text{-}i)}\\ =o_i^{(\text{G})}[n]-A_i^{(\text{G})}[n]{\left(A_i^{(\tilde{\text{M}}\text{-}i)}[n]\right)}^{-1}{o}_i^{(\tilde{\text{M}}\text{-}i)}[n]\\ +A_i^{(\text{G})}[n]{\left(A_i^{(\tilde{\text{M}}\text{-}i)}[n]\right)}^{-1}{\tilde{j}}_i^{(\tilde{\text{M}}\text{-}i)}.\end{array}& (22)\end{array}$

Reminding that we have ${\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]$ from the captured marker positions as in Equations (6) and (7), and ${\text{o}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ and ${\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ expressed by the Fourier series as in Equations (12) and (13), we now have the landmark-based joint position ${\tilde{\text{j}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ for the link-i in the global coordinate system represented as the function of unknown Fourier coefficients of STA ${\text{q}}_{ij,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{ij,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, the distance between the landmarks C_i, as well as ${\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$.

The joint constraint between the link-A and the link-B is satisfied by an appropriate set of the Fourier coefficients of STA, the distance between landmarks, and the constant joint position vectors in the landmark-coordinate system, if the following function $\tilde{F}$ becomes zero:

$\begin{array}{cc}\begin{array}{r}\tilde{F}={\displaystyle \sum _n}|{\tilde{\text{j}}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]-{\tilde{\text{j}}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]|.\end{array}& (23)\end{array}$

In the proposed algorithm, the function $\tilde{F}$ is used as the cost function to be minimized for obtaining the optimal Fourier coefficients of STA, the distance between landmarks, and the joint positions in the landmark-coordinate systems. However, it is not easy to find the optimal parameters that minimize the cost function $\tilde{F}$, because $\tilde{F}$ includes nonlinear terms of the unknown parameters to be solved, such as ${\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ times ${\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, both of which are the functions of the unknown parameters. Interestingly, a certain preprocessing, associated with the fact that the distance between two markers on each link should change periodically in time under the periodicity assumption for the STA, will reduce the difficulty in solving the optimization problem. This is a maneuver that makes the proposed algorithm practically useful. We illustrate the maneuver prior to minimizing $\tilde{F}$.

2.4.3. The Maneuver for Solving the Distance Between Two Landmarks for Each Link and the Posture of the Marker-Coordinate System Relative to the Landmark-Coordinate System

Here we show that C_i (the distance between the landmark-i1 and the landmark-i2 for the link-i) and ${\text{q}}_{i2,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}-{\text{q}}_{i1,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{i2,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}-{\text{q}}_{i1,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ (differences between the Fourier coefficients of STA in the landmark-coordinate system of the link-i) can be obtained without optimizing the cost function $\tilde{F}$. Then, substituting them into Equation (14), we can obtain the posture of the marker-coordinate system ${\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ relative to the landmark-coordinate system.

To this end, we utilize the following identity,

$\begin{array}{cc}\begin{array}{r}|{\text{m}}_{i2}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]-{\text{m}}_{i1}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]|=|{\text{m}}_{i2}^{(\mathrm{\text{G}})}\left[n\right]-{\text{m}}_{i1}^{(\mathrm{\text{G}})}\left[n\right]|,\end{array}& (24)\end{array}$

which expresses the fact that the distance between two motion-captured markers on the link-i in the landmark-coordinate system (the left-hand-side of the identity) and the distance between those in the global coordinate system (the right-hand-side of the identity) are independent of the coordinate system, and they are identical for any instant of time n, although the distance may change periodically in time due to the periodicity of STA.

For the left-hand-side of Equation (24), we rewrite the square of it using Equations (10) and (11), and differences between vectors of Fourier coefficients of STA and C_i, according to the Pythagorean theorem, as follows.

$\begin{array}{cc}\begin{array}{r}|{\text{m}}_{i2}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]-{\text{m}}_{i1}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]{|}^2={\left(C_i+\text{P}\left[n\right]{\xi }_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\right)}^2+{\left(\text{P}\left[n\right]{\xi }_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\right)}^2,\end{array}& (25)\end{array}$

where

$\begin{array}{l}{\text{ξ}}_{i,x}^{(\tilde{\text{M}}\text{-}i)}=q_{i2,x}^{(\tilde{\text{M}}\text{-}i)}-q_{i1,x}^{(\tilde{\text{M}}\text{-}i)},\\ {\text{ξ}}_{i,y}^{(\tilde{\text{M}}\text{-}i)}=q_{i2,y}^{(\tilde{\text{M}}\text{-}i)}-q_{i1,y}^{(\tilde{\text{M}}\text{-}i)}.\end{array}$

That is, ${\text{ξ}}_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{ξ}}_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ represent the differences between the Fourier coefficients of STA in the landmark-coordinate system. Note that, by expanding the square terms in the right-hand-side of Equation (25) concretely, and utilizing the product-to-sum formulae for sine and cosine functions of the Fourier series, the right-hand-side of Equation (25) can be expressed as the sum of terms in the form of v_{i, k}cos(2πkn/N) and w_{i, k}sin(2πkn/N), where k runs from 1 to 2K, with an additional constant term $C_i^{\prime }$. Note also that, in this expanded form of Equation (25), v_{i, k}, w_{i, k}, and $C_i^{\prime }$ are the functions of ${\text{ξ}}_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, ${\text{ξ}}_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, and C_i.

For the right-hand-side of Equation (24), we expand the square of it by another Fourier series as

$\begin{array}{cc}{\left|m_{i2}^{(\text{G})}[n]-m_{i1}^{(\text{G})}[n]\right|}^2={\gamma }_i+{\displaystyle \sum _{k=1}^{2K}\left\{{\alpha }_{i,k}\cos \frac{2\pi kn}{N}+{\beta }_{i,k}\sin \frac{2\pi kn}{N}\right\}},& (26)\end{array}$

which is possible because temporal changes in the distance between two markers in the global coordinate system are also periodic by the periodicity assumption for STA. It is important to note that the left-hand-side of Equation (27) can be obtained directly from the motion-captured data ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$, which can easily be Fourier-expanded to obtain the coefficients α_{i, k}, β_{i, k}, and γ_i.

Then the identity between the right-hand-side of Equation (25) and that of Equation (26) can be guaranteed by term-wise equalization for the coefficients of cos(2πkn/N), (i.e., v_{i, k} = α_{i, k}) and those of sin(2πkn/N), (i.e., w_{i, k} = β_{i, k}) for k = 1, 2, ⋯ , 2K. The unknown parameters ${\xi }_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, ${\xi }_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, and C_i can be solved from the 4K+1 set of equations, since v_{i, k} and w_{i, k} are the quadratic functions of those unknown. Substituting the solved parameter values of the difference in the Fourier coefficients of STA, namely ${\xi }_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}={\text{q}}_{i2,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}-{\text{q}}_{i1,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, ${\xi }_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}={\text{q}}_{i2,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}-{\text{q}}_{i1,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, and C_i into Equation (14), we obtain the posture of the marker-coordinate system ${\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$ relative to the landmark-coordinate system.

2.4.4. Solving the Positions of the Landmarks and the Joint in the Global Coordinate System by Minimizing the Cost Function

So far, we have obtained the differences in the Fourier coefficients of STA (${\xi }_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\xi }_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$) and the posture of the marker-coordinate system (${\text{A}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}\left[n\right]$) relative to the landmark-coordinate system. Those solutions make the optimization of the cost function $\tilde{F}$ defined by Equation (23) easy, as shown in this sequel. Indeed, one can confirm that the cost function can be rewritten as follows.

$\begin{array}{cc}\begin{array}{l}\tilde{F}={\displaystyle \sum _n|\{o_{\text{A}}^{(\text{G})}[n]-\frac{1}{2}{B}_{\text{A}}[n]\left[\begin{array}{c}{{\displaystyle \text{ξ}}}_{A,x}^{(\tilde{\text{M}}\text{-A})}\\ {{\displaystyle \text{ξ}}}_{A,y}^{(\tilde{\text{M}}\text{-A})}\end{array}\right]-B_{\text{A}}[n]\left[\begin{array}{c}{q}_{\text{A}1,x}^{(\tilde{\text{M}}\text{-A})}\\ q_{\text{A}1,y}^{(\tilde{\text{M}}\text{-A})}\end{array}\right]}\\ +{\tilde{A}}_{\text{A}}^{(\text{G})}[n]{\tilde{j}}_{\text{A}}^{(\tilde{\text{M}}\text{-A})}\}-\{o_{\text{B}}^{(\text{G})}[n]-\frac{1}{2}{B}_{\text{B}}[n]\left[\begin{array}{c}{{\displaystyle \text{ξ}}}_{B,x}^{{}^{(\tilde{\text{M}}\text{-B})}}\\ {{\displaystyle \text{ξ}}}_{B,y}^{(\tilde{\text{M}}\text{-B})}\end{array}\right]\\ -B_{\text{B}}[n]\left[\begin{array}{c}{q}_{\text{B}1,x}^{(\tilde{\text{M}}\text{-B})}\\ q_{\text{B}1,y}^{(\tilde{\text{M}}\text{-B})}\end{array}\right]+{\tilde{A}}_{\text{B}}^{(\text{G})}[n]{\tilde{j}}_{\text{B}}^{(\tilde{\text{M}}\text{-B})}\}|,\end{array}& (27)\end{array}$

where

$B_i[n]={\tilde{A}}_i^{(\text{G})}[n]\left[\begin{array}{cc}P[n]& 0\\ 0& P[n]\end{array}\right].$

In Equation (27), the unknown parameters that should minimize the cost are reduced only to the Fourier coefficients of STA for the marker-i1 (${\text{q}}_{i1,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{i1,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$) and the time-invariant constant vector of the joint position in the landmark-coordinate system (${\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$). Since $\tilde{F}$ is now a linear function of those unknown parameters, the optimal values of those parameters can be obtained easily as

$\begin{array}{cc}\begin{array}{r}{\displaystyle \underset{{\text{q}}_{i1,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)},{\text{q}}_{i1,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)},{\tilde{\text{j}}}_i^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}}{argmin}}\tilde{F}.\end{array}& (28)\end{array}$

Once we have ${\text{q}}_{i1,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{i1,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, it follows simply that ${\text{q}}_{i2,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}={\xi }_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}+{\text{q}}_{i1,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\text{q}}_{i2,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}={\xi }_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}+{\text{q}}_{i1,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$. In this way, we could successfully determine the STA in the landmark-coordinate system as the Fourier series of Equation (11), as well as the joint position in each of the landmark-coordinate systems for the link-A and the link-B.

Using the Fourier coefficients of STA in the landmark-coordinate system and the distance between two landmarks for each link, the position and the posture of the landmark-coordinate system in the global coordinate system (${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$) can also be obtained directly from Equations (20) and (21). Moreover, the position of the landmarks and the joint in the global coordinate system (${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{j}}}_i^{(\mathrm{\text{G}})}\left[n\right]$) can be calculated using Equations (10) and (22). As in the case with the naive algorithm, the joint position is estimated as the average of ${\tilde{\text{j}}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}$ and ${\tilde{\text{j}}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}$ as follows.

$\begin{array}{cc}\begin{array}{r}{\tilde{\text{j}}}^{(\mathrm{\text{G}})}\left[n\right]=\frac{{\tilde{\text{o}}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]+{\tilde{\text{A}}}_{\mathrm{\text{A}}}^{(\mathrm{\text{G}})}\left[n\right]{\tilde{\text{j}}}_{\mathrm{\text{A}}}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-A}})}+{\tilde{\text{o}}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]+{\tilde{\text{A}}}_{\mathrm{\text{B}}}^{(\mathrm{\text{G}})}\left[n\right]{\tilde{\text{j}}}_{\mathrm{\text{B}}}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-B}})}}{2}\end{array}& (29)\end{array}$

From these positions, we can obtain the exact position and posture of the link-coordinate system for each link.

3. Methods for Evaluating the Proposed Algorithm

To assess the efficiency of the proposed algorithm, we conducted a numerical experiment using a rigid seven-link model that moves in the sagittal plane (Figure 2), as utilized in Yamasaki et al. (2003). The model consists of a head-arm-trunk (HAT) link, left and right thigh links (l/r-T), left and right shank links (l/r-S), and left and right foot links (l/r-F), which are connected by the pin joints (Figure 2A). See Table 1 for dynamic variables and parameters of the model. Two landmarks (and two landmark-markers), i.e., the landmark-i1 and the landmark-i2 for i={HAT,l/r-T,l/r-S,l/r-F} fixed on each link are defined as in Figure 2A. In this experiment, we set the link-coordinate system for each link such that the origin is located at the CoM of the link and the $x_i^{(\mathrm{\text{L}})}$-axis directs from the CoM to the distal end of the link. For a given set of two landmarks for each link, position and posture of the link-coordinate system can be determined by the centroid of the triangle having the vertices at the joint and two landmarks as in section 2.2, since the CoM position may be available from the statistics of body-parameters, and thus the relative position of the CoM and the centroid can be obtained.

FIGURE 2

Figure 2. Rigid seven-link model of human walking. (A) Positions of landmarks and rigid seven-link model of human body. Rigid seven-link model consists of Head-Arm-Trunk link (HAT), left and right Thigh links (l/r-T), left and right Shank links (l/r-S), and left and right Foot links (l/r-F). Blue circles represent landmarks of each link, and each landmark corresponds to anatomical landmark of human body as follows: landmark-HAT1 is seventh cervical vertebra, landmark-HAT2 is xiphoid process, landmark-l/r-T1 is greater trochanter, landmark-l/r-T2 is lateral condyle of femur, landmark-l/r-S1 is lateral condyle of tibia, landmark-l/r-S2 is lateral malleolus, landmark-l/r-F1 is calcaneus, and landmark-l/r-F2 is first Metatarsal bone. ${\theta }_{\mathrm{\text{HAT}}}^{(\mathrm{\text{G}})}$ is posture of HAT link, which is represented as tilt angle from vertical axis. Joint angles at left/right-Hip (θ_l/r-H), left/right-Knee (θ_l/r-K), and left/right-Ankle (θ_l/r-A), are defined as relative angle for the proximal link. These angles take positive value for clockwise rotation, and negative value for counterclockwise rotation. See Table 1 for details of variables and parameters of the model. (B) Conceptual diagram of methods for evaluating the proposed algorithm. Black lines show walking dynamics of rigid link model, and blue circles are landmarks which move with the rigid link model. Red circles are captured markers which excurse due to STA. See main text for details.

TABLE 1

Table 1. Variables and parameters of the rigid seven-link model.

First, we performed dynamic simulations (forward-dynamics) of the model to obtain stable steady-state periodic walking and the corresponding kinematics of each link, as well as two landmark-markers fixed on each link, i.e., 14 markers in total (Figure 2B). See Yamasaki et al. (2003) for details of the simulation. Specifically, from the simulated kinematics, we obtained time-series of the position and the posture of HAT link (${\tilde{\text{o}}}_{\mathrm{\text{HAT}}}^{(\mathrm{\text{G}})}\left[n\right]$ and ${\theta }_{\mathrm{\text{HAT}}}^{(\mathrm{\text{G}})}\left[n\right]$) and six joint angles of θ_m[n] for m={l/r-H,l/r-K,l/r-A} with the subscript m specifying the joint (not the link) as the left/right hip joint (l/r-H), the left/right knee joint (l/r-K), and the left/right ankle joint (l/r-A). We then calculated the position and the posture of the landmark-coordinate system of each link in the global coordinate system, i.e., ${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ for i={HAT,l/r-T,l/r-S,l/r-F}, referred to as the true kinematics, which we would like to retrieve from motion data of STA-affected markers in this experiment. Note that the link-i for i={l/r-T,l/r-S} is connected to two adjacent links at proximal and distal joints in this model, unlike in the simple two-link model used in section 2, and thus the symbol ${\text{j}}_i^{(\mathrm{\text{M-}}i)}\left[n\right]$ used for the joint position cannot specify the joint uniquely, whether it is proximal or distal joint of the link-i. To resolve this situation, the position of each joint in a given coordinate system is specified with the joint-name subscript m, such as ${\text{j}}_{i,m}^{(\mathrm{\text{M-}}i)}\left[n\right]$ for m={l/r-H,l/r-K,l/r-A}, meaning that, in this case, the position of the m-th joint in the marker-coordinate system of the link-i. For example, ${\text{j}}_{\mathrm{\text{r-T,r-H}}}^{(\mathrm{\text{M-T}})}\left[n\right]$ and ${\text{j}}_{\mathrm{\text{r-T,r-K}}}^{(\mathrm{\text{M-T}})}\left[n\right]$ can be distinguished, respectively, as the position of the right hip joint in the marker-coordinate system of the right thigh and that of the right knee joint also in the marker-coordinate system of the right thigh.

Subsequently, we introduced STA-affected marker positions in the same dynamic simulation as above (Figure 2B). In our previous study, we characterized temporal changes in STA-affected markers that were experimentally motion-captured from seven subjects during periodic walking, in which the markers were attached on the landmarks corresponding to the landmark-ij in Figure 2A (Inoue et al., 2016). Moreover, we obtained a rough estimation of the STA profile for each of seven links in the marker-coordinate system using the naive algorithm. We showed that the STA in the marker-coordinate system for each link might be periodic, and could be fitted by the fourth order Fourier series of the gait cycle. In this study, we considered those STA estimates [i.e., the Fourier coefficients (a_{ij, x, k}, b_{ij, x, k}) and (a_{ij, y, k}, b_{ij, y, k}) as in Equation (11) for k = 1, ⋯ , 4 (i.e., K = 4)] obtained in the previous study as simulated STA profiles in the landmark-coordinate system. Figure 3 shows the Fourier coefficients (a_{ij, x, k}, b_{ij, x, k}) and (a_{ij, y, k}, b_{ij, y, k}) of those STA profiles. As confirmed in Figure 3, the Fourier coefficients of each marker are qualitatively similar across subjects, but also exhibit some individual characteristics quantitatively. It should be noted that those simulated STA profiles used in this experiment were not the “true STA,” because they were estimated by the naive algorithm in the marker-coordinate system, not in the link-coordinate, nor the landmark-coordinate systems. However, we assumed that the STAs estimated by the naive algorithm might be qualitatively similar to the true STA, and they might be good enough to examine performance of the proposed algorithm.

FIGURE 3

Figure 3. Fourier coefficients of STA of markers, used in the numerical experiments for assessment of efficiency of the new posture estimation algorithm (Equation 31). These coefficients were obtained from experimental measurement data in our previous work, using old algorithm (“naive algorithm,” Inoue et al., 2016). We obtained the coefficients from seven subjects, and different color indicates different subject (purple: sub 1, blue: sub 2, aqua: sub 3, green: sub 4, yellow: sub 5, orange: sub 6, red: sub 7). Circles connected with solid lines are Fourier coefficients of sine terms STA of marker-i1 (a_{i1, x, k}, a_{i1, y, k}), triangles with dashed lines are those of cosine terms of STA of marker-i1 (b_{i1, x, k}, b_{i1, y, k}), tetragons with chain lines are those of sine terms of marker-i2 (a_{i2, x, k}, a_{i2, y, k}), and diamonds with dotted lines are those of cosine terms of marker-i2 (b_{i2, x, k}, b_{i2, y, k}). (A,B) Fourier coefficients of x-elements and y-elements of STA of markers, which are attached on HAT. (C,D) Those of markers on left Thigh. (E,F) Those of markers on left Shank. (G,H) Those of markers on left Foot.

The simulated marker positions ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ in the global coordinate system (corresponding to motion-captured positions of the STA-affected markers) for the j-th marker of the link-i were represented by the Fourier-expanded STA and the kinematic data of the model (${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$) as follows.

$\begin{array}{cc}\begin{array}{l}{m}_{ij}^{(\text{G})}[n]\\ ={\tilde{o}}_i^{(\text{G})}[n]\\ +{\tilde{A}}_i^{(\text{G})}[n]\left\{{\tilde{m}}_{ij}^{(\tilde{\text{M}}\text{-}i)}+\left[\begin{array}{c}{\displaystyle {\sum }_{k=1}^4\left\{a_{ij,x,k}\cos \frac{2\pi kn}{N}+b_{ij,x,k}\sin \frac{2\pi kn}{N}\right\}}\\ {\displaystyle {\sum }_{k=1}^4\left\{a_{ij,y,k}\cos \frac{2\pi kn}{N}+b_{ij,y,k}\sin \frac{2\pi kn}{N}\right\}}\end{array}\text{}\right]\right\},\end{array}& (30)\end{array}$

where ${\tilde{\text{m}}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ is the landmark marker position in the landmark-coordinate system of the link-i. Since this is a simulated examination of the algorithm, we know all of the parameters in this equation, but we consider a situation that we have a set of time-series data for the STA-affected ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ only, and the others, i.e., ${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$, ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$, ${\tilde{\text{m}}}_{ij}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, and the STA, are all unknown. The result section shows that we could estimate all of those unknowns from the simulated STA-affected marker positions ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ using the proposed algorithm.

3.1. Gait Data Assimilation by the Naive and Proposed Algorithms

3.1.1. Assimilation by the Naive Algorithm for Comparison

In the naive algorithm, position and posture of the marker-coordinate system in the global coordinate system ${\text{o}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and ${\text{A}}_i^{(\mathrm{\text{G}})}\left[n\right]$ for the link-i were calculated from the marker positions by Equations (6) and (7). Position of the joint-m, ${\text{j}}_{i_1,m}^{(\mathrm{\text{M-}}{i}_1)}\left[n\right]$, in the marker-coordinate system of the link-i₁ and ${\text{j}}_{i_2,m}^{(\mathrm{\text{M-}}{i}_2)}\left[n\right]$ in the marker-coordinate system of the link-i₂ were estimated by Equation (8) for two links i₁ and i₂ that are connected at the joint m. The corresponding joint position j^(G)[n] for the joint m in the global coordinate system was then calculated from ${\text{o}}_{i_1}^{(\mathrm{\text{G}})}\left[n\right]$, ${\text{o}}_{i_2}^{(\mathrm{\text{G}})}\left[n\right]$, ${\text{A}}_{i_1}^{(\mathrm{\text{G}})}\left[n\right]$, ${\text{A}}_{i_2}^{(\mathrm{\text{G}})}\left[n\right]$, ${\text{j}}_{i_1,m}^{(\mathrm{\text{M-}}{i}_1)}$, and ${\text{j}}_{i_2,m}^{(\mathrm{\text{M-}}{i}_2)}$ by Equation (9). Finally, we retrieved the link-coordinate systems ${\bar{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ located at the CoM of the link-i and ${\bar{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ for all i, and then obtained the estimates of the joint angles θ_m[n] of the joint-m for all m.

3.1.2. Assimilation by the Proposed Algorithm

In the proposed algorithm, we first determined the distance between two landmarks C_i and the differences between Fourier coefficients of STA ${\xi }_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\xi }_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, defined by Equations (10), (11), and (25), for each of the link-i from the STA-affected marker positions ${\text{m}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$. In this process, we used function fsolve in MATLAB to find a solution of the quadratic simultaneous equations derived from Equations (25) and (26). An initial condition for the searching difference between the Fourier coefficients in fsolve was set to the Fourier coefficients of STA obtained by the naive algorithm (see Inoue et al., 2016 for details).

Subsequently, we estimated the Fourier coefficients ${\text{q}}_{ij,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, ${\text{q}}_{ij,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, and the joint position ${\tilde{\text{j}}}_{i,m}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ in the landmark-coordinate system using C_i, ${\xi }_{i,x}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ and ${\xi }_{i,y}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$ by minimizing the cost function $\tilde{F}$ of Equation (28) as Equation (29). We then determined the position ${\tilde{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ and the posture ${\tilde{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ of the landmark-coordinate system for each link using Equations (20) and (21), from which we calculated the positions of the landmarks ${\tilde{\text{m}}}_{ij}^{(\mathrm{\text{G}})}\left[n\right]$ and the joints ${\tilde{\text{j}}}_{i,m}^{(\mathrm{\text{G}})}\left[n\right]$ in the global coordinate system using ${\tilde{\text{j}}}_{i,m}^{(\widetilde{\mathrm{\text{M}}}\mathrm{\text{-}}i)}$, Equations (10) and (22), and then ${\tilde{\text{j}}}_m^{(\mathrm{\text{G}})}\left[n\right]$ by Equation (30). Finally, we retrieved the link-coordinate systems ${\bar{\text{o}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ located at the CoM of the link-i and ${\bar{\text{A}}}_i^{(\mathrm{\text{G}})}\left[n\right]$ for all i, and then obtained the assimilated joint angles θ_m[n] of the joint-m for all m.

3.2. Comparison Between the Naive and the Proposed Algorithms

We compared the assimilated posture of HAT link ${\theta }_{\mathrm{\text{HAT}}}^{(\mathrm{\text{G}})}\left[n\right]$ and the joint angles θ_m[n] by the naive and the proposed algorithms with the true kinematics. Moreover, we examined length of each link, i.e., the distance between proximal and distal joints that were estimated by each of two algorithms. Due to the STA, the origin of the link-coordinate system and the joint position in the global coordinate system for each link assimilated by the naive and the proposed algorithms could include errors, which might induce errors in the link-length. To examine those errors, we calculated the lengths of the thigh and the shank links as the distances between the associated joints, and for those of HAT and the foot links, the distances between the origin of the link-coordinate system (the CoM position) and the distal or the proximal joint were calculated as the functions of time for one gait cycle. The obtained mean length of each link was used in the following inverse dynamics analysis.

3.2.1. Comparison in the Inverse Dynamics Analysis

We performed the inverse dynamics analysis to compare errors in the joint torques based on the assimilated HAT and the six joint angles by the naive and the proposed algorithms with the true joint torques. The true joint torques ${\tilde{\tau }}_m$ for the joint-m were obtained from the true kinematics, the corresponding ground reaction forces that were obtained in the forward dynamics simulation, and the model parameters shown in Table 1 (see Appendix A for details of the inverse dynamics analysis).

Estimations of the joint torques, denoted by τ_m, were also calculated for two assimilated motions by the naive and the proposed algorithms, in which we used the estimated joint angles and the estimated length of each link that had been obtained in the assimilation process as described above. We used the true values for the mass and the inertia moment for each link, and for the ground reaction forces. We then compared the joint torque τ_m with the true joint torque ${\tilde{\tau }}_m$ for all m to assess the efficiency of the proposed algorithm.

4. Results of the Evaluation

In the numerical experiment, we successfully simulated a stable periodic sequence of gait and the corresponding true kinematics, from which we could obtain the STA-affected marker positions for each of seven motion-captured subjects shown in Figure 3. As shown in Figure 4, the assimilated kinematics of the model by the proposed algorithm, and thus the resultant joint torques τ_m by the inverse dynamics were the exactly the same as those of the true kinematics and joint torques ${\tilde{\tau }}_m$ for m = {l/r-H,l/r-K,l/r-A}, while those by the naive algorithm exhibit substantial errors from the true kinematics and joint torques.

FIGURE 4

Figure 4. Kinematics and joint torques estimated from spatio-temporal data of markers including STA, which correspond to Fourier coefficients shown in Figure 3. In each panel, there are seven colored curves and one black curve. Colored curves in (A–D) are estimated kinematics by using naive algorithm from simulated marker positions with STA, of which Fourier coefficients are plotted with same color in Figure 3. Those curves in (E–G) are joint torques calculated from the estimated kinematics (τ_l-H, τ_l-K, τ_l-A), and those in (H–J) are difference between the calculated joint torques and joint torques of original kinematics (${\tilde{\tau }}_{\mathrm{\text{l-H}}}$, ${\tilde{\tau }}_{\mathrm{\text{l-K}}}$, ${\tilde{\tau }}_{\mathrm{\text{l-A}}}$). Black curve in each panel shows kinematics and joint torques, which were estimated by using new algorithm proposed in this study. Since the proposed algorithm estimated model kinematics that was exactly same with original kinematics for all (=seven) case of STA, only one curve are described in each panel, and black straight lines at zero are drawn in (H–J). Gray shaded region in each panel represents duration when left Foot contacts ground (stance phase of left leg). (A) Postures of HAT link (${\theta }_{\mathrm{\text{HAT}}}^{(\mathrm{\text{G}})}$). (B) Joint angles of left-Hip (θ_l-H). (C) Those of left-Knee (θ_l-K). (D) Those of left Ankle (θ_l-A). (E) Joint torques at left-Hip (τ_l-H). (F) Those at left-Knee (τ_l-K). (G) Those at left-Ankle (τ_l-A). (H) Difference between joint torques at left-Hip (${\tau }_{\mathrm{\text{l-H}}}-{\tilde{\tau }}_{\mathrm{\text{l-H}}}$). (I) Those at left-Knee (${\tau }_{\mathrm{\text{l-K}}}-{\tilde{\tau }}_{\mathrm{\text{l-K}}}$). (J) Those at left-Ankle (${\tau }_{\mathrm{\text{l-A}}}-{\tilde{\tau }}_{\mathrm{\text{l-A}}}$).

More specifically, Figures 4A–D show, respectively, the assimilated joint angles of the left leg for ${\theta }_{\mathrm{\text{HAT}}}^{\mathrm{\text{(G)}}}$, θ_l-H, θ_l-K, and θ_l-A by the naive and the proposed algorithms. In each panel, seven colored curves represent ${\theta }_{\mathrm{\text{HAT}}}^{\mathrm{\text{(G)}}}$ or θ_m that were estimated by the naive algorithm from each of the STA profiles for seven subjects shown in Figure 3, for which the colors in each panel of Figures 4A–D correspond to those in each panel of Figure 3. The black curves in Figures 4A–D, but visible only in Figures 4A,D, are those estimated by the proposed algorithm, and they were exactly the same with the true kinematics. In Figures 4B,C, the back curves for the proposed algorithm and the colored joint angles θ_l-H and θ_l-K for the naive algorithm are closely overlapped, meaning that the hip and knee joint angles are less-influenced by STA, regardless of the assimilation algorithms.

The left ankle joint angles (Figure 4D) estimated by the naive algorithm were largely different from the corresponding true joint kinematics, in comparison with the angles of the other joints (Figures 4B,C). It is noteworthy that the ankle joint angles estimated by the naive algorithm tended to be smaller than the true joint angles. Such large deviation from the true joint ankle angle was especially large for the subject 7 (red curve), although the STA Fourier coefficients for the subject 7 were not particularly deviated from those for the other subjects, except the cosine terms of y-element of marker-2 on the left-shank link (Figure 3F, red diamonds). Moreover, interestingly, despite the fact that the STA Fourier coefficients of the subject 4 for the cosine terms of y-element of marker-1 on the left-thigh link (Figure 3D green triangle) was particularly large, the assimilated kinematics for the subject 4 (Figure 4, green curves) were not particularly deviated from those for the other subjects. These observations imply non-trivial (nonlinear) relationships between the amounts of STA and the errors in the assimilated joint kinematics.

Figures 4E–G are joint torques of the left-hip (τ_l-H), the left-knee (τ_l-K), and the left-ankle (τ_l-A), and Figures 4H–J show differences between the estimated joint torques (τ_m) and the true joint torques ${\tilde{\tau }}_{\mathrm{\text{l-H}}}$, ${\tilde{\tau }}_{\mathrm{\text{l-K}}}$, and ${\tilde{\tau }}_{\mathrm{\text{l-A}}}$ as the functions of time. Since the assimilated kinematics by the proposed algorithm was exactly the same with the true kinematics, the estimated joint torques were also exactly the same with the true joint torques, leading to the zero differences between them (the black horizontal lines in Figures 4H–J).

The joint torques estimated by the naive algorithm exhibited unignorable errors during stance phase, which can be confirmed in Figures 4H–J, where the gray shaded region in each panel indicates the stance phase of the left foot of the left leg. In all of seven STA profiles, the ankle joint angles and torques estimated by the naive algorithm showed similar error profiles and took smaller values than the true joint angles and torques. Such similarity between Figures 4D,G implies that the errors in the joint torques have its origin in the misestimation of the joint angles of the left ankle (Figure 4G), which propagates to the joint torques at the knee (Figure 4I) and the hip (Figure 4H), even though the estimated posture of HAT link and the hip and knee joint angles did not include large errors (Figures 4A–C). In case of the subject 7 (red curve), joint torques at the left-ankle (Figure 4G) took negative values during most of the left-leg stance phase, while the true joint torque of the left-ankle (the black curve in Figure 4G) took positive values for more than half of the stance phase. This error in the naive algorithm for the inverse dynamics analysis was absent in the proposed algorithm that removes the effect of STA.

5. Summary and Discussion

5.1. Summary

In this paper, we proposed a simple but efficient algorithm to assimilate data of motion-captured marker positions affected by soft tissue artifact (STA) into models of multi-rigid-body systems. In the proposed algorithm, an unknown STA profile for each marker was Fourier-expanded, and its Fourier coefficients were determined so that the captured motion was optimally represented by a model of the rigid multi-link system. The determination of the Fourier coefficients of the STA profiles, i.e., the estimation of the non-measurable STA, allows to estimate the temporal changes in the global positions of the markers that are firmly fixed to the skeletal links (the landmark markers in this paper), leading to the STA-free assimilated motion that is consistent with the multi-rigid-link model. The key idea to determine the unknown STA was simply the periodicity assumption for the STA and kinematic constraints requiring that any two adjacent rigid-links are connected by a rotary joint.

To assess the efficiency of the proposed algorithm, we conducted a numerical experiment using a rigid seven-link model of human gait in the sagittal plane. In the experiment, we estimated kinematics of the model from seven sets of STA-affected marker positions, and showed that the proposed algorithm could determine the true kinematics of the model accurately for all sets of the STA profiles.

5.2. Difference Between the Proposed Algorithm and Similar Algorithms

Since the proposed algorithm estimate kinematics of multi-rigid-body model using the joint constraints for the whole body, the proposed algorithm can be considered as a kind of the “global optimization”(Leardini et al., 2005) including multi-body kinematics optimization (Duprey et al., 2010; Richard et al., 2017) and kinematics estimation based on Kalman filter (Cerveri et al., 2003, 2005; Bonnet et al., 2017a). Here, let us clarify the novel aspect of the proposed algorithm complementary to the existing global optimization methods. In most studies of the global optimization, multi-rigid link models with virtual markers that are firmly fixed to the skeletal links are used. Kinematics of the multi-rigid-body system is optimized such that total distance between positions of captured markers and the corresponding virtual markers is minimized, under so-called “hard constraint” such as the joint constraint (Richard et al., 2017) or a physical constraint based on kinetic data (Cerveri et al., 2003, 2005). Since the cost function is the total distance between the captured and the virtual markers (Cerveri et al., 2003, 2005; Duprey et al., 2010), kinematics of the model is estimated such that motions of the virtual markers are similar to that of the corresponding captured markers, the STA of which are not treated in this process. Consequently, the estimated kinematics obtained by the global optimization can be strongly affected by the STA. That is, in this process, the larger STA amplitude, the stronger the effect on the kinematics estimation is. In contrast to these algorithms of the global optimization, in the proposed algorithm, an unknown STA profile is determined, and STA-free marker positions and the corresponding kinematics of a multi-rigid link model are estimated. This process is independent of the amplitude of STA, and can be used as long as the STA is periodic.

The major difference between the proposed algorithm and the other algorithms of the global optimization is whether STA profiles are explicitly determined or not. Since STA profiles are explicitly determined in the proposed algorithm, it may be claimed that we should use database of STA (Cereatti et al., 2017). Nevertheless, we believe that the proposed algorithm has some advantages over the algorithms that use the databases. Since STA is task-dependent and not reproducible for different subjects, effect of STA remains in kinematics estimated by the databases, unless a captured trial exactly matches the data in the database. In contrast, when we use the proposed algorithm, STA profiles are determined for each trial from the captured data, therefore, we do not need worry about task-dependent or subject-dependent STA, as long as the STA is periodic.

The proposed algorithm and the existing global optimization methods should be used complementary and collaboratively with other, rather than compared in their performance. Although we showed that the proposed algorithm for removing the STA can improve the naive algorithm as a very primitive global optimization method in this paper, a data assimilation can be performed better in general, if the proposed algorithm is utilized jointly with more sophisticated global optimization methods such as used in AnyBody (Andersen et al., 2009, 2010; Lund et al., 2015), among others.

5.3. Spurious Solutions

Although the proposed algorithm succeeded to find the true kinematics of the human gait model for all examined STA profiles in the numerical experiment, it could fail depending on some other STA profiles. Since the minimization of the cost function $\tilde{F}$ of Equation (27) is a linear optimization problem, possible failures might be caused by the quadratic simultaneous equations arising from Equations (24–26) for solving the distance between two landmarks and the difference between the Fourier coefficients of STA. This is because the quadratic simultaneous equations have multiple solutions including the true and spurious solutions, and the spurious ones may be obtained if an initial condition for the solution search is not appropriate.

A promising way to set a good initial condition for solving the quadratic simultaneous equation is to use the naive algorithm based on the marker-coordinate system defined by Equations (6) and (7), which provide a rough estimation of the STA profiles without assuming the periodicity of the STA. We then expand the rough estimation of the STA into Fourier series, and utilize the obtained coefficients as the initial values of the STA Fourier coefficients for the proposed algorithm. Moreover, it is expected that the proposed algorithm would become more robust by finding appropriate positions of the landmarks for reducing a risk of being trapped by spurious solutions.

5.4. Extension to Three-Dimensional Motions

An obvious issue that should be addressed in the future is to extend the proposed algorithm to three-dimensional motions and models. It would be able to achieve, because every processing performed in the proposed algorithm is independent of the dimensionality, but application of the algorithm to three-dimensional motions and models would just increase the number of markers (from two markers to three for each link) and unknown Fourier coefficients. That is, to specify a posture of a single rigid body in the three-dimensional space, at least three markers (marker-ij for j = {1,2,3}) should be attached on landmarks of the link-i (the landmark-ij) for a motion capture experiment.

However, there might be difficulty again in solving the quadratic simultaneous equations arising from Equations (24–27), particularly in the three-dimensional space. This is because the number of unknown parameters to be solved becomes larger than the number of the quadratic simultaneous equations if we consider distances between two markers on each link as in the two-dimensional case, leading to an ill-posed set of the quadratic simultaneous equations. Nevertheless, we still might be able to derive additional quadratic equations about distances between other pairs of markers among three markers on each link. Thus, it is expected that we can assimilate motion-captured kinematics of each link in the three-dimensional space into three-dimensional multi-link models using the proposed algorithm extended to the three-dimensional space.

Author Contributions

YS and TN designed the study, and wrote manuscript. YS and TI developed the proposed algorithm and the various computer programs.

Funding

This work was supported in part by JSPS grants-in-aid 17K13016.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

1. ^The origin and axes of the link-coordinate system can be defined in different ways. For example, the origin can also be at the center of mass (CoM) of the link-i, and the $x_i^{(\mathrm{\text{L}})}$-axis can be directed from the origin to distal end of the link as in Inoue et al. (2016).

References

Alexander, E., and Andriacchi, T. (2001). Correcting for deformation in skin-based marker systems. J. Biomech. 34, 355–361. doi: 10.1016/S0021-9290(00)00192-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Alexander, N., and Schwameder, H. (2016). Lower limb joint forces during walking on the level and slopes at different inclinations. Gait Posture 45, 137–142. doi: 10.1016/j.gaitpost.2016.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersen, M., Damsgaard, M., MacWilliams, B., and Rasmussen, J. (2010). A computationally efficient optimisation-based method for parameter identification of kinematically determinate and over-determinate biomechanical systems. Comput. Methods Biomech. Biomed. Eng. 13, 171–183. doi: 10.1080/10255840903067080

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersen, M., Damsgaard, M., and Rasmussen, J. (2009). Kinematic analysis of over-determinate biomechanical systems. Comput. Methods Biomech. Biomed. Eng. 12, 371–384. doi: 10.1080/10255840802459412

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersen, M., Damsgaard, M., Rasmussen, J., Ramsey, D., and Benoit, D. (2012). A linear soft tissue artefact model for human movement analysis: proof of concept using in vivo data. Gait Posture 35, 606–611. doi: 10.1016/j.gaitpost.2011.11.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Andriacchi, T., Alexander, E., Toney, M., Dyrby, C., and Sum, J. (1998). A point cluster method for in vivo motion analysis: applied to a study of knee kinematics. J. Biomech. Eng. 120, 743–749. doi: 10.1115/1.2834888

PubMed Abstract | CrossRef Full Text | Google Scholar

Benoit, D., Ramsey, D., Lamontagne, M., Xu, L., Wretenberg, P., and Renström, P. (2006). Effect of skin movement artifact on knee kinematics during gait and cutting motions measured in vivo. Gait Posture 24, 152–164. doi: 10.1016/j.gaitpost.2005.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Blache, Y., Dumas, R., Lundberg, A., and Begon, M. (2017). Main component of soft tissue artifact of the upper-limbs with respect to different functional, daily life and sports movements. J. Biomech. 62, 39–46. doi: 10.1016/j.jbiomech.2016.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonnet, V., Dumas, R., Cappozzo, A., Joukov, V., Daune, G., Kulić, D., et al. (2017a). A constrained extended kalman filter for the optimal estimate of kinematics and kinetics of a sagittal symmetric exercise. J. Biomech. 62, 140–147. doi: 10.1016/j.jbiomech.2016.12.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonnet, V., Richard, V., Camomilla, V., Venture, G., Cappozzo, A., and Dumas, R. (2017b). Joint kinematics estimation using a multi-body kinematics optimisation and an extended kalman filter, and embedding a soft tissue artefact model. J. Biomech. 62, 148–155. doi: 10.1016/j.jbiomech.2017.04.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Camomilla, V., Dumas, R., and Cappozzo, A. (2017). Human movement analysis: the soft tissue artefact issue. J. Biomech. 62, 1–4. doi: 10.1016/j.jbiomech.2017.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappello, A., Cappozzo, A., La Palombara, P., Lucchetti, L., and Leardini, A. (1997). Multiple anatomical landmark calibration for optimal bone pose estimation. Hum. Mov. Sci. 16, 259–274. doi: 10.1016/S0167-9457(96)00055-3

CrossRef Full Text | Google Scholar

Cappello, A., Stagni, R., Fantozzi, S., and Leardini, A. (2005). Soft tissue artifact compensation in knee kinematics by double anatomical landmark calibration: performance of a novel method during selected motor tasks. IEEE Trans. Biomed. Eng. 52, 992–998. doi: 10.1109/TBME.2005.846728

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappozzo, A., Catani, F., Leardini, A., Benedetti, M., and Della Croce, U. (1996). Position and orientation in space of bones during movement: experimental artefacts. Clin. Biomech. 11, 90–100. doi: 10.1016/0268-0033(95)00046-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappozzo, A., Croce, U. D., Leardini, A., and Chiari, L. (2005). Human movement analysis using stereophotogrammetry: Part 1: theoretical background. Gait Posture 21, 186–196. doi: 10.1016/j.gaitpost.2004.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cereatti, A., Bonci, T., Akbarshahi, M., Aminian, K., Barré, A., Begon, M., et al. (2017). Standardization proposal of soft tissue artefact description for data sharing in human motion measurements. J. Biomech. 62, 5–13. doi: 10.1016/j.jbiomech.2017.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerveri, P., Pedotti, A., and Ferrigno, G. (2005). Kinematical models to reduce the effect of skin artifacts on marker-based human motion estimation. J. Biomech. 38, 2228–2236. doi: 10.1016/j.jbiomech.2004.09.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerveri, P., Rabuffetti, M., Pedotti, A., and Ferrigno, G. (2003). Real-time human motion estimation using biomechanical models and non-linear state-space filters. Med. Biol. Eng. Comput. 41, 109–123. doi: 10.1007/BF02344878

PubMed Abstract | CrossRef Full Text | Google Scholar

Clément, J., Dumas, R., Hagemeister, N., and de Guise, J. (2015). Soft tissue artifact compensation in knee kinematics by multi-body optimization: performance of subject-specific knee joint models. J. Biomech. 48, 3796–3802. doi: 10.1016/j.jbiomech.2015.09.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Dicharry, J. (2010). Kinematics and kinetics of gait: from lab to clinic. Clin. Sports Med. 29, 347–364. doi: 10.1016/j.csm.2010.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumas, R., and Cheze, L. (2009). Soft tissue artifact compensation by linear 3d interpolation and approximation methods. J. Biomech. 42, 2214–2217. doi: 10.1016/j.jbiomech.2009.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Duprey, S., Cheze, L., and Dumas, R. (2010). Influence of joint constraints on lower limb kinematics estimation from skin markers using global optimization. J. Biomech. 43, 2858–2862. doi: 10.1016/j.jbiomech.2010.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, C., Suzuki, Y., Kiyono, K., Morasso, P., and Nomura, T. (2014). An intermittent control model of flexible human gait using a stable manifold of saddle-type unstable limit cycle dynamics. J. R. Soc. Interface 11:101. doi: 10.1098/rsif.2014.0958

PubMed Abstract | CrossRef Full Text | Google Scholar

Halvorsen, K., Johnston, C., Back, W., Stokes, V., and Lanshammar, H. (2008). Tracking the motion of hidden segments using kinematic constraints and kalman filtering. J. Biomech. Eng. 130:011012. doi: 10.1115/1.2838035

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, G. F., and Smith, P. A. editors (1996). Human Motion Analysis: Current Applications and Future Directions. New York, NY:IEEE Press.

Google Scholar

Holden, J., Orsini, J., Siegel, K., Kepple, T., Gerber, L., and Stanhope, S. (1997). Surface movement errors in shank kinematics and knee kinetics during gait. Gait Posture 5, 217–227. doi: 10.1016/S0966-6362(96)01088-0

CrossRef Full Text | Google Scholar

Inoue, T., Suzuki, Y., Kiyono, K., and Nomura, T. (2016). “Skin motion artifact in motion capturing of human bipedal gait: characterization and influence on joint torque estimation,” in Proceedings of the 8th International Workshop on Biosignal Interpretation, BSI2016 (Osaka), 213–216.

Google Scholar

Lafortune, M., Cavanagh, P., Sommer, III H., and Kalenak, A. (1992). Three-dimensional kinematics of the human knee during walking. J. Biomech. 25, 347–357. doi: 10.1016/0021-9290(92)90254-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Leardini, A., Chiari, A., Della Croce, U., and Cappozzo, A. (2005). Human movement analysis using stereophotogrammetry part 3. Soft tissue artifact assessment and compensation. Gait Posture 21, 212–225. doi: 10.1016/j.gaitpost.2004.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J.-D., Lu, T.-W., Lin, C.-C., Kuo, M.-Y., Hsu, H.-C., and Shen, W.-C. (2017). Soft tissue artefacts of skin markers on the lower limb during cycling: effects of joint angles and pedal resistance. J. Biomech. 62, 27–38. doi: 10.1016/j.jbiomech.2017.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, T.-W., and Chang, C.-F. (2012). Biomechanics of human movement and its clinical applications. Kaohsiung J. Med. Sci. 28, S13–S25. doi: 10.1016/j.kjms.2011.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, T.-W., and O'Connor, J. (1999). Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints. J. Biomech. 32, 129–134. doi: 10.1016/S0021-9290(98)00158-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Lucchetti, L., Cappozzo, A., Cappello, A., and Della Croce, U. (1998). Skin movement artefact assessment and compensation in the estimation of knee-joint kinematics. J. Biomech. 31, 977–984. doi: 10.1016/S0021-9290(98)00083-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lund, M. E., Andersen, M. S., de Zee, M., and Rasmussen, J. (2015). Scaling of musculoskeletal models from static and dynamic trials. Int. Biomech. 2, 1–11. doi: 10.1080/23335432.2014.993706

CrossRef Full Text | Google Scholar

Peters, A., Galna, B., Sangeux, M., Morris, M., and Baker, R. (2010). Quantification of soft tissue artifact in lower limb human motion analysis: a systematic review. Gait Posture 31, 1–8. doi: 10.1016/j.gaitpost.2009.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, L., Jones, R., and Howard, D. (2008). Whole body inverse dynamics over a complete gait cycle based only on measured kinematics. J. Biomech. 41, 2750–2759. doi: 10.1016/j.jbiomech.2008.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Richard, V., Cappozzo, A., and Dumas, R. (2017). Comparative assessment of knee joint models used in multi-body kinematics optimisation for soft tissue artefact compensation. J. Biomech. 62, 95–101. doi: 10.1016/j.jbiomech.2017.01.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Riemer, R., Hsiao-Wecksler, E., and Zhang, X. (2008). Uncertainties in inverse dynamics solutions: a comprehensive analysis and an application to gait. Gait Posture 27, 578–588. doi: 10.1016/j.gaitpost.2007.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Rusaw, D., and Ramstrand, N. (2011). Motion-analysis studies of transtibial prosthesis users: a systematic review. Prosthet. Orthot. Int. 35, 8–19. doi: 10.1177/0309364610393060

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryu, J., Miyata, N., Kouchi, M., Mochimaru, M., and Lee, K. (2006). Analysis of skin movement with respect to flexional bone motion using mr images of a hand. J. Biomech. 39, 844–852. doi: 10.1016/j.jbiomech.2005.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Sati, M., De Guise, J., Larouche, S., and Drouin, G. (1996). Quantitative assessment of skin-bone movement at the knee. Knee 3, 121–138. doi: 10.1016/0968-0160(96)00210-4

CrossRef Full Text | Google Scholar

Sibella, F., Galli, M., Romei, M., Montesano, A., and Crivellini, M. (2003). Biomechanical analysis of sit-to-stand movement in normal and obese subjects. Clin. Biomech. 18, 745–750. doi: 10.1016/S0268-0033(03)00144-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, S. R. (2004). Quantification of human motion: gait analysis—benefits and limitations to its application to clinical problems. J. Biomech. 37, 1869–1880. doi: 10.1016/j.jbiomech.2004.02.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Winter, D. A. (2009). Biomechanics and Motor Control of Human Movement, 4th Edition. Hoboken, NJ: Wiley.

Google Scholar

Yamasaki, T., Nomura, T., and Sato, S. (2003). Possible functional roles of phase resetting during walking. Biol. Cybern. 88, 468–496. doi: 10.1007/s00422-003-0402-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshikawa, N., Suzuki, Y., Kiyono, K., and Nomura, T. (2013). A theoretical study on a computational algorithm for human posture estimation based on motion capture of a small number of markers. Adv. Biomed. Eng. 2, 107–116. doi: 10.14326/abe.2.107

CrossRef Full Text | Google Scholar

Yoshikawa, N., Suzuki, Y., Ozaki, W., Yamamoto, T., and Nomura, T. (2012). “4d human body posture estimation based on a motion capture system and a multi-rigid link model,” in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS, 4847–4850. doi: 10.1109/EMBC.2012.6347079

CrossRef Full Text | Google Scholar

Appendix A

Inverse Dynamics Analysis of Rigid Seven-Link Model

Motion equation of rigid seven-link model in the sagittal plane is described as follows

$\begin{array}{cc}\begin{array}{r}\text{M}(\text{θ})\ddot{\text{θ}}+\text{B}(\text{θ},\dot{\text{θ}})+\text{K}(\text{θ})=\text{F}+\tau ,\end{array}& (\mathrm{A1})\end{array}$

where θ is the vector of position and postures of all segments, M is the inertia matrix, B is the centrifugal and Coriolis force, K is the gravitational force and torque, F is the ground reaction force, and τ is the joint torque vector. See Yamasaki et al. (2003) or Fu et al. (2014) for details of elements of each matrix and vectors.

In the numerical experiment, we performed forward dynamics simulation, and obtained marker position data and ground reaction force. Subsequently, model kinematics (θ) were estimated by using old algorithm (“naive algorithm”) or new algorithm proposed in this study. From the estimated kinematics, joint angular velocities ($\dot{\theta }$) and joint angular accelerations ($\ddot{\theta }$) were calculated, and then, inertia matrix (M), vector of centrifugal and Coriolis force (B), and vector of gravitational force and torque (K) at each time instance were calculated. Finally, joint torques (τ) were calculated according to Equation (A1).