Experimenting with IK solver

Happy new year everyone!

I have spent most of my free time between Christmas and New Year by moving the preliminary IK solver code from animation merge to generic pose editor. The results are mostly working, so expect a release sometimes soon.

Inverse kinematics is a set of algorithms everyone is talking about, but there is actually very little useful information found in web. Other than studying the source code of existing software, that is. There are actually two complexities - first the algorithm itself and then specifying meaningful constraints. Plus all kinds of small technical things.

The algorithm is an easier part. Just went with Downhill Simplex, mainly because it is conceptually easy and reasonably fast. The function is simply the distance between IK target and the predefined endpoint of moving limb.

Constraints are trickier. We have to add the properly scaled tensions of joints to distance function, so joints try to find the easiest (most relaxed) position for given target. As much as I hate Euler angles, I cannot imagine quick and easy way to specify joint tension using only quaternions and/or matrices. It should be possible to create 3-dimensional grid of values and interpolate normalized quaternion in that space. But I do not want to experiment with it before I can get everything working simpler way.
So what had to be done was to decompose all bone relative movements (compared to rest position) into Euler angles. The order of axes is configurable (XYZ, XZY and so on).

Once we have Euler angles, it is quite logical to specify the maximum positive and negative rotations allowed to given bone. But there is one caveat - minimization algorithms normally require the tension function to be smooth. So we cannot simply clip the tension function at threshold. Also it is normally not a good idea to keep the tension value constant below threshold, because joints always prefer to be as close to rest position as possible.
At moment I made a mix of two power functions for tension calculation:

T = R * (angle / P) ** 1.25 + (angle / P) ** 16

T - calculated tension value          
R - rigidity (0-1) of the movement    
P - threshold value (maximum rotation)

There are two different set of threshold and rigidity values - one for positive, one for negative values. Rigidity determines, how fast the tension increases below the threshold value (it has small exponent 1.25 to make the overall function differentiable at zero).

The example of tension function

Here is a graph of example tension function with parameters -1.0;1.0 (negative) and 1.5;0.1 (positive). It describes a joint, that moves backwards up to about one radian with considerable tension and forward up to 1.5 radians relatively freely. Although it does not seem so looking at graph, the function is actually smooth at zero.
Certain joints have some rotations blocked - like twisting and sideways bending for knees. This is better done by turning those axes off completely instead of setting very small movement thresholds, because we want to keep the degrees of freedom as low as possible.

And it works - at least wee enough for static positioning. For some reason the simplex solver diverges a lot for small changes in target position if the overall movement is relatively big.

Positioning of limbs with IK targets

The latest version of Khayyam has a set of meaningful IK chains and rotation constraints built into XPP importer (only for SBZ and compatible skeletons).