I just had a pm related to my previous post and the reply got a bit long. I’ll post it here in case it helps anyone else:
That little test took me about a year of thinking about it, getting confused, putting it to one side and then going back to it many times.
I’m kinda working on something related right now and this is not the kind of thing that can be explained in just a few words. I’ll probably post an in-depth explanation (pdf maybe) when I’ve got my current project working.
Until then, I’d advise looking at http://www.euclideanspace.com/ to find a bit about matrices and quaternions and how they apply to rotations, and then set up some simple tests in Blender.
Some tests that I did are:
Rotate the basic cube using matrices.
90 degrees at a time so you can see the axes changing.
45 degrees at a time so you can see the axes changing.
Rotate the basic cube using quaternions.
90 degrees so it is easy to see.
45 degrees at a time so it is easy to see.
Apply rotations from a cube to a single bone by using action channels. This will be the same as the bone in a chain that has no parents.
Apply rotations from 2 cubes to a 2-bone armature using two cubes and action channels.
Finally include the third cube and bone. If you can deal with the second and third bones in the same way then you can easily scale your script up to as many bones as you like.
Keep it really simple, it’s very easy to get lost when playing with longer bone chains.
The octahedrons in my .blend are not just pretty. Using the convention x, y, z = r, g, b I have applied vertex colours so that the points of the octahedron represent the axes with light red = positive x, dark red = negative x, light blue = positive z, etc. Feel free to append them to your experiments, I found they made it easier to visualise orientations.
In my final implementation of this, I needed the rest orientation of the driver object, the current orientation of the driver object (from which I could obtain the rotation that it had undergone from rest position using Quaternion.difference(Quaternion) see http://www.blender.org/documentation/250PythonDoc/mathutils.html) I also needed the rest orientation of the bone and the amount that the parent bone had been rotated from rest position.
I know that puting all the info in a dict is cool, but its throwing me off.
The dict just made it easy for me to keep track of the values. Try placing them each in their own variable, but be careful of variable names. It’s difficult to find a short but descriptive variable name. Maybe base_bone_rest_orientation_relative_to_world_coordinates, etc. 
#The composite matrix represents the world orientation of the bone in rest position.
obj['bone_rest_orientations']['COG'] = COGmat.copy().invert().to_quat()
Bone1CompositeMat = COGmat * Bone1mat
obj['bone_rest_orientations']['Bone.001'] = Bone1CompositeMat.copy().invert().to_quat()
Bone2CompositeMat = Bone1CompositeMat * Bone2mat
obj['bone_rest_orientations']['Bone.002'] = Bone2CompositeMat.copy().invert().to_quat()
Assuming the armature has not been rotated (armature rotation will mess up this script). All these matrix calculations are to find the rest orientation relative to world coordinates.
The first rest orientation is simply taken from the bone matrix since it has no parent bones to rotate it.
(first bone orientation) = (first bone matrix)
The second rest orientation is the first rest orientation rotated by the second bone matrix since each bone’s orientation is given relative to it’s parent and the parent rotation is required to turn the bone matrix into one that is relative to world space.
(second bone orientation) = (first bone matrix)*(second bone matrix)
The third rest orientation is the second rest orientation rotated by the third bone matrix to allow for the parent rotation.
(third bone orientation) = (second bone orientation)(third bone matrix)
This is the same as
(third bone orientation) = (first bone matrix)(second bone matrix)*(third bone matrix).
Note the invert() functions to invert the matrices.
I’m not entirely sure if I’ve done the inverts at the correct time or used the Quaternion.difference() the right way around, but this works so I’m gonna leave it for now.
In that .blend is an old text file UseRagDoll_a. Ignore this file. In fact, you can unlink it from the .blend. It is not used. To understand how I drive armatures with objects, the only script you need to understand is Control_Armature.py. The dictionaries are necessary only for storing the bone matrix in rest position since the game does not enable you to read this information. The script is about as simple as I can make it. It’s written in a simple linear programming style with the first section (defining the rest position matrices of the bones) only running at the start.
The second section calculates rotations for each of the 3 object/bone pairs in turn, each in it’s own code block. Notice that Bone.001 and Bone.002 are both treated in exactly the same way. The base bone (“Bone”) is treated slightly differently since it doesn’t have a parent rotation to take into account.
I hope this helps a little.
Good luck. 