Rotation of object using Euler Matrix in python

Question

I am trying to rotate a roll ( or cylinder) using Euler matrix. For that purpose I use the following function.

def roll( R, zi, zf, Euler):

    # R is the radius of the cylinder
    # t is the angle which is running from 0 to 2*pi
    # zi is the lower z co-ordinate of cylinder
    # zf is the upper z co-ordinate of cylinder
    t = np.arange( 0, 2* np.pi + 0.1, 0.1)
    z = np.array([zi, zf])    
    t, z = np.meshgrid(t, z)
    p, q = t.shape
    r = R* np.ones([p,q], float)
    # polar co-ordinates to Cartesian co-ordinate
    x, y, z = pol2cart(r,t,z)

    # Euler rotation
    rot0 = np.array([x[0,:], y[0,:], z[0,:]])
    rot1 = np.array([x[1,:], y[1,:], z[1,:]])
    # mult is the matrix multiplication
    mat0 = mult( Euler, rot0)
    mat1 = mult( Euler, rot1)
    #
    x[0,:] = mat0[0,:]
    y[0,:] = mat0[1,:]
    z[0,:] = mat0[2,:]
    #
    x[1,:] = mat1[0,:]
    y[1,:] = mat1[1,:]
    z[1,:] = mat1[2,:]
    #
    return x, y, z

the function works well when Euler rotation matrix is Euler = np.array([[1,0,0],[0,1,0],[0,0,1]]) and the inputs for function are x, y, z = roll(1, -2, 2, np.array([[1,0,0],[0,1,0],[0,0,1]]) ). Using ax.plot_surface(x,y,z) I got the following figure.

But when I try to rotate the object by Euler matrix Euler = np.array([[1,0,0],[0,1/np.sqrt(2),-1/np.sqrt(2)],[0,1/np.sqrt(2),1/np.sqrt(2)]]) i got the unexpected result.

Here the rotation is 45 degree which is correct but the shape of object is not proper.

normally, np.dot is the matrix multiplication. why do you use mult, from which module do you get it ? — B. M., Dec 10 '17 at 10:40
I didn't know that. so, I wrote my own function for matrix multiplication @B.M. — Divyang Prajapati, Dec 10 '17 at 11:16
I defined mult for matrix multiplication. As discussed in above comment np.dot can be used instead of mult. — Divyang Prajapati, Dec 10 '17 at 11:45

j-i-l · Accepted Answer · 2017-12-10T16:27:33.723

You were almost there. A few things:

You are actually using cylindrical coordinates not spherical ones. I did not check if numpy has a cyl2cat but this is also not really hard to write yourself:

def cyl2cat(r, theta, z):
    return (r*np.cos(theta), r*np.sin(theta), z)

For the rotation I do not quite understand why you make it in two steps. You can use numpy's ravel to do the rotation of a meshgrid:

# ...
rot = np.dot(Euler,np.array([x.ravel(), y.ravel(), z.ravel()]))

and reshape the rotated coordinates:

x_rot = rot[0,:].reshape(x.shape)
# ...

Putting it together:

import numpy as np

def cyl2cart(r,theta,z):
    return (r*np.cos(theta), r*np.sin(theta), z)

def roll( R, zi, zf, Euler):               
    t = np.arange( 0, 2* np.pi + 0.1, 0.1)          
    z = np.array([zi, zf])                          
    t, z = np.meshgrid(t, z)                        
    p, q = t.shape                                  
    r = R* np.ones([p,q], float)                    
    # cylindrical coordinates to Cartesian coordinate   
    x, y, z = cyl2cart(r,t,z)                       

    # Euler rotation                                
    rot = np.dot(                                                
        Euler,                                            
        np.array([x.ravel(), y.ravel(), z.ravel()]) 
    )                                               
    x_rot = rot[0,:].reshape(x.shape)               
    y_rot = rot[1,:].reshape(y.shape)               
    z_rot = rot[2,:].reshape(z.shape)               
    return x_rot, y_rot, z_rot

Now roll does what you want:

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax=fig.add_subplot(111, projection='3d')
x,y,z=roll(1,-2,2,np.array([[1,0,0],[0,1/np.sqrt(2),-1/np.sqrt(2)],[0,1/np.sqrt(2),1/np.sqrt(2)]]))
ax.plot_surface(x,y,z)
plt.show()

Et voilà:

Note that the aspect ratio of the axes is not the same which is why the cylinder does appear with an elliptic curvature. Getting equal axis in a Axes3D is not straightforward but can be achieved with a workaround by plotting a cubic bounding box (almost copy/pasted from this SO answer):

ax.set_aspect('equal')    
max_range = np.array([x.max()-x.min(), y.max()-y.min(), z.max()-z.min()]).max()
Xb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][0].flatten() + 0.5*(x.max()+x.min())
Yb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][1].flatten() + 0.5*(y.max()+y.min())
Zb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][2].flatten() + 0.5*(z.max()+z.min())
# Comment or uncomment following both lines to test the fake bounding box:
for xb, yb, zb in zip(Xb, Yb, Zb):
   ax.plot([xb], [yb], [zb], 'w')

Simply add this after the ax.plot_surface(... and the cylinder appear with circular curvature.

Rotation of object using Euler Matrix in python

1 Answers1