How to transfer 3d joints position to Rotation of limbs/torso

I have a set of 3d positions of joints.

I’d like to convert them to Euler/Quaternion which can be applited to bones of GLB/FBX model.

My previous approach was to simply calculate a vector by minus joints positions. eg. elbow - shoulder. Then get the quaternion from this vector and the initial vector (0,-1,0).
But the problem is, it doesn’t consider the orientation of the limbs. It works if the limbs are just simple Cylinder geometry, but with real 3d model, I assume there is a more accurate approach.