Rapid bilevel optimization to concurrently solve musculoskeletal scaling, marker registration, and inverse kinematic problems for human motion reconstruction

Creating large-scale public datasets of human motion biomechanics could unlock data-driven breakthroughs in our understanding of human motion, neuromuscular diseases, and assistive devices. However, the manual effort currently required to process motion capture data is costly and limits the collection and sharing of large-scale biomechanical datasets. We present a method to automate and standardize motion capture data processing: bilevel optimization that is able to scale the body segments of a musculoskeletal model, register the locations of optical markers placed on an experimental subject to the markers on a musculoskeletal model, and compute body segment kinematics given trajectories of experimental markers during a motion. The optimization requires less than five minutes of computation to process a subject’s motion capture data, compared with about one day of manual work for a human expert. On a sample of 34 trials of experimental data, the root-mean-square marker reconstruction error (RMSE) was 1.38 cm, approximately 40% lower than the 2.58 cm achieved manually by 3 experts. Optimization solutions reconstructed known joint angle trajectories from four diverse motion trials of synthetic data to an average of 0.79 degrees RMSE. We have published an open source cloud service at AddBiomechanics.org to process experimental motion capture data, which is available at no cost and asks that users agree to share processed and de-identified data with the community. Reducing the barriers to processing and sharing high-quality human motion biomechanics data will enable more people to engage in state-of-the-art biomechanical analysis in their work, do so at lower cost, and share larger and more accurate datasets. Author summary Creating large-scale public datasets of human motion could unlock data-driven breakthroughs in our understanding of neuromuscular diseases, assistive devices, and human motion more broadly. The manual effort currently required to process these motion datasets is costly and limits the collection and sharing of large-scale datasets. Our cloud-based software tool, called AddBiomechanics, uses state-of-the-art optimization techniques to automatically scale the body segments of a musculoskeletal model to match the subject of interest, and then compute body segment kinematics during a motion. The optimization requires less than five minutes of computation to process a subject’s motion capture data, compared with about one day of manual work for a human expert. The accuracy of the approach in quantifying the body segment kinematics is as good or better than the results achieved manually by experts. Reducing the barriers to processing and sharing high-quality human motion biomechanics data will enable more people to engage in state-of-the-art biomechanical analysis, do so at lower cost, and share larger and more accurate datasets.

Quantitative analysis of human movement dynamics is a powerful tool that has been 2 widely used to solve problems such as assessing muscle and joint function e.g. [1][2][3][4][5], 3 better understanding joint loading and consequent adverse health outcomes e.g. [6][7][8][9][10][11][12], 4 analyzing the performance of assistive devices for improving human movement 5 e.g. [13][14][15][16], quantifying the effects of complex conditions like Parkinson's disease 6 e.g. [17,18], and even generating more realistic computer graphics e.g. [19,20]. The 7 benchmark data acquisition technique for quantifying human biomechanics remains 8 optical motion capture [21,22]. This process involves placing a number of optical 9 markers on a subject's body segments and having the subject perform actions in a 10 laboratory space surrounded by specialized cameras. These camera systems are able to 11 reconstruct the three-dimensional locations of the optical markers in the lab, and given 12 the marker trajectories over time, one can use proprietary, open, or custom software to 13 reconstruct the kinematics of the subject's body segments. 14 The current state-of-the-art software for reconstructing the motion of a human 15 musculoskeletal model from optical marker trajectories requires substantial iterative 16 "guess-and-check" refinement, and can take a large amount of expert time to achieve a 17 high quality result on each new subject. This manual interaction increases costs, limits 18 scalability, and reduces the reproducibility of motion capture studies [23,24]. 19 To reconstruct movement kinematics from optical motion capture data, software 20 must address several sources of noise, ambiguity, and model error. Initially the optical 21 markers all appear as undifferentiated bright spots to the cameras, so the first challenge 22 is identifying which bright spots correspond to which locations of markers placed on the 23 body. This is called marker labeling, and is generally handled by having a user manually 24 label each marker using commercial software from motion capture vendors [25,26]. Once 25 the markers have been labeled, software must reconstruct a digital twin of the captured 26 subject, with segment dimensions that match the physical subject as closely as possible. 27 This process is called scaling, and a variety of approaches have been described [23,[27][28][29][30][31][32][33][34]. 28 Finding accurate scaling is especially important when using motion capture data to 29 create muscle-driven simulations because the muscle-tendon parameters are scaled by 30 the body segment dimensions [35]. To achieve accurate kinematic results, the locations 31 of the markers on the scaled digital twin must be adjusted to account for variations 32 caused by human error in attaching the markers to the body and the variations in the 33 dimensions of human subjects [24]. This is called marker registration. Finally, the 34 positions and orientations of the body segments over time must be determined, which is 35 typically done using an optimization process called inverse kinematics [36][37][38][39][40]. Inverse 36 kinematics algorithms generally produce more accurate results when the solutions are 37 constrained by an underlying skeletal model [14,24,41]. 38 In practice, the interdependence between scaling, marker registration, and inverse 39 kinematics means that experts must follow an iterative guess-and-check procedure, 40 where they refine each of the steps several times, making small adjustments to each 41 value until a desired accuracy is achieved [42,43]. For example, increasing the length of 42 the upper arm segment in a subject's digital twin will require also changing the marker 43 registrations for any markers on the forearm and the hands, because otherwise those 44 August 17, 2022 2/18 markers would move as a result of the longer upper arm. A longer upper arm will also, 45 all else being equal, change the resulting motion found by inverse kinematics. While 46 there are best practices for conducting validation at each step [35], the process typically 47 requires subjective input from an expert and time-consuming interaction. 48 Automating the scaling and registration process has been studied before, in the 49 pioneering work of Reinbolt et. al. [44] and Charlton et. al. [45]. These papers found it 50 is possible to use gradient-free optimization methods to automatically guess body 51 segment scales and marker registrations while solving gradient-based inverse-kinematics 52 problems repeatedly in an inner-loop to evaluate optimization progress. These methods 53 require large amounts of compute time because every iterative guess the outer optimizer 54 makes about body segment scaling and marker offset must solve its own inner 55 optimization problem (inverse kinematics) to evaluate the quality of the guess, which is 56 itself computationally costly. The method of Reinbolt and colleagues [44] also finds best 57 results using a particle-based optimizer for their outer optimization problem, to combat 58 the non-convexity of the problem, but this comes at a further multiple of computational 59 cost. 60 Given the interconnected nature of body segment scaling, marker registration, and 61 inverse kinematics, one might also consider posing all three problems as a single 62 optimization problem. However, such a formulation leads to a nonconvex optimization 63 in which a global solution is not guaranteed [46]. Instead, we can only guarantee to find 64 a local optimum close to an initial guess, so providing a high quality initial guess for 65 parameters is crucial. Andersen and colleagues [47] have formulated such nonconvex 66 optimization problems, but did not address the problem of reliably finding an initial 67 guess for the non-convex optimization problem proposed.

68
This paper introduces an automated method (Fig 1) that uses a combination of 69 traditional kinematic solvers and modern bilevel optimization techniques to achieve high 70 quality results in reasonable computation time. Our method requires neither a 71 user-provided initial guess for the decision variables nor large computational resources. 72 We first apply a sequence of simple optimizations to individually approximate the initial 73 values for the body segment scales, marker registrations, and inverse kinematics [44,45]. 74 Then, rather than iteratively repeat those optimization problems hundreds of times as 75 in previous work, we apply bilevel optimization techniques to simultaneously optimize 76 body scaling, marker registration, and inverse kinematics. This process takes 77 approximately 3-5 minutes on a consumer grade laptop computer or low-end server.

78
Our approach is appropriate for use by non-experts to process large amounts of 79 motion capture data automatically. We have released the software, which we call 80 "AddBiomechanics," as open source. We also provide a cloud-based service hosting 81 AddBiomechanics, available at addbiomechanics.org, where researchers can process data 82 without downloading or installing any software. We provide a drag-and-drop interface 83 that has the potential to save thousands of hours and enable quantitative movement   Once raw marker data has been collected and uploaded (steps 1-3), addBiomechanics (step 4), replaces the time-consuming and error-prone manual step in previous workflows. Our method processes input marker data through several steps automatically. First, it finds the functional joint centers from the data (step 4.1), and then it uses the marker data and those joint centers to make an initial guess for body segment scaling and marker registrations (step 4.2). The initial guess then serves as the starting point for a bilevel optimization problem that matches the model to the experimental data as closely as possible (step 4.3). The final output is a musculoskeletal model scaled to the subject with registered markers and joint angles over time.
minimize the deviation of estimated marker positions fromx 1:T : where M is the number of markers, and f F K (q t , s, p (i) ) is the forward kinematic individually initializes each type of variable utilizing independent sources of information. 108 Specifically, we use kinematic constraints to initialize q 1:T , a geometric invariant to initialized individually, the final bilevel optimization simply makes sure that they agree 111 with one another, given the observed data and model priors.

112
Input Marker Data

113
The output of a commercial motion capture system is a series of frames, often at capture volume at the moment in time corresponding to this frame. Each 3D coordinate 117 is "labeled" with a tag corresponding to an experimental marker location on the subject 118 (e.g. "C7" for the optical marker placed on the C7 spinal segment). A full list of marker 119 tags, and their location on an idealized digital twin of the experimental subject, is 120 known as the "marker set." In practice, markers are almost never placed exactly at 121 their ideal locations, and these small deviations in experimental marker placement must 122 be accounted for during the marker registration step. Not every marker from the marker 123 set is observed on every frame, because they may be occasionally obstructed during a 124 motion capture experiment. The only assumption we make is that the tag associated 125 with an observed marker is correct.
The first term is a conditional probability of the observed data given the estimated 148 parameters. This formulation is equivalent to the standard least-squares inverse 149 kinematics objective term if we assume Gaussian noise in our optical marker multivariate Gaussian skeleton scaling prior is conditioned on that information before 154 any optimization. The third term is a zero-mean Gaussian distribution that regularizes 155 the deviation of the marker locations from their intended locationsp provided by the 156 experimenter, encoding that markers are generally placed close to their intended 157 locations, even if they do not perfectly align.

158
Note also that this is a bilevel optimization problem, because in order to evaluate 159 the quality of given skeleton scaling s and marker locations p, we need to optimize over 160 the possible joint positions q t . To efficiently solve the bilevel optimization problem, we 161 observe that at the optimal values of q t for max qt Px(x t |q t , s, p)), the gradient of the 162 inner optimization problem will be zero. Using this observation, we reformulate the 163 bilevel optimization problem as a single-level nonconvex optimization problem with 164 nonconvex constraints. For numerical stability, we minimize the negative log of the 165 above objective function: 166 min s,p,qt − ln(Px(x t |q t , s, p)) − ln(P s (s)) − ln(P p (p|p)) subject to ∂ ∂qt ln(Px(x t |q t , s, p)) = 0 At a locally-optimal point, the gradient of the objective term with respect to any of 167 the decision variables is zero, so it must be zero with respect to q t : Thus, at a locally-optimal point for the objective function our constraint must hold 169 regardless, and so we could theoretically omit it from the optimization problem without 170 loss of correctness. However, we found that including the constraint explicitly allows the 171 optimizer to converge to a high-quality solution much more quickly. 172 We could use any nonlinear optimization solver to solve Equation 3. In practice, we 173 use IPOPT [51], which is a high-quality and open source solver. However, due to our 174 problem's non-convexity, a good initial guess for the decision variables is needed to 175 produce reasonable results.

177
Prior to solving the optimization problem in Equation 3, we need to get "close-enough" 178 initial guesses for the decision variables. We do this through a sequence of optimization 179 problems shown as a flow-chart in Fig 2. We obtain initial guesses for q t , s, and p 180 individually based on independent sources of information such that the cascading errors 181 can be mitigated.
where r i is the estimated distance between the i-th markerx  The sphere-fitting approach to finding functional joint centers can yield ambiguities 201 when marker motion adjacent to a joint is primarily confined to the sagittal plane, as 202 commonly happens in locomotion. In such cases, we could move our joint center  To address these ambiguities, we formulate another optimization problem to 208 simultaneously find the joint axis and the joint center, building on Equation 5. This 209 problem is similar in spirit to the axis-of-rotation problem described in [55], but can be 210 implemented without any matrix factorizations. The goal of the axis fit problem is to 211 identify not only a joint center c, but also the direction of axis a at each frame. We also 212 estimate a fixed distance from the center for each marker, parameterized by a distance 213 u i along the axis a and a distance v i perpendicular to the axis a. The result of a 214 successful axis fit is that we capture a line at each frame, where the functional joint For each markerx  We run both sphere fitting and axis fitting at each joint. Because the axis fit is a 221 strictly more demanding problem, if we succeed in axis fitting, then we pass the axis on 222 as a constraint for subsequent problems. If we fail at axis fitting, then it must be 223 because there is out-of-plane marker motion, which means that our sphere fit is not 224 ambiguous, so then we pass along the exact joint center to subsequent problems.

225
Since Equation 5 and 6 are nonconvex, we initialize c 1:T using the joint trajectory 226 q 1:T solved by Step 2 and use any nonconvex optimization solver to optimize them.

227
Once we determine the joint center and the joint axis, we formulate another 228 optimization to initialize the scaling parameters s: where the zero vector 0 (j) indicates the local coordinate of the joint j and N is the does not exist for the joint j, we set a j to zero and remove α from the optimization.

234
Open Source Implementation

235
To facilitate adoption, we provide the algorithm as an open source cloud-based tool that 236 allows researchers to automate scaling and marker registration on their motion capture 237 data without downloading or installing any software, available at addbiomechanics.org. 238 Users can drag and drop files for automated processing, and then visualize on the web 239 or download results for analysis in OpenSim ( Figure 3). The cloud tool also allows 240 researchers to automatically generate comparisons of their own hand-scaled data versus 241 the output of the automated system. To evaluate our method on experimental data, we recruited three human experts who 254 had recently collected and processed experimental motion capture data as part of three 255 separate research studies. Expert 1 collected data for a gait perturbation-recovery study 256 using a marker set based on [56] and analyzed it using the model from [48], Expert 2 257 analyzed the sprinting data from [57] using the model from [49], and Expert 3 collected 258 data for sit-to-stand, squatting, jumping, and walking using a custom marker set and 259 analyzed it using the model from [49]. Each expert scaled and registered between 4 and 260 17 trials, for a total of 34 trials. We then used the AddBiomechanics software to automatically compute body segment scales, marker registration, and inverse kinematics 262 from the same raw experimental marker data.

263
Quantitative comparison of the solved joint angles with the ground truth joint angles 264 is another critical test of our method. However, ground truth joint angles cannot be 265 directly measured from the real life lab experiments. As a result, we generated synthetic 266 marker trajectory data from a subject model with known body segment scalings, marker 267 offsets, and joint angle trajectories and then attempted to recover the known joint angle 268 trajectories from the synthetic marker data.

269
To generate synthetic data, we used several scaled OpenSim models, registered 270 marker sets, and corresponding joint angle trajectories selected from the data created by 271 our experts. We calculated marker positions at each frame of motion, assuming all 272 markers were completely determined by the model and marker set (e.g. no non-rigid 273 effects). We then ran AddBiomechanics on the synthetic marker data, giving it the 274 generic unscaled version of the corresponding OpenSim model, and an unregistered 275 version of the appropriate marker set to optimize, and compared the recovered motion 276 to the known joint angles.

278
Our bilevel optimization algorithm to find body segment scales, marker offsets, and 279 joint angle trajectories reduced the average root-mean-square marker reconstruction 280 error by approximately 40% compared to a group of experts manually processing the 281 same motion capture data (Fig 4), and faithfully reconstructed known joint trajectories 282 from synthetic motion data (Fig 5, 6). The automated approach required no expert 283 intervention.

285
The automated pipeline had significantly lower errors than manual processing by 286 experts and a smaller standard deviation (Fig 4), where errors were measured as RMSE 287 discrepancy between experimentally measured markers and simulated markers. This 288 indicates that the automated system is able to produce combinations of scaling, marker 289 registration, and inverse kinematics that more closely replicate experimental data.

290
Our experts perceived the manual body segment scaling and marker registration to 291 be labor intensive: our three experts reported that manual processing required 292 approximately one working day per subject. Average computation time for a trial 293 processed with AddBiomechanics was 3 minutes on a consumer laptop.

294
Low marker reconstruction error is a necessary but not sufficient condition to 295 establish the correctness of our system: it is possible to score well on marker RMSE and 296 still have a model of the subject with hands, feet, and head that are scaled in an 297 unrealistic way. This is possible because unrealistic bone scaling can be compensated for 298 by adjusting marker registrations. We thus improved our confidence in the correctness 299 of the automated scaling by fitting a multivariate gaussian over many individuals' body 300 measurements on the publicly available ANSUR II anthropometric dataset [50], and 301 measuring how automated scaling compared with expert scaling on that gaussian. For 302 each subject, we measured P (automated)/P (expert), and found that across our 303 comparison dataset the ratio was 1.65 ± 0.43. The value greater than 1 indicates that 304 autoscaling consistently produced a more probable skeleton than expert scaling did. We found that AddBiomechanics was able to recover the ground truth joint angles from 308 synthetic marker data to an average of 0.79 deg RMSE and a maximum of 1.0 deg 309 RSME (computed over all joints in a trial together) for motions including a standing 310 calibration pose, walking, sprinting (Fig 5), and jumping (Fig 6). Representative joint 311 angles are shown for sprinting and jumping, with the full results included as Table 1 in 312 the supplementary data.

314
In addition to being computationally efficient, our method is able to recover joint angles 315 on realistic synthetic data at state-of-the-art accuracy and outperform human experts in 316 marker reconstruction error on experimental data. For comparison, the method in [24] 317 assumed that all the body segment scalings were known to the algorithm and only 318 attempted to find the marker offsets and the joint angles, and it reported 1.21 degree 319 joint angle RMSE. Our method solves a harder problem, because it must recover 320 segment scaling information from the data, and yet achieves similar results, with joint 321 angle RMSE ranging from 0.46 degrees to 1.05 degrees. Our findings on marker 322 reconstruction error on experimental data are also consistent with previous work on 323 optimization based approaches, which all outperform human experts when fitting a 324 model to the same data [24,44,45,[58][59][60]. However, previous approaches required large 325 amounts of compute time, were limited to one specific skeleton, or only addressed part 326 of the body segment scaling and marker registration problem.

327
Our method goes from labeled marker trajectories to a scaled and registered 328 musculoskeletal model and corresponding human motion in a few minutes on a standard 329 consumer laptop or a low-end server. We also provide a web version where processing and if the angles of the hips and spine are appropriately adjusted then the markers will 348 still closely match the target data. Fig 6 provides an example of this happening on the 349 jumping subject in our synthetic data experiments. If this effect is observed in practice, 350 AddBiomechanics users can leverage the fact that the optimizer will prioritize solutions 351 that move the markers as little as possible, and adjust the marker starting locations on 352 the bones to more closely match the experimental placement. Secondly, the optimizer 353 applies a statistical prior to body segment scales to bring them more in-line with 354 population statistics as represented by the ANSUR II anthropometric dataset [50]. If 355 the optimizer can find a way to fit the marker data with a skeleton that is more likely 356 to exist in the ANSUR II population (such as by tilting the pelvis forward 2 degrees), it 357 will choose that one, even if the "true" underlying skeleton was slightly different. The 358 data in ANSUR II is large and detailed, but was collected from active duty military 359 personnel, and so is not reflective of many patient populations. A broader 360 anthropometric dataset could help address this limitation.

361
Our method also does not incorporate any dynamic quantities, ignoring ground 362 reaction force data even when it is available. Ground reaction forces could provide very 363 August 17, 2022 13/18 useful constraints on the acceleration of the center of mass of the subject, but exploiting 364 these constraints to produce more accurate models is left to future work.

365
The open source cloud-based tool addbiomechanics.org reduces the barrier for 366 researchers to integrate the optimization approach into their workflow. Researchers can 367 get their data processed for free, if they agree to share the resulting anonymized motion 368 data with the scientific community under a creative commons license. While large-scale 369 public datasets of human motion already exist (e.g. [20]), they are focused on computer 370 graphics applications and do not provide biomechanically accurate models of the 371 experimental subjects. Because AddBiomechanics is focused on human motion 372 biomechanics, the cloud application also encourages sharing ground reaction force data 373 along with kinematics, even though the optimizer does not currently leverage this data 374 to increase the accuracy of its solutions. We hope AddBiomechanics will enable creation 375 of a large-scale public dataset of accurately modeled human motion biomechanics,