Presuming all objects have initially been rotated from their origin and there is no other operation on them (no translation for instance).
We need to find two vertices that are not aligned with the mesh origin.
If these two exist, we can calculate two rotations in order to find the wanted result.
Below, green/blue is the object and yellow/red is the reference.
First rotation is the one that goes from one of the vectors to the other:
Once this rotation applied, to the green and blue, we have green over yellow and we have blue rotated:
Now we want this rotation, with the green/yellow as axis:
The following code does it. Select all the objects you want to "unrotate" with the reference one as active object and run the code.
import bpy
from math import degrees
from mathutils import Matrix
# loop over the vertices in order to find
# 2 of them that are not aligned with the
# origin of the mesh
def pick_points(verts, epsilon = 0.01):
for i1 in range(len(verts)):
v1 = verts[i1].co.normalized()
for i2 in range(i1 + 1, len(verts)):
v2 = verts[i2].co.normalized()
if v1.cross(v2).length > epsilon:
return i1, i2, v1, v2
return None, None, None, None
# Take the reference object, the active one
reference = bpy.context.active_object
# Find reference vertices for it
i1, i2, v1, v2 = pick_points(reference.data.vertices)
if i1 is None:
print("All vertices are aligned")
else:
# The vector formed by the projection of v2 on v1 to v2
pv2 = v2 - v2.project(v1)
# Loop over the selected object
for o in [o for o in bpy.context.selected_objects if o != reference]:
# Get the two vertices at same indices
ov1 = o.data.vertices[i1].co.normalized()
ov2 = o.data.vertices[i2].co.normalized()
# Get the rotation separating ov1 and the reference v1
q1 = ov1.rotation_difference(v1)
# Rotate ov2 with this rotation
ov2 = q1 @ ov2
# Get the vector formed by the projection of ov2 on v1 to ov2
pov2 = ov2 - ov2.project(v1)
# Get the rotation separating these to "p" vectors
q2 = pov2.rotation_difference(pv2)
# Combine the rotations = how to go from o to the reference
q = q2 @ q1
# Extract current world matrix components of o
t, _, s = o.matrix_world.decompose()
# Assign it the same component, but replacing the rotation part
# which is calculated above inverted
o.matrix_world = Matrix.LocRotScale(t, q.inverted(), s)
# Link objects to the ref/active one
bpy.ops.object.make_links_data(type='OBDATA')
Here is your file (I've used V4.0 here), with a copy of the screws and the script applied.
All works as expected:
Select them all with the original screw as active.
Run the script.
And the screws are instances with the wanted rotations.