How to get local rotation from Matrix but taking parent transforms into account beI have this simplified blend file were I need to match the bone 'shoulder' to this other bone called 'shoulder_extra_snap':

https://pasteall.org/blend/9f468cbe96f145eb9d77553f3ebf1b89

So, both bones have constraints that rotate them. In order to make the rotations match, I'm calculating a rotational difference with matrix.to_euler. and then I'm applying that rotation to 'shoulder'.

This is the script you will find inside of the blend file:

import bpy
from mathutils import Vector

Original_rot = bpy.data.objects['Armature'].pose.bones['shoulder'].matrix_basis.copy()

bpy.data.objects['Armature'].pose.bones['shoulder'].rotation_euler = Vector((0,0,0))
bpy.context.evaluated_depsgraph_get()


a = bpy.data.objects['Armature'].pose.bones['shoulder'].matrix.copy()
b = bpy.data.objects['Armature'].pose.bones['shoulder_extra_snap'].matrix.copy()

bpy.data.objects['Armature'].pose.bones['shoulder'].matrix = b

difEuler = Vector(b.to_euler()[:]) - Vector(a.to_euler()[:]) + Vector(Original_rot.to_euler()[:])


bpy.data.objects['Armature'].pose.bones['shoulder'].rotation_euler = difEuler


This calculation works as long as I don't rotate the 'master' bone, which is the parent of the other bones. Therefore I need to fix that.

I'm pretty sure the solution must be to multiply some matrix by some other matrix so that I can get the local rotation difference but also take the transforms of the hierarchy into account, but I don't know anything about it! (Although I did read some answers regarding similar topics around here, but still I don't get it and it's not my area of expertise, haha)

Any help will be appreciated! thanks!
So quaternions worked! Thanks BatFinger!

import bpy
from mathutils import Vector

Original_rot = bpy.data.objects['Armature'].pose.bones['shoulder'].matrix_basis.copy()

bpy.data.objects['Armature'].pose.bones['shoulder'].rotation_euler = Vector((0,0,0))
bpy.context.evaluated_depsgraph_get()


a = bpy.data.objects['Armature'].pose.bones['shoulder'].matrix.copy()
b = bpy.data.objects['Armature'].pose.bones['shoulder_extra_snap'].matrix.copy()

bpy.data.objects['Armature'].pose.bones['shoulder'].matrix = b

q1 = a.to_quaternion()

q2 = b.to_quaternion()

difQuat = q1.rotation_difference(q2)

difEuler = difQuat.to_euler()

bpy.data.objects['Armature'].pose.bones['shoulder'].rotation_euler = Vector(difEuler) + Vector(Original_rot.to_euler())


Might not be the best code but it works!