~/imallett (Ian Mallett)

CS 6665.001 "Character Animation" Projects, Fall 2014

by Ian Mallett (ian AT geometrian.com)

This is my project 1 for the course CS 6665.001, being taught at the University of Utah in Fall 2014 by Adam Bargteil. Project 1 is to implement an inverse kinematics solver.


The inverse kinematics algorithm essentially casts the degrees of freedom for each joint as being "directions" one can go toward a solution. Formally, each degree of freedom works out to be an element of a Jacobian matrix describing how much some feature on the linkage moves toward the goal.

A simple example could be an arm with one hinge joint (so it can swing around in a circle) and you're trying to move one endpoint to a point on that circle. The Jacobian in this case describes how much moving the joint \(\Delta\theta_0\) radians moves you toward/away the point in each of \(\vec{i}\), \(\vec{j}\), and \(\vec{k}\).

The Jacobian forms a system that looks like this:\[ \mathbf{J} \cdot \Delta\vec{\Theta} = \Delta\vec{x} \]In the above, \(\mathbf{J}\) is the Jacobian, \(\Delta\vec{\Theta}=\{\Delta\theta_0,\Delta\theta_1,\cdots,\Delta\theta_{n-1}\}\) is a list of changes to parameters of the joints (usually angles, but possibly, for example, scalar translations), and \(\Delta\vec{x}\) is the difference between the feature's current value (for example, its current position) and the feature's goal value.

The system is solved by using SVD to compute the pseudoinverse, thus obtaining \(\Delta\vec{\Theta}\). These changes are made to all the joints and the process is repeated. Hopefully, it converges to a configuration where the linkage satisfies the features.

If you want to keep some particular parameter \(\theta_i\) as close to some value \(\hat{\theta}_i\) as possible (relative to the null-space) while still solving the system, you can instead use the difference \(\Delta\theta_i+(\mathbf{J}^T\mathbf{J}-I)\hat{\theta}_i\) instead of \(\Delta\theta_i\).


My inverse kinematics solver implements arbitrary rigid bodies with up to two connections. Bodies have an arbitrary center of mass (which for my solver is, unfortunately, meaningless) and can have arbitrary joints attached at any location. To make things simpler though, my examples below keep joint locations axis-aligned. All joints have degrees of freedom of rotation around their local axes as well as a sliding translational value. However, all joints have at least one degree of freedom disabled. Disabled degrees of freedom are excluded from the calculations above. To enforce joint constraints, the values are clamped after each step in the iteration. The algorithm proceeds until a certain tolerance is reached.

I had a lot of trouble with the SVD implementation. Everywhere I looked, I saw claims that the algorithm was "well-known". However, if you Google for it, you'll find about a dozen papers talking about ludicrously intricate techniques. One paper I found spent eleven pages discussing previous work. If you look for code, you'll find a few libraries, like LAPACK, GLS, and Eigen, and not much else. I did find a small Python script ported directly from ALGOL which gets kindof the right answers, sortof, for a subset of the problem. My math degree is in abstract math, not applied. Ask me about fundamental group of the octonions, analytic continuations, operator-valued Banach algebras, sheesh whatever. But ask me about factoring a matrix? No. This became one of the very rare occasions when I defer implementation to the experts; I gave up and used Eigen.

The inverse kinematics problem can not have a solution in particular cases (for example if the linkage tries to overextend itself), but sometimes the iteration can get caught in local minima, in the sense that a solution is possible but the algorithm is trying to converge to it in a manner that is "through" a constraint. My solver implements a basic global algorithm (which can be disabled): if convergence is not reached after a given number of iterations, the joint parameters are randomized and the solver tries again. This repeats indefinitely. In practice this helps resolve most of these cases.

My implementation is quite configurable. In the examples below, \(10\) iterations are done per-frame until the solution reaches a tolerance of a Euclidean, world space distance of \(\frac{1}{1000}\)th unit. If the global solver is enabled, the global randomization kicks in after \(100\) iterations without a solution being reached.


In the following examples, the images are the expected result when the provided (below) binaries are run. The images have been resized here to better fit on the screen; view them individually to get the full effect.

The white points are the centers of mass of the rigid bodies (cyan). Rigid bodies are connected by joints (red points) and are constrained (lying within green portions of red gimbals or else inside red cones). The effector (yellow) tries to reach the goal (orange). The controls are left mouse to drag the goal around (hold shift to raise it up), right click and drag and scrollwheel to move the camera.

I apologize that the binaries are Windows-only. Mac makes it literally impossible to write GUI programs in C++ without horrific compiling cruft, and Linux platforms vary too much. Requires GL 4.0 or higher and support for 32x multisampling FBOs. No installation or DLL horror required. These are standalone, statically linked applications. The only known problem is that the orange point may slowly drift when holding shift. I know why this happens, and I don't care about fixing it.

Example 1: Simple Arm
In this (fundamentally 2D) example, the arm has two unconstrained hinge joints. The global solver is not enabled; since the target point can move off the plane that the arm can reach, convergence would never be reached and the arm would start all over, causing the arm to jump. As it is, the arm will continuously keep trying, although it will probably never reach the required tolerance.

Example 2: Translational Arm
In this example, the arm has a single constrained translational degree of freedom and a ball joint to maneuver it. The global solver is not enabled since it is unnecessary for this configuration.

Example 3: Complex Arm
This example tries to mimic a human arm as much as possible down through the wrist. I measured approximate degrees of freedom on my own right arm and used them for the simulation. The global solver is enabled. In the video below, you can see this in several places where the arm will reconfigure itself and try again.



Ian Mallett - Contact -
- 2021 - Creative Commons License