Cogs and Levers A blog full of technical stuff

numpy

numpy is an excellent library for dealing with array-based math problems. According to the website:

NumPy is the fundamental package for scientific computing with Python

In today’s post, I’ll go through some very basic operations of the library while piecing together the transformation segment that you might see in the back-end of a 3D rendering pipeline. We’ll focus on some vector and matrix math, but most importantly we’ll look at the creation of specialised matricies and matrix multiplication.

Vectors

All points in our 3D world is going to be re-presented by a 4D Homogenous co-ordinate. numpy provide the perfect abstraction here with the array function.

p = np.array([1, 0, 0, 1])

There are many other construction functions available for this single dimension value series that you can find in the documentation.

Once we’ve created the array, we can start to perform arithmetic on it.

>>> np.array([1, 0, 0, 1]) + np.array([0, 1, 0, 1])
array([1, 1, 0, 2])

>>> np.linalg.norm(np.array([1, 0, 0, 1]))
1.4142135623730951

>>> np.array([1, 0, 0, 1]).dot(np.array([0, 1, 0, 1]))
1

Matrices

Graduating from the vector, we now need a tabular representation for our linear algebraic operations that we’ll perform. To do this, we’ll use a matrix which numpy has a direct analog for.

First of all, every math library needs a shortcut to the identity matrix:

>>> np.identity(4)
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

Getting back into our 3D engine, we’re going to need a way to translate, scale and rotate the array points.

Translation

A translation matrix in 3-space using a 4x4 matrix looks like this:

[ 1 0 0 tx ]
[ 0 1 0 ty ]
[ 0 0 1 tz ]
[ 0 0 0 1  ]

So, when we apply this matrix to a column vector; we’ll end up with a simple addition of each of the components:

[ 1 0 0 tx ] [ vx ]   [ tx + vx ]
[ 0 1 0 ty ] [ vy ]   [ ty + vy ]
[ 0 0 1 tz ] [ vz ] = [ tz + vz ]
[ 0 0 0 1  ] [ 1  ]   [    1    ]

In python using numpy, this is greatly simplified:

def create_translation(x, y, z):
  return np.matrix([
      [1, 0, 0, x], 
      [0, 1, 0, y], 
      [0, 0, 1, z], 
      [0, 0, 0, 1]
  ])

With this function in hand, we can now start to move our vectors around in space:

# v sits at the origin
v = np.array([0, 0, 0, 1])

# translation matrix to move vectors along
# the x-axis by 1
t = create_translation(1, 0, 0)

# move v
t.dot(v)

With the output of this function:

matrix([[1, 0, 0, 1]])

Scale

We scale points geometry to adjust proportionality of objects. A scale matrix will do just this on any point around the origin 0, 0, 0, 1.

[ sx 0  0  0 ]
[ 0  sy 0  0 ]
[ 0  0  sz 0 ]
[ 0  0  0  1 ]

Applying a matrix of this shape to a column vector, we now get a multiplication effect:

[ sx 0  0  0 ] [ vx ]   [ sx * vx ]
[ 0  sy 0  0 ] [ vy ]   [ sy * vy ]
[ 0  0  sz 0 ] [ vz ] = [ sz * vz ]
[ 0  0  0  1 ] [ 1  ]   [    1    ]

The python for this process is pretty simple, again:

def create_scale(x, y, z):
  return np.matrix([
    [x, 0, 0, 0], 
    [0, y, 0, 0], 
    [0, 0, z, 0], 
    [0, 0, 0, 1]
  ])

Rotation

We’ll treat rotation on each axis separately.

x-axis rotation performs spherical movement between the y and z axis:

[  1   0     0    0  ]
[  0  cos∅ -sin∅  0  ]
[  0  sin∅  cos∅  0  ]
[  0   0     0    1  ]

y-axis rotation performs spherical movement between the x and z axis:

[  cos∅   0   sin∅  0 ]
[   0     1    0    0 ]
[ -sin∅   0   cos∅  0 ]
[   0     0    0    1 ]

z-axis rotation performs spherical movement between the x and y axis:

[  cos∅ -sin∅   0    0 ]
[  sin∅  cos∅   0    0 ]
[   0     0     1    0 ]
[   0     0     0    1 ]

These translate really well into python code, too!

def create_rotate_x(theta):
  ct = math.cos(theta)
  st = math.sin(theta)
  
  return np.matrix([
    [1, 0, 0, 0], 
    [0, ct, -st, 0], 
    [0, st, ct, 0], 
    [0, 0, 0, 1]
  ])

def create_rotate_y(theta):
  ct = math.cos(theta)
  st = math.sin(theta)
  
  return np.matrix([
    [ct, 0, st, 0], 
    [0, 1, 0, 0], 
    [-st, 0, ct, 0], 
    [0, 0, 0, 1]
  ])

def create_rotate_z(theta):
  ct = math.cos(theta)
  st = math.sin(theta)
  
  return np.matrix([
    [ct, -st, 0, 0], 
    [st, ct, 0, 0], 
    [0, 0, 1, 0], 
    [0, 0, 0, 1]
  ])

Armed with all of these functions, we have a basic (local) co-ordinate transform pipeline. We can start to make a little more sense out of these mathematical constructs by binding the results to variables that are named by their side effects.

# the cartesian plane moves positive, to the right
move_right_by_2 = create_translation(2, 0, 0)

# make any point twice the distance away from the origin
make_twice_bigger = create_scale(2, 2, 2)

# π radians is 180° 
turn_the_other_way = create_rotate_x(math.pi)

What is interesting about the matricies that we’ve just created, is that they’re reusable for as many points as we want. We can uniformly transform a group of vectors with the same matrix.

So, if I wanted to not only move right by 2, but also make twice bigger and turn the other way; I could multiply these matricies together to form a new transformation:

do_all_the_things = (turn_the_other_way * make_twice_bigger * move_right_by_2)

Now we can use do_all_the_things to transform our vector objects to perform all of these transformations at once:

>>> do_all_the_things.dot(np.array([0, 0, 1, 1]))
matrix([[  4.00000000e+00,  -2.44929360e-16,  -2.00000000e+00,
           1.00000000e+00]])

We did start, pointing towards the camera 0, 0, 1 and then after moving to the right by 2 and doubling our size, we’re now at 4 on the x-axis. We turn around about the x axis and we’re now facing away from the camera; twice the distance (-2 on the z-axis). Note that the y-axis has some epsilon garbage after our multiplies.

Well, it was a quick tour; but this has been numpy in (a little bit of) action.