13.2.8. Homogeneous Transformation Matrices and Quaternions — MDAnalysis.lib.transformations

A library for calculating 4x4 matrices for translating, rotating, reflecting, scaling, shearing, projecting, orthogonalizing, and superimposing arrays of 3D homogeneous coordinates as well as for converting between rotation matrices, Euler angles, and quaternions. Also includes an Arcball control object and functions to decompose transformation matrices.

Authors

Christoph Gohlke, Laboratory for Fluorescence Dynamics, University of California, Irvine

Version

2010.05.10

Licence

BSD 3-clause

13.2.8.1. Requirements

Notes

The API is not stable yet and is expected to change between revisions.

This Python code is not optimized for speed. Refer to the transformations.c module for a faster implementation of some functions.

Documentation in HTML format can be generated with epydoc.

Matrices (M) can be inverted using numpy.linalg.inv(M), concatenated using numpy.dot(M0, M1), or used to transform homogeneous coordinates (v) using numpy.dot(M, v) for shape (4, *) “point of arrays”, respectively numpy.dot(v, M.T) for shape (*, 4) “array of points”.

Use the transpose of transformation matrices for OpenGL glMultMatrixd().

Calculations are carried out with numpy.float64 precision.

Vector, point, quaternion, and matrix function arguments are expected to be “array like”, i.e. tuple, list, or numpy arrays.

Return types are numpy arrays unless specified otherwise.

Angles are in radians unless specified otherwise.

Quaternions w+ix+jy+kz are represented as [w, x, y, z].

A triple of Euler angles can be applied/interpreted in 24 ways, which can be specified using a 4 character string or encoded 4-tuple:

  • Axes 4-string: e.g. ‘sxyz’ or ‘ryxy’

    • first character : rotations are applied to ‘s’tatic or ‘r’otating frame

    • remaining characters : successive rotation axis ‘x’, ‘y’, or ‘z’

  • Axes 4-tuple: e.g. (0, 0, 0, 0) or (1, 1, 1, 1)

    • inner axis: code of axis (‘x’:0, ‘y’:1, ‘z’:2) of rightmost matrix.

    • parity : even (0) if inner axis ‘x’ is followed by ‘y’, ‘y’ is followed by ‘z’, or ‘z’ is followed by ‘x’. Otherwise odd (1).

    • repetition : first and last axis are same (1) or different (0).

    • frame : rotations are applied to static (0) or rotating (1) frame.

References

  1. Ronald Goldman. Matrices and transformations. In Graphics Gems I, pages 472–475. Morgan Kaufmann, USA, 1990.

  2. Ken Shoemake. Uniform random rotations. In David Kirk, editor, Graphics Gems III, pages 124–132. Academic Press, 1992.

  3. Charles F.F. Karney. Quaternions in molecular modeling. Journal of Molecular Graphics and Modelling, 25(5):595–604, 2007. doi:https://doi.org/10.1016/j.jmgm.2006.04.002.

  4. Itzhack Y. Bar-Itzhack. New method for extracting the quaternion from a rotation matrix. Journal of Guidance, Control, and Dynamics, 23(6):1085–1087, 2000. doi:10.2514/2.4654.

  5. Ronald N. Goldman. Vii.4 - more matrices and transformations: shear and pseudo-perspective. In JAMES ARVO, editor, Graphics Gems II, pages 338–341. Morgan Kaufmann, San Diego, 1991. URL: https://www.sciencedirect.com/science/article/pii/B9780080507545500724, doi:https://doi.org/10.1016/B978-0-08-050754-5.50072-4.

  6. Spencer W. Thomas. Vii.1 - decomposing a matrix into simple transformations. In JAMES ARVO, editor, Graphics Gems II, pages 320–323. Morgan Kaufmann, San Diego, 1991. URL: https://www.sciencedirect.com/science/article/pii/B9780080507545500694, doi:https://doi.org/10.1016/B978-0-08-050754-5.50069-4.

  7. Ronald N. Goldman. Vii.2 - recovering the data from the transformation matrix. In James Arvo, editor, Graphics Gems II, pages 324–331. Morgan Kaufmann, San Diego, 1991. URL: https://www.sciencedirect.com/science/article/pii/B9780080507545500700, doi:https://doi.org/10.1016/B978-0-08-050754-5.50070-0.

  8. Ken Shoemake. Euler angle conversion. In Paul S. Heckbert, editor, Graphics Gems IV, pages 222–229. Morgan Kaufmann, 1994.

  9. Ken Shoemake. Arcball rotation control. In Paul S. Heckbert, editor, Graphics Gems IV, pages 175–192. Morgan Kaufmann, 1994.

  10. Ken Shoemake. Quaternions. In URL: https://www.ljll.math.upmc.fr/~frey/papers/scientific%20visualisation/Shoemake%20K.,%20Quaternions.pdf.

  11. James Diebel. Representing attitude: euler angles, unit quaternions, and rotation vectors. In 2006.

  12. W. Kabsch. A discussion of the solution for the best rotation to relate two sets of vectors. Acta Crystallographica Section A, 34(5):827–828, Sep 1978. URL: https://doi.org/10.1107/S0567739478001680, doi:10.1107/S0567739478001680.

  13. Berthold Horn. Closed-form solution of absolute orientation using unit quaternions. Journal of the Optical Society A, 4:629–642, 04 1987. doi:10.1364/JOSAA.4.000629.

  14. JMP van Waveren. From quaternion to matrix and back. In 2005. URL: http://www.intel.com/cd/ids/developer/asmo-na/eng/293748.htm.

Examples

>>> alpha, beta, gamma = 0.123, -1.234, 2.345
>>> origin, xaxis, yaxis, zaxis = (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)
>>> I = identity_matrix()
>>> Rx = rotation_matrix(alpha, xaxis)
>>> Ry = rotation_matrix(beta, yaxis)
>>> Rz = rotation_matrix(gamma, zaxis)
>>> R = concatenate_matrices(Rx, Ry, Rz)
>>> euler = euler_from_matrix(R, 'rxyz')
>>> numpy.allclose([alpha, beta, gamma], euler)
True
>>> Re = euler_matrix(alpha, beta, gamma, 'rxyz')
>>> is_same_transform(R, Re)
True
>>> al, be, ga = euler_from_matrix(Re, 'rxyz')
>>> is_same_transform(Re, euler_matrix(al, be, ga, 'rxyz'))
True
>>> qx = quaternion_about_axis(alpha, xaxis)
>>> qy = quaternion_about_axis(beta, yaxis)
>>> qz = quaternion_about_axis(gamma, zaxis)
>>> q = quaternion_multiply(qx, qy)
>>> q = quaternion_multiply(q, qz)
>>> Rq = quaternion_matrix(q)
>>> is_same_transform(R, Rq)
True
>>> S = scale_matrix(1.23, origin)
>>> T = translation_matrix((1, 2, 3))
>>> Z = shear_matrix(beta, xaxis, origin, zaxis)
>>> R = random_rotation_matrix(numpy.random.rand(3))
>>> M = concatenate_matrices(T, R, Z, S)
>>> scale, shear, angles, trans, persp = decompose_matrix(M)
>>> numpy.allclose(scale, 1.23)
True
>>> numpy.allclose(trans, (1, 2, 3))
True
>>> numpy.allclose(shear, (0, math.tan(beta), 0))
True
>>> is_same_transform(R, euler_matrix(axes='sxyz', *angles))
True
>>> M1 = compose_matrix(scale, shear, angles, trans, persp)
>>> is_same_transform(M, M1)
True

13.2.8.2. Functions

Changed in version 0.11.0: Transformations library moved from MDAnalysis.core.transformations to MDAnalysis.lib.transformations

class MDAnalysis.lib.transformations.Arcball(initial=None)[source]

Virtual Trackball Control.

>>> ball = Arcball()
>>> ball = Arcball(initial=np.identity(4))
>>> ball.place([320, 320], 320)
>>> ball.down([500, 250])
>>> ball.drag([475, 275])
>>> R = ball.matrix()
>>> np.allclose(np.sum(R), 3.90583455)
True
>>> ball = Arcball(initial=[1, 0, 0, 0])
>>> ball.place([320, 320], 320)
>>> ball.setaxes([1,1,0], [-1, 1, 0])
>>> ball.setconstrain(True)
>>> ball.down([400, 200])
>>> ball.drag([200, 400])
>>> R = ball.matrix()
>>> np.allclose(np.sum(R), 0.2055924)
True
>>> ball.next()

Initialize virtual trackball control.

initial : quaternion or rotation matrix

down(point)[source]

Set initial cursor window coordinates and pick constrain-axis.

drag(point)[source]

Update current cursor window coordinates.

getconstrain()[source]

Return state of constrain to axis mode.

matrix()[source]

Return homogeneous rotation matrix.

next(acceleration=0.0)[source]

Continue rotation in direction of last drag.

place(center, radius)[source]

Place Arcball, e.g. when window size changes.

centersequence[2]

Window coordinates of trackball center.

radiusfloat

Radius of trackball in window coordinates.

setaxes(*axes)[source]

Set axes to constrain rotations.

setconstrain(constrain)[source]

Set state of constrain to axis mode.

MDAnalysis.lib.transformations.arcball_nearest_axis(point, axes)[source]

Return axis, which arc is nearest to point.

MDAnalysis.lib.transformations.compose_matrix(scale=None, shear=None, angles=None, translate=None, perspective=None)[source]

Return transformation matrix from sequence of transformations.

This is the inverse of the decompose_matrix function.

Sequence of transformations:

scale : vector of 3 scaling factors shear : list of shear factors for x-y, x-z, y-z axes angles : list of Euler angles about static x, y, z axes translate : translation vector along x, y, z axes perspective : perspective partition of matrix

>>> scale = np.random.random(3) - 0.5
>>> shear = np.random.random(3) - 0.5
>>> angles = (np.random.random(3) - 0.5) * (2*math.pi)
>>> trans = np.random.random(3) - 0.5
>>> persp = np.random.random(4) - 0.5
>>> M0 = compose_matrix(scale, shear, angles, trans, persp)
>>> result = decompose_matrix(M0)
>>> M1 = compose_matrix(*result)
>>> is_same_transform(M0, M1)
True
MDAnalysis.lib.transformations.concatenate_matrices(*matrices)[source]

Return concatenation of series of transformation matrices.

>>> M = np.random.rand(16).reshape((4, 4)) - 0.5
>>> np.allclose(M, concatenate_matrices(M))
True
>>> np.allclose(np.dot(M, M.T), concatenate_matrices(M, M.T))
True
MDAnalysis.lib.transformations.decompose_matrix(matrix)[source]

Return sequence of transformations from transformation matrix.

matrixarray_like

Non-degenerative homogeneous transformation matrix

Return tuple of:

scale : vector of 3 scaling factors shear : list of shear factors for x-y, x-z, y-z axes angles : list of Euler angles about static x, y, z axes translate : translation vector along x, y, z axes perspective : perspective partition of matrix

Raise ValueError if matrix is of wrong type or degenerative.

>>> T0 = translation_matrix((1, 2, 3))
>>> scale, shear, angles, trans, persp = decompose_matrix(T0)
>>> T1 = translation_matrix(trans)
>>> np.allclose(T0, T1)
True
>>> S = scale_matrix(0.123)
>>> scale, shear, angles, trans, persp = decompose_matrix(S)
>>> scale[0]
0.123
>>> R0 = euler_matrix(1, 2, 3)
>>> scale, shear, angles, trans, persp = decompose_matrix(R0)
>>> R1 = euler_matrix(*angles)
>>> np.allclose(R0, R1)
True
MDAnalysis.lib.transformations.euler_from_quaternion(quaternion, axes='sxyz')[source]

Return Euler angles from quaternion for specified axis sequence.

>>> angles = euler_from_quaternion([0.99810947, 0.06146124, 0, 0])
>>> np.allclose(angles, [0.123, 0, 0])
True
MDAnalysis.lib.transformations.projection_from_matrix(matrix, pseudo=False)[source]

Return projection plane and perspective point from projection matrix.

Return values are same as arguments for projection_matrix function: point, normal, direction, perspective, and pseudo.

>>> point = np.random.random(3) - 0.5
>>> normal = np.random.random(3) - 0.5
>>> direct = np.random.random(3) - 0.5
>>> persp = np.random.random(3) - 0.5
>>> P0 = projection_matrix(point, normal)
>>> result = projection_from_matrix(P0)
>>> P1 = projection_matrix(*result)
>>> is_same_transform(P0, P1)
True
>>> P0 = projection_matrix(point, normal, direct)
>>> result = projection_from_matrix(P0)
>>> P1 = projection_matrix(*result)
>>> is_same_transform(P0, P1)
True
>>> P0 = projection_matrix(point, normal, perspective=persp, pseudo=False)
>>> result = projection_from_matrix(P0, pseudo=False)
>>> P1 = projection_matrix(*result)
>>> is_same_transform(P0, P1)
True
>>> P0 = projection_matrix(point, normal, perspective=persp, pseudo=True)
>>> result = projection_from_matrix(P0, pseudo=True)
>>> P1 = projection_matrix(*result)
>>> is_same_transform(P0, P1)
True
MDAnalysis.lib.transformations.quaternion_imag(quaternion)[source]

Return imaginary part of quaternion.

>>> quaternion_imag([3.0, 0.0, 1.0, 2.0])
[0.0, 1.0, 2.0]
MDAnalysis.lib.transformations.quaternion_real(quaternion)[source]

Return real part of quaternion.

>>> quaternion_real([3.0, 0.0, 1.0, 2.0])
3.0
MDAnalysis.lib.transformations.reflection_from_matrix(matrix)[source]

Return mirror plane point and normal vector from reflection matrix.

>>> v0 = np.random.random(3) - 0.5
>>> v1 = np.random.random(3) - 0.5
>>> M0 = reflection_matrix(v0, v1)
>>> point, normal = reflection_from_matrix(M0)
>>> M1 = reflection_matrix(point, normal)
>>> is_same_transform(M0, M1)
True
MDAnalysis.lib.transformations.rotation_from_matrix(matrix)[source]

Return rotation angle and axis from rotation matrix.

>>> angle = (random.random() - 0.5) * (2*math.pi)
>>> direc = np.random.random(3) - 0.5
>>> point = np.random.random(3) - 0.5
>>> R0 = rotation_matrix(angle, direc, point)
>>> angle, direc, point = rotation_from_matrix(R0)
>>> R1 = rotation_matrix(angle, direc, point)
>>> is_same_transform(R0, R1)
True
MDAnalysis.lib.transformations.rotaxis(a, b)[source]

Return the rotation axis to rotate vector a into b.

Parameters

  • a (array_like) – two vectors

  • b (array_like) – two vectors

Returns

c – vector to rotate a into b

Return type

np.ndarray

Note

If a == b this will always return [1, 0, 0]

MDAnalysis.lib.transformations.scale_from_matrix(matrix)[source]

Return scaling factor, origin and direction from scaling matrix.

>>> factor = random.random() * 10 - 5
>>> origin = np.random.random(3) - 0.5
>>> direct = np.random.random(3) - 0.5
>>> S0 = scale_matrix(factor, origin)
>>> factor, origin, direction = scale_from_matrix(S0)
>>> S1 = scale_matrix(factor, origin, direction)
>>> is_same_transform(S0, S1)
True
>>> S0 = scale_matrix(factor, origin, direct)
>>> factor, origin, direction = scale_from_matrix(S0)
>>> S1 = scale_matrix(factor, origin, direction)
>>> is_same_transform(S0, S1)
True
MDAnalysis.lib.transformations.shear_from_matrix(matrix)[source]

Return shear angle, direction and plane from shear matrix.

>>> angle = (random.random() - 0.5) * 4*math.pi
>>> direct = np.random.random(3) - 0.5
>>> point = np.random.random(3) - 0.5
>>> normal = np.cross(direct, np.random.random(3))
>>> S0 = shear_matrix(angle, direct, point, normal)
>>> angle, direct, point, normal = shear_from_matrix(S0)
>>> S1 = shear_matrix(angle, direct, point, normal)
>>> is_same_transform(S0, S1)
True
MDAnalysis.lib.transformations.translation_from_matrix(matrix)[source]

Return translation vector from translation matrix.

>>> v0 = np.random.random(3) - 0.5
>>> v1 = translation_from_matrix(translation_matrix(v0))
>>> np.allclose(v0, v1)
True