Monday, August 29, 2011

GLSL water reflections

I have been mostly working on sehle (Khayyam rendering engine) based game framework lately. There have been some new things in Khayyam - that will be the scene/level builder for it - as well. The latest addition is reflective water material.

Esk standing on near water - with ripples and reflection
You can convert all static objects (OBJ, 3DS and Collada) to water bodies - although to get realistic results, they should be bowl- or plane shaped with surface at Z zero (in object coordinates). Everything is fake, of course - there are no real waves or refraction, but only distorted reflections and transparency. But it still looks reasonably good in my opinion.

At moment there are following water properties implemented:
  • Full reflection (everything that can be rendered can be reflected, except other reflective surfaces)
  • Viewing angle based reflectivity - i.e. if looking directly down, water is mostly transparent, looking along water surface it becomes almost perfect mirror
  • Depth-based transparency - shallow water is almost transparent, deeper water becomes opaque
  • Up to 8 simultaneous ripple generators, that can be assigned to "shore" vertices of water body
At the image below you can see how the water color changes from almost fully transparent near bank to completely opaque brown at the middle of river. As the viewing angle is high, overall reflectivity is low and the colors of river bottom and water can be seen. Ripples change the apparent angle between viewer and water surface, so although sky is uniform color, ripples still show as they change the reflectivity of surface.

A river from above showing water and bottom colors
 At the next image the viewing angle is low and thus the surface of water is almost completely reflective.
 
A river from low angle, showing almost perfect mirror
Ripples are circular and their wavelengths, amplitudes and generation points are updated randomly over time. Thus the actual pattern almost never repeats - although sometimes it is not as nice as on the image above.


When I find some free time, I'll write the technical details about the actual mathematics in shaders. And hopefully there will be new Khayyam release too.

Have fun!

Monday, July 18, 2011

Encoding normals for Gbuffer

Deferred shading has nowadays become nearly ubiquitous in real-time rendering due to the separation of geometry and lighting passes. The general idea of this technique is to render all important properties (color, normal, reflectivity and so on) of all visible surfaces to a set of screen-sized textures. Later rendering passes (lighting) get surface parameters from those textures and thus complex geometries do not have to be rendered multiple times.

One design problem of Gbuffer implementation is the efficient encoding of surface normals. Normals in three-dimensional space have three components:

N = Nx,Ny,Nz

The most intuitive way would be to store those in three texture channels. But this would be wasteful, because we know, that:

Nx2 + Ny2 + Nz2 = 1

Gbuffer requires many texture channels for other data, thus we want to cut off as many unnecessary channels as possible to conserve memory and fillrate.

Another simple alternative would be to store only two components (Nx and Ny) and reconstruct third (Nz) in pixel shader. Unfortunately we will lose the sign of Nz in that process. While in ideal scene there should be no normals with negative Z (facing away from camera), they happen in real scenes because of realistic polygon counts.

There exist several techniques for encoding normals while (at least partially) preserving Z component. Unfortunately no one of those is ideal. Either they are not able to encode full range of Z (various projections) or waste precision and calculation (polar coordinates).

While implementing foliage shader for Khayyam I needed a way to encode full range of possible normals into two-component vector (that will later be stored in 4-component 8-bit RGBA texture of Gbuffer). The following is the brief introduction of the technique used. It solves all the problems I am aware of - with the price of slightly more shader instructions than simple projection encoders.

Efficient projection of normal space

We forget for awhile the need of storing the sign of Z component and look, how we can store X and Y components most efficently.
The trivial way is to simply forget Z and store X and Y. In geometric terms it means projecting the hemisphere of normal space to an unit circle on XY plane along Z vector as in following figure:

The trivial projection of normal space to XY plane

The problem with this projection is, that the uniform angular difference between normals near equator gives much smaller difference in X and Y coordinates than near poles. If the precision of normal texture is limited to 16 bits (4x8 bit or 2x16 bit integer buffers) or even less (FP16 buffers) , the loss of precision may become visible. For example, for 16 bit buffers, the minimum representable angle near equator is:

acos (1 - 2e-15) = 1.44°

Near poles it is much better:

asin (2e-15) = 0.002°

Clearly we are wasting precision near poles and introduce errors near equator.

To get better distribution of precision, we can switch projection from parallel to radial, with the projection center at the point (0,0,-1) as on the following figure:

The radial projection of normal space

The distribution of angles is still not uniform but is clearly much better (it differs only by the factor of 2 at poles and equator).
Actually, having the center of projection at (0,0,-sqrt(2)) would result in even better distribution, but I wanted to keep it at (0,0,-1) to simplify calculations a bit.

The hemisphere is still projected at unit circle at XY plane, although now we cannot project the normals on negative hemisphere directly but have to take the absolute value of Z coordinate before projecting.

This solves the problem of effecient use of normal texture precision. Now we can go on to encoding the sign of Z into the same buffer.

Encoding the sign of normal Z coordinate

Our positive hemisphere of normal space got projected to unit circle on XY plane as in following figure:


The normal space projection to XY plane and compaction schema


As we can see, we are actually wasting the numerical space of normal texture. The texture is able to encode values in unit square, but only unit circle is used - the possible values on dashed areas can never happen in our projected normals. We could store normals with negative Z in those areas, but unfortunately there is no easy algorithm to do that.

Now we do another conversion and pack all the values from unit circle into diamond-shaped area, given by the following formula:

|X| + |Y| <= 1

This can be done with the following algorithm (let the normal be Nx', Ny' and the projection be Nx'', Ny''):

r = sqrt (Nx' * Nx' + Ny' * Ny')
D = |Nx'| + |Ny'|
Nx" = Nx' * r / D
Ny" = Ny' * r / D

We lose some precision at  the normals close to 45° angles, but the maximum loss is by the factor of sqrt(2).

The resulting projection of normal space is shown in following figure:

Compacted normal space and Z flipping directions

Now the whole projection of the positive hemisphere of normal space takes exactly the half of the size of unit square. The extra space can be used to encode normals with negative Z.

We do this by mirroring the final normal projections by the closest one of diagonal lines (|x| + |y| = 1).

Thus we have encoded the full normal space to unit square with minimal loss of precision.

Decoding normals

To decode normals we have first to find the sign of Z coordinate. This is done with the following formula:

Zsign =  1 if (|x| + |y|) <= 1
Zsign = -1 if (|x| + |y|) >  1

Now, if Z is negative, we have to mirror our encoded normal by the closest  diagonal line (|x| + |y| = 1) so the resulting normal is inside the diamond shaped area.

We unpack normals to unit circle with reverse packing algorithm:


d = sqrt (Nx" * Nx" + Ny" * Ny")
r = |Nx"| + |Ny"|
Nx' = Nx" * d / r
Ny' = Ny" * d / r

From normal projection on unit sphere (N') unit sphere we unproject to normal on positive hemisphere (N).

Finally we set the sign of the Z component of N from previously found Zsign.

And that is it.

Encoding quality

I ran simulated tests of this algorithm on CPU with 100 000 000 randomly generated normals, including normals with any 1 or 2 components zero. Encoded vector was written and then read back from 4-component unsigned char buffer to simulate storing of 2-component normal in RGBA8 texture.

The maximum angular error between original and retrieved normal was 0.028°.

GLSL shader code

For those interested, here is current GLSL shader code of libsehle implementation. It is a bit heavy, but if your lighting is fill-rate bound, as is mine, you can spare some GPU cycles.

Encoding.

vec2 encodeNormal (vec3 normal)
{
    // Project normal positive hemisphere to unit circle
    // We project from point (0,0,-1) to the plane [0,(0,0,-1)]
    // den = dot (l.d, p.n)
    // t = -(dot (p.n, l.p) + p.d) / den
    vec2 p = normal.xy / (abs (normal.z) + 1.0);

    // Convert unit circle to square
    // We add epsilon to avoid division by zero
    float d = abs (p.x) + abs (p.y) + EPSILON;
    float r = length (p);
    vec2 q = p * r / d;

    // Mirror triangles to outer edge if z is negative
    float z_is_negative = max (-sign (normal.z), 0.0);
    vec2 q_sign = sign (q);
    q_sign = sign (q_sign + vec2 (0.5, 0.5));
    // Reflection
    // qr = q - 2 * n * (dot (q, n) - d) / dot (n, n)
    q -= z_is_negative * (dot (q, q_sign) - 1.0) * q_sign;

    return q;
}

Decoding.

vec3 decodeNormal (vec2 encodedNormal)
{
    vec2 p = encodedNormal;
   
    // Find z sign
    float zsign = sign (1.0 - abs (p.x) - abs (p.y));
    // Map outer triangles to center if encoded z is negative
    float z_is_negative = max (-zsign, 0.0);
    vec2 p_sign = sign (p);
    p_sign = sign (p_sign + vec2 (0.5, 0.5));
    // Reflection
    // qr = q - 2 * n * (dot (q, n) - d) / dot (n, n)
    p -= z_is_negative * (dot (p, p_sign) - 1.0) * p_sign;

    // Convert square to unit circle
    // We add epsilon to avoid division by zero
    float r = abs (p.x) + abs (p.y);
    float d = length (p) + EPSILON;
    vec2 q = p * r / d;

    // Deproject unit circle to sphere
    float den = 2.0 / (dot (q, q) + 1.0);
    vec3 v = vec3(den * q, zsign * (den - 1.0));

    return v;
}

Monday, June 13, 2011

More about Matrix to Euler angles conversion

I previously posted the algorithm to convert rotation matrix to arbitrarily ordered Euler angles.
It solved two problems of the trivial conversion method in Matrix and Quaternion FAQ:
  • Use arbitrarily ordered axes (like ZXY) in addition to XYZ
  • Remove the discontinuity happening at Ax = +/- 90 degrees
Meanwhile I discovered two remaining problems with my implementation:
  • The algorithm always return the minimum set of absolute Euler angles (closest to {0;0;0}, as defined by the sum of squares of the angles) from the two alternatives. Often it is desirable to get the set of angles that are closest not to identity {0;0;0} but to some predefined rotation {Rx;Ry;Rz}.
  • When gimbal lock occured, the X rotation (or corresponding rotation if axes were reordered) was set to zero. This may not be the best solution, if we need the Euler angles to be as close to some specific set of values as possible.
The problem surfaces if we want to describe animations by Euler angles to allow easy adjustments by IK solvers (and have joint constraints defined) and interpolate between keyframes. In that case we absolutely have to keep the angles of adjacent keyframes as similar as possible - because although the alternate sets of angles make identical matrices, interpolations between them generally do not.

Stable conversio of matrix to Euler angles closest to given set

The matrices have row-major order, so they have to be transposed for OpenGL.

Base is a set of three angles {Bx;By;Bz} that will be used as reference while determining the best set of resulting angles.

void
Matrix3x4f::getEulerAngles (f32 angles[], const f32 base[]) const
{
 float aX, aY, aZ;
 float aX2, aY2, aZ2;

 aY = asin (c[2]);
 aY2 = (aY >= 0) ? M_PI_F - aY : -M_PI_F - aY;
 if (fabs (c[2]) < 0.99999f) {
  // Over 0.25 degrees from right angle
  float C = cos (aY);
  float C2 = cos (aY2);
  float trX = c[10] / C;
  float trY = -c[6] / C;
  aX = atan2 (trY, trX);
  aX2 = atan2 (-c[6] / C2, c[10] / C2);
  trX = c[0] / C;
  trY = -c[1] / C;
  aZ = atan2 (trY, trX);
  aZ2 = atan2 (-c[1] / C2, c[0] / C2);
 } else {
  // Less than 0.25 degrees from right angle
  // Gimbal lock occured
  aX = aX2 = base[0];
  // Keep X rotation the same as base
  float Sx = sin (aX);
  float Cx = cos (aX);
  if (fabs (Cx) > fabs (Sx)) {
   float Sz = (c[4] + (Sx/Cx) * c[8]) / (Sx*Sx/Cx + Cx);
   aZ = aZ2 = asin (Sz);
  } else {
   float Sz = (c[8] + (Cx/Sx) * c[4]) / (Cx*Cx/Sx + Sx);
   aZ = aZ2 = asin (Sz);
  }
 }

 float d1X = clamp_pi (aX - base[0]);
 float d1Y = clamp_pi (aY - base[1]);
 float d1Z = clamp_pi (aZ - base[2]);
 float d2X = clamp_pi (aX2 - base[0]);
 float d2Y = clamp_pi (aY2 - base[1]);
 float d2Z = clamp_pi (aZ2 - base[2]);

 if ((d1X * d1X + d1Y * d1Y + d1Z * d1Z) < (d2X * d2X + d2Y * d2Y + d2Z * d2Z)) {
  angles[X] = aX;
  angles[Y] = aY;
  angles[Z] = aZ;
 } else {
  angles[X] = aX2;
  angles[Y] = aY2;
  angles[Z] = aZ2;
 }
}

If aY is close to +/- 90 degrees (gimbal lock occured) we want to set aX to the base value instead of setting it to zero.
Also we calculate the differences of both alternate sets of angles from base and choose the final set based on the squared sum of those differences.
Clamp_pi simply converts angle to -PI...PI range.

Conversion of matrix to arbitrarily ordered Eular angles

The matrices have row-major order, so they have to be transposed for OpenGL.

Order is array of three integers, determining the order of rotations ({0;1;2} is XYZ, {2;0;1} is ZXY).
Base values should be ordered according to order.

void
Matrix3x4f::getEulerAngles (f32 angles[], const int order[], const f32 base[]) const
{
 float sign = ((order[1] == (order[0] + 1)) || (order[2] == (order[1] + 1))) ? 1.0f : -1.0f;
 // Permute matrix to switched axes space
 Elea::Matrix3x4f rs;
 for (int row = 0; row < 3; row++) {
  for (int col = 0; col < 3; col++) {
   rs[4 * row + col] = c[4 * order[row] + order[col]];
  }
 }
 rs[2] *= sign;
 rs[6] *= sign;
 rs[8] *= sign;
 rs[9] *= sign;
 // New base
 float nbase[3];
 for (int i = 0; i < 3; i++) nbase[i] = base[i];
 nbase[2] *= sign;
 // Decompose to Euler angles
 rs.getEulerAngles (angles, nbase);
 angles[2] *= sign;
}

The trick is to permute the rotation matrix in such way, that the XYZ angles of new matrix are properly ordered decomposition angles of the original matrix and then do the trivial decomposition.
Hopefully this time the algorithm is final ;-)

Have fun!

Tuesday, April 5, 2011

The future of Khayyam

Khayyam 3D did not have very clear beginning. Few years ago I was experimenting with 3D rendering and animation under projects codenamed Floyd and SDLMu. The initial impulse was to write an engine and scene editor for my Ant-Attack inspred game.
The unorganized codebase soon grew out of control, so I started from scratch. Or more precisely ported the program structure from my former project Sodipodi. When choosing name for the new project, I discovered Omar Khayyám - a Persian mathematician, who sort of discovered the inherent link between geometry and algebra. And thus the program was born.

As of now, Khayyam is evolving in several directions based on what I am interested in any given moment.
  • Scene builder. My game is still in design stage, but at least I have the main character - Sara.
  • Virtual photography tool. It was partially inspired by HongFire hentai game modding community and partially by Poser and DAZ applications.
  • Animation editor. This came from the need for some decent animations for Sara.
  • Model importer and converter. Implemented step-by step as needed.
  • 3D engine. Mostly for learning, but also because I think I have some ideas, how it can be done right. 
As of now (April 2011) the program has at last become sort of usable for some simple tasks, like converting model formats, building POV-Ray scenes and so on. To go further, there should be a somewhat clearer working plan - which I have tried to lay out for myself.

1. Scene builder

Khayyam itself is not game engine and never will - the two-level DOM structure is too heavy for gameplay. Instead there will eventually be an engine based on libsehle.
There will be custom object formats. At least an octree-based map node, skeletal figures and static objects. Those can be imported from Blender Collada files, optimized and written to custom format.
Objects need JavaScript bindings. I plan to make JavaScript the scripting language for game engine, so the actions should be programmable and testable in Khayyam.

2. Posing application 

The Poser figure loading and manipulating is still very flaky. This will get better eventually, so expect that in some future date Khayyam will be able to render (with the help of POV-Ray) everything Poser and DAZ studio can. Virtual photography is becoming more and more used, and Linux definitely needs a free application for that. The fist milestones are unlit objects (like backgrounds) and smoothscale channel implementation for full body morphs.

3. Animation editor

Animations will be integrated globally into Khayyam document structure, so you can render not only still pictures, but full animation sequences.
I am halfway implementing animation graph and biped walking controller for figures. Eventually you could import animations from other formats (like Biovision BVH), define sequences, merge, split, edit and join them, and build fully controllable animation graph from these. The controller will have IK solvers for fine limb/finger adjustments and will be controllable via JavaScript.

4. Model converter

While several popular models are supported, the actual features vary a lot. There will be more model and animation formats (like MD5 animations) and more features. I also plan to support more game model formats - definitely all future Illusion games, but probably some others too.

5. 3D engine

The big missing parts are proper shadow map algorithms, particle effects and post processing.

Other than that, Kahhaym will have some UI overhauls in future. Tools will have their own dynamic property pages. Some extra tools, like bone and object pickers are needed. And there have to be UI controls for constraint placement. Also I am particularly not happy with Gtk+ default style - which wastes enormous amount of screen space. At least the animation editor needs more compact style

Sara and african elephant (DAZ 3D model)

Have fun!

Monday, April 4, 2011

How to convert animations in Khayyam

Starting from the 20110402 release Khayyam can import, convert and export Biovision BVH animations from one figure/skeleton to another. Following is a mini-tutorial guiding you through the process.
You may take a look to the overview of pose fitting algorithm to better understand which joints will be moved and how.

0. The example problem

Say, you have some good animations in BVH format and want to apply these to Poser characters. Unfortunately, skeleton topologies  and/or bone orientations do not match, so they are not usable as-is.
Khayyam can actually export BVH files from supported animation types (currently only Illusion game formats). How to do this is outside of this tutorial - but as a hint look at skinAnimation and sequence objects in document tree editor.

1. Setting up target

Before we can convert animation to target figure, we need it's skeleton imported into Khayyam. Skeleton comes together with model, so at first we have to import Poser cr2 figure.
It is probably good idea to use reduced resolution figures for this process, because animating the whole iteration 4 (V4, M4) models can be quite slow. I will use Chibibel for current example - she has much lower polycount and can show animations without too much frame-skipping.

Chibibel imported to Khayyam

At moment Khayyam does not set up joint constraints automatically for cr2 objects, so you have to define at least elbow and knee constraints. Otherwise you risk, that legs or arms may suddenly bend sideways during animation.
Select imported figure, switch to pose mode and select left shinbone.

Chibibel skeleton with shinbone selected

The colors of rotation arrows mark the bone-relative axes of rotations (red - X, green - Y, blue Z). Normally knee bends only in one direction - in current case this is green rotation, i.e. Y axis. So you have to disable X and Z rotations.

Open document tree dialog. It will automatically show property page for selected bone - i.e. left shin.

The property page of left shinbone

In constraints frame select first X, then Z channel and uncheck the Enabled box for both. This turns off the rotations around those axes.
On the same page you can set more fine-grained behavior by defining tension parameters for bones. Angles are basically maximal allowed rotations in either direction, rigidity values control, how "difficult" is given rotation. For example, knee joint moves very easily (in allowed range), but twisting shoulder is difficult, so character tries to minimize it. But if the pose is not achievable otherwise, the twisting still can happen - unlike the hard cutoff angles, that cannot be exceeded.
In current example we fortunately do not need to set up tensions.
Repeat the disabling X and Z channel rotations for right shinbone. Then select left forearm and disable X and Y rotations. The same with right forearm.
Now we have completed the minimum setup needed for successful pose merge.

2. Converting animation

In document tree select the main figure node. From the property page click on Merge animations... to open the converter dialog.
The first step is to load source animation or pose. Click Load animation... and select your BVH file.

Running animation applied to Chibibel skeleton

Once the animation is loaded, Khayyam tries to adjust the destination figure to starting pose. This is iterative process and may take more than one try. So if the pose is not close to the original, you may try to click Adjust pose few times.
The red armature is the original animation skeleton. The blue one is the skeleton of target figure. These can be toggled on and off to see clearer picture.
The other controls you may tweak are:
  • Scale - changes the size of source skeleton. This should be used (together with ZScale) to make the skeletons approximately the same height - otherwise shoulders and hips can not be posed correctly.
  • ZScale - changes the height of source skeleton. Thus Scale can be used to make the widths of hips and shoulders approximately right, and ZScale then to make the height of character right.
  • Guess scale - Khayyam tries to determine scale from the distances between hips and shoulders. This may or may not work correctly.
  • EMax - tha maximum error tolerance - if the combined distance between poseable chains (body, limbs...) is less than EMax, iterative pose refinement will stop. Smaller value may give better adjusting but makes the process slower.
  • TScale - the scaling factor of bone tensions. Usually this can be 1 for default tensions or some smaller value (0.01) if all tensions have been set up in skeleton.
  • Runs - number of distinct runs of IK solver
  • Iterations - the maximum number of tries in each run
Playing with these may give better results. If animation is jerky, increasing the number of runs often helps. If limbs do not track the originals well, lowering TScale may help. Experiment!
Also, you should keep in mind, that while hips and shoulders are placed using the absolute positions, other joints are adjusted by directions. Previous blog entry describes the process in more detail.
Khayyam knows the bone names of some common skeletons and can determine the corresponding standard parts automatically. But sometimes it may not know, which standard humanoid bone certain name corresponds to. In that case you have to define bone names manually, using the bone mapper tool. This can be done by clicking Map buttons in source and destination skeleton frames.
Not all skeletal bones are 1:1 mappable and this is not even needed. What is important is, that standard joints of movable body parts are mapped correctly.

Bone mapper dialog for Poser skeleton

In the example above, you can see the bone mapper dialog for standard Poser biped skeleton (V4, M4, Chibibel). In the left-hand tree there are all bone names that appear in skeleton file, at right are the corresponding Khayyam standardized body parts. Once set up the mapping can be saved and loaded (actually the last saved mapping will be loaded automatically).
Once the bone mappings are set up (and saved), click Close and the conversion solver will update itself.
You can play the animation with current conversion parameters using play button.

3. Importing animation to document

Once you are satisfied with the conversion results, click Import. Khayyam will then generate a tree of animation objects as child nodes of parent Figure node.
Hint: You can set NRuns to bigger value and EMax to lower value for import so the actual conversion will be done with maximal accuracy.
The object of interest for us is a sequence node, that will be created as the last child of figure. It defines the imported animation sequence and can be played from play button on its property page.

The imported animation sequence object

If the imported sequence works as intended, you can export it to new BVH file. This new file uses already the proper skeleton of the target figure.


 And that's it. I have many ideas how to make this process smoother, so stay tuned!

Sunday, April 3, 2011

Snapshot release 20110402

A new snapshot release 20110310 is available from SourceForge page as source tarball and precompiled Win32 binary.

The most important new feature is refactored biped pose/animation converter:
  • Both animations (from BVH files) and static poses (from other objects) can be converted
  • Animations can be saved to sequences (and exported to BVH again)
  • Almost all body parts of biped are adjusted (body, limbs, head, fingers)
  • IK solver parameters can be adjusted
  • Bone definitions can be set by user (and saved for future use)
Some problems remain - collar bones do not move, you cannot set T-pose and so on.
The better the constraint system of the destination figure, the better pose matching will be. At minimum you should set knee and elbow constraints to avoid wrong bendings.


A demonstration video of making Sara-chan to run. Red armature is the original BVH animation file, blue is the Sara skeleton. It is not ideal, because the BVH does not have finger poses.

Have fun!

Saturday, March 26, 2011

Converting animations between skeletons

I have written one more piece about the skeleton merge about half year ago - you can read it here. What follows is the a much better version of the original algoritm.

The problem

Say you have some good animations (there are many free BVH sites in web) and want to use these for your own models. Unfortunately - unless you made your models using exactly the same skeleton layout, the bones do not match and copying animations is not trivial problem.

There are many ways, how skeletons may differ:
  • Bone names do not match. This is easy to overcome - either build mapping or rename bones.
  • Rest pose orientations use different systems (like XYZ euler angles vs. ZXY euler angles vs. quaternions). This is also minor problem, as normally you want to convert everything to quaternions or matrices anyways.
  • Bone lengths do not match. Thus copying animation directly can make movement unnatural. 
  • Rest pose bone positions (directions) do not match. Thus animations can not be immediately copied between skeletons, but have to be transformed to absolute space (object space) before.
  • Skeleton topologies do not match. This is often a showstopper, because if you do not even have similar bones, there is no way to copy animations. 
Example of different common skeleton topologies

The topology mismatch is the hardest problem to overcome. In the example above, if you, say have animation for the right skeleton, the 4 parts of backbone should move the 2 backbone parts and pelvis of the left skeleton. While in theory it should be possible to solve this problem analytically, it is not easy - and has to be done by hand for each combination of skeletons.
Instead of that, using IK solvers gives us easier way to map those skeletons. It is still not universal - the mapping algorithms have to be adjusted to object itself - in current case a humanoid. But normally this is not big problem, because in 99% cases humanoid is precisely what you want to animate.

In any case you need good IK solver, that can do at least 6-7 bones long chains, adjust both positions and rotations and use both absolute positions and bone directions for targets. Also, the solver has to support comparing up to 4 target points or directions simultaneously. I am currently using simplex method in Khayyam (the code is in helpers/ik.cpp), but there are other options as well.
Also the skeleton should have bone constraints defined - otherwise you may discover that elbows or knees will be bent to completely wrong direction.

We start implementing the algorithm by a small observation. Namely - while there are different ways to construct skeletons, the joints have to be still placed to precise spots - otherwise the movement of joints will not be natural. Thus we can use the joints as anchor spots for the algorithm.

The body

The first and hardest part is moving the main bodies of skeletons to the same position. First - the root bones may be arranged differently. Second - root bones may be joined with body either below or above pelvis. Third - there may be different number of backbone elements. Fourth - there may or may not be extra bones, like hips or collar bones.
To solve this problem, we can do the following:
  • Set 4 IK target points to the corner joints of body - both hips and both shoulders.
  • Scale the source skeleton so, that the distances between hips and shoulders are approximately the same, as in target skeleton.
  • Find the topmost common parents of both hips and shoulders Ah and As.
  • Find the topmost common parent of Ah and As B.
  • Build rotating IK chain of bones from  B to Ah and from B to As, not including B.
  • Add both the rotations and translations of the rootmost bone R to chain. 
The schematics of body IK solver (green: translate + rotate, red: rotate, blue: not moving)
After solving this IK chain, the corner joins of body will be approximately aligned to the original. Normally you do not achieve perfect alignment, because the width of shoulders and hips is different in skeletons - but this is not big problem, because the direction of hip-hip and shoulder-shoulder lines should be aligned.
It is important to do this alignment in one step - otherwise the alignment of hips may involve extensive rotations and translations of root bone R - and you cannot later align shoulders at all.

Limbs 

To align limbs, there are two strategies - depending on what we want to achieve.
  • We can align the endpoints - if we want to model some specific action - like taking some item from table. In thus case the hand will be in correct place - but because the lengths of arms and forearm probably differ, elbow will be off.
  • Or we can align the directions of bones - for example if modeling hand movement during walking - in that case both elbow and wrist will be off, but the overall look of hand will be more similar.
At the first case it is just ordinary single-target IK solver. Or if you want to keep the elbow/knee direction more similar to original, there should be additional directional target at midpoint.
At the second case, the IK chain has to be solved with two directional targets (shoulder-elbow + elbow-wrist or hip-knee + knee-ankle) simultaneously. Otherwise the topmost bone may turn it to such direction, that 1 DOF joint at midpoint (elbow or knee) cannot move other bone to the right direction.

Palm

There are probably different strategies for adjusting palms. In khayyam I did the following:
  • Set 4 directional targets from wrist to the topmost joints of fingers
  • Add the 3 DOF rotations of palm to IK chain
  • If forearm includes extra rotational bone (like at the skeleton at left), include the axial rotation of it to chain 
By solving this chain, the positions of finger joints will be approximately aligned to the original.

Fingers

Fingers can be aligned in the same way as limbs - but here we need three targets instead of two. Also, as the fingertips are normally the end sites, the last target does not coincide with any joint, but has to be placed about 2-3 cm from the last joint along the bone direction.

Head

Head should be placed by adjusting three positional targets - eyes and skulltop. This can become tricky, because many skeletons may not include eye positions at all.

Sara chan in skeleton merge editor (blue: figure skeleton, red: animation skeleton)
At the above picture you can see the results of adjustment of Sara-chan skeleton to Illusion Characolle animated skeleton. Although the skeleton topologies and bone lengths differ, the pose is pretty close to the original.

At moment the skeleton merge works only in CVS version of Khayyam, but expect it in the next release.