How can I manually calculate bpy.types.PoseBone.matrix using Blender's Python API? I have this .blend file which has a simple scene, and a python script. There's a readme section in the script, but to give a general idea:

There are two identical meshes in the scene, one skinned, and the other unskinned. The skinned mesh has a 20 frame animation. If you move to any frame of the animation, and then run the script, the script will modify the position of the unskinned mesh's vertices, to match the vertex positions of the skinned mesh, for the current frame of animation.

The formula I use to calculate the new vertex positions is:

nv = nv + ((boneWeight / TW) * ( bms[boneName] * (ibps[boneName] * v.co) ))


Where:


boneWeight: weight of the current bone that influences this vertex
TW: total weight of all bones that affect this vertex
bms[boneName]: bind pose of the current bone (bpy.data.objects[ARMATURE_NAME].pose.bones[BONE_NAME].matrix)
ibps[boneName]: inverted bind pose of the current bone (bpy.data.objects[ARMATURE_NAME].data.bones[BONE_NAME].matrix_local.inverted())
v.co: the original vertex coordinates


It's worth noting that, in my Python script, I do not traverse the hierarchy of bones to calculate the correct transform matrices; simply accessing the pose matrix and inverted bind pose matrix for the current bone is enough. 

So with this figured out, I now want to export the bind pose and inverted bind pose of every bone, for every keyframe in this animation. In order to keep my exported file size to a minimum, I only want to export the keyframes of an animation, and have my game engine generate the bind pose and inverted bind pose of all bones for all non-keyframes in the animation. 

I have found the interpolation functions for Matrices in Blender's source code, and I have re-implemented them in Java. For the root bone of an armature, the matrices that my engine generates for the bind pose and inverted bind pose are correct for all frames in my animation. 

However, for non-root bones, the bind pose and inverted bind pose are incorrect, i.e. they do not match the pose_bone.matrix and data_bone.matrix_local.inverted() readouts that Blender gives in the python console.

So, it seems to me that, it isn't enough to export the bind pose and inverted bind pose matrices, and interpolate between them for given keyframes.

Instead, I either need to: 


Export and interpolate between the matrices/vectors from which the bind pose and inverted bind pose are derived from, and then multiply them all together to give the bind pose and inverted bind pose of the bone on the engine side


OR


I continue to export the bind pose and inverted bind pose from blender, but after I have generated the matrices on the engine side, for non-root bones, I multiply the generated matrix with a matrix/inverted matrix belonging to that bone's parent (but what parent bone matrix should I multiply against?)


Essentially, I need to know how, using the Python API in Blender, I can generate a matrix equivalent to: 

bpy.data.objects['Armature'].data.pose.bones['Bone.001'].matrix


Without directly accessing bpy.data.objects['Armature'].data.pose.bones['Bone.001'].matrix, if that makes sense?

So, I need to know, step-by-step, what child/parent bone matrices (available via the Python API) to multiply together, to give a matrix equal to the bind pose (bpy.data.objects['Armature'].data.pose.bones['Bone.001'].matrix) of a given bone.

I've been linked to this diagram before: 

But I wasn't able to follow along with it using only the Python API, as I was only able to find the equivalent functions for each step in the C code, not the Python API.
armature.data.bones["Bone"].matrix_local


Is the bone's world un-posed matrix (the pose_bone.matrix in rest state).

 armature.pose.bones["Bone"].matrix_basis


Is the bone's local pose matrix.

From there we can calculate any bone's world matrix by recursively traversing it's parents. Note: This will not include any constraints or IK (but it will include drivers or actions):

def matrix_world(armature_ob, bone_name):
    local = armature_ob.data.bones[bone_name].matrix_local
    basis = armature_ob.pose.bones[bone_name].matrix_basis

    parent = armature_ob.pose.bones[bone_name].parent
    if parent == None:
        return  local * basis
    else:
        parent_local = armature_ob.data.bones[parent.name].matrix_local
        return matrix_world(armature_ob, parent.name) * (parent_local.inverted() * local) * basis


You can test it on an empty and some armature with some active pose_bone. Also the armature_object's transformation is accounted here:

import bpy

empty_ob = bpy.data.objects["Empty"]
armature_ob = bpy.data.objects["Armature"]

bone_name = bpy.context.active_pose_bone.name
empty_ob.matrix_world = armature_ob.matrix_world * matrix_world(armature_ob, bone_name )