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

Human movement analysis is often performed with a model of multi-rigid-body system, whereby reflective-marker-based motion capture data are assimilated into the model for characterizing kinematics and kinetics of the movements quantitatively. Accuracy of such analysis is limited, due to motions of the markers on the skin relative to the underlying skeletal system, referred to as the soft tissue artifact (STA). Here we propose a simple algorithm for assimilating motion capture data during periodic human movements, such as bipedal walking, into models of multi-rigid-body systems in a way that the assimilated motions are not affected by STA. The proposed algorithm assumes that STA time-profiles during periodic movements are also periodic. We then express unknown STA profiles using Fourier series, and show that the Fourier coefficients can be determined optimally based solely on the periodicity assumption for the STA and kinematic constraints requiring that any two adjacent rigid-links are connected by a rotary joint, leading to the STA-free assimilated motion that is consistent with the multi-rigid-link model. To assess the efficiency of the algorithm, we performed a numerical experiment using a dynamic model of human gait composed of seven rigid links, on which we placed STA-affected markers, and showed that the algorithm can estimate the STA accurately and retrieve the non-STA-affected true motion of the model. We also confirmed that our STA-removal processing improves accuracy of the inverse dynamics analysis, suggesting the usability of the proposed algorithm for gait analysis.


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-linkrigid-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 mechanodynamic 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 . 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 multirigid-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(Cappello et al., , 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., 2012Yoshikawa et al., , 2013Clément et al., 2015;Richard et al., 2017), approaches with Kalman filters (Cerveri et al., 2003(Cerveri et al., , 2005Halvorsen 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 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 . 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.

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. We denote by m (G) ij [n] 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 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. i={A,B} and j={1,2}. The integers n = 1, 2, · · · indicate the data numbers or the sampling times. m (G) ij [n] corresponds to an experimentally motion-captured positions of the markerij contaminated by unknown STA. Now, the positions of the markers located exactly on the landmarks are denoted asm (G) ij [n] (Figure 1B), which cannot be measured experimentally, but can be estimated if we can estimate the unknown STA time-profiles asm (G) ij The problem here can be stated as follows: the motion capture experiments provide time-series data of m (G) ij [n], 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 thenm (G) ij [n], as described in this sequel.

Definitions of Marker-Positions and the Multi-Link System
We define a local coordinate system, referred to as the linkcoordinate system, which is fixed on the link-i as in Figure 1C to describe the position and posture of the link-i. The linkcoordinate system of the link-i is denoted by x (L) i -y (L) i , 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ō (G) i [n] represents the global position of the link-i, and the coordinate transformation matrixĀ (G) i [n], consisting of normalized orthogonal basis vectors in the directions of x (L) iand y (L) i -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 (m (G) i1 [n] andm (G) i2 [n]). As in Figure 1C, the y (L) i -axis of the link-coordinate system of the linki is defined so that it is parallel to the vector directing from , and then the x (L) i -axis is defined so that it becomes orthogonal to the y (L) i -axis 1 . It is essential for the proposed algorithm to notice thatō (G) i [n] andĀ (G) i [n] of the linkcoordinate system defined here cannot be determined easily from the motion-captured marker data m (G) ij [n], since three vertices at the joint,m (G) i1 [n] andm (G) i2 [n] are all unknown. That is, the maker positions m (G) ij [n] may show temporal deviations from the landmark positionsm (G) ij [n] due to unknown STA, and thus the landmark positionsm (G) i1 [n] andm (G) i2 [n] 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 positionsm (G) ij [n] are exactly the same as the motion-captured marker positions m (G) ij [n], and the position of the joint could be calculated as described below.

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 positionsm (G) ij [n] (Figure 1C), referred to as the landmarkcoordinate system and denoted by x (M) i -y (M) i . The origin of the landmark-coordinate system in the global coordinate system, o (G) i [n] for the link-i, is set as the midpoint of the landmark-i1 and the landmark-i2, and the x (M) i -axis of the landmarkcoordinate system relative to the global coordinate system, which determines the posture matrixÃ (G) i [n], is set as the normalized vector directing from the landmark-i1 to the landmark-i2: The joint position in the landmark-coordinate system of the linki 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 vectorj (M-i) i , where the superscript (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 vectorj (G) i [n], can be described by the position and the posture of the landmark-coordinate system of the link-i as follows: 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 thẽ j (G) [n], 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., This is the joint constraint that assures the two adjacent links are connected at a single joint. The two unknown constant , 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 Then the joint position in the global coordinate system,j (G) [n], can be obtained by Equation (3). In this way, using the joint positionj (G) [n] in the global coordinate system and the positions of two landmarks for each linkm (G) i1 [n] andm (G) i2 [n] that are experimentally available in the special case withm (G) ij [n] = m (G) ij [n], we can determine the positionõ (G) i [n] and the posturẽ A (G) i [n] 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.

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(Yoshikawa et al., , 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 where 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 landmarkcoordinate system, and the relative position between the markercoordinate system of link-i (o (G) i [n] and A (G) i [n]) and the actual joint position (j (G) ) changes in time. Therefore, unlike in the landmark-coordinate system, the joint position vector in the marker-coordinate system j (M-i) i inevitably changes in time. Nevertheless, in the naive algorithm, the joint position in the Frontiers in Bioengineering and Biotechnology | www.frontiersin.org marker-coordinate system of the link-i is estimated as the timeindependent constant vectors j (M-i) i (i=A and B) that minimizes the cost function as arg min ) 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 (9) The position and the posture of each link are then estimated by using j (G) . 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.

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 [n] that is the position of non-STA-affected landmark marker in the landmark coordinate system. m (M-i)

ij
[n] can be described formally as with where C i is the distance between two landmark markers on the link-i, and e (M-i)

ij
[n] 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 i -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 i -axis of the coordinate system. We expand the STA for the marker-ij in the landmarkcoordinate system into the following Fourier series: with 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 [n] be the position and the posture of the marker-coordinate system for the link-i, respectively ( Figure 1D), relative to the landmark-coordinate system (õ ). 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 can be expressed by using the Fourier-expanded STA as follows. where Frontiers in Bioengineering and Biotechnology | www.frontiersin.org In this sequel, we show that C i and the Fourier coefficients [n] (i.e., STA in the landmark-coordinate system), q (M-i) ij,x and q (M-i) ij,y , 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, [n] in the landmark-coordinate system, which can be transformed into the landmark-marker positionsm (G) ij [n] in the global coordinate system. Once we determinem (G) ij [n] successfully, the position and the posture of each link (õ (G) in the landmarkcoordinate system, and the joint positionj (G) [n] in the global coordinate system can be obtained using Equations (1, 2, 5), and Equation (3), respectively.

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 in the landmark-coordinate system, which is also represented in the global coordinate as by the coordinate transformation directly from the markercoordinate system (with its origin at o (G) i [n] and the posture 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 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 By rearranging Equations (18) and (19), the landmark-coordinate system in the global coordinate system can be related to the marker-coordinate system as Equations (20) and (21) imply that we can obtain the landmarkcoordinate system in the global coordinate system, if, somehow, we can obtain the position and the posture of the markercoordinate system relative to the landmark-coordinate system, for the marker-coordinate system are available from the motion-captured marker positions using Equations (6) and (7).

The Cost Function for Ensuring the Joint Constraint
To estimate the Fourier coefficient of STA, q (M-i) ij,x and q (M-i) ij,y , a cost function to ensure the joint constraint, such as in Equation (8) (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 Reminding that we have o (G) i [n] and A (G) i [n] from the captured marker positions as in Equations (6) and (7) (12) and (13), we now have the landmark-based joint positionj (G) i [n] for the link-i in the global coordinate system represented as the function of unknown Fourier coefficients of STA q (M-i) ij,x and q (M-i) ij,y , the distance between the landmarks C i , as well asj (M-i) 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 functionF becomes zero: In the proposed algorithm, the functionF 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 functionF, becauseF includes nonlinear terms of the unknown parameters to be solved, such as timesj (M-i) 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 minimizingF. Here we show that C i (the distance between the landmark-i1 and the landmark-i2 for the link-i) and q (M-i) i2,x − q (M-i) i1,x and q (M-i) i2,y −q (M-i) i1,y (differences between the Fourier coefficients of STA in the landmark-coordinate system of the link-i) can be obtained without optimizing the cost functionF. Then, substituting them into Equation (14), we can obtain the posture of the marker- [n] relative to the landmark-coordinate system.
To this end, we utilize the following identity, which expresses the fact that the distance between two motioncaptured 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-handside 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.
That is, ξ (M-i) i,x and ξ (M-i) i,y 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 .
Note also that, in this expanded form of Equation (25), v i,k , w i,k , and C ′ i are the functions of ξ (M-i) i,x , ξ (M-i) i,y , and C i . For the right-hand-side of Equation (24), we expand the square of it by another Fourier series as 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 (26) can be obtained directly from the motion-captured data m (G) ij [n], 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 ξ (M-i) i,x , ξ (M-i) i,y , 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 ξ (M-i) ,y , and C i into Equation (14), we obtain the posture of the marker-coordinate system A (M-i) i [n] relative to the landmark-coordinate system. [n]) relative to the landmarkcoordinate system. Those solutions make the optimization of the cost functionF defined by Equation (23) easy, as shown in this sequel. Indeed, one can confirm that the cost function can be rewritten as follows. where
In Equation (27), the unknown parameters that should minimize the cost are reduced only to the Fourier coefficients of STA for the marker-i1 (q (M-i) i1,x and q (M-i) i1,y ) and the time-invariant constant vector of the joint position in the landmark-coordinate system (j (M-i) i ). SinceF is now a linear function of those unknown parameters, the optimal values of those parameters can be obtained easily as arg min Once we have q (M-i) i1,x and q (M-i) i1,y , it follows simply that 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 landmarkcoordinate system and the distance between two landmarks for each link, the position and the posture of the landmarkcoordinate system in the global coordinate system (õ (G) i [n] and A (G) i [n]) can also be obtained directly from Equations (20)

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 (L) i -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.
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 (õ (G) HAT [n] and θ (G) HAT [n]) 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.,õ (G) i [n] andÃ (G) i [n] 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 j (M-i) i [n] 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 j (M-i) i,m [n] 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 markercoordinate system of the link-i. For example, j (M-T) r-T,r-H [n] and j (M-T) r-T,r-K [n] 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 STAaffected markers that were experimentally motion-captured from seven subjects during periodic walking, in which the markers were attached on the landmarks corresponding to the landmarkij 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 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. θ (G) HAT 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 wherem (M-i) ij is the landmark marker position in the landmarkcoordinate 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 m (G) ij [n] only, and the others, i.e.,õ (G) , 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 m (G) ij [n] using the proposed algorithm.

Assimilation by the Naive Algorithm for Comparison
In the naive algorithm, position and posture of the markercoordinate system in the global coordinate system o (G) i [n] and A (G) i [n] for the link-i were calculated from the marker positions by Equations (6) and (7). Position of the joint-m, j (M-i 1 ) i 1 ,m [n], in the marker-coordinate system of the link-i 1 and j (M-i 2 ) i 2 ,m [n] in the marker-coordinate system of the link-i 2 were estimated by Equation (8) for two links i 1 and i 2 that are connected at the joint ). 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 ).

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 ξ (M-i) i,x and ξ (M-i) i,y , defined by Equations (10), (11), and (25), for each of the link-i from the STA-affected marker positions m (G) ij [n]. 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 q (M-i) ij,x , q (M-i) ij,y , and the joint positionj (M-i) i,m in the landmark-coordinate system using C i , ξ (M-i)

Comparison Between the Naive and the Proposed Algorithms
We compared the assimilated posture of HAT link θ (G) HAT [n] 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.

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τ m for the jointm 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τ m for all m to assess the efficiency of the proposed algorithm.

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τ 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.
More specifically, Figures 4A-D show, respectively, the assimilated joint angles of the left leg for θ (G) HAT , θ l-H , θ l-K , and θ l-A by the naive and the proposed algorithms. In each panel, seven colored curves represent θ (G) HAT 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, 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 nontrivial (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 leftknee (τ 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τ l-H ,τ l-K , andτ 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).
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 (τ l-H ,τ l-K ,τ 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 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 leftankle ( 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.

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.

Difference Between the Proposed Algorithm and Similar Algorithms
Since the proposed algorithm estimate kinematics of multi-rigidbody model using the joint constraints for the whole body, the proposed algorithm can be considered as a kind of the "global optimization"  including multi-body kinematics optimization (Duprey et al., 2010;Richard et al., 2017) and kinematics estimation based on Kalman filter (Cerveri et al., 2003(Cerveri et al., , 2005Bonnet 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(Cerveri et al., , 2005. Since the cost function is the total distance between the captured and the virtual markers (Cerveri et al., 2003(Cerveri et al., , 2005Duprey 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 multirigid 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(Andersen et al., , 2010Lund et al., 2015), among others.

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 functionF 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.

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 threedimensional 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-26), 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.