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()



The Cult of Software Engineering Academia

This morning I was online checking-out a software architecture conference that I had attended before and was considering to attend this year. I was looking over the keynote speakers and found the bio of one particular speaker to be interesting.

The speaker, a PhD and college professor is a leading authority on modern software development methodologies.  Her focus is guiding organizations to institute the cultural changes required to adopt modern software development methods.

Her bio included a quotable-quote or catch-phrase that went like this: "You can take a man out of the Stone Age, but you can't take the Stone Age out of the man".

I was rocked-back by the arrogance of that statement.  I am sure that she and I would be on the same page regarding software methods and the need for cultural change, but to publicly blast those who are not ready to accept your opinions as "less-evolved" is arrogant and absurd.

Perhaps we are seeing the deficiencies in her culture. The culture of academia. A culture accustomed to feeding bright, fresh, hungry minds. A culture accustomed to unquestionable omnipotence over its audience.

This culture of education has attempted to recast itself into a "culture of change", targeting the software industry. But in the process of recasting itself, it has done little in the way of introspection.

Perhaps the most important skills fostered by working software professionals include communication, collaboration and negotiation. Influence does not come easy and neither does experience. Both have to be earned.

Veteran software professionals have learned to listen as well as to speak, to accept criticism and to consider the ideas and opinions of others. Veteran software professionals would feel negatively about a one-way monologue that targets topics so close to the practice of their art.

I am not suggesting that "change agents" engage in conversation as a way to "handle" an audience, using dialog merely as a soft-skill. Instead I would hope that any discussion would be an honest, open, two-way exchange of ideas.

Of course, this may require the "change agent" to closely engage with her audience. In effect, to lose some of her elite status in order to gain the acceptance of her ideas.

To be a participant instead of a prescriber.

Write Software Requirements For The Love Of It

Yeah sure. Right. Who am I kidding, writing software requirements kind of stinks. But right now, I cannot imagine a better use of my time.

Writing Good Software Requirements is Hard

Writing good requirements is a least as hard as writing good code. One key attribute of a good software requirement is that it must be unambiguous. The English language is notoriously ambiguous. The legal profession, political debate and religious holy wars are based on differing interpretations of the written word. Conversely, software programming languages are designed to be unambiguous because they have to be interpreted or parsed by a machine. Perhaps this is why programmers love to code - the machine detects and points out errors, or better yet, validates our efforts by correctly performing the operations that we have specified.

This brings takes us to another key attribute of a good software requirement, the requirement must be correct. This is obvious but the critical point is that the requirement must be perceived as correct by all of the software system's stakeholders. The point of the SRS is to capture a commonly agreed upon view of how the software system will function. The act of generating the SRS can end up being a process of socialization and team building. (Right now most programmers are saying "Have fun, I will be in the lab writing code".)

Looking at some job postings, they all seem to seek programming language expertise and experience. Bah! Programming is natural, joyous and simple compared to writing requirements. Writing requirements is like being a salmon, battling its way up stream to spawn and die. Well, maybe its not that bad, but its not as gratifying as writing code.

The Greatest Software Requirements Pitfall

Specifying software requirements that includes design and implementation details will lead you and your project to your doom.

It is actually difficult, unnatural, and odd to specify system requirements in way that is bereft of design and implementation. For one thing, as software developers, we are primarily problem solvers and solution creators. We are chomping on the bit to design, and implement. Secondly, we are comfortable talking in terms of the software system itself. Specifying software functions without including implementation details is not easy.

Writing software requirements that includes design and implementation is a trap. Essentially, if the implementation is used to define requirements, the implementation itself must be controlled. This leads to an explosion of requirements documentation. It becomes difficult to make even simple changes to the system without stakeholder consensus and documentation maintenance.

A Software Requirements Fallacy

There is an old software engineering joke, one developer says to the other, "You start coding while I go and get the requirements". At this point all learned software engineers slap their knees and laugh because its funny to think that some developers would start building a system before they knew what the system was required to do.

I believe this is line of thought is incorrect. Software requirements are not necessary to start software system development. Software requirements are paramount when trying to figure out if software development is complete.

Initial development can begin with some basic knowledge of the intended system. Initial development, proof-of-concepts, and prototypes will foster an understanding of the domain, encourage communications and team-building.

Essentially, I am suggesting that the act of software requirements specification does not have to imply the use of the waterfall methodology. I am suggesting that software requirements specification works well with iterative and incremental methods.

Unfortunately I cannot speak to how Agile methodologies factor in software requirements. I would bet that most include some form of functional description document. Perhaps the system is developed in increments that leads to some collective consensus of the intended systems functions.

A Software Requirements Truth

An SRS is simply a checklist of functions that a software system must perform.

Ultimately, the effectiveness and usability of the software system, as perceived by its end-users, is primarily driven by the talent, dedication and expertise of the software development team. This is because the one, unwritten requirement levied on all software systems should be:

The software system shall not stink.