1 Introduction
The human body in computer graphics has evolved from a stick figure with torque actuators to a muscledriven skeleton. The musculoskeletal model for physicsbased simulation includes many parts and parameters that require careful orchestration to make it anatomically plausible. Designing such a musculoskeletal model is laborintensive and handcrafted models are often inaccurate to actuate an anatomical body in physically simulated environments. We present a new method that retargets a reference musculoskeletal model to new bodies of different sizes, body proportions, muscle capability, and joint Range of Motion (ROM).
The musculature is a system of muscles and tendons actuating the skeleton. The motion of the skeleton is the result of the harmonious coordination of many muscles. Muscles are correlated with each other in various ways. A group of muscles works together to actuate a single joint and the ROM around the joint is also determined by multiple adjacent muscles. Some adjacent muscles share a tendon and thus the contraction of one muscle may affect the activation of the others. A single (multiarticular) muscle may actuate two or more joints simultaneously.
The functional role of each individual skeletal muscle is thoroughly studied and well documented in anatomy textbooks [gray2005]. Each muscle is responsible for particular movement (e.g., flexion/extension, adduction/abduction, internal/external rotation) of body parts. For example, the active contraction (shortening) of the rectus femoris (a muscle in the thigh) results in knee extension and hip flexion. Even if the muscle is inactive, its background elasticity generates passive force to prevent the adjacent joints from excessive movement, which influences the ROM around the joints.
Our focus is on adapting the musculature of a reference anatomical model to new bodies while preserving the functionality of the musculature. The muscles are expected to perform the same functional roles and have the same ROMs in the new body after retargeting. The key parameters that affect muscle functions and joint ROMs are the length of musculotendons and their geometric routing paths. The core of our retargeting algorithm is a numerical solver that optimizes the key parameters throughout the entire body. The algorithm allows us to create a variety of musculoskeletal models including exotic creatures.
We also implemented a human body editing system that allows the user to interactively edit the height, width, bone length, bending and torsion angles, and joint ROMs of the skeleton. The musculature of the anatomical model is automatically retargeted to fit the new skeletal body. Additionally, the user is allowed to edit physics parameters, such as the mass and inertia of individual body parts, and make the body stronger or weaker by adjusting Hilltype muscle parameters, such as forcelength curves and forcevelocity curves. The physics and muscle parameters are exploited in a physically based simulation of muscleactuated characters.
The simulation and control of underactuated biped locomotion under gravity have been studied for decades in computer graphics and robotics. The recent progress in deep reinforcement learning research has been very successful in continuous control problems and made it possible to simulate dynamic human activities without compromising fundamental physics laws
[peng2018deepmimic, lee2019scalable]. We will demonstrate that our muscleactuated anatomical model is compatible with stateoftheart simulation algorithms. Our retargeted characters with extreme body proportions can be physically simulated and controlled to walk, run, kick, and jump while maintaining balance under gravity.2 Related Work
There has been a stream of studies for reproducing natural human motion by incorporating the increasingly more accurate models of human anatomy and the mechanics of anatomical structures. Musclebased anatomical modeling and simulation have been extensively explored and exploited in Computer Graphics and Biomechanics.
Muscle Modeling
The anatomical model includes an articulated skeleton driven by its musculature. The Hilltype muscle model [hill1938heat, zajac1989muscle, millard2013flexing] has been broadly adopted to encode the nonlinear contraction dynamics of muscles. The geometry of a muscle is often simplified into a sequence of line segments for the efficiency of computation. There have been continuous efforts in Computer Animation and Biomechanics to demonstrate the computational plausibility of volumetric FEM muscles [lee2009comprehensive, si2014realistic, fan2014active, berranen2014real]. In vivo estimation of muscle parameters relies primarily on medical 3D imaging. Matias et al. [matias2009transformation] studied the estimation of musclebone attachments based on bony landmarks. Levin et al. [levin2011extracting]
studied the estimation of muscle fiber directions from diffusion tensor images. Arnold et al.
[arnold2010model] and Holzbaur et al. [holzbaur2005model] provided a comprehensive reference of muscle modeling parameters for lower and upper limbs. Anatomical modeling for specific body parts, such as face [sifakis2005automatic, ichim2017phace], feet [park2018multi], hands [Sueda2008Musculotendon, sachdeva2015biomechanical], shoulders [van1994finite, kaptein2004estimating, maurel2000human], tongue [stavness2012automatic], and jaw [zoss2018empirical], has been explored for last two decades. Comprehensive fullbody musculoskeletal modeling and simulation systems are also available freely [Delp2007opensim] and commercially [Damsgaard2006anybody]. Muscle routing plays an important role in muscleactuated fullbody simulation. Many modeling approaches, such as conditional waypoints [delp1990interactive], obstacle sets [garner2000obstacle], surface SSD [murai2016anatomographic], and wrapping surfaces [rajagopal2016full], have been exploited to improve the estimate of muscle lengths during joint motion.Joint Modeling
The joints of a skeletal model are often simplified as either revolute (1 degree of freedom) or ballandsocket (3 degrees of freedom) joints, though anatomical joints are structurally more complicated than simple mechanical joints
[lee2008spline]. Maurel and Thalman [maurel2000human] designed a skeletal model of human shoulders and represented the flexibility of the ballandsocket joint by a combination of the sinus cone and an interval of the arm twist angle. Lee [lee2000hierarchical] studied a general approach of describing the range of motion around a joint based on quaternion halfspaces and their boolean combination. Akhter and Black [akhter2015pose] collected a motion capture dataset that explores a wide range of human poses, and learned the posedependent range of motion around each joint. Jiang and Liu [jiang2018data]represented the boundary of valid human joint configurations by using a fullyconnected neural network. Measuring the joint range of motion and muscle lengths are of great clinical interest. A variety of clinical examination techniques are performed frequently in clinical practice
[Reese2017].MuscleDriven Simulation and Control
Recent accomplishments in Computer Animation made it possible to reproduce realistic human motion in physicsbased simulation [liu2012terrain, coros2010generalized, lee2010data, sok2007simulating]. The emergence of deep reinforcement learning accelerated the progress and successfully simulated skillful actions, such as jump, flip, cartwheel, basketball dribbling, and even aerobatic flapping flight [peng2018deepmimic, liu2018learning, yu2018learning, won2018aerobatics, lee2019scalable]. Muscledriven anatomical simulation poses further challenges of dealing with the complexity of anatomical modeling and scalability in physicsbased simulation. Wang et al. [wang2012optimizing] demonstrated a muscleactuated biped with eight musculotendon units on each leg. The biped was essentially twodimensional because all muscles are aligned in the sagittal plane. Geijtenbeek et al. [geijtenbeek2013flexible] improved the flexibility in biped design by allowing offsagittal muscles and optimizing muscleattachments sites for locomotion control. Comprehensive 3D musculoskeletal models with up to 120 musculotendon units were successfully simulated and controlled using a twolevel control architecture that consists of lowlevel muscle coordination based on quadratic programming and higherlevel gait modulation based on stochastic optimization [lee2014locomotion]. Lee et al. [lee2018dexterous] demonstrated that volumetric FEM muscles are viable for dextrous manipulation, such as juggling, that requires a high level of control precision. Nakada et al. [Nakada2018biomimetic] built a comprehensive fullbody neuromuscular system to learn biomimetic sensorimotor control of human animation. The muscleactuated control system proposed by Lee et al. [lee2019scalable] employed a hierarchical deep reinforcement learning algorithm to cope with the full details (46 joints and 346 muscles) of the musculoskeletal model and successfully simulated dynamic motor skills. They also demonstrated preoperative and postoperative simulation of pathologic gaits. We built our musculoskeletal model on top of their simulation system, which provides favorable features such as analytic differentiation of muscleactuated dynamic systems, muscle routing based on LBS (linear blend skinning), and efficient handling of joint ROMs based on LCP (linear complementary problem).
Anatomy Transfer and Retargeting
Constructing a highfidelity anatomical model is timeconsuming and laborintensive. Therefore, there has been a series of studies that retarget a reference anatomical model to a range of target bodies. Dicko et al. [ali2013anatomy] transferred the geometry of internal anatomy from an input template to a target body while maintaining anatomic constraints. Saito et al. [saito2015computational] generated a spectrum of human body shapes with various degrees of muscle growth. Kadleček et al. [kadlevcek2016reconstructing] generated a personalized anatomic model by retargeting a template model to fit fullbody 3D scans while accounting for body shape variations and pose variations. Our work is on the line of these approaches with emphasis that our retargeted model is optimized to function as closely as possible to its reference model in the realm of muscledriven simulation and control.
3 Musculoskeletal System
Our musculoskeletal model has an articulated skeleton with 5 revolute joints, 13 ballandsocket joints, and 6 DoFs (Degrees of Freedom) at the skeletal root yielding 50 DoFs in total (see Figure 2). Our work is focusing on simulating human motion driven by four limbs and spinal joints. The model includes 282 musculotendon units corresponding to skeletal muscles, which are attached to bones on each end by tendons. The attachment site on the proximal end is called the origin of the muscle, while the attachment site on the opposite end is called its insertion. Muscle contraction pulls the attached bones and moves the joint between them. Our model comprehensively includes skeletal muscles that contribute to the motion in any limb joints. Since our model has no fingers or toes, the muscles that originate and insert within hands or feet are omitted.
3.1 Muscle Model
The Hilltype muscle is a threeelement scheme that models the muscle contractionforce relation (see Figure 3). A musculotendon unit is defined by an active contractile element (CE), a passive parallel element (PE), and a passive serial element (SE). The active force of the contractile element comes from muscular contraction. All skeletal muscles have a rest length. When the muscle is stretched to its ideal length, it can maximize muscular contraction. The serial element represents the tendons on either side of the contractile segment of the muscle. The tendon is a tough band of fibrous connective tissue that connects muscle to bone. The parallel element models the background elasticity of the muscle. The parallel element is responsible for muscle passive behavior when it is stretched.
Let , and be muscle fiber length, tendon length, and the total muscletendon length, respectively. Optimal fiber length is the length of muscle fibers when it develops maximum isometric force. Let be tendon slack length. Pennation is the angle between the muscle fibers and the tendon axis. Assuming a quasistatic setting, the tension in muscle fibers is a function of muscle length and its activation , while the tension in the tendon is a function of only muscle length. The tension in the three elements holds
(1) 
where and are normalized muscle and tendon lengths, respectively. Here, we assume that the pennation angle is constant regardless of muscle contraction.
Actual muscles wrap around bones and soft tissues during joint motion and they are also constrained by surrounding ligaments and skin tissue layers. Therefore, the straight line segment between the origin and insertion sites of a muscle is a lousy approximator of its actual length. Better approximators can be achieved by exploiting geometric proxies, such as waypoints [delp1990interactive], obstacle sets [garner2000obstacle], and wrapping surfaces [rajagopal2016full], that approximate muscle deformation during joint motion. We use waypoints as geometric proxies and expressed the location of a waypoint relative to nearby bones using the idea of linear blend skinning (LBS) to improve the accuracy and flexibility over boneattached waypoints (see Figure 4).
(2) 
where is the transformation matrix of th bone, is skinning weight, and is the coordinates relative to th bone. Since the derivatives of LBSbased waypoints can be derived in an analytic form, they are wellsuited for efficient dynamics simulation with inequality constraints [lee2019scalable].
3.2 Parametric Skeleton Model
The skeletal body is a parametric model that can generate a wide variety of human body shapes (see Figure
5). The model has a reference shape and a set of shape parameters. The trunk and four limbs are parameterized to have variations in their length, size, and alignment relative to the reference shape. The limb bones including Femur, Tibia, Humerus, Ulna exhibit an elongated shape with its proximal head, long shaft, and distal head. The limb bone has four parameters: proximal head scale, distal head scale, shaft elongation, and torsion. Even though the trunk consists of many bones, the trunk is considered as a lumped body and parameterized by three parameters: elongation, expand, and bend. The hands and feet have one parameter for scaling their size. We use a geometric deformation method [kadlevcek2016reconstructing] to produce parametricallyvaried shapes from the reference model.We scaled physics and muscle parameters relative to geometric parameters inspired by the work of Hodgins and Pollard [hodgins1997adapting]. Assuming uniform scaling by factor
in all dimensions, they suggested to scale time, force, mass, moment of inertia, velocity, stiffness, and damping by factors
, , , , , , and , respectively. This scaling rule provides a rough guideline as to how physics and muscle parameters should be adjusted relative to body scaling. The user can edit the mass and inertia of each individual body part and scale the forcelength curve of each muscle following the guideline. The guideline is not strict since stable dynamics simulation can be often accomplished for a wide range of parameters.4 Joint Range of Motion Modeling
Many articulated characters in computer graphics have the ROM around joints described by an interval of min/max angles for each DoF. This approach is inherently limited because the posedependence of joint ROMs cannot be captured in perDoF intervals. An alternative approach is to estimate a highdimensional function from a motion capture dataset, which fits the boundary of valid poses. The isValid function
takes a fullbody pose as input and outputs a boolean value to answer whether the pose is valid or not. The function estimation requires nonlinear regression over a large collection of fullbody poses
[akhter2015pose, jiang2018data]. The musculoskeletal system readily provides an isValid function without requiring nonlinear regression because the musculoskeleton represents the fundamental mechanism of joint movements and their limits.4.1 MuscleInduced ROM
The joint ROM is affected by many anatomic factors, such as bony structures, surrounding ligaments, joint capsules, and muscular tension. We are particularly interested in the factors induced by muscles. The elasticity of muscle fibers represented by a parallel element determines the distance a joint can move to its full potential in the direction of the muscle. This extension of muscle stretching is described by an implicit function , where is the generalized coordinates of the fullbody skeletal pose, and is its dimension. The constraint induced by passive fiber tension is
(3) 
where muscle length is a function of skeletal pose. Note that there exists a direct mapping between skeletal poses and muscle lengths. Given a pose, the lengths of all muscles are uniquely determined and vice versa. We assume that passive elongation reaches its limit when muscle tension is larger than maximum isometric tension in muscle fibers by a certain margin (see Figure 3(right)). Deriving from muscle contraction dynamics [thelen2003adjustment], we set the maximal muscletendon length by
(4) 
where and . We can also express other factors induced by bones and ligaments as implicit functions (e.g., the constraint by the kneecap can be denoted by a function of skeletal pose). The collection of all constraints defines an isValid function.
(5) 
It is often the case that multiple muscles cross a joint and all those muscles affect the ROM in the joint. Among those muscles, some are multiarticular meaning that it crosses multiple joints. The presence of multiarticular muscles makes joint ROM posedependent. For example, we have a wider range of motion in the hip when the knee is flexed than it is fully extended. Our isValid function based on the musculoskeletal model inherently captures the posedependence of ROM.
4.2 Muscle Length Estimation
Even though carefully specifying muscle parameters would yield appropriate joint ROMs, parameter tuning is laborious and timeconsuming. We estimate muscle parameters (in particular, muscle fiber and tendon lengths) automatically such that the isValid function fits the boundary of fullbody poses captured in a motion dataset. Specifically, we used a data set collected by Akhter and Black [akhter2015pose], which includes an extensive variety of stretching and yoga poses performed by trained athletes and gymnasts. In the dataset, we can find a fullbody pose that maximizes the length for each muscle (see Figure 6). Assuming that muscle fiber and tendon length ratios are invariant, solving Equation (1) and Equation (4) gives us optimal fiber length and tendon slack length at the maximal length . This estimation of muscle lengths guarantees that all poses in the dataset are valid.
The dataset actually provides us with the lower bound for the maximal length of individual muscles and therefore it tends to underestimate the maximal lengths. For example, the range of knee extension is bounded by the bone and ligament structures around the knee rather than muscle passive force. It is highly likely that knee flexor muscles (e.g., Hamstring and Gastrocnemius) do not extend to their maximal length when the knee is straight. Hence, the dataset does not provide a precise estimation of the maximal length of the knee flexors. We address this underestimation issue by measuring joint torques at several keyposes, for which muscles should not develop passive force. The keyposes include standing, zerogravity, and Tposes. We gradually increase and of muscles if they develop passive force until the force magnitude is reduced below a certain threshold.
4.3 ROM Editing
Muscle lengths thus estimated yield an isValid function , which serves as a reference ROM model. The reference model represents the level of flexibility of athletes, which is far higher than the flexibility level of average, nonathletic individuals. We would like to be able to edit or modify the reference model to create various individualized ROM functions. ROM editing involves improving or reducing the ROM in a particular joint and shifting the range. Note that ROM editing is a supplementary step that provides flexible and convenient user control over the modeling and retargeting procedures. The ROM editing step can be skipped if it is not necessary.
Since is a highdimensional function defined implicitly depending on many factors and complex anatomical structures, it is not practical to edit the function directly. Alternatively, we define a new, modified function indirectly with transformation over the input pose. The action of affects the ROM inversely. For example, if the transformation exaggerates poses, the ROM with the modified function becomes narrower than the reference model because the test pose will be exaggerated first and then passed to the reference isValid test .
We can define the transformation for each joint or each DoF independently. For a revolute joint, the ROM is simply an interval , where is the average of the joint angles in the dataset. The transformation of joint angle in pose results in scaling of the interval followed by shifting. Note that we do not need to know the value of to define the transformation. Similarly, we can define ROMediting transformations for ballandsocket joints. The configuration of a ballandsocket joint is a threedimensional rotation, which can be decomposed into a twist about the principle (bone shaft) axis in the reference pose followed by another rotation about its orthogonal axis (see Figure 7). The decomposed configuration is , where is a twist angle and
is a unit vector indicating the direction of the bone shaft. The ROM of
is a simple interval, while the ROM of is a cone with its center direction estimated from the dataset. The scaling and shifting transformation on the cone is also defined in a straightforward manner and can be exploited for ROM editing. Through simple user interfaces, the user can edit ROMs interactively by modulating scale and shift factors.5 Musculature Retargeting
In this section, we describe a retargeting algorithm that adapts the fullbody musculature to the change of the skeleton (see Figure 8). Our parametric skeletal model can generate a rich variety of human body shapes. The goal of retargeting is to make the varied models viable for physicsbased simulation and control.
The functional roles of skeletal muscles are classified by the plot of their length changes and the direction of muscle forces acting on the attached bones. For example, a muscle is called a flexor (or extensor) of the joint if the muscle shortens as the joint flexes (or extends). Table
1 shows examples of muscle actions during various joint motions. If the muscle is functionally significant for causing a particular joint motion (either as an agonist or an antagonist), its lengthangle curve is monotonically increasing, monotonically decreasing, or close to unimodal (see Figure 9). Therefore, the lengthangle curves characterize the functional roles of muscles. Naïve scaling of the musculoskeletal model could alter the monotonicity and modality of the lengthangle curves, muscle force directions at attachment sites and joint ROMs, eventually affecting the functionality of muscles in the new body.Joint Motion  Active Muscles 

Hip Flexor  Psoas major, Iliacus, Rectus femoris*, 
Tensor fasciae latae*, Sartorius*  
Hip Extensor  Gluteus maximus, Biceps femoris*, 
Semitendinosus*, Semimembranosus*  
Hip  Gluteus medius, Pectineus, Tensor fasciae latae* 
Medial Rotator  
Hip  Quadratus femoris, Obtrator externus, 
Lateral Rotator  Biceps femoris* 
Hip Abductor  Gluteus minimus & medius, 
Tensor fasciae latae*  
Hip Adductor  Adductor longus & brevis & magnus 
Knee Extensor  Vastus medialis & intermedius & laterlias, 
Rectus femoris*  
Knee Flexor  Sartorius*, Biceps femoris*, Semitendinosus*, 
Semimembranosus*, Gastrocnemius*, Plantaris*  
Foot  Gastrocnemius*, Plantaris*, Soleus, 
Plantarflexor  Tibialis posterior 
Foot  Tibialis anterior, Extensor hallucis longus, 
Dorsiflexor  Extensor digitorum longus, Fibularis tertius 
Musculature retargeting involves three categories of parameters. The parameters in each category have the desired impacts on muscle functionality. Our algorithm solves for three categories of parameters sequentially.

Optimize muscle routing to modulate lengthangle curves and force directions, where the muscle routing is expressed by LBS coordinates of muscle waypoints.

Optimize muscle fiber and tendon length ratios to modulate the angle of peak torque at joints.
Since we already discussed the second step in the previous section, this section focuses on the first and third steps which are formulated as nonlinear optimization.
5.1 Muscle Routing Optimization
Waypoints play an important role in estimating muscle lengths during joint motion. While the use of LBS coordinates improves the estimation accuracy substantially, manual tuning of LBS coordinates is laborious. Our algorithm decides LBS coordinates of waypoints such that the energy function is minimized:
(6) 
The first term encourages muscle force directions through the polyline approximations of muscles to be preserved
(7) 
where is the index of muscles, is the index of waypoints of the muscle (including the origin and insertion), is the type of joint motion in which the muscle participates, and is the normalized joint angle. and are normalized force directions in the reference and retargeted bodies, respectively, at point to .
The second term encourages the qualitative characteristics of lengthangle curves to be preserved in the retargeted body. Assuming that the lengthangle curve is monotonic or unimodal on interval , the curve can be characterized by three parameters: , , . The curve has its maximum value at , minimum value at , and the difference between the maximum and minimum value is (see Figure 10). The second term is defined using the characteristics parameters.
(8) 
We use a gradient descent algorithm to minimize the energy function. The gradient of the energy function is estimated by finite difference and line search along gradient direction accelerates the algorithm. In our example, and .
The performance and quality of gradientbased optimization depend on the initial guess of the optimized parameters. During the scaling of the skeleton, each waypoint moves together with the skeleton anchored to the nearest point on the bone. The anchor position relative to the skeleton serves as an initial guess of the optimization.
5.2 Angle of Peak Torque
Joint torque is the sum of muscle torques crossing the joint. Peak torque usually occurs in the joint ROM when the muscle fibers are close to their optimal length. For example, knee and elbow joints generate their maximum bending torque when they are moderately flexed. The angle of peak torque is related to muscle fiber and tendon length ratios.
Peak torque angle adjustment is a constrained optimization process based on gradient descent. The optimization objective is that all joints in the retargeted model have their peak torque at the same angles as in the reference model. The optimization parameters are fibertendon length ratios at all muscles. The ratios are constrained to vary within from their original values.
(9) 
where is the angle of peak torque in the reference model and is the net torque in joint motion when its normalized angle is .
6 Experimental Results
We can generate a wide spectrum of body proportions through interactive user interfaces. The user interface provides a set of slide bars to manipulate body parameters and allows the user to edit joint ROMs. In this section, we show experiments conducted with four representative characters: Human, Hulk, Dwarf, and Alien (see Figure FunctionalityDriven Musculature Retargeting and Table 2). The Human is a reference model we created based on a human skeleton geometry. We annotated origins, insertions, and waypoints of all muscles and tuned their parameters to create a model that can be used for physicsbased simulation. The Hulk has a large muscular body with long arms and short legs and the trunk leans slightly forward. The Alien has a short torso and extremely long limbs, while the Dwarf is the opposite with a long torso and short limbs. The musculature of the Human is retargeted to the others to build musculoskeletal dynamic systems.
Reference  Alien  Hulk  Dwarf  

Trunk  Elongate  1.0  1.0  1.0  2.0 
Expand  1.0  1.0  2.0  0.6  
Mass  1.0  1.0  1.0  2.0  
Femur  Elongate  1.0  2.5  0.8  0.6 
Torsion  
Mass  1.0  2.5  0.8  0.6  
Tibia  Elongate  1.0  2.5  0.8  0.6 
Torsion  
Mass  1.0  2.5  0.8  0.6  
Humerus  Elongate  1.0  2.5  1.8  0.6 
Torsion  
Mass  1.0  2.5  1.8  0.6  
Ulna  Elongate  1.0  2.5  1.4  0.6 
Torsion  
Mass  1.0  2.5  1.4  0.6  
Neck  Elongate  1.0  2.0  1.0  1.0 
Mass  1.0  2.0  1.0  1.0 
6.1 ROM Construction
We constructed the ROMs of our reference model from a motion capture data set [akhter2015pose]
. We excluded corrupted outliers from the dataset, doubled the data by mirror reflection, and subsampled at a 1/10 ratio to take approximately 70,000 fullbody poses. The mirror reflection makes the ROMs symmetric. The maximum lengths of all muscles are estimated from the pose set and further relaxed to address the underestimation issue as explained in Section
4.2. Muscle lengths thus estimated represent posedependent joint limits. Figure 11 shows that the hip ROM is wider when the knee is flexed than it is straight. This observation is due to the presence of biarticular muscles. The hamstring muscles (Biceps femoris, Semitendinosus, and Semimembranosus) are all biarticular, originating at the bottom of the Pelvis, crossing both the hip and knee joints, and inserting at the head of the Tibia. The hamstring muscles therefore involved in both hip extension and knee flexion. The hamstring muscles are extended when the knee is straight and thus the potential for the hip flexion is limited.The ROM around a joint is strongly influenced by the change of the body shape. Figure 12 shows how the longitudinal elongation (in the range of 70% to 130% of the original length) and torsion (in the range of 30 to 30) of the femur affect the hip ROM. Naïve linear scaling of the muscles attached to the femur results in undesirable changes in the hip ROM, which is discretized into cells. Each cell has a boolean value. We measured the error rate by the percentage of false positives and true negatives in the cells. Without retargeting, the error rates are in the range of 5.8% to 75.4%. The maximum error drops down to 1.7% after retargeting.
Our system provides a simple user interface that allows the user to edit the ROM of any joint interactively. The user can make the ROM wider, narrower, or shifted on the space of conic and twist rotations. The editing operations are defined by transformations as explained in Section 4.3. Our system recalculates the muscles at relevant cites such that the muscleinduced ROM matches the user’s intention approximately. In Figure 13, the hip ROM tilted forward by 30 and shrank in its width on the cone by a factor of 0.63 by userspecified transformations. Recalculating muscles lengths approximates the user editing in the error rate of 7.3%.
6.2 Muscle Routing Evaluation
We found that LBSbased waypoints work noticeably better than fixedanchor type waypoints in the previous studies [geijtenbeek2013flexible, lee2014locomotion]. Nonetheless, transplanting the musculature of the reference model to varied skeletons without retargeting often results in undesirable distortion in lengthangle curves and disturbance in muscle force directions. Figure 14 shows the lengthangle curves of the reference and varied models with and without retargeting. The qualitative and quantitative changes of muscle curves can be observed in many major muscles (e.g., Iliacus for the Hulk and the Alien, and Gluteus Medius of the Dwarf in the lower limbs). The influence on force direction is broader. The force direction at almost all muscle origins and insertions are disturbed by the change of body proportions(see Figure 16). Both curve distortions and force direction disturbances affect the maneuverability and ROM of joints. Our muscle routing optimization addresses both problems simultaneously to build a retargeted model that functions as closely as possible to the reference model.
Figure 15 shows that many muscles change their functional role as the limbs lengthen or shorten in the range of 60% to 250% of their original length. Each value indicates the percentage of functional disorders (increasingtodecreasing and vice versa) in the lengthangle curve during joint motion with and without retargeting. Some muscles (e.g., Tibialis Anterior, Tibialis Posterior, Pectineus, and Semimembranosus) change their function with shorter limbs, while functional disorders occurs for some muscles (e.g., Rectus Femoris and Vastus muscles) with longer limbs. For example, the functional disorders for Pectineus muscles occur with shorter (60% of its original length) legs in the half of the range of hip extension and flexion. Our retargeting algorithm alleviate this issue substantially with mild body shape variations, though the problem is not completely fixed with extreme body proportions (e.g., longer limbs by 250%).
6.3 MuscleCoordination using QP
We tested the usability of our retargeted models in our muscleactuated simulation system. We briefly summarize the dynamics formulation and its control methods based on Quadratic Programming (QP). The Lagrangian dynamics of the musculoskeletal model is defined by
(10) 
where , , and are generalized coordinates, the mass matrix, and Coriolis/centrifugal force, respectively. and are muscle and constraint forces, respectively. Jacobian matrices and map the forces to generalized coordinates. is external forces. We compute the constraint forces by solving linear complementary conditions:
(11) 
where is the constraint velocity representing the rate of change of the constraint.
Exercising a particular joint is an underspecified control problem, since there are more muscles crossing the joint than minimally required to actuate the joint. Solving the control problem requires muscle coordination to decide activation levels at the muscles that achieve desired joint accelerations
. A typical solution method based on QP formulates the problem such that(12)  
subject to  
where is a vector concatenating muscle activations. This formulation solves for muscle activations while satisfying the Lagrangian equation of motion and regularizing large activations.
The formulation allows us to move the joint in an arbitrary direction by specifying desired joint accelerations if the musculoskeletal model is carefully calibrated. The QP with a varied model without retargeting often fails to solve for muscle activations because force directions are degenerate and unable to span the solution space. Our retargeting algorithm recovers the force directions to better condition the problem (see Figure 16). The results may be best viewed in the supplementary video.
6.4 Learning motor skills using DRL
Our characters also learned to perform challenging fullbody motor skills. Deep Reinforcement Learning has recently shown its promise in the control of physicallysimulated bipeds and muscleactuated skeletal figures. The goal of trajectory tracking control is to learn a control policy that imitates reference motion data in physics simulation. It has been shown in previous studies that the control policy represented by deep networks can reproduce highlydetailed human movements [peng2018deepmimic, lee2019scalable].
Let be the dynamic state of all musculotendon units at time . The state includes muscle lengths and the rate of their length changes. Let be the level of muscle activations and be a reward of taking action at state . The agent interacts with the environment by taking actions according to its control policy . The objective of policy learning is to find a control policy that maximizes the expected cumulative rewards
(13) 
where is a discount factor.
The reward evaluates how well joint angles and endeffector positions match the reference trajectory. The reward is defined by a multiplication of a joint angle match reward and an endeffector match reward, where
(14) 
where are the joint configurations in generalized coordinates and are endeffectors positions. and are the indices of joints and endeffectors, respectively. Our implementation of the endeffector reward includes the positions of both hands, both feet, and the head relative to the moving coordinate frame of the skeletal root (pelvis). The hat symbol indicates the positions and joint configurations taken from the reference trajectory. The joint configurations are represented by unit quaternions. The geodesic distance measures the shortest rotation angle between two joint configurations [lee2008geometric]. We employed a hierarchical RL algorithm proposed by Lee et al. [lee2019scalable] to learn muscleactuated control policies mimicking reference motion data.
We collected motion capture data available on the web and retargeted each motion clip to our kinematic models using Autodesk MotionBuilder. The DRL algorithm is able to learn muscleactuated control policies successfully with mild body shape variations. The DRL algorithm can also cope with larger variations to a certain extent. The highlyvaried characters were able to successfully imitate motor skills, including walking, running, jumping, kicking, and dancing while maintaining their balance (see Figure 17 and the results may be best viewed in the accompanied video).
6.5 Moment Arm
The moment arm of a muscle to a joint is a quantity of interest often referred to by biomechanics literature. The moment arm is the perpendicular distance from the joint to the line of muscle tension force and serves as a measure of effectiveness that the muscle has on joint torque. The moment arm about a particular axis of action at joint angle can be calculated by [lieber1993skeletal]
(15) 
where is the effective torque applied to joint by muscle and is the tension force at the attachment site generated by muscle activation.
Moment arm is not an invariant factor we want to preserve across different body sizes and proportions. We rather wish to have moment arms to be adjusted appropriately during the retargeting process. In Figure 18, the anglemoment arm curves in the retargeted models differ from thes curves in the reference model. The Alien has larger moment arms for hip flexors and extensors to compensate for longer Femurs. As the femurs gets longer, larger moment arms are preferred to use flexor and extensor muscles more efficiently. The retargeting algorithm reduces the moment arms for the Hulk and the Dwarf on the contrary to compensate for their short Femurs.
6.6 Individualized Anatomy Modeling
We would like to create individualized anatomic models from medical images and physical examination data. There are standard physical examination procedures of measuring joint ROMs in physical practice, such as popliteal angle, Staheli, and Silfverskiöld tests. In this example, we used biplanar Xray images of an anonymous person to create a skeletal model and demonstrated how ROM editing can be performed according to physical examination data.
The biplanar image data include lateral and anteroposterior Xray images acquired simultaneous at precisely calibrated view angles (see Figure 19). The calibrated orthographic images allow us to reconstruct 3D points from a pair of 2D feature points on the images. We selected joint points on the Xray images and scaled our parametric skeleton model to match the reconstructed 3D joint points. Musculature retargeting to the scaled skeleton creates an individualized musculoskeletal model.
The ROM measurements in physical examination data can be incorporated into our anatomic model. For example, the unilateral popliteal angle is a test for measuring hamstring tightness, which is related to the ROM in the knee. The test measures the range of a knee angle when the subject is in the supine position with its leg flexed to 90 at the hip and the knee is extended passively. We simulated this procedure to have our musculoskeletal model take the specified pose (see Figure 20). Interactive ROM editing allows us to change the knee ROM to match any target range by deciding the scaling and shifting factors as explained in Section 4.3.
7 Discussion
Our attempt to build a musculoskeletal model for a variety of body shapes is an important step toward the goal of expanding the applicability of anatomic human modeling and simulation. Nonetheless, our current framework still has several limitations in its scope of applicability and biomechanical accuracy. Human anatomy is an extremely complicated system. There has been a trail of research works that incrementally added details to the computer model of human anatomy. Currently, our model ignores a lot of important anatomic features and phenomena, such as ligaments and contact coupling. Joint ROMs are largely influenced by the collision between bony features and tensions generated by ligaments around the joints. Muscle paths wrapping around bones and other muscles will best be approximated by contact coupling between anatomic features if their highresolution geometry and volumetric FEM simulation are provided [lee2018dexterous]. The use of waypoints is a crude, yet computationallyfeasible approximation.
Our demonstration of reconstructing a musculoskeletal model from biplanar Xrays and clinical examination shows the potential applicability of our study in clinical practice. Ideally, we wish to be able to construct individualized musculoskeletal models for any human subjects with minimal effort. Although the potential is apparent, we have not tried to validate or verify its effectiveness in medical applications, which may require rigorous experiments with human subjects [lee2015push]. This ambitious goal remains the subject of future research beyond the scope of this paper.
Our algorithm assumes the presence of a reference model, which can be retargeted only to structurally equivalent bodies. If we want to build a musculoskeletal model of nonanthropomorphic characters, there is no highquality reference model to begin with. An interesting direction for future research is to build anatomical models from scratch, potentially with help of various types of measurements including 3D scans, medical images, EMG, motion capture, and videos [kadlevcek2016reconstructing].