The answer is this:

1. Jacobian needs to be calculated globally. For each joint in global coordinates rotation axis needs to be calculated as a cross product between joint-to-end and joint-to-target. The vector that goes into the Jacobian position field is a cross product between joint-to-end and rotation axis - it is a vector along which movement needs to be made (wikipedia names it: cross-radial component).

    // For each joint (except end-effector) JointPosition = ToWorld(jointChainID[j]).Translation(); JointToEnd = EndPosition - JointPosition; JointToTarget = TargetPosition - JointPosition; // Calculate Rotation Axis RotationAxis = JointToEnd.Cross(JointToTarget); RotationAxis.Normalize(); // Cross-Radial Vector Result = JointToEnd.Cross(RotationAxis); // Jacobian column entry jacobian(0, j) = Result.x; jacobian(1, j) = Result.y; jacobian(2, j) = Result.z; jacobian(3, j) = RotationAxis.x; jacobian(4, j) = RotationAxis.y; jacobian(5, j) = RotationAxis.z; 

2. The result dQ contains rotation angles for each joint along RotationAxis (in global coordinates). Then I apply changes locally (this is weird, as I apply global changes to local joint coordinates, but somehow this works):

    // Apply changes to local skinning matrices // For each joint (except end-effector) jointID = JointChainID[j]; gTransform = Matrix::CreateFromAxisAngle(cache.RotationAxes[j], ToRadians(dQ(j))); SkinningLocal[jointID] = gTransform * SkinningLocal[jointID]; 

The answer is this:

1. Jacobian needs to be calculated globally. For each joint in global coordinates rotation axis needs to be calculated as a cross product between joint-to-end and joint-to-target. The vector that goes into the Jacobian position field is a cross product between joint-to-end and rotation axis - it is a vector along which movement needs to be made (wikipedia under "angular velocity" names it: cross-radial component).

    // For each joint (except end-effector) JointPosition = ToWorld(jointChainID[j]).Translation(); JointToEnd = EndPosition - JointPosition; JointToTarget = TargetPosition - JointPosition; // Calculate Rotation Axis RotationAxis = JointToEnd.Cross(JointToTarget); RotationAxis.Normalize(); // Cross-Radial Vector Result = JointToEnd.Cross(RotationAxis); // Jacobian column entry jacobian(0, j) = Result.x; jacobian(1, j) = Result.y; jacobian(2, j) = Result.z; jacobian(3, j) = RotationAxis.x; jacobian(4, j) = RotationAxis.y; jacobian(5, j) = RotationAxis.z; 

2. The result dQ contains rotation angles for each joint along RotationAxis (in global coordinates). Then I apply changes locally (this is weird, as I apply global changes to local joint coordinates, but somehow this works):

    // Apply changes to local skinning matrices // For each joint (except end-effector) jointID = JointChainID[j]; gTransform = Matrix::CreateFromAxisAngle(cache.RotationAxes[j], ToRadians(dQ(j))); SkinningLocal[jointID] = gTransform * SkinningLocal[jointID];