Advanced Computer Graphics, Fall 2020

Final Project: Line-based Musculoskeletal Simulation and Control

Click here for source codes!

Date updated: 5 Nov 2020

Ying Wang

Project Description

Project Proposal

Summary

Description:

The goal of this project is to simulate and control the muscle line system with complex routing constraints in a more efficient and accurate way.

Importance:

Most simulators in animation and biomechanics use algorithms for chains of rigid bodies connected by joints. It is very important to account for muscle mass in those models. A widely used method is "Segment Lumping" which combines muscles and bones in a segment. This method is easy but it also introduces large error during dynamic simulation. The proposed project will generate motions which are more natural by accounting for the underlying muscles and bones.

Previous works:

Biomechanical Simulation and Control of Hands and Tendinous Systems
Large-Scale Dynamic Simulation of Highly Constrained Strands
Musculotendon Simulation for Hand Animation:
This paper proposed a new method for simulating thin strands such as tendons and muscles with complex routing constraints. However, the rigid bones are simulated using maximal coordinate and they need a lot of constraints (joint constraints, surface constraints, fixed constraints etc). We want to further reduce the system using new Jacobian matrices.
Muscle Path Wrapping on Arbitrary Surfaces:
This paper can generate the massless and elastic strand for arbitrary wrapping surfaces. However, in our case, we want the muscles to have mass, and the mass of muscle lines should be comparable to the mass of bones.

What I am planning to do:

I plan to simulate the cases where we approximate the muscle geometry as a polyline (via-point model) and also different obstacles such as spheres, cylinders, and eclipsoids (wrapping surface model) in the reduced system, by deriving different Jacobian matrices based on the constrained quasistatic condition.

List of Goals

Intermediate goals:

I plan to finish muscle line simulation with different wrapping surfaces (obstacles), with the correct energy plot by the update point. During this stage, everything is implemented in MATLAB for easy testing.

Final goals:

For the final goal, I want to port everything from MATLAB to C++ and generate some interesting demos which can showcase the advantage of this method.

Additional goals:

If I have more time, I will implement this with a higher order implicit integration such as Backward differentiation formula BDF1 and BDF2 to improve the accuracy.

Project Update

Summary of Work to Date

I have finished the via-point muscle model integrated with the REDMAX articulated rigid body framework. And I also tested all the Jacobian matrices in our method. Now both the simulation and energy plot are correct. I have also finished the muscle model with stationary wrapping surfaces and did all the testing with that model. However, I have some energy drifting problem with the moving wrapping surfaces. Therefore, I am currently working on solving those energy issues.



Analysis of Work

My original plan for the intermediate goal is to finish all the wrapping surfaces with the correct energy plot. Up until now, the goal is partially met. But I still need to find out what is causing the energy drifting problem in my code (or math). My hunch is that I missed some terms in the equations or my Jacobian time derivative is not correct.

Plan for Completion

I think I am slightly behind my schedule because I spent more time than I expected in figuring out why the energy is drifting in the via point model. Before I assumed that the energy was correct in the via point model, but the problem became worse when I moved on to the next part of the project. It took me sometime to go back and check it again using the simplest test case that I can possibly use. The good thing is, after all the checking and testing, now I have created a lot of test functions for my code and did more writeups in different methods. Therefore, I think I can finish my goals during the rest of the time.

Final Report

Problem Summary

The goal of this project is to simulate and control the muscle line system with complex routing constraints in a more efficient and accurate way.

Most simulators in animation and biomechanics use algorithms for chains of rigid bodies connected by joints. It is very important to account for muscle mass in those models. A widely used method is "Segment Lumping" which combines muscles and bones in a segment. This method is easy but it also introduces large error during dynamic simulation. The proposed project will generate motions which are more natural by accounting for the underlying muscles and bones.

Previous Work

  • Musculotendon Simulation for Hand Animation
    This paper proposed a new method of generating the motion of tendons and muscles under the skin of an animated character.
  • Large-Scale Dynamic Simulation of Highly Constrained Strands
    This paper proposed a new method for simulating thin strands such as tendons and muscles with complex routing constraints. However, the rigid bones are simulated using maxmal coordinate and they need a lot of constraints (joint constraints, surface constraints, fixed constraints etc). We want to further reduce the system using new Jacobian matrices. Also, the strands in this paper can only be inextensible.
  • Biomechanical Simulation and Control of Hands and Tendinous Systems
    This paper proposed a new method for eliminating the Eulerian DOFs by introducing a new quasistatics Jacobian matrix for via point model. The muscle can be flexible in this work. However, there are some underlying energy problems due to the lack of Jacobian time derivatives and quadratic forces.
Description of Work

Assume that we have the full equations of motion of the system, including all degrees of freedom (DoFs) from all the rigid bodies and strand nodes that need to be solved. The goal is to simplify this full problem into a reduced problem which only has DoFs of the rigid bodies, while taking into account the motion (intertia) of the strands. We will achieve this by forming different Jacobian terms.

Static Condensation

We used static condensation to eliminate the internal DoFs of the muscle by assuming these DoFs are in static equilibrium. This method will involve the computation of stiffness matrix, displacements and forces such as gravity. We will partition the DoFs of the muscle into external DoFs which are the end points of the muscle, and the internal DoFs which are the nodes along the muscle path (straight or curved). By using block Gaussian elimination, we can arrive at the Jacobian which maps the DoFs of the internal nodes to those of the external nodes.

Constrained Quasistatics

We can incorporate velocity-level constraints between the internal DoFs of the muscle and rigid body DoFs. For example, we can constrain an internal muscle node to lie on a line that is defined with respect to a rigid body. Or we can constrain an internal node to lie on a wrapping surface that is attached to a rigid body. This is mainly done by forming a linear KKT system, rearrange the row and columns to obtain the final forms of the Jacobian which expresses the internal muscle node velocity as an affine function of the joint velocity.

Obstacles I Met

It is really helpful to examine whether my simulation is correct or not by looking at the energy plot of the system. However, there are so many things that can make the energy look bad. It is getting more and more difficult for me to debug once the system gets bigger. For example, I did double check all the Jacobians and Jacobian time derivatives using finite differencing method. And I also implemented some toy example such as a single wrapping surface with four nodes (the first one is fixed, and the second and the third nodes are tangential to surface, and the forth node is moving in a circle). In this toy example, the only DOF is the rotational velocity of the last node. The purpose of this test is that I want to see what's going to happen if I use some alternative methods. I can compute all the inertia matrices and quadratic forces by intergrating along the straight line and the arc line. Also I can compute the Jacobian that maps the only DOF to the full DOFs of the system by implicit function theorem. However, even with the toy example, I found the same energy problem. The bright side is that instead of finding the bug in the bulky code I had before, now I had a much simpler example to debug. Now I feel more promising in solving this lingering issue.

Results

This is a toy example of the via point model.


This is a toy example of the wrapping surface model.


I did a simple hand simulation using my current work. I used five chains of rigid cuboids to represent the hand bone structures and some cylinders attached to the bones to represent the wrapping surfaces near the joints. I approximate the sizes of the bones and joint surfaces by measuring my own hands. The bones and muscles have equal mass. Although I have some energy issues, visually the simulation doesn't have obvious artifacts or problems. The goal of this example is to show that we can combine the two methods we mentioned above, and simulate a much complicated human muscle model only using the DOFs of the joint configuration. The total DOFs of the system are much less than the DOFs in the previous works (maximal coordinates and lots of constraints).

Analysis
New Results
I think my project has some new contributions in this subject, for example, the muscle models are incorporated in the joint space, and we include some new Jacobian terms and some quadratics velocity vectors terms which were neglected in the previous works. Also, if the energy problem is solved in future, we can obtain the correct motions of the muscle wrapping model, which will be a new feature of our work.

Meeting Goals

I think I made some progress in my research project, such as solving the energy problem in the via point model, and solving the position problem of some collision points in the wrapping surface model. But I didn't finish all the original goals in my plan, such as solving the energy problem in the wrapping surface model. And another thing is that I didn't port my current MATLAB code into C++ OpenGL program. Instead, I did the simulation in MATLAB which looks much coarser than it will be in a C++ program that probably would have a much better bone models, lighting, rendering, etc. There are still a lot of things that I need to do in the future.

Future Work

I definitely think there are more future works to be done in the short and long term. My short term goal is to find out what is causing the energy problem in the muscle wrapping surface model. In the long term, I need to add the muscle forces with activation level into the system. Another task is to use a higher order implicit integration method, for example, Backward differentiation formula BDF1 and BDF2. Those methods can increase the accuracy of the simulation, especially for the solution of stiff differential equations. Finally, I think incorporating into the OpenSim framework would be a interesting future work.

Reference

  • David Baraff and Andrew Witkin. 1998. Large Steps in Cloth Simulation. In Proc. SIGGRAPH 98, Annual Conference Series. 43–54.
  • Scott L Delp, Frank C Anderson, Allison S Arnold, Peter Loan, Ayman Habib, Chand T John, Eran Guendelman, and Darryl G Thelen. 2007. OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE transactions on biomedical engineering 54, 11 (2007), 1940–1950.
  • Shinjiro Sueda, Andrew Kaufman, and Dinesh K. Pai. 2008. Musculotendon Simulation for Hand Animation. ACM Transactions on Graphics 27, 3 (Aug 2008), 83:1–83:8.
  • Shinjiro Sueda, Garrett L. Jones, David I.W. Levin, and Dinesh K. Pai. 2011. Large-Scale Dynamic Simulation of Highly Constrained Strands. ACM Transactions on Graphics 30, 4 (Aug 2011), 39:1–39:9.
  • Prashant Sachdeva, Shinjiro Sueda, Susanne Bradley, Mikhail Fain, and Dinesh K. Pai. 2015. Biomechanical Simulation and Control of Hands and Tendinous Systems. ACM Transactions on Graphics 34, 4 (Jul 2015), 42:1–42:10.