Extrinsic/intrinsic convention swapped in `scipy.spatial.transform.Rotation.from_euler()` ?

scipy-from-euler-convention.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "39f88dae-9a36-489b-be06-d6727dc391b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from scipy.spatial.transform import Rotation\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01482be0-aac2-491f-83f7-2914e58b6db2",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Convention of SciPy's `Rotation.from_euler()`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88397ae7-62a2-4409-a67a-fdde9c924f0a",
   "metadata": {},
   "source": [
    "I think the declared convention of SciPy's `Rotation.from_euler()` may be wrong or ambiguously documented."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d5fe781-0660-4447-a17f-9b3b2498bea4",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Basic Rotations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82165157-e8ea-4d6f-97b6-fc6816ebe454",
   "metadata": {},
   "source": [
    "As a preliminary for things to follow, let's' first define the basic rotation matrices about X, Y, and Z. These are matrices which, when pre-multiplying vectors, perform an active, right-handed rotation of the vector by the angle $\\theta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2869a175-07e1-4ec6-a5c0-83c4d329cbb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = symbols('theta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "df82e9be-742b-4b23-828f-5370ffd52aa0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0\\\\0 & \\cos{\\left(\\theta \\right)} & - \\sin{\\left(\\theta \\right)}\\\\0 & \\sin{\\left(\\theta \\right)} & \\cos{\\left(\\theta \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,          0,           0],\n",
       "[0, cos(theta), -sin(theta)],\n",
       "[0, sin(theta),  cos(theta)]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rot_x = Matrix([\n",
    "    [1, 0, 0],\n",
    "    [0, cos(theta), -sin(theta)],\n",
    "    [0, sin(theta), cos(theta)]\n",
    "])\n",
    "rot_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ecf02fe2-3ae5-4fb2-a37e-83d5b2bc3357",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\cos{\\left(\\theta \\right)} & 0 & \\sin{\\left(\\theta \\right)}\\\\0 & 1 & 0\\\\- \\sin{\\left(\\theta \\right)} & 0 & \\cos{\\left(\\theta \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[ cos(theta), 0, sin(theta)],\n",
       "[          0, 1,          0],\n",
       "[-sin(theta), 0, cos(theta)]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rot_y = Matrix(\n",
    "    [\n",
    "        [cos(theta), 0, sin(theta)],\n",
    "        [0, 1, 0],\n",
    "        [-sin(theta), 0, cos(theta)]\n",
    "    ])\n",
    "rot_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5bbadf3b-338c-420b-98b3-b0b7219ce413",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\cos{\\left(\\theta \\right)} & - \\sin{\\left(\\theta \\right)} & 0\\\\\\sin{\\left(\\theta \\right)} & \\cos{\\left(\\theta \\right)} & 0\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[cos(theta), -sin(theta), 0],\n",
       "[sin(theta),  cos(theta), 0],\n",
       "[         0,           0, 1]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rot_z = Matrix(\n",
    "    [\n",
    "        [cos(theta), -sin(theta), 0],\n",
    "        [sin(theta), cos(theta), 0],\n",
    "        [0, 0, 1]\n",
    "    ])\n",
    "rot_z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0ffa32a-a463-47b0-a034-d8f6f137baff",
   "metadata": {},
   "source": [
    "### Verify active right-handed rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fdbbddd-e3a2-4227-ad7b-091b57bb4e9c",
   "metadata": {},
   "source": [
    "Just to be sure, let's verify that these matrices do indeed perform an active, right-handed rotation of a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2f51ed85-442b-4a9c-88ab-e5cf5499a548",
   "metadata": {},
   "outputs": [],
   "source": [
    "R = np.array(rot_z.subs(theta, np.radians(45.0)), dtype=float)\n",
    "vec = np.array([1, 0, 0])\n",
    "res = R @ vec\n",
    "\n",
    "assert np.allclose(res, (np.sqrt(2)/2, np.sqrt(2)/2, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ee52aab-90e7-497d-bacf-3703e12a502f",
   "metadata": {
    "tags": []
   },
   "source": [
    "### An intrinsic rotation (_x-y'-z''_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2240c1da-3045-4f7a-9535-554cc3b8277b",
   "metadata": {},
   "source": [
    "Let's consider the case of a matrix representing three intrinsic rotations, in the following order:\n",
    "- _x_ by $\\gamma$ ($R_1$)\n",
    "- _y'_ by $\\beta$ ($R_2$)\n",
    "- _z''_ by $\\alpha$ ($R_3$)\n",
    "\n",
    "The rotation matrix $R$ can be created by multiplying the basic rotations in the **reverse order**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c925f69a-1792-4a36-b9b1-e6637bea466c",
   "metadata": {},
   "source": [
    "To understand why the order of multiplication is reversed from the transformation order, consider the rotation of a point $p$.\n",
    "\n",
    "We wish to apply the transformations in sequence.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    p' &= R_1 \\cdot p \\\\\n",
    "    p'' &= R_2 \\cdot p' \\\\\n",
    "    p''' &= R_3 \\cdot p''\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "It can be seen that this is equivalent to\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    p''' &= R_3 \\cdot \\left( R_2 \\cdot \\left( R_1 \\cdot p \\right)\\right) \\\\\n",
    "         &= R_3 \\cdot R_2 \\cdot R_1 \\cdot p\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Giving $ R = R_3 \\cdot R_2 \\cdot R_1$, noting that [matrix multiplication is associative](https://en.wikipedia.org/w/index.php?title=Matrix_multiplication&oldid=1119617061#Associativity).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fa7a606-4ac8-443e-b936-daeea834e4ad",
   "metadata": {},
   "source": [
    "Multiplying these matrices to give $R$,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b15cc122-3457-4b1f-88b6-8b59be17326d",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha, beta, gamma = symbols('alpha,beta,gamma')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "aed29345-3b4a-49ca-8a53-114bd6b6f505",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\cos{\\left(\\alpha \\right)} \\cos{\\left(\\beta \\right)} & - \\sin{\\left(\\alpha \\right)} \\cos{\\left(\\gamma \\right)} + \\sin{\\left(\\beta \\right)} \\sin{\\left(\\gamma \\right)} \\cos{\\left(\\alpha \\right)} & \\sin{\\left(\\alpha \\right)} \\sin{\\left(\\gamma \\right)} + \\sin{\\left(\\beta \\right)} \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\gamma \\right)}\\\\\\sin{\\left(\\alpha \\right)} \\cos{\\left(\\beta \\right)} & \\sin{\\left(\\alpha \\right)} \\sin{\\left(\\beta \\right)} \\sin{\\left(\\gamma \\right)} + \\cos{\\left(\\alpha \\right)} \\cos{\\left(\\gamma \\right)} & \\sin{\\left(\\alpha \\right)} \\sin{\\left(\\beta \\right)} \\cos{\\left(\\gamma \\right)} - \\sin{\\left(\\gamma \\right)} \\cos{\\left(\\alpha \\right)}\\\\- \\sin{\\left(\\beta \\right)} & \\sin{\\left(\\gamma \\right)} \\cos{\\left(\\beta \\right)} & \\cos{\\left(\\beta \\right)} \\cos{\\left(\\gamma \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[cos(alpha)*cos(beta), -sin(alpha)*cos(gamma) + sin(beta)*sin(gamma)*cos(alpha), sin(alpha)*sin(gamma) + sin(beta)*cos(alpha)*cos(gamma)],\n",
       "[sin(alpha)*cos(beta),  sin(alpha)*sin(beta)*sin(gamma) + cos(alpha)*cos(gamma), sin(alpha)*sin(beta)*cos(gamma) - sin(gamma)*cos(alpha)],\n",
       "[          -sin(beta),                                     sin(gamma)*cos(beta),                                    cos(beta)*cos(gamma)]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = rot_z.subs(theta, alpha) * rot_y.subs(theta, beta) * rot_x.subs(theta, gamma)\n",
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb459efc-cd9d-4282-b93a-84698212f35e",
   "metadata": {},
   "source": [
    "We see that this is the same as the first matrix in [this section](https://en.wikipedia.org/w/index.php?title=Rotation_matrix&oldid=1122062121#General_rotations) from Wikipedia."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82154b5a-0d4a-41d4-ad82-bf759db98eab",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Comparison with `Rotation.from_euler()`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdb84888-2099-490f-b26a-0190048185b6",
   "metadata": {},
   "source": [
    "First we evaluate the \"truth,\" at some discrete $\\alpha$, $\\beta$, $\\gamma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "54c0ff36-d4c4-431b-af9b-049aa3abf09b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.981060262190407 & -0.0394135507189039 & 0.189650557527834\\\\0.0858316511774313 & 0.966167267147208 & -0.24321541799288\\\\-0.17364817766693 & 0.254887002244179 & 0.951251242564198\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[ 0.981060262190407, -0.0394135507189039, 0.189650557527834],\n",
       "[0.0858316511774313,   0.966167267147208, -0.24321541799288],\n",
       "[ -0.17364817766693,   0.254887002244179, 0.951251242564198]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha_val = np.radians(5)\n",
    "beta_val = np.radians(10)\n",
    "gamma_val = np.radians(15)\n",
    "\n",
    "mat_symbolic_truth = R.subs(alpha, alpha_val).subs(beta, beta_val).subs(gamma, gamma_val)\n",
    "mat_symbolic_truth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d64bc0-e30c-48bc-94a2-1cea04109a48",
   "metadata": {},
   "source": [
    "As a NumPy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "179aad64-908c-4d22-a3ac-8e50d793abdf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.98106026, -0.03941355,  0.18965056],\n",
       "       [ 0.08583165,  0.96616727, -0.24321542],\n",
       "       [-0.17364818,  0.254887  ,  0.95125124]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mat_truth = np.array(mat_symbolic_truth, dtype=float)\n",
    "mat_truth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "453249d1-e051-42ec-ab7f-67ed72035d5e",
   "metadata": {},
   "source": [
    "Now try to build the same matrix using SciPy's `Rotation.from_euler()`.\n",
    "\n",
    "The [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.from_euler.html) defines the input as follows:\n",
    "\n",
    "> `seq`: `string`\n",
    ">\n",
    "> Specifies sequence of axes for rotations. Up to 3 characters belonging to the set {‘X’, ‘Y’, ‘Z’} for intrinsic rotations, or {‘x’, ‘y’, ‘z’} for extrinsic rotations. Extrinsic and intrinsic rotations cannot be mixed in one function call.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bbdb3ab-1d42-4c1a-98b2-c75f0cbb6247",
   "metadata": {},
   "source": [
    "Recall the sequence of transformations from above:\n",
    "\n",
    "> - _x_ by $\\gamma$ ($R_1$)\n",
    "> - _y'_ by $\\beta$ ($R_2$)\n",
    "> - _z''_ by $\\alpha$ ($R_3$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ea7100-37b6-4a79-9c31-91311df70b3f",
   "metadata": {},
   "source": [
    "Thus, the following should reproduce the truth matrix - but it doesn't!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a4d7a906-e039-4b10-af05-caddc710f824",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.98106026 -0.08583165  0.17364818]\n",
      " [ 0.12895841  0.95833311 -0.254887  ]\n",
      " [-0.14453543  0.2724529   0.95125124]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/gt/2y3njwd96374rxc4x5ls706c0000gp/T/ipykernel_14546/2075304721.py:8: UserWarning: \n",
      "BAD: Matrix does not match expected result!\n",
      "  warnings.warn('\\nBAD: Matrix does not match expected result!')\n"
     ]
    }
   ],
   "source": [
    "mat = Rotation.from_euler('XYZ', [gamma_val, beta_val, alpha_val]).as_matrix()\n",
    "print(mat)\n",
    "\n",
    "if np.allclose(mat, mat_truth):\n",
    "    print('\\nOK: Matrix matches expected result!')\n",
    "else:\n",
    "    import warnings\n",
    "    warnings.warn('\\nBAD: Matrix does not match expected result!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47320089-1a2c-4249-90f1-c93bdb9c1d66",
   "metadata": {},
   "source": [
    "But using the lowercase `xyz`, which should specify **an extrinsic** rotation, does!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "72f885f4-ea3d-48a9-b21a-4c65ad1fc07f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.98106026 -0.03941355  0.18965056]\n",
      " [ 0.08583165  0.96616727 -0.24321542]\n",
      " [-0.17364818  0.254887    0.95125124]]\n",
      "\n",
      "OK: Matrix matches expected result!\n"
     ]
    }
   ],
   "source": [
    "mat = Rotation.from_euler('xyz', [gamma_val, beta_val, alpha_val]).as_matrix()\n",
    "print(mat)\n",
    "\n",
    "if np.allclose(mat, mat_truth):\n",
    "    print('\\nOK: Matrix matches expected result!')\n",
    "else:\n",
    "    import warnings\n",
    "    warnings.warn('\\nBAD: Matrix does not match expected result!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96710c80-62e5-49b4-bd1d-943c071e2833",
   "metadata": {},
   "source": [
    "## Verify some other conventions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb8cf258-4f72-4425-89e3-a7cac9b7aa36",
   "metadata": {},
   "source": [
    "### SciPy's `Rotation.as_matrix()` produces a matrix which should be pre-multiplied to perform an active rotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "53b7c003-d54d-4f30-8df5-9d12c8a32f75",
   "metadata": {},
   "outputs": [],
   "source": [
    "rot_mat = Rotation.from_euler('z', np.radians(45)).as_matrix()\n",
    "assert np.allclose(rot_mat @ [1, 0, 0], [np.sqrt(2)/2, np.sqrt(2)/2, 0])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}