0

Objective: Apply an Euler rotation to an ellipsoid, then plot it using matplotlib and mplot3d.

I found a function which applies an Euler rotation to a vector or array of vectors:

import numpy as np
from scipy.linalg import expm

def rot_euler(v, xyz):
''' Rotate vector v (or array of vectors) by the euler angles xyz '''
# https://stackoverflow.com/questions/6802577/python-rotation-of-3d-vector
for theta, axis in zip(xyz, np.eye(3)):
    v = np.dot(np.array(v), expm(np.cross(np.eye(3), axis*-theta)))
return v

However, the ellipsoid I need to rotate is represented as a set of three coordinate matrices (in order to be plotable using ax.plot_surface()):

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import pi,sin,cos

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.set_aspect('equal','box')
ax.set_xlim3d(-1,1)
ax.set_ylim3d(-1,1)
ax.set_zlim3d(-1,1)
ax.view_init(90,90)

ellipseSteps= 100
diamCoef = 500
widthCoef = 25
coefs = (widthCoef, diamCoef, diamCoef)  # Coefficients in a0/c x**2 + a1/c y**2 + a2/c z**2 = 1 
# Radii corresponding to the coefficients:
rx, ry, rz = 1/np.sqrt(coefs)

# Set of all spherical angles:
u = np.linspace(0, 2 * pi, ellipseSteps)
v = np.linspace(0, pi, ellipseSteps)

# Cartesian coordinates that correspond to the spherical angles:
# (this is the equation of an ellipsoid):
ex = rx * np.outer(cos(u), sin(v))
ey = ry * np.outer(sin(u), sin(v))
ez = rz * np.outer(np.ones_like(u), cos(v))

# Plot:
ax.plot_surface(ex, ey, ez,  rstride=4, cstride=4, color='blue')

plt.show()

How can I apply the Euler rotation to this object? I thought of just converting the object from three coordinate matrices to vectors and then feeding that into the existing formula, but it occurs to me that doing so might be computationally inefficient... which makes me wonder if the rotation function can be modified to work on the coordinate matrices?

I realize it's probably a rather trivial ask, but it's been years since I last did any linear algebra, and would very much appreciate the advice of an expert here.

Thanks!

TraxusIV
  • 1,388
  • 1
  • 19
  • 34

1 Answers1

2

While I don't think you can do much better than applying the rotation point-wise there still is room for significant economy.

(1) Using the matrix exponential to compute simple rotation matrices is ridiculously wasteful. Much better to use the scalar exponential or sine and cosine

(2) To a lesser extent the same applies to using the cross product for shuffling. Indexing is preferable here.

(3) The order of matrix multiplication matters. When bulk rotating more than three vectors, the left most multiplication should be done last.

Together these measures speed up the computation by a factor of six:

(insert before last two lines of original script)

from scipy.linalg import block_diag
from timeit import timeit

def rot_euler_better(v, xyz):
    TD = np.multiply.outer(np.exp(1j * np.asanyarray(xyz)), [[1], [1j]]).view(float)
    x, y, z = (block_diag(1, TD[i])[np.ix_(*2*(np.arange(-i, 3-i),))] for i in range(3))
    return v @ (x @ y @ z)

# example
xyz = np.pi * np.array((1/6, -2/3, 1/4))

print("Same result:",
      np.allclose(rot_euler(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz),
                  rot_euler_better(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz))

print("OP:       ", timeit(lambda: rot_euler(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz), number=1000), "ms")
print("optimized:", timeit(lambda: rot_euler_better(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz), number=1000), "ms")

ex, ey, ez = map(np.reshape, rot_euler_better(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz).T, map(np.shape, (ex, ey, ez)))

Output:

Same result: True
OP:        2.1019406360574067 ms
optimized: 0.3485010238364339 ms
Paul Panzer
  • 47,318
  • 2
  • 37
  • 82