Understanding GLTF 2.0 Skin Animation Joint Matrix Calculations

Question

I am a novice in skin animation so I've been spending a good month trying to figure out why my animation is not currently outputting the proper matrix palette, after performing calculations on my joint matrices. In the GLTF 2.0 skin tutorials, it states how to calculate the joint matrix for skinning:

jointMatrix(j) =
 globalTransformOfNodeThatTheMeshIsAttachedTo^-1 *
 globalTransformOfJointNode(j) *
 inverseBindMatrixForJoint(j);

and so I proceed to do this for my animation engine:

void Animation::DoSampleJob(AnimJobSubmitInfo& job, r32 gt)
{
  if (!job._output->_currState._bEnabled) { return; }
  // This part is just calculating the local time progression, no issue here.
  r32 tau = job._output->_currState._tau;
  r32 rate = job._output->_currState._fPlaybackRate;
  r32 lt = job._output->_currState._fCurrLocalTime + gt * rate;
  if (lt > job._pBaseClip->_fDuration) {
    lt -= job._pBaseClip->_fDuration;
    job._output->_currState._tau = gt;
  }
  if (lt < 0.0f) {
    lt = job._pBaseClip->_fDuration + lt;
    if (lt < 0.0f) {
      lt += job._pBaseClip->_fDuration;
      job._output->_currState._tau = gt;
    }
  }
  job._output->_currState._fCurrLocalTime = lt;
  Skeleton* pSkeleton = Skeleton::GetSkeleton(job._pBaseClip->_skeletonId);

  u32 currPoseIdx = 0;
  u32 nextPoseIdx = 0;

  GetCurrentAndNextPoseIdx(&currPoseIdx, &nextPoseIdx, job._pBaseClip, lt);

  ApplyMorphTargets(job._output, job._pBaseClip, currPoseIdx, nextPoseIdx, lt);

  if (EmptyPoseSamples(job._pBaseClip, currPoseIdx, nextPoseIdx)) { return; }

  AnimPose* currAnimPose = &job._pBaseClip->_aAnimPoseSamples[currPoseIdx];
  AnimPose* nextAnimPose = &job._pBaseClip->_aAnimPoseSamples[nextPoseIdx];

  for (size_t i = 0; i < job._pBaseClip->_aAnimPoseSamples[currPoseIdx]._aLocalPoses.size(); ++i) {
    JointPose* currJoint = &currAnimPose->_aLocalPoses[i];
    JointPose* nextJoint = &nextAnimPose->_aLocalPoses[i];
    Matrix4 localTransform = LinearInterpolate(currJoint, nextJoint, currAnimPose->_time, nextAnimPose->_time, lt);
    job._output->_currentPoses[i] = localTransform;
  }

  ApplySkeletonPose(job._output->_finalPalette, job._output->_currentPoses, pSkeleton);
}


void Animation::ApplySkeletonPose(Matrix4* pOutput, Matrix4* pLocalPoses, Skeleton* pSkeleton)
{
  if (!pSkeleton) return;
  // This is where the issue is at, somewhere...
  for (size_t i = 0; i < pSkeleton->_joints.size(); ++i) {
    Matrix4 parentTransform;
    Matrix4 currentPose;
    u8 parentId = pSkeleton->_joints[i]._iParent;
    if (parentId != Joint::kNoParentId) {
      parentTransform = pLocalPoses[parentId];
    }
    // Now become world space joint matrices
    currentPose = pLocalPoses[i] * parentTransform;
    pLocalPoses[i] = currentPose;
  }

  for (size_t i = 0; i < pSkeleton->_joints.size(); ++i) {
    pOutput[i] = pSkeleton->_joints[i]._InvBindPose * pLocalPoses[i] * pSkeleton->_joints[i]._invGlobalTransform;

  }
}

Unfortunately, the outcome was not expected:

Investing it, I proceeded to remove the calculations of the animation, and just calculate the joint matrix as the inverse bind pose along with the global joint transform:

void Animation::ApplySkeletonPose(Matrix4* pOutput, Matrix4* pLocalPoses, Skeleton* pSkeleton)
{
  if (!pSkeleton) return;

  for (size_t i = 0; i < pSkeleton->_joints.size(); ++i) {
    Matrix4 parentTransform;
    Matrix4 currentPose;
    u8 parentId = pSkeleton->_joints[i]._iParent;
    if (parentId != Joint::kNoParentId) {
      parentTransform = pLocalPoses[parentId];
    }
    // Now become work space joint matrices
    currentPose = pLocalPoses[i] * parentTransform;
    pLocalPoses[i] = currentPose;
  }

  for (size_t i = 0; i < pSkeleton->_joints.size(); ++i) {
    // Just calculating only the inverse bind pose, and global joint transform, removing the current pose.
    pOutput[i] = pSkeleton->_joints[i]._InvBindPose * pSkeleton->_joints[i]._invGlobalTransform.Inverse();

  }
}

And the outcome is as expected, when not animating:

And here is where I am at my wits end. Unsure why, or how, to fix this issue with calculating my joints. Could it be something I am not doing correctly with calculating the current world joint matrices, prior to applying the inverse bind pose and global joint transform? Or something even prior to that? I don't specialize in animation, only graphics, but it's a neat thing to understand how it all works from behind the scenes :). Unfortunately, there are many ways to do animation skinning, so I hope to find some help in this, as I have been spending a good month on this particular issue in gltf. Much appreciated on the assistance!

Additionally you may also want to look at how I parse the skin, and global joint transforms too:

static skeleton_uuid_t LoadSkin(const tinygltf::Node& node, const tinygltf::Model& model, Model* engineModel, const Matrix4& parentMatrix)
{
  if (node.skin == -1) return Skeleton::kNoSkeletonId;

  Skeleton skeleton;
  tinygltf::Skin skin = model.skins[node.skin];
  b32 rootInJoints = false;
  for (size_t i = 0; i < skin.joints.size(); ++i) {
    if (skin.joints[i] == skin.skeleton) {
      rootInJoints = true; break;
    }
  }
  skeleton._joints.resize(skin.joints.size());
  skeleton._name = skin.name;
  skeleton._rootInJoints = rootInJoints;

  const tinygltf::Accessor& accessor = model.accessors[skin.inverseBindMatrices];
  const tinygltf::BufferView& bufView = model.bufferViews[accessor.bufferView];
  const tinygltf::Buffer& buf = model.buffers[bufView.buffer];

  struct NodeTag {
    i32               _gltfParent;
    u8                _parent;
    Matrix4           _parentTransform;
  };

  std::map<i32, NodeTag> nodeMap;
  for (size_t i = 0; i < skin.joints.size(); ++i) {
    size_t idx = i;
    Joint& joint = skeleton._joints[idx];
    i32 skinJointIdx = skin.joints[i];
    const tinygltf::Node& node = model.nodes[skinJointIdx];
    NodeTransform localTransform;

    auto it = nodeMap.find(skinJointIdx);
    if (it != nodeMap.end()) {
      NodeTag& tag = it->second;
      localTransform = CalculateGlobalTransform(node, tag._parentTransform);
      joint._iParent = tag._parent;
      joint._invGlobalTransform = localTransform._globalMatrix.Inverse();
    } else {
      localTransform = CalculateGlobalTransform(node, Matrix4());
      joint._iParent = 0xff;
      joint._invGlobalTransform = localTransform._globalMatrix.Inverse();
    }

    DEBUG_OP(joint._id = static_cast<u8>(skinJointIdx));
    for (size_t child = 0; child < node.children.size(); ++child) {
      NodeTag tag = { static_cast<u8>(skinJointIdx), i, localTransform._globalMatrix };
      nodeMap[node.children[child]] = tag;
    }
  }

  const r32* bindMatrices = reinterpret_cast<const r32*>(&buf.data[bufView.byteOffset + accessor.byteOffset]);

  for (size_t i = 0; i < accessor.count; ++i) {
    Matrix4 invBindMat(&bindMatrices[i * 16]);
    skeleton._joints[i]._InvBindPose = invBindMat;
  }

  Skeleton::PushSkeleton(skeleton);

  engineModel->skeletons.push_back(Skeleton::GetSkeleton(skeleton._uuid));
  return skeleton._uuid;
}


static void LoadNode(const tinygltf::Node& node, const tinygltf::Model& model, Model* engineModel, const Matrix4& parentMatrix, const r32 scale)
{
  NodeTransform transform = CalculateGlobalTransform(node, parentMatrix);
  if (!node.children.empty()) {
    for (size_t i = 0; i < node.children.size(); ++i) {
      LoadNode(model.nodes[node.children[i]], model, engineModel, transform._globalMatrix, scale);
    }
  }

  if (node.skin != -1) {
    skeleton_uuid_t skeleId = LoadSkin(node, model, engineModel, transform._globalMatrix);
    Mesh* pMesh = LoadSkinnedMesh(node, model, engineModel, transform._globalMatrix);
    pMesh->SetSkeletonReference(skeleId);
  }
  else {
    LoadMesh(node, model, engineModel, transform._globalMatrix);
  }
}

If more information is required, please feel free to see the source code on my github.

Many thanks!