kabsch alignment algorithm swapped mobile and fixed when creating rotation matrix

[jonathanking-absci]

Hello,

I'm working on some code that suggests that when creating the rotation matrix via superimpose, the kabsch algorithm implementation simply has fixed and mobile swapped as in superimpose.py:_superimpose(). Do you think this is an accurate assessment? I think the definitions of x and y in that function should be flipped. Or, in other words, I think the kabsch algorithm should have mobile.T dot reference, when this implementation is doing reference.T dot mobile.

My method:

• create a set of points
• create a rotation matrix R
• rotate points by R to get points*
• use superimpose(points*, points) to calculate R2
• use superimpose(points, points*) to calculate R3
• R2 != R
• R3 == R

Just wanted to bring that to your attention! Or perhaps I'm mistaken, and that would be nice to know as well 🙂

Best,
Jonathan

[padix-key]

Hi, I will look into it. Intuitively, I would also presume, that R2 = R should be the case.

[padix-key]

I took this issue to perform some refactoring on superimpose() in #526, including some new tests. One test also tests simply if the rotation matrix is as expected for a 90 degrees rotation around the z-axis:

biotite/tests/structure/test_superimpose.py

Lines 46 to 69 in 29ecc1e

	def test_rotation_matrix():
	"""
	Create randomly generated coordinates *a* and rotate them via a
	rotation matrix to obtain new coordinates *b*.
	``superimpose(b, a)`` should give the rotation matrix.
	"""
	N_COORD = 100
	# A rotation matrix that rotates 90 degrees around the z-axis
	ref_rotation = np.array([
	[0, -1, 0],
	[1, 0, 0],
	[0, 0, 1]
	])
	np.random.seed(0)
	original_coord = np.random.rand(N_COORD, 3)
	# Rotate about 90 degrees around z-axis
	rotated_coord = struc.rotate(original_coord, angles=(0, 0, np.pi/2))
	_, transform = struc.superimpose(rotated_coord, original_coord)
	test_rotation = transform.rotation
	assert test_rotation.flatten().tolist() \
	== pytest.approx(ref_rotation.flatten().tolist(), abs=1e-6)

However, the rotation matrix is as expected, but it might also be fixed by the refactoring.
As soon as #526 is merged or released, could you check again?

If the issue persists, could it be possible the error is this step:

rotate points by R to get points*

which, I think, could easily happen as the inverse of a rotation matrix is simply its transpose.

[jonathanking-absci]

I went back and looked at this. You were right. I was rotating points via points @ rotation_matrix and not points @ rotation_matrix.T. Sorry for the confusion, and thanks for the help!