[Numpy-discussion] composing Euler rotation matrices

[Robert McLeod]

Instead of trying to decipher what someone wrote on a Wikipedia, why don't
you look at a working piece of source code?

e.g.

https://github.com/3dem/relion/blob/master/src/euler.cpp

Robert


On Wed, Feb 1, 2017 at 4:27 AM, Seb <spluque at gmail.com> wrote:

> On Tue, 31 Jan 2017 21:23:55 -0500,
> Joseph Fox-Rabinovitz <jfoxrabinovitz at gmail.com> wrote:
>
> > Could you show what you are doing to get the statement "However, I
> > cannot reproduce this matrix via composition; i.e. by multiplying the
> > underlying rotation matrices.". I would guess something involving the
> > `*` operator instead of `@`, but guessing probably won't help you
> > solve your issue.
>
> Sure, although composition is not something I can take credit for, as
> it's a well-described operation for generating linear transformations.
> It is the matrix multiplication of two or more transformation matrices.
> In the case of Euler transformations, it's matrices specifying rotations
> around 3 orthogonal axes by 3 given angles.  I'm using `numpy.dot' to
> perform matrix multiplication on 2D arrays representing matrices.
>
> However, it's not obvious from the link I provided what particular
> rotation matrices are multiplied and in what order (i.e. what
> composition) is used to arrive at the Z1Y2X3 rotation matrix shown.
> Perhaps I'm not understanding the conventions used therein.  This is one
> of my attempts at reproducing that rotation matrix via composition:
>
> ---<--------------------cut here---------------start------
> ------------->---
> import numpy as np
>
> angles = np.radians(np.array([30, 20, 10]))
>
> def z1y2x3(alpha, beta, gamma):
>     """Z1Y2X3 rotation matrix given Euler angles"""
>     return np.array([[np.cos(alpha) * np.cos(beta),
>                       np.cos(alpha) * np.sin(beta) * np.sin(gamma) -
>                       np.cos(gamma) * np.sin(alpha),
>                       np.sin(alpha) * np.sin(gamma) +
>                       np.cos(alpha) * np.cos(gamma) * np.sin(beta)],
>                      [np.cos(beta) * np.sin(alpha),
>                       np.cos(alpha) * np.cos(gamma) +
>                       np.sin(alpha) * np.sin(beta) * np.sin(gamma),
>                       np.cos(gamma) * np.sin(alpha) * np.sin(beta) -
>                       np.cos(alpha) * np.sin(gamma)],
>                      [-np.sin(beta), np.cos(beta) * np.sin(gamma),
>                       np.cos(beta) * np.cos(gamma)]])
>
> euler_mat = z1y2x3(angles[0], angles[1], angles[2])
>
> ## Now via composition
>
> def rotation_matrix(theta, axis, active=False):
>     """Generate rotation matrix for a given axis
>
>     Parameters
>     ----------
>
>     theta: numeric, optional
>         The angle (degrees) by which to perform the rotation.  Default is
>         0, which means return the coordinates of the vector in the rotated
>         coordinate system, when rotate_vectors=False.
>     axis: int, optional
>         Axis around which to perform the rotation (x=0; y=1; z=2)
>     active: bool, optional
>         Whether to return active transformation matrix.
>
>     Returns
>     -------
>     numpy.ndarray
>     3x3 rotation matrix
>     """
>     theta = np.radians(theta)
>     if axis == 0:
>         R_theta = np.array([[1, 0, 0],
>                             [0, np.cos(theta), -np.sin(theta)],
>                             [0, np.sin(theta), np.cos(theta)]])
>     elif axis == 1:
>         R_theta = np.array([[np.cos(theta), 0, np.sin(theta)],
>                             [0, 1, 0],
>                             [-np.sin(theta), 0, np.cos(theta)]])
>     else:
>         R_theta = np.array([[np.cos(theta), -np.sin(theta), 0],
>                             [np.sin(theta), np.cos(theta), 0],
>                             [0, 0, 1]])
>     if active:
>         R_theta = np.transpose(R_theta)
>     return R_theta
>
> ## The rotations are given as active
> xmat = rotation_matrix(angles[2], 0, active=True)
> ymat = rotation_matrix(angles[1], 1, active=True)
> zmat = rotation_matrix(angles[0], 2, active=True)
> ## The operation seems to imply this composition
> euler_comp_mat = np.dot(xmat, np.dot(ymat, zmat))
> ---<--------------------cut here---------------end--------
> ------------->---
>
> I believe the matrices `euler_mat' and `euler_comp_mat' should be the
> same, but they aren't, so it's unclear to me what particular composition
> is meant to produce the matrix specified by this Z1Y2X3 transformation.
> What am I missing?
>
> --
> Seb
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>



-- 
Robert McLeod, Ph.D.
Center for Cellular Imaging and Nano Analytics (C-CINA)
Biozentrum der Universität Basel
Mattenstrasse 26, 4058 Basel
Work: +41.061.387.3225
robert.mcleod at unibas.ch
robert.mcleod at bsse.ethz.ch <robert.mcleod at ethz.ch>
robbmcleod at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20170201/7a3d743c/attachment.html>