Force Shadows: An Online Method to Estimate and Distribute Vertical Ground Reaction Forces from Kinematic Data

Kinetic models of human motion rely on boundary conditions which are defined by the interaction of the body with its environment. In the simplest case, this interaction is limited to the foot contact with the ground and is given by the so called ground reaction force (GRF). A major challenge in the reconstruction of GRF from kinematic data is the double support phase, referring to the state with multiple ground contacts. In this case, the GRF prediction is not well defined. In this work we present an approach to reconstruct and distribute vertical GRF (vGRF) to each foot separately, using only kinematic data. We propose the biomechanically inspired force shadow method (FSM) to obtain a unique solution for any contact phase, including double support, of an arbitrary motion. We create a kinematic based function, model an anatomical foot shape and mimic the effect of hip muscle activations. We compare our estimations with the measurements of a Zebris pressure plate and obtain correlations of 0.39≤r≤0.94 for double support motions and 0.83≤r≤0.87 for a walking motion. The presented data is based on inertial human motion capture, showing the applicability for scenarios outside the laboratory. The proposed approach has low computational complexity and allows for online vGRF estimation.


Introduction
The analysis of forces driving the human motion and the forces exerted to the surrounding via contacts with the human body are part of kinetic human motion analysis, which is an important tool for instrumented human motion analysis [1,2]. According to Newton's third law, the ground exerts a force on the human body by mirroring the acceleration of the body's center of mass (CoM). This force is called ground reaction force (GRF) and carries important information to analyze gait parameters, e.g., in rehabilitation, sports, robotics, etc. [3]. The peak vertical GRF is a reliable parameter to evaluate lower extremity functional strength during sport-specific movements and can be used, together with the vertical impulse, which is the product of force multiplied by time during the propulsion phase of a vertical jump, for prognoses after sports-related injuries [4]. The landing phase in running sports produce high-frequency peaks in the vertical GRF trajectories and are assumed to play a crucial role in the development of injuries [5]. The vertical GRF can be used for online estimation of (unknown) carried loads, e.g., in working environments, in order to assess physical stress [6][7][8]. In laboratories, e.g., when human gait is measured and analyzed using force plates, subjects often need to strike a locally fixed measuring device during the measurement which can lead to disturbances in their natural walking behavior [9,10]. To avoid the restriction induced by a force plate, inverse dynamics can be used to derive kinetic parameters from kinematic data in order to estimate the current stress on joints during human motion and on an external contact point. Since kinematic data can be measured from mobile inertial motion capture, this allows for kinetic analysis outside the laboratory. The starting point for inverse dynamics calculations is the knowledge about the external forces. In this work, we focus on GRF, especially the vertical component. In the case of only one connected contact region, e.g., during single support phase in walking, the GRF can be uniquely derived from motion data, cf. [11][12][13]. However, in double support phases, more than one contact region is present and the human body forms a closed kinetic chain which makes the estimation of GRFs indeterminate (indeterminancy problem) [14][15][16]. In that case, we need to provide additional information to obtain a unique solution for the GRF. This was done via a smooth transition assumption and an empirically tuned function for gait in [11,13]. Static optimization and a detailed musculoskeletal model with additional muscle actuators under the feet were used in [14], which yielded good results for different types of motion; however, it has a very high computational cost, preventing online usage to date. A neural network was used in [10], which has low computational cost but is inherently data dependent. A further approach is to employ additional measurements for the vGRF via pressure insoles, cf. [17].
In contrast to the mentioned work, we propose a method that is independent of the motion and generalizes well to different motions, thus is not restricted to walking. Additionally, it has a low computational complexity which makes it real-time capable. The overall idea is to compute a conservative load distribution on each contact region (here, the feet) via a suitable projection of each segment along the negative GRF vector to the next surface (here, the ground) and utilize anatomical shape information of the foot as well as rough modeling of hip flexion as additional information. This allows for a conservative and consistent distribution of the vGRF in the double stance phase. The method is evaluated with three subjects that performed weight shifting movements with stationary placements of the feet (double support) as well as normal gait cycles on a single pressure plate.
The paper is structured as follows: Section 2 introduces the force shadow method (FSM), Section 3 describes the experimental results and Section 4 presents our conclusions as well as future work.

Methods
In this section, we formalize the underlying idea behind the FSM. According to Newton's third law of motion, each body that exerts a force on a second one receives a force equal in magnitude and opposite direction. With this law in mind, and assuming that a standing person is in contact with the ground only via the feet, we observe an exerted force on the ground caused by the subject's body weight, which in turn causes the ground to react with a force in opposite direction. This reaction force which acts on the human body (likewise in the dynamic case) is called the ground reaction force One key question motivating this work is the following: How can the GRF be consistently distributed over multiple contact regions, if only the kinematics and mass distributions of the body are present? The approach should be simple and consistent, i.e., match direct measurements at the contact regions.
One simple idea for the GRF estimation is to use only the center of mass (CoM) and compute the GRF with respect to one contact point [18]. This simplifies the kinematics of the complete human body to the kinematics of the CoM and reduces the interaction surface to one contact point, but does not allow to deal with multiple contact points.
The basic idea behind the FSM is to approximate a body mass distribution and project it onto a possible contact or interaction surface along the overall force vector, similar to a "shadow". In a further step, this projected distribution is considered as weighting for the contact regions.
In Section 2.1, we formally explain how to determine F GRF and how we use it to construct a force shadow (surface). As a direct application of the force shadow, Section 2.1.4 offers a general approach to obtain a plantar distribution of the total mass. In Section 2.2, we adjust the approach and make it more suitable to the human foot structure by incorporating the curvature of the foot arch and hip flexor muscles and propose the adjusted formula in Section 2.2.3. Then, we complete the method chapter presenting the algorithm in Section 2.4. The notation is explained in the abbreviations section, at the end of this article.

Force Shadow Method
Let us first start with formalizing the GRF with respect to the CoM of a multi-body system. In the following, we assume the positions x of the CoM and all segments to be given with respect to one global coordinate system G, i.e., x = x G , but we leave out the superscript for better readability.

Ground Reaction Force
Newton's third law tells us that the ground reaction force mirrors every acceleration of the body's CoM, based on the assumption that the feet define the only contact regions between the body and its surrounding. Let S be the set of all body segments, i.e., we assume that the body is a multi-body system consisting of |S| rigid body segments with masses m s , then the body's CoM, x CoM ∈ R 3 , can be written as the weighted sum of the segments' CoMs which can be derived as follows (cf. [19]): where m total denotes the subject's total mass. Provided that we know the accelerations of each segment, e.g., by kinematics, Equation (2) allows to compute the acceleration of the body's CoM: The ground reaction force can be derived as follows:

Construction of Force Shadows for a Human Body
The force shadow method (FSM) constructs a surface by summing up multiple bivariate elliptical normal distributions using projected CoMs as mean value and the segments' orientations to construct the corresponding covariances.
Let x CoM s ,t ∈ R 3 be the CoM for an arbitrary segment s ∈ S at time t which is assumed to be in the middle of the segment. The time stamp t ∈ R >0 denotes the current time at which the measurement was taken. Then, x CoM s ,t is projected onto the bottom (we assume z = 0) along the ground reaction force vector in order to acquire the respective projected Gaussian center, denoted as µ s,t ∈ R 3 , i.e., where n = [0, 0, 1] T is the unit vector which is orthogonal to the ground. The factor k ∈ R is chosen, such that the vertical component of µ s,t is zero, i.e., it corresponds to a point on the interaction surface. This allows to consider variations in the 2-dimensional surface since only the projection onto the bottom plane is of interest, hence, we assume x, µ s,t ∈ R 2 hereafter. The covariance matrix Σ s,t ∈ R 2×2 is constructed in such a way that its principal components have the same orientation as the corresponding segment.
Let a s , b s ∈ R be the parameters representing the dispersion of the ellipsoid of the Gaussian function and R s,t ∈ R 2×2 the rotation matrix responsible to rotate this ellipsoid around the vertical axis into the direction of the corresponding segment's orientation, then the covariance is determined via See Section 2.3 for the estimation of a s , b s . Equation (6) ensures the symmetry and the positive definiteness of the covariance matrix. Thus, for segment s ∈ S we obtain the function The constant d s ∈ [0, 1] represents the contribution of segment s to the total body mass in terms of body mass proportion as proposed in [20], they are given in percentage. We repeat the projection for all segments. The shadow function is then obtained via superposition of all projected contributions of the segments For an illustration, see Figure 1.

Contact Region Estimation
In order to obtain contact regions for a general motion, we need a model for the respective parts of the body that incorporates possible contact with the environment. In this work, we consider only the foot regions, depicted in Figure 1, and a flat ground as interaction surface. The foot region model is registered with a 2D foot model that is tailored to the person at hand.
The pose of the region model is computed as follows: Note, in the following, we consider the left foot only; the right foot can be treated analogously. We denote x i ∈ R 2 for i ∈ LM := {FM, V M, CA} as the anatomic landmarks taken from the current pose and p i ∈ R 2 the respective landmarks of the foot model; then we minimize the energy function where t ∈ R 2 , M ∈ R 2×2 , in order to obtain a transformation of the foot model, see Figure 2. We consider the transformation as a matrix multiplication using homogeneous coordinates, i.e., where 0 = [0, 0] T is a column vector consisting of zeros. In order to apply the transformation to the points, these are transformed into homogeneous coordinates by adding 1 to the third component, Let X := x FMxV MxCA and P := p FMpV MpCA , X, P ∈ R 3×3 , be the homogeneous marker points in matrix form. Then, we can rewrite Equation (12) as The linear least-squares problem can be solved in closed-form via the pseudo-inverse, e.g., utilizing the singular value decomposition (SVD). In each time step, the regions Ω L , Ω R are transformed this way. The (sub-)regions in Ω L , Ω R are defined to be in contact with the ground if the registered center falls below a distance h min ≈ 0.03 m from the ground and has a velocity below v min ≈ 0.8 m s . These values are chosen similar to [14]. If, for a (sub-)region, at least one condition (height below h min or velocity below v min ) is not fulfilled, this particular region is assumed to have no ground contact and, therefore, set to zero. With these conditions, the method is independent of stance times or gait event patterns.

Force Distribution
Assuming unconnected sets of regions Ω L ⊆ R 2 for the left and Ω R ⊆ R 2 for the right foot obtained from the registration described in the previous section: Consider an arbitrary region Ω k ⊆ Ω L of the left foot model (Ω k ⊆ Ω R analogously). Let F GRF,z be the vertical component of the ground reaction force from Equation (3), then, the estimated load on this region is generally given by where the load distribution factor I Ω k ,t ∈ [0, 1] for the corresponding region is computed by integration and normalization, i.e, Here, f t corresponds to the force shadow defined in (11). The integral of an arbitrary region Ω k ⊆ (Ω L ∪ Ω R ) is approximated via a triangulation of this region (see Figure 1 on the right) using the midpoint rule, that is, wheref i,t is the mean value of the points evaluated at the nodes and |∆ i | is the area of the triangle ∆ i . The number of triangles in region Ω k is denoted as N ∆ k .

Anatomical Models
In this section, we describe additional, biomechanically motivated modifications of the load distribution.

Arch of Foot
The above described force distribution method puts equal weight on all points of the foot region model. This leads to a similar load, e.g., on the heel and the arch, during a neutral stance which would be only the case for a pathologic foot anatomy. This would have an impact on the center of pressure (CoP) for each foot, i.e., on the vGRFs' application points. To compensate for this effect, we take additional weights into account. To model the curvature of the foot, we fit a surface w(x) = a 0 + a 1 x + a 2 y + a 3 xy + a 4 x 2 + a 5 y 2 , for x := [x, y] ∈ R 2 , such that w(x) = 1 for the contact points of the metatarsal heads and the calcaneus bone and w(x) = h a dist CP for the midpoints between the front and rear contact points, where dist CP := x CA − x V M 2 is the distance between the the subject's fifth metatarsal head and calcaneus (not to be confused with the fixed scalar d CP from Figure 3). The scalar factor h a ∈ R is a parameter modeling the anatomical height of the cuboid bone relative to the distance between the support points of calcaneus and the head of the fifth metatarsal bone, cf. Figure 3 on the left. To guarantee the conservation of integrals, we rescale the weight function pointwise with a time-dependent scalarw Note that the scalars n t belonging to the left and right foot differ from each other, thus, we have to do this twice to obtainw L,t (x) andw R,t (x).

Hip Flexion Muscles
So far, the proposed approach cannot detect decreasing loads resulting from lifting off the feet slightly, for example, in double support movements in which the whole CoM is located within the region of one foot. This can hardly be measured by human motion capturing systems, especially when the feet do not lift off, but influences the vGRF. Therefore, we model the effect of the muscles responsible for hip flexion, see Figure 3 on the right, as follows: We assume that the muscles of the unloaded leg are increasingly activated, the longer the mass is balanced on the other foot, which is interpretable as raising the foot gradually, the longer the person stands on the other foot.
In order to model the described effect, the total loads are distributed accordingly to the left and right foot. Therefore, for each time step t ∈ R >0 we compute the current velocity,ẋ CoM,t , of the CoM to predict its position for the nextt ∈ N time instances, i.e., x CoM,t+k := x CoM,t + k∆t ·ẋ CoM,t , for k ∈ {0, . . . ,t} ⊆ N.
Let s L,t , s R,t ∈ [0, 1], initialized as s L,t 0 , s R,t 0 = 1, be the distribution parameters for the left and right foot. Then, if one of the predicted CoMs is contained in the convex hull of the left foot, that is, x CoM,t+k ∈ Conv(Ω L ), for at least one k, we increase s L,t (if s L,t < 1) and decrease s R,t (if s R,t > 0) using a ramp function with slope m (see Section 2.3 for the estimation of m). There are three cases in total which can occur: Note that the increment and decrement should provide that s L,t , s R,t ∈ [0, 1], e.g., we do not increase s L,t 0 , s R,t 0 as they already are at their maximum. The weightsw t (x) from (20) are then multiplied by a linear combination of s L,t and s R,t , which leads to the modified weight functionŝ ).

Anatomy Aware Load Estimation
Incorporating the anatomy of a normal foot and mimicking the effect of hip flexion muscles on the FSM in (16) and (17) for an arbitrary region Ω k we obtain the final vGRF distribution formula as follows

Utilized Hyper Parameters
For optimizing the hyper parameters introduced in the previous sections, i.e., the dispersion parameters a s , b s constructing the Gaussian bell functions in Equation (7), the time parametert, and the slope m of the ramp in Section 2.2.2, we used the load trajectories of a representative dataset (we choose sway all around (AR) of one representative person, see Section 3.2 for the description of the performed motions) for each of the six depicted areas of the foot model, see Figure 1 on the right; that is, we define the vector L Ω L ,i t , L Ω R ,i t ∈ R 6 as the estimated load vector for the time frame i t and P Ω L ,i t , P Ω R ,i t ∈ R 6 regarding the corresponding loads measured by the pressure plate. Let N f denote the number of frames of the trajectories and P the set of all parameters, i.e., then, we minimize the following energy function with respect to P using the interior-point method with a step tolerance of 10 −8 , that is, the optimization ends if P k+1 − P k < 10 −8 (1 + P k ). Note the dependency of the FSM on the parameter set, i.e., L Ω L ,i t = L Ω L ,i t (P ) and L Ω R ,i t = L Ω R ,i t (P ).

FSM Algorithm
The subsequent algorithm summarizes the proposed method and requirements. We denote the locations of the foot landmarks as x FM , x V M , x CA and p FM , p V M , p CA regarding the kinematic pose and the foot model, respectively (cf. Figure 1).

Materials
For the validation of the proposed method, we used two measurement systems. The Xsens full body inertial measurement system with 17 IMUs at a sampling rate of 60 Hz to capture human motion; MVN studio 4.97.1 for data acquisition. The validated Zebris pressure plate type FDM 1.5 [21] acted as reference measuring system for pressure values. The pressure plate has an update rate of 100 Hz, a sensor area with dimensions 149 cm × 54.2 cm (length × width) and a spatial pressure value resolution of 0.8469 cm 2 , i.e., per time frame, the pressure plate data is stored in a 176 × 64 pressure matrix. The methods were implemented and evaluated in MATLAB, version 9.4.0.813654 (R2018a).

Experiments
Three subjects (body weight 80.3 ± 7.7 kg, i.e. between 73 and 87 kg, age between 27 and 29 years) walked onto the pressure plate and performed double support movements (weight shifting), i.e., sway all around (AR), sway from side to side (S2S) and sway back and forth (BF), with eighteen cycles (three cycles per measurement) in total for each movement. Further, more dynamic motions like squatting (SQ) and walking (NW) were performed. As in the double support case, eighteen squats (three squats per measurement) in total were measured. Regarding NW, the subjects walked over the plate in normal speed measuring nine non-consecutive gait cycles (heel-strike to toe-off) per measurement; in total, 54 gait cycles were measured. Due to poor inertial tracking data, 12 cycles were omitted.
After receiving all relevant study information, the participants signed an informed consent to the study including a permission to publish the data.
The focus is on the evaluation of the predicted vGRFs and the forces for the three regions (the front part is summarized into one) per foot shown in Figure 1 using the proposed method in Algorithm 1 as compared to the pressure plate data.

Registration of the Pressure Plate Data
In order to obtain reference data from the pressure plate, the locations of the foot prints on the pressure plate need to be known, such that we can detect which region Ω k is activated and which force is acting on this respective region. For that reason, we extract the foot prints on the pressure plate by summing up the 176 × 64 pressure matrices for every time frame (given a static pose of the feet, i.e., we know that in the measured range the subject's foot pose was fixed) and use an indicator function for each activated cell, i.e., we set the cells with non-vanishing pressure values to one; that way, we obtain a binary matrix storing the spatial information about the 2D foot prints. The connected non-zero entries of the matrix represent the Cartesian coordinates of the foot prints. These coordinates are registered to the foot models via a 2D point registration, using the rigid extended coherent point drift (rECPD), cf. [22]. To determine the vGRF for the region Ω k the pressure values that are active at the respective region and time are multiplied pointwise by the cell's area (0.8469 cm 2 ) and summed up in that region. For the total vGRF the total foot prints need to be considered.

Error Metrics
To evaluate the performance of the FSM, the estimations y(k), for k = 1, . . . , N, were compared with the measurementsỹ(k) of the pressure plate, both normalized by bodyweight, e.g., y Heel (k) = L Ω Heel ,k m total regarding the estimation for the heel at time step k. In order to enable the comparison to recent studies, we include the metrics defined in [13], i.e., the root-mean-square error (RMSE), where e(k) = y(k) −ỹ(k) is defined as the error at the k th sample, and the relative root-mean-square error (rRMSE), which constitutes the normalization of RMSE by the peak-to-peak amplitude between both discrete trajectories. Further metrics included in the evaluation are the mean absolute error (MAE) of the sample data the standard deviation of the absolute error (SD) and the Pearson correlation coefficient (r)

Subregions
The results of the subregions in Ω L and Ω R , according to the foot model shown in Figure 1 on the right, are listed in Table 1. The subregions comprise heel, arch of the foot and front foot, which consists of Ω Metatarsum and Ω Toes . Table 1. Error metrics for the subregions, i.e., heel (vGRF H ), foot arch (vGRF A ), front foot (vGRF F ) and the total foot (vGRF), per movement for each foot. Walking with normal speed (NW), squatting (SQ), swaying all around (AR), from side to side (S2S) and back and forth (BF). The chosen metrics, i.e., root-mean-square error (RMSE), relative root-mean-square error (rRMSE), mean-absolute error (MAE), standard deviation (SD) and the correlation coefficient (r), are determined including all measured cycles per motion. The smallest deviations can be found in the mid foot with a RMSE of at most 0.76 N kg for weight shifting movements and 1.03 N kg for walking. The RMSEs of estimated loads on heel and front foot do not exceed 2.36 N kg and 2.13 N kg , respectively, which can be found in S2S and NW movements. Regarding the entire foot area, the RMSE lies in between 0.60 N kg and 1.64 N kg for all movements. Small deviations in the CoM can already falsely detect the activation of hip flexion muscles leading to lower correlations, see the correlation coefficients in SQ. Hence, overall, the results show that the FSM can approximate measured pressure values and vGRF reasonably well, based on the kinematics of the human body, its CoM and the force shadow function.     Figure 7 shows the vGRF trajectories obtained from the FSM as compared to the reference data during walking with normal speed (normalized with the corresponding total weight). Shown are the gait cycles from heel-strike to toe-off. According to [13,23,24], the trajectories follow the typical structural pattern of gait cycles.

Execution Time
In Table 2, the execution times for the substeps of the proposed FSM algorithm are listed; some substeps are taken together due to interdependencies. As the computation time is independent of the movement, one AR measurement with 1308 time frames was chosen for the statistics. Note, the execution time was produced with (non-optimized) MATLAB code, which can be much slower as compared to an optimized C++ implementation. A single iteration of the algorithm takes 114.1 ms (average over 1308 frames) making the method usable for online applications. The algorithm was executed using an Intel Core i7-7700 processor at 3.6 GHz. Table 2. Execution times regarding one full iteration (last line) and the substeps during a sway all around measurement. The third column shows the mean values and one standard deviation of the required time for 1308 executions.

Code Line
Procedure Execution Time

Discussion on Application Scenarios and Accuracy
The FSM predicts the vGRF for arbitrary contact regions online and is designed to be robust to different stance times or gait events, e.g., initial contacts; via biomechanically motivated models and few hyper-parameters. However, the method has not been validated for patients with severe gait abnormalities, so far. Possible differences to the presented results could arise, e.g., for different single stance times, due to abnormalities of gait or balance that could have an influence on the hip flexion model. The parameterst and m for the hip flexion model were treated as hyper-parameters in this work and were only tested for the presented motion types for healthy participants.
Regarding accuracy comparisons with existing methods for GRF prediction, the presented approach has comparable to inferior accuracy regarding vGRF prediction as compared to the approach presented in [14], with respect to the squatting motion and an increased error (by about a factor of 2) with respect to normal walking (cf. Table 3 of [14] with Table 1). Note that the approach [14] is an offline approach that considers an optimization over a detailed musculoskeletal model and is, thus, considerably more expensive regarding computation time and only applicable with the detailed musculoskeletal model. Moreover, our approach was not yet optimized for walking scenarios, e.g., the hyper-parameters were only tuned on one person for the AR motion and used for all other motions, including NW.
One application scenario, where the FSM can be applied, is related to the Ovako Working Posture Analysing System (OWAS) [6][7][8] to assess work postures. This system is partitioned in four parts, the first three parts are concerned with the worker's pose, in particular, the pose of the head, back, arms and legs. The fourth OWAS category classifies the carried load in the following classes, from 0 kg to 10 kg, 10 kg to 20 kg and above 20 kg with the respective scores. With this specific application in mind and assuming a target accuracy of 5 kg to 10 kg, we tolerate errors of about 0.61 N kg to 1.22 N kg (using the subjects' average body weight of 80.3 kg for weight normalization). Comparing with the results in Table 1, this target accuracy can be reached. However, to estimate unknown carried loads, additional data is needed, which can be provided by pressure insoles. In this scenario, the proposed method can act as a predictor for vGRF (in comparison to the the insoles' measurements), especially when ambiguities due to increased accelerations (caused by motion) occur, as this would lead to ambiguities with increased loads, if measured with insoles alone.

Main Limitations
We used full body kinematics, based on inertial motion tracking, and a mass distribution table [20] to obtain the CoM kinematics. Inertial motion capture can have sensor-to-segment calibration errors with severe impacts on the tracked kinematics [25]. Figure 8 demonstrates how errors in tracking the CoM can result in vGRF deviations. In this situation, the FSM is not able to distribute the vGRF accordingly, especially in movements where the CoM is predominantly located in the center of the body, see the correlation coefficients for BF or SQ in Table 1. Further, the pressure estimates of the different regions can be ambiguous, in the case where the pressure distribution under the foot is changed due to muscle contraction, but the kinematics do not change. Thus, the distributions under each foot have to be interpreted with care. However, the total vGRF, predicted for each independent contact region, e.g., left or right foot, is invariant to this due to integration over the region. The only variation here can happen in the application point of the vGRF, i.e., the CoP under the foot, which is derived from the distribution of the respective contact region.
Finally, the results have been, up to this point, only tested on healthy male subjects of rather similar weight.

Conclusions
In this work, we presented the FSM which is an approach to estimate and distribute vGRF for multiple contact regions, using only kinematic data of a human body, approximate segment weight distributions and the person's overall weight. This has been addressed in literature so far only for specific motion patterns [11,13], completely data-based, i.e., with neural networks [10], or with a high computational burden, utilizing a detailed musculoskeletal model [14,26]. In contrast, the proposed approach is generally applicable, model-based and has a small computational complexity which allows for online (real-time) usage. The experimental results also show that the proposed FSM can approximate measured pressure values and vGRF reasonably well with kinematic data from inertial body motion capture, hence, allowing for a mobile setup. Therefore, this work is meant as another step in the direction of consistent (online) mobile kinematic and kinetic analysis of human motion which is currently lacking in literature.

Future Work
The FSM is not limited to the shown movements and should be evaluated for other, more dynamic, movements, e.g., running and cutting maneuvers, which are relevant in the field of sports medicine.
In the current work, only healthy foot anatomy was considered. However, the anatomical model can be extended to pathologic feet.
The projected distribution parameters (a s , b s ) introduced in (7) are currently optimized to fit a small set of subjects. According to the constructive nature of the approach, these could be replaced by projected segment volumes of an adapted human body model, e.g., automatically generated based on a portable RGB-D sensor as proposed in [27][28][29].
The influence of sensor-to-segment calibration errors already described in Section 3.10 could be addressed by including automatic sensor alignment in the inertial tracking in future studies, similar to [30,31], or by applying other corrections of the body's CoM. One possibility is the exploitation of sensor insoles to provide additional data. In this scenario, the FSM could be used as a predictor (in the sense of statistical inference) to improve the insoles' measurements via bias estimation or estimate additionally carried loads, see also Section 3.9. The combination of full-body kinematics and carried load estimation is particularly interesting in the context of ergonomic assessment, see e.g., [6][7][8]32,33].
In order to get access to complete gait analyzes, it is inevitable to estimate the full 3D GRF and the ground reaction torques for each foot. Based on the latter, one can also obtain the center of pressure, which can be used as additional information to correct the CoM.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

FSM
proposed force shadow method R set of real numbers N set of natural numbers F GRF three dimensional ground reaction force F GRF,z third, vertical component of F GRF S set of all segments of a human body |S| number of segments m mass in kilogram m s mass of segment s ∈ S in kilogram m total total mass in kilogram x position vector, x = [x, y, z] or x = [x, y], respectivelÿ x acceleration vector CoM center of mass x CoM location of total mass center x CoM,s location of segment's center of mass, s ∈ S g acceleration of gravity R rotation matrix N (µ, Σ, x) bell-shaped probability density function representing a normal distribution with mean value µ, covariance Σ, evaluated at x Ω general region specifier (Ω L and Ω R indicate the regions for the left and right foot) p position vector regarding the landmarks of the foot model |·| absolute value function L Ω estimated load on region Ω ∆t difference between two measured time frames Conv(Ω) convex hull of region Ω w(x),w(x),ŵ(x) weighting functions evaluated at x P set of hyper parameters · p p-norm or l p -norm RMSE root mean square error rRMSE relative root mean square error MAE mean absolute error SD standard deviation r Pearson correlation coefficient; r ∈ [0, 1] AR motion indicator: shifting the weight all around S2S motion indicator: shifting the weight from side to side BF motion indicator: shifting the weight back and forth SQ motion indicator: performing squats NW motion indicator: walking with normal speed