How to calculate the angles between axis and bounding box normals? I would like to determine the torsion of an object bounding box in reference to the global coordinate system (such that for a axis aligned object all angles will be 0). 

My idea was to create a bounding box object bb via bpy.ops.mesh.primitive_cube_add() and copy its parameters from the object obj:

centre = sum((Vector(b) for b in obj.bound_box), Vector())
bb.dimensions = obj.dimensions
bb.rotation_euler = obj.rotation_euler
bb.location = (centre[0], centre[1], centre[2])
bb.location = obj.matrix_world * bb.location


Than I tried to calculate the angles:

for poly in obj.data.polygons:
    # angle of bb normal against z axis
    angle = Vector((0,0,1)).angle(obj.matrix_world * poly.normal)
    print("%6i  %.3f°" % (poly.index, degrees(angle)))


But the values don't make any sense. I am not sure whether this is a good approach at all.

Any suggestion to solve this problem?
The orientation of an object is obj.rotation_euler or obj.rotation_quaternion (and whichever one does not match the object's rotation_mode can be invalid or inaccurate).

The orientation can also be recovered from obj.matrix_world.decompose()[1], but that's a little bit of overkill.

For a slightly deeper dive down the rabbit hole, check out Relationship between global vertex coordinates and mesh object&#39;s matrix attributes .