CODEWITHSUNDEEP
pythonnumpymatplotlibrotationeuler-angles
Divyang Prajapati
Divyang Prajapati

Reputation: 158

Rotation of object using Euler Matrix in python

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. enter image description here

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.

enter image description here

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

Upvotes: 2

Views: 3074

Answers (1)

j-i-l
j-i-l

Reputation: 10957

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à:

enter image description here

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.

Upvotes: 4

Related Questions