How to get quaternion from rotation matrix

For now probably the whole generation of programmers has learned the basics of projective geometry from the famous Matrix and Quaternion FAQ. I myself have consulted it numerous times while writing Khayyam (or more precisely Elea support library). But experimenting with the 3D transformations I found that while theoretically correct, the two most esoteric methods of the FAQ give suboptimal results. Namely - converting the rotation matrix to quaternion or Euler angles becomes unstable for certain rotations and - at least for Euler angles - does not give the most intuitive solution. So decided to write slightly enhanced variants.

The original code from FAQ:

T = 1 + mat[0] + mat[5] + mat[10]
if ( T > 0.00000001 ) {
    S = sqrt(T) * 2;
    X = ( mat[9] - mat[6] ) / S;
    Y = ( mat[2] - mat[8] ) / S;
    Z = ( mat[4] - mat[1] ) / S;
    W = 0.25 * S;
} if ( mat[0] > mat[5] && mat[0] > mat[10] ) {
    S = sqrt( 1.0 + mat[0] - mat[5] - mat[10] ) * 2;
    X = 0.25 * S;
    Y = (mat[4] + mat[1] ) / S;
    Z = (mat[2] + mat[8] ) / S;
    W = (mat[9] - mat[6] ) / S;
} else if ( mat[5] > mat[10] ) {
    S = sqrt( 1.0 + mat[5] - mat[0] - mat[10] ) * 2;
    X = (mat[4] + mat[1] ) / S;
    Y = 0.25 * S;
    Z = (mat[9] + mat[6] ) / S;
    W = (mat[2] - mat[8] ) / S;
} else {
    S = sqrt( 1.0 + mat[10] - mat[0] - mat[5] ) * 2;
    X = (mat[2] + mat[8] ) / S;
    Y = (mat[9] + mat[6] ) / S;
    Z = 0.25 * S;
    W = (mat[4] - mat[1] ) / S;
}

The problem lies in how the calculation method is switched, depending on the comparison of certain matrix elements. If, for two sequential rotations, the matrix values happen to be only slightly different, but fall into different calculation paths, the resulting quaternions will be very different and thus interpolating between them results in sudden twists. This should not happen, if your matrices are perfectly orthogonal - but if they come into existence as the result of many multiplications and (especially) inversions chances are good that they are not. Now add into the mix rotation angles close to 180 degrees, and the there will be distortions.

Stable variant from libelea (please notice, that while usually libelea has column-major order, packed 3x4 matrices, as used here, have row-major order. So to port the code to OpenGL matrices you have to transpose matrix):

/* Get normalized transformed unit vectors of rotated space */
Vector3f axx(Vector3f(m[0], m[4], m[8]).normalize ());
Vector3f axy(Vector3f(m[1], m[5], m[9]).normalize ());
Vector3f axz(Vector3f(m[2], m[6], m[10]).normalize ());

/* Calculate transposition vectors */
/* I.e. how the unit vector endpoints will be transformed */
Vector3f tx = axx - Vector3fX;
Vector3f ty = axy - Vector3fY;
Vector3f tz = axz - Vector3fZ;

/* Find two biggest transpositions */
f32 txl2 = tx.length2 ();
f32 tyl2 = ty.length2 ();
f32 tzl2 = tz.length2 ();
Vector3f *a = &tx;
Vector3f *b = &ty;
if ((txl2 < tyl2) && (txl2 < tzl2)) {
    a = &ty;
    b = &tz;
} else if ((tyl2 < txl2) && (tyl2 < tzl2)) {
    a = &tz;
    b = &tx;
}

/* Get axis (cross product of two biggest transpositions)*/
Vector3f ax(*a * *b);
if (!ax.length2 ()) {
    /* Zero rotation */
    return Quaternionf0;
}
ax.normalizeSelf ();

/* Get vector perpendicular to axis and transform it */
Vector3f s(a->normalize ());
Vector3f t(m.transformVector3 (s));
t.normalizeSelf ();

/* Find rotation angle (between vector and its transform) */
Vector3f s_t(s * t);
f32 e = Vector3f::scalarProduct (s_t, ax);
if (e < 0) {
    ax = -ax;
}
f32 theta = Vector3f::angle (s, t);
f32 sint2 = sinf (theta / 2);
f32 cost2 = cosf (theta / 2);
return Quaternionf(ax[X] * sint2, ax[Y] * sint2, ax[Z] * sint2, cost2);

The idea is, that any vector perpendicular to axis will be transformed to vector also perpendicular to axis. So we can get axis as the cross product between such vector and its transform.
Now, as you can see, this is much slower, involving the vector normalization, cross products and trigonometric functions. But hopefully you do not have to do this very often - I mostly only needed it for importing foreign animation formats and compacting rotation matrices for serializing.
I have used this method quite extensively and it has excellent stability.

Next time I'll show, how to avoid sudden instabilities if converting rotation matrix to Euler angles.

DAZ Victoria 4

Have been trying to get DAZ 3D base model Victoria 4 to load and render in Khayyam. Things are slowly getting better, although it is not in very usable state yet.

In-engine screenshot of Victoria 4 with Jamie hair

• Body is loaded normally and most joints behave well
• Morphs do not work yet 
• Hair can be loaded only as static object (no morphing, not attached to skeleton)
• Textures can be applied (only diffuse, specular and transparency)

With some patience it can be made looking pretty good. Like most Poser models, the textures are either missing from base model, or so ugly they are pretty useless. So one has to apply textures one-by-one, using material overrides.

POV-Ray rendering of Victoria 4 with Jamie hair

Getting Poser stuff to work is pretty high-priority now. Especially because Illusion does not seem to be planning to release anything noteworthy for some time. Meanwhile at places like Renderosity there is almost infinite amount of very high quality content available - although most of it is a bit too realistic for my taste.

Sara and Cyberdemon (and lighting)

I have been ill (sick if you are American) for the last week, so there was very little interesting things happening with Khayyam. But at last I am feeling better now, so hopefully more will come.

I took few evenings some weeks ago and fixed diffuse and specular lighting (and tangent map calculation). As as result, there is again a reason to read normal- and specular map textures for objects.

Sara does not surrender easily (of course AK-47 itself is also legendary)

At least Doom creatures will look much better this way. Unfortunately Illusion games do not use normal maps, so there are no visible enhancements - unless you want to emulate the in-game plastic oily look.

But - it is possible to replace textures on Illusion (and other formats too) models with materialDNS node (look at tutorial for how exactly) - and this material node supports both normal and specular maps. So if someone wants to experiment, surface effects can be done with some effort.

There are still some big things to do for lighting - like cascading shadow maps, shadows for point and parallel lights and so on. Although immediately I want to play with skin shader instead...

Snapshot release 20110130

I have made a snapshot release of the latest development tree as of 2011/01/30. This time there are both source tarball and precompiled Win32 binaries.
It can be found from SourceForge project page, or use direct download links:

• Win32 binary (9.2 MB) 
• Source package (3,7 MB)

What is new (compared to 4th of January release)?

• OBJ and 3DS importers (static objects)
• OBJ exporter
• Fixed Collada importer/exporter (works with Blender 2.5)
• Preliminary IK solver for skeleton posing
• Bone constraints for IK solver
• Object library with thumbnails to save useful models

Basically the available list of models increased dramatically. There are literally hundreds of different sources of free or cheap 3D models in OBJ and 3DS formats.

Bonus model

As a extra bonus, Sara tween girl Collada model is included in release (in both versions).

Sara tween girl, Blender model

The Blender source file can be downloaded from TurboSquid.

Mini-tutorials

I have previously written two mini tutorials, both available from this site. In case you missed these, they can be found at:

• How to import SBZ figures (2010/12/22)
• How to export composed scenes for POV-Ray rendering (2011/01/30)

Have fun!
Lauris Kaplinski

Rendering with POV-Ray

POV-Ray (The Persistence of Vision Raytracer) is high quality free raytracing renderer. Khayyam can export most models and scenes into POV-Ray file format (.pov) to allow special effects, like focal blur, global illumination and so on. Following is a mini-tutorial for scene composing and exporting.

1. Build your scene

This is really outside of this tutorial so only few hints.

As the first thing, open Document Tree editor (Dialogs -> Document Tree). It has property pages for the most of object types in scenes and most fine-tuning has to be done with it.

Add lights. It can be done later, but a good idea is to add at least one directional light as the first step (Objects -> Add -> Lights -> Directional Light) because otherwise model will look plain and it is hard to understand their poses exactly.

Import objects and figures (Objects -> Import). Currently there is support for most Illusion model formats (althoughthe ones from cell-shaded games do no look nice), MD5 (Doom3 and Quake4), OBJ, 3DS, Collada (partial) and Poser (partial) file formats.

Scale models. Default scales are mostly OK for Illusion, MD5 and Poser files, but the ones imported from other formats may be of very different scales. For scaling, select object (click on it when it has yellow highlight box) - the highlight box should turn orange. Now look at Document Tree dialog - either object of type Figure (poseable object) or StaticMesh (non-poseable) should be selected.

Document tree showing selected figure

Figure and StaticMesh are toplevel containers that can include many geometries (body, hair, clothes...). Select one geometry node in tree and the content of the rightmost notebook pane should change. Each notebook page corresponds to one level in selected object type hierarchy so switching pages gives you different levels of controls.

Click on Geometry tab and there are controls applicable to all Geometry types, among others a scale spinner.

Document tree with SkinnedGeometry selected Geometry property page open

Place models. For this select model and transform control should appear. If it is not visible, try pressing Q key for few times - it cycles the control placement between object origin, object bounding box and the center of screen. Pressing SPACE key toggles between move and rotation controls. The same can be achieved with clicking on already selected model.

Take a look at hotkeys list (Help -> Hotkeys) - there are some useful keyboard and mouse hotkeys, like hiding controls temporarily and so on.

Pose figures - to do this, select figure (poseable object, like Illusion character, MD5 or Collada model) and choose skeleton tool from toolbar. If the selected object is of figure type, bluish skeleton should appear. Nothing will happen, if the object is not of posebale type.

Selected object with skeleton control

Posing individual bones is quite similar to placing objects - although normally you do not want to move bones (it breaks body very unnaturally) but only rotate them.

Save scene. Save often, as Khayyam is still alpha and crashes a lot.

2. Set up lighting

Special lighting is the main reason you want to raytrace your scene. So take some time and try to set up lighting as well as possible - it really makes difference in final images.
Add spotlights (Objects -> Add -> Lights -> Spotlight). There can be practically unlimited number of spotlights, but they make rendering slower. Also semitransparent objects (hair) cannot currently be lighted by more than four spotlights.
Set light direction for spotlights. The easiest method to do this, is copying light direction from camera.

• Move camera to show the scene light you want it to be seen from the spotlight.
• From document tree select SpotLight node
• From right pane select Light tab
• Click From Camera button
• From the same place you can set the light intensity and color

Document tree with Spotlight selected and Light property page open

Ambient lighting adds undirectional light to the whole scene. It can be used to make shadow regions to have less contrast - but it tends to make scene plain.
Diffuse lighting determines the actual spotlight intensity and color.
Direct lighting if for specular highlights, but at moment it does not work reliably.


If you want to add light with soft shadows, use Area Lights. They are similar to spotlights, but are composed of many small lights on rectangular grid. Consequently they can slow down rendering quite a lot, so by default they only have 2 x 2 nodes. To achieve softer effect, you can add more sublights from document tree dialog.

3. Export scene to POV-Ray format

To export POV-Ray scene, select File -> Export POV-Ray Scene...

POV-Ray export dialog

The exporter dialog has three tabs. Normally you want to start from default values, as these render reasonably fast and you can check how the scene looks like. Khayyam and POV-Ray use slightly different lighting formula (actually very different if you use radiosity) so you may want to adjust light intensities and colors before doing final render.

If you have exported scene once and have not added any new objects to it, you can uncheck Export Textures as PNG files checkbox - this saves some exporter time, as it does not have to write texture files over and over.

Lighting properties from POV-Ray export dialog

If you are otherwise content with the scene lookout, you can enable Radiosity. This add indirect lighting to objects - i.e. light bouncing back and forth from surfaces and thus makes the sene much more natural. But it makes rendering considerably slower and can mess up lighting, as in Khayyam you cannot estimate easily, how much reflecting light changes the overall scene brightness - so experiment.

Camera properties from POV-Ray export dialog

Adding Aperture makes objects out of focus (nearer and further) slightly blurred - how much, depends on aperture value (smaller value reduces blur). The focus point is determined by the rotation center in Khayyam (fading target image that appears if you move camera) so place it close to the objects or features you want to make sharp.
Aperture makes rendered images much more photo-like and can add a lot to the mood. But it makes rendering slower.


Click Export and save file with the extension .pov. If you have already installed POV-Ray, it has registered itself as default handler for those file types.


 4. Render it with POV-Ray

Using POV-Ray is outside this tutorial. Some hints though:

• POV-Ray opens all files in builtin editor, so to make small changes you do not have to re-export scene again.
• Rendering starts with button Run
• You can select render size from top-left menu box
• Render result will automatically be saved as png file to the same directory as pov file.
• POV-Ray uses custom scene description language that can be quite complex to understand.

Fully set up scene in Khayyam

The same scene rendered with POV-Ray (using radiosity and aperture). As you can see, using radiosity added lot of extra light to environment.

And this is it! Have fun!

Sara

I have for long time planned to bundle some freeware models with Khayyam. Now, at last, the first iteration of Sara - a young female - is nearing completion. She is quite low-poly, so she will be better suited for real-time animation than for renders.

Sara rigged in Blender

She still needs hair and some more clothes before making a public debut.

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).