Python, NumPy and matplotlib - Visualize Matrix Rotations

I am using Python, NumPy and matplotlib to experiment with, and visualize 3-D graphics techniques.  The first step was figuring out how to generate polar plots with matplotlib and to determine the domain and range of spherical to Cartesian coordinate conversions.

The next step is to rotate a point in space using matrices, and then to explore the concept of gimbal lock.  I will probably tackle gimbal lock in a future posting.  Before delving into 3-D rotations, I needed to figure out how to extend my plotting knowledge.

I cleaned up the code from the previous posting, moving the spherical to Cartesian coordinates conversion code to a separate module - sphcoords. sphcoords also includes a function that takes spherical coordinates and returns a 1x3 NumPy matrix that represents a point in space.

import numpy as np
from numpy import matrix

def x_cart(incl_sph_coords, azim_sph_coords):
   return np.sin(incl_sph_coords) * np.cos(azim_sph_coords)

def y_cart(incl_sph_coords, azim_sph_coords):
   return np.sin(incl_sph_coords) * np.sin(azim_sph_coords)

def z_cart(incl_sph_coords):
   return np.cos(incl_sph_coords)

def sph2cart(incl_sph_coord, azim_sph_coord):
   return np.matrix([
   x_cart(incl_sph_coord, azim_sph_coord),
   y_cart(incl_sph_coord, azim_sph_coord),
   z_cart(incl_sph_coord)])

Below is an example of sphcoords function, sph2cart which generates a 1x3 matrix containing Cartesian coordinates. In this case, spherical coordinates of 30 degrees inclination and 180 degrees azimuth are converted.

poi = sphcoords.sph2cart(np.pi/6, -np.pi+(10*azi_incr))

By using matplotlib's subplot capability along with spherical to Cartesian coordinate conversions, I was able to depict a 3-D spherical model of space.  A single red "point of interest" will be rotated using rotational matrices.

The initial position of the red point is arbitrary.  The red point will be initially rotated about the z axis.  The rotated position of the red point will remain along the point's original circle of inclination.  This is what I was hoping  to illustrate with the cyan-colored annotations.

Before going further, let me apologize for my inability to generate circles. The elliptical, squashed figures in my plots will be circular after I figure out how to set fixed aspect ratios.

The Initial Position of the Red Point of Interest

Here are the rotational matrices that were implemented in module rotmat using the NumPy Matrix. 

import numpy as np
from numpy import matrix

def x_rot(rads):
  return np.matrix([
  [1,0,0], 
  [0, np.cos(rads), -np.sin(rads)],
  [0, np.sin(rads), np.cos(rads)]])

def y_rot(rads):
  return np.matrix([
  [np.cos(rads), 0, np.sin(rads)],
  [0, 1, 0], 
  [-np.sin(rads), 0, np.cos(rads)]])

def z_rot(rads):
  return np.matrix([
  [np.cos(rads), -np.sin(rads), 0],
  [np.sin(rads), np.cos(rads), 0],
  [0,0,1]]) 

The 1x3 matrix containing the point of interest's (POI) Cartesian coordinates is rotated about the z axis thusly:

poi = poi*rotmat.z_rot(-5*azi_incr)

The plot below shows the POI rotated 90 degrees about the z axis.  The cyan annotations show the range of the next rotation which is about the x axis.  The POI will be rotated 30 degrees about the x axis.  The most interesting update will be on the upper right plot that will show the the POI jumping from the yellow 30 degrees circle to the black circle at 60 degrees.

After Rotating the POI about the Z axis.

After Rotating the POI about the X axis.

 

The final position of the POI after rotating about each Euler Axis.

 


import numpy as np
import matplotlib.pylab as plt
from numpy import matrix
import rotmat
import sphcoords


def generate_plot(inclin, azi, poi):
# plot x-y plane
plt.subplot(221)
plt.plot(
   sphcoords.x_cart(inclin[0],azi[0]),
   sphcoords.y_cart(inclin[0],azi[0]), "g-o",
   sphcoords.x_cart(inclin[1],azi[1]),
   sphcoords.y_cart(inclin[1],azi[1]), "y-o",
   sphcoords.x_cart(inclin[2],azi[2]),
   sphcoords.y_cart(inclin[2],azi[2]), "y-o",
   sphcoords.x_cart(inclin[3],azi[3]),
   sphcoords.y_cart(inclin[3],azi[3]), "k-o",
   sphcoords.x_cart(inclin[4],azi[4]),
   sphcoords.y_cart(inclin[4],azi[4]), "k-o",
   sphcoords.x_cart(inclin[5],azi[5]),
   sphcoords.y_cart(inclin[5],azi[5]), "b-o",
   sphcoords.x_cart(inclin[6],azi[6]),
   sphcoords.y_cart(inclin[6],azi[6]), "b-o",
   poi[0,0], poi[0,1], "r-o")
plt.grid(True)
plt.ylabel('Y Axis')
plt.xlabel('X Axis')

# plot y-z plane
plt.subplot(222)
plt.plot(
   sphcoords.z_cart(inclin[0]),
   sphcoords.y_cart(inclin[0],azi[0]), "g-o",
   sphcoords.z_cart(inclin[1]),
   sphcoords.y_cart(inclin[1],azi[1]), "y-o",
   sphcoords.z_cart(inclin[2]),
   sphcoords.y_cart(inclin[2],azi[2]), "y-o",
   sphcoords.z_cart(inclin[3]),
   sphcoords.y_cart(inclin[3],azi[3]), "k-o",
   sphcoords.z_cart(inclin[4]),
   sphcoords.y_cart(inclin[4],azi[4]), "k-o",
   sphcoords.z_cart(inclin[5]),
   sphcoords.y_cart(inclin[5],azi[6]), "b-o",
   sphcoords.z_cart(inclin[6]),
   sphcoords.y_cart(inclin[6],azi[6]), "b-o",
   poi[0,2], poi[0,1], "r-o")
plt.xlabel('Z Axis')
plt.grid(True)

#plot x-z plane
plt.subplot(223)
plt.plot(
   sphcoords.x_cart(inclin[0],azi[0]),
   sphcoords.z_cart(inclin[0]),"g-o",
   sphcoords.x_cart(inclin[1],azi[1]),
   sphcoords.z_cart(inclin[1]),"y-o",
   sphcoords.x_cart(inclin[2],azi[2]),
   sphcoords.z_cart(inclin[2]),"y-o",
   sphcoords.x_cart(inclin[3],azi[3]),
   sphcoords.z_cart(inclin[3]),"k-o",
   sphcoords.x_cart(inclin[4],azi[4]),
   sphcoords.z_cart(inclin[4]),"k-o",
   sphcoords.x_cart(inclin[5],azi[5]),
   sphcoords.z_cart(inclin[5]),"b-o",
   sphcoords.x_cart(inclin[6],azi[6]),
   sphcoords.z_cart(inclin[6]),"b-o",
   poi[0,0], poi[0,2], "r-o")
plt.ylabel('Z Axis')
plt.grid(True)
plt.show()

if __name__ == '__main__':
   inclin = []
   azi = []

   inclin_points = range(0,20)
   azi_incr = (2 * np.pi)/float(20)
   azi_points = np.arange(-np.pi, np.pi, azi_incr)

   #90 deg inclination
   inclin.append([np.pi/2 for i in inclin_points])
   #30,150 deg inclinations
   inclin.append([np.pi/6 for i in inclin_points])
   inclin.append([np.pi-np.pi/6 for i in inclin_points])
   #60,120 deg inclinations
   inclin.append([np.pi/3 for i in inclin_points])
   inclin.append([np.pi-np.pi/3 for i in inclin_points])
   # poles
   inclin.append([0])
   inclin.append([np.pi])

   azi.append(azi_points)
   azi.append(azi_points)
   azi.append(azi_points)
   azi.append(azi_points)
   azi.append(azi_points)
   azi.append([np.pi])
   azi.append([np.pi])

   poi = sphcoords.sph2cart(np.pi/6, -np.pi+(10*azi_incr))
   poi = poi*rotmat.z_rot(-5*azi_incr)
   poi = poi*rotmat.x_rot(np.pi/6)
   poi = poi*rotmat.y_rot(np.pi/2)
   generate_plot(inclin, azi, poi)

Spherical to Cartesian Coords - Thanks Python!

Game and graphics programming requires the modeling of three-dimensional geometries in software.  This stuff is relatively simple in theory but easy to screw up and difficult to debug. 

The math used in 3-D geometrical programming has its limitations and special cases.  One example is the gimbal lock problem that occurs when rotating a body using matrices and Euler angles.

Another problem that I recently stumbled across occurred when trying to find the intersection point of a line and plane.  If the line is pointing away from the plane, the intersection point ends up being behind the line's start point.  Essentially, if we shoot an arrow directly away from a target, it will still hit the target, but with the back of the arrow.

The problem I describe above with the intersection point, line and plane is not really a problem at all.  The special case is easily detected while computing the result and the code can respond as required.

The real issue is that we must have a complete understanding of the mathematical techniques we use in our software.

There are two ways to understand math, one is to use analysis, the other is intuition/visualization.  Analysis is the stronger of the two methods and includes things such as proofs.  Intuition, visualization and experimentation is the approach that I take.

My latest interest is converting spherical coordinates into Cartesian coordinates.  I can easily plug and chug using the standard formulas, but after learning a few hard lessons, I know I had better gain some understanding.

When it comes to three-dimensional programming, it is not as easy to visualize 3-D space as one would imagine.  I needed some help to visualize this stuff.

Locations on earth are identified using spherical coordinates, using the intersection of two angles, an angle of inclination (latitude) and an angle of azimuth (longitude).  In my investigation, I would like to consider inclination angles from 0 to 180 degrees, azimuths from 0 to 360 degrees and a spherical radius of one.

The plot simplifies the problem by reducing three-dimensions to a two-dimensional display with circles of common inclination.  This is essentially a polar-plot. Each circle of common inclination (latitude) represents two inclinations, the inclination (northern hemisphere) and 180 - inclination (southern hemisphere).

As a two-dimensional polar plot, the problem is solved using the azimuth angle, and the distance from the origin to the proper circle of common inclination.

It appears from this work, that spherical to cartesian coordinate conversions is a simple and safe process (with the exception of when the inclination angle is 0 or 180 degrees, of course).


import numpy as np
import matplotlib.pylab as plt

def x(inclin, azi):
""" compute the x cartesian coord """
   return np.sin(inclin) * np.cos(azi)

def y(inclin, azi):
""" compute the y cartesian coord """
   return np.sin(inclin) * np.sin(azi)

def z(inclin):
""" compute the z cartesian coord """
   return np.cos(inclin)


inclin = []
inclin.append([np.pi/2 for i in range(0,20)])
inclin.append([np.pi/3 for i in range(0,20)])
inclin.append([np.pi/4 for i in range(0,20)])
inclin.append([np.pi/6 for i in range(0,20)])
inclin.append([np.pi/8 for i in range(0,20)])
inclin.append([np.pi/12 for i in range(0,20)])
inclin.append([np.pi/24 for i in range(0,20)])
inclin.append([0 for i in range(0,20)])
azi = np.arange(-np.pi, np.pi, (2 * np.pi)/float(20))

plt.plot(
   x(inclin[0],azi), y(inclin[0],azi), "g-o",
   x(inclin[1],azi), y(inclin[1],azi), "r-o",
   x(inclin[2],azi), y(inclin[2],azi), "c-o",
   x(inclin[3],azi), y(inclin[3],azi), "m-o",
   x(inclin[4],azi), y(inclin[4],azi), "y-o",
   x(inclin[5],azi), y(inclin[5],azi), "k-o",
   x(inclin[6],azi), y(inclin[6],azi), "b-o",
   x(inclin[7],azi), y(inclin[7],azi), "w-o")
plt.legend(
   ('90', '60, 120', '45, 135',
   '30, 150', '23, 157',
   '15, 165', '8, 172', '0, 180'))
plt.grid(True)
plt.ylabel('Y Cartesian Coord')
plt.title('Cartesian from Spherical')
plt.xlabel('X Cartesian Coord')
plt.show()