Matrix Transformations

Just in case it wasn't obvious yet, I've always been a dunce and never listened in school. Recently I had an issue I knew matrices could solve but I had no idea how those doodads did their thing. Most documentation I found online was either impossible to understand ("draw those squiggly shapes between the brackets and it'll solve itself" or "here's a library, just RTFM," nothing short of cruel jokes,) too scholar (and as boring as the real thing,) or downright wrong (looking at you, "pay-for-more-great-articles" website - you know who you are, you dirty slut kept edging me through the gloryhole in your paywall.) So here's something (I hope) will clear things up a bit. Matrices are simple in concept but a tiny bit confusing when first approached. Despair not m'lady! Uncle Mike learned how to do it five minutes ago so you can too!

Matrices are an easy way to transform (translate, scale, rotate, shear or mirror) 2D and 3D points without using dedicated functions (sort of,) and they can be merged together so that multiple transformations can be done in one step!

I'm no teacher, I don't speak math, but more or less, matrices can be seen as double-entry tables where each row is a function, each column is an argument and each cell contains a value. For each row, you multiply the content of each cell with one of the point's coordinate, then return the sum of those products.

For instance, if we want to process a 3D point, we'd need at least a 3x3 matrix. Let's take point (x, y, z) and throw it into the following matrix to obtain point2 (x2, y2, z2):

X goes in column 1, Y in column 2, Z in column 3, X2 is the result of row 1, Y2 of row 2 and Z2 of row 3... and presto:

x2 = x * 11 + y * 12 + z * 13
y2 = x * 21 + y * 22 + z * 23
z2 = x * 31 + y * 32 + z * 33

Now as you can see, you can put any arbitrary value you want in your matrix, but that would be completely counter-productive. Serious people worked seriously hard on matrices, so let's not be that nonchalant about it. Let's write a matrix to double the coordinates of point (x, y, z):

And the result:

x2 = x * 2 + y * 0 + z * 0
y2 = x * 0 + y * 2 + z * 0
z2 = x * 0 + y * 0 + z * 2

That's it. Matrices are confusing to look at, but really simple stuff in theory. At this point you got 99% of it down, but let's hack into it a little bit. What if we want to translate (move) a point? If you wanted point2 (x2, y2, z2) to be point (x, y, z) moved by (1, 2, 3), you would just do this:

x2 = x + 1
y2 = y + 2
z2 = z + 3

"But Mike!" I hear you say (because now I hear voices too,) "how can you add stuff with matrices? They only multiply things!". Except they don't and that's the trick: they first multiply and then add stuff together. What if we could send a fourth dummy argument to our functions so that we can have our offset unchanged? What if we hacked that extra dummy parameter into our 3D point so it isn't (x, y, z) but rather (x, y, z, w)? What if "w" would always be 1, regardless of the situation? What if we had a 4x4 matrix instead? What if I had cold pizza and cola for breakfast every morning for ten years?

x2 = x * 1 + y * 0 + z * 0 + w * 0
y2 = x * 0 + y * 1 + z * 0 + w * 0
z2 = x * 0 + y * 0 + z * 1 + w * 0
w2 = x * 0 + y * 0 + z * 0 + w * 1

I like to call this dummy coordinate a dirty hack in the matrix, but professionals seem to call this "homogenous coordinates." By the way, that previous matrix is called "identity matrix" and it returns the point you provided without any change whatsoever.

We're going to need a bunch of things. First, we need a type for our 3D points which will include the mysterious "w" dummy:

TYPE vec3
  x AS SINGLE
  y AS SINGLE
  z AS SINGLE
  w AS SINGLE
END TYPE

Then we need a way to store the current matrix. A simple array of 16 single-precision floating-point will do:

DIM m(15) AS SINGLE

And finally we need some routines... let's start with the one that will apply the matrix to the point (separating matrix configuration, merging and application will be super useful later on when we'll apply multiple transformations at once!) and just in case, let's write another routine to initialize/reset the matrix.

'
' Apply matrix
'
SUB mApply (m() AS SINGLE, v AS VEC3)
  DIM t AS VEC3 ' We need a temporary point

  t.x = v.x * m(0)  + v.y * m(1)  + v.z * m(2)  + v.w * m(3)
  t.y = v.x * m(4)  + v.y * m(5)  + v.z * m(6)  + v.w * m(7)
  t.z = v.x * m(8)  + v.y * m(9)  + v.z * m(10) + v.w * m(11)
  t.w = v.x * m(12) + v.y * m(13) + v.z * m(14) + v.w * m(15)
  v = t
END SUB

'
' Set identity
'
SUB mIdentity (m() AS SINGLE)
  m(0)  = 1: m(1)  = 0: m(2)  = 0: m(3)  = 0
  m(4)  = 0: m(5)  = 1: m(6)  = 0: m(7)  = 0
  m(8)  = 0: m(9)  = 0: m(10) = 1: m(11) = 0
  m(12) = 0: m(13) = 0: m(14) = 0: m(15) = 1
END SUB

Alright, now we can start doing stuff or whatever.

Like said previously, you can put any value you want into your matrices, so why not put variables so we can control the outcome better? Let's start with a simple translation matrix using tX, tY and tZ:

In our code, we're going to achieve this via a custom routine that will create a brand new matrix:

'
' Translation matrix, m() must be a 16-entry array!
'
SUB mOffset (m() AS SINGLE, tX AS SINGLE, tY AS SINGLE, tZ AS SINGLE)
  m(0)  = 1: m(1)  = 0: m(2)  = 0: m(3)  = tX
  m(4)  = 0: m(5)  = 1: m(6)  = 0: m(7)  = tY
  m(8)  = 0: m(9)  = 0: m(10) = 1: m(11) = tZ
  m(12) = 0: m(13) = 0: m(14) = 0: m(15) = 1
END SUB

And to test it:

DIM v AS VEC3

v.x = 10: v.y = 15: v.z = 30 ' Coordinates
v.w = 1 ' Always remember to set W to 1

PRINT "Before:"; TAB(10); v.x; v.y; v.z
mOffset m(), -5, 0, -15
mApply m(), v
PRINT " After:"; TAB(10); v.x; v.y; v.z

Another easy one to configure, we've done that earlier. This one will use sX, sY and sZ:

And the custom routine:

'
' Scaling matrix, m() must be a 16-entry array!
'
SUB mScale (m() AS SINGLE, sX AS SINGLE, sY AS SINGLE, sZ AS SINGLE)
  m(0)  = sX: m(1)  = 0:  m(2)  = 0:  m(3)  = 0
  m(4)  = 0:  m(5)  = sY: m(6)  = 0:  m(7)  = 0
  m(8)  = 0:  m(9)  = 0:  m(10) = sZ: m(11) = 0
  m(12) = 0:  m(13) = 0:  m(14) = 0:  m(15) = 1
END SUB

We'll need three matrices for that, one per axis. Here's rotation around X (the X value remains intact:)

'
' X-rotation matrix, m() must be a 16-entry array, a is in degrees
'
SUB mRotateX (m() AS SINGLE, a AS SINGLE)
  DIM s AS SINGLE, c AS SINGLE

  s = SIN(a * 0.01745329#) ' convert angle, get sine
  c = COS(a * 0.01745329#) ' convert angle, get cosine

  m(0)  = 1: m(1)  = 0: m(2)  = 0: m(3)  = 0
  m(4)  = 0: m(5)  = c: m(6)  =-s: m(7)  = 0
  m(8)  = 0: m(9)  = s: m(10) = c: m(11) = 0
  m(12) = 0: m(13) = 0: m(14) = 0: m(15) = 1
END SUB

Rotation around Y (Y is untouched this time:)

'
' Y-rotation matrix, m() must be a 16-entry array, a is in degrees
'
SUB mRotateY (m() AS SINGLE, a AS SINGLE)
  DIM s AS SINGLE, c AS SINGLE

  s = SIN(a * 0.01745329#) ' convert angle, get sine
  c = COS(a * 0.01745329#) ' convert angle, get cosine

  m(0)  = c: m(1)  = 0: m(2)  = s: m(3)  = 0
  m(4)  = 0: m(5)  = 1: m(6)  = 0: m(7)  = 0
  m(8)  =-s: m(9)  = 0: m(10) = c: m(11) = 0
  m(12) = 0: m(13) = 0: m(14) = 0: m(15) = 1
END SUB

And rotation around Z (guess what doesn't happen to Z:)

'
' Z-rotation matrix, m() must be a 16-entry array, a is in degrees
'
SUB mRotateZ (m() AS SINGLE, a AS SINGLE)
  DIM s AS SINGLE, c AS SINGLE

  s = SIN(a * 0.01745329#) ' convert angle, get sine
  c = COS(a * 0.01745329#) ' convert angle, get cosine

  m(0)  = c: m(1)  =-s: m(2)  = 0: m(3)  = 0
  m(4)  = s: m(5)  = c: m(6)  = 0: m(7)  = 0
  m(8)  = 0: m(9)  = 0: m(10) = 1: m(11) = 0
  m(12) = 0: m(13) = 0: m(14) = 0: m(15) = 1
END SUB

Rotation can only be done one axis at a time and will always happen around the origin (0, 0, 0). It is possible to change the center of rotation by moving the points with mOffset (subtract the value of the coordinates of the new center,) rotate, and finally return the points to their original location (add the value back.) Now, what if I told you it's possible to combine several transformations into one matrix so we only need to call mApply once? How neat would that be? Very neat. That's how neat. Let's do that.

This is the part where things get messy. It's not hard, just super confusing to look at. Basically, what's going on is that we're taking a 4x4 matrix and "applying" it to another 4x4 matrix as if it were four (x, y, z, w) points all going in at once. Get ready...

'
' Combine matrix m2 on top of matrix m1 (m1 is modified)
'
SUB mProduct (m1() AS SINGLE, m2() AS SINGLE)
  DIM t(15) AS SINGLE ' Temporary matrix

  t( 0) = m2( 0)*m1( 0) + m2( 1)*m1( 4) + m2( 2)*m1( 8) + m2( 3)*m1(12)
  t( 1) = m2( 0)*m1( 1) + m2( 1)*m1( 5) + m2( 2)*m1( 9) + m2( 3)*m1(13)
  t( 2) = m2( 0)*m1( 2) + m2( 1)*m1( 6) + m2( 2)*m1(10) + m2( 3)*m1(14)
  t( 3) = m2( 0)*m1( 3) + m2( 1)*m1( 7) + m2( 2)*m1(11) + m2( 3)*m1(15)
  t( 4) = m2( 4)*m1( 0) + m2( 5)*m1( 4) + m2( 6)*m1( 8) + m2( 7)*m1(12)
  t( 5) = m2( 4)*m1( 1) + m2( 5)*m1( 5) + m2( 6)*m1( 9) + m2( 7)*m1(13)
  t( 6) = m2( 4)*m1( 2) + m2( 5)*m1( 6) + m2( 6)*m1(10) + m2( 7)*m1(14)
  t( 7) = m2( 4)*m1( 3) + m2( 5)*m1( 7) + m2( 6)*m1(11) + m2( 7)*m1(15)
  t( 8) = m2( 8)*m1( 0) + m2( 9)*m1( 4) + m2(10)*m1( 8) + m2(11)*m1(12)
  t( 9) = m2( 8)*m1( 1) + m2( 9)*m1( 5) + m2(10)*m1( 9) + m2(11)*m1(13)
  t(10) = m2( 8)*m1( 2) + m2( 9)*m1( 6) + m2(10)*m1(10) + m2(11)*m1(14)
  t(11) = m2( 8)*m1( 3) + m2( 9)*m1( 7) + m2(10)*m1(11) + m2(11)*m1(15)
  t(12) = m2(12)*m1( 0) + m2(13)*m1( 4) + m2(14)*m1( 8) + m2(15)*m1(12)
  t(13) = m2(12)*m1( 1) + m2(13)*m1( 5) + m2(14)*m1( 9) + m2(15)*m1(13)
  t(14) = m2(12)*m1( 2) + m2(13)*m1( 6) + m2(14)*m1(10) + m2(15)*m1(14)
  t(15) = m2(12)*m1( 3) + m2(13)*m1( 7) + m2(14)*m1(11) + m2(15)*m1(15)

  ' Replace target matrix
  FOR i% = 0 TO 15
    m1(i%) = t(i%)
  NEXT i%
END SUB

There is one thing you need to know about combining matrices: product of matrices is non-commutative. It means that multiplying m1() by m2() will not provide the same result as multiplying m2() by m1(). This also means that the order in which transformations are made is important (scaling and then translating will not provide the same result as translating and then scaling.)

Now, with all this good stuff, we can operate complex transformations on multiple points by configuring a single matrix! This code includes mProduct at the end of each transformation routine so the matrix is automatically updated:

DEFINT A-Z

DECLARE SUB mIdentity (m() AS SINGLE)
DECLARE SUB mProduct (m1() AS SINGLE, m2() AS SINGLE)
DECLARE SUB mApply (m() AS SINGLE, v AS ANY)
DECLARE SUB mOffset (m() AS SINGLE, tX AS SINGLE, tY AS SINGLE, tZ AS SINGLE)
DECLARE SUB mScale (m() AS SINGLE, sX AS SINGLE, sY AS SINGLE, sZ AS SINGLE)
DECLARE SUB mRotate (m() AS SINGLE, x AS INTEGER, d AS SINGLE)
DECLARE SUB vSet (v AS ANY, x AS SINGLE, y AS SINGLE, z AS SINGLE)

'
' 3D vector
'
TYPE VEC3
  x AS SINGLE
  y AS SINGLE
  z AS SINGLE
  w AS SINGLE
END TYPE

DIM v AS VEC3       ' Point
DIM m(15) AS SINGLE ' Matrix

' Initialize matrix and point
mIdentity m()
vSet v, 10, 15, 30

' Stack multiple transformations
mOffset m(), -5, 0, -15
mScale m(), 2, 2, 2
mRotate m(), 3, 90

' Apply to point
PRINT "  Init:"; TAB(10); v.x; v.y; v.z
mApply m(), v
PRINT " Final:"; TAB(10); v.x; v.y; v.z

SLEEP
END


'
' Apply matrix
'
SUB mApply (m() AS SINGLE, v AS VEC3)
  DIM t AS VEC3

  t.x = v.x * m(0)  + v.y * m(1)  + v.z * m(2)  + v.w * m(3)
  t.y = v.x * m(4)  + v.y * m(5)  + v.z * m(6)  + v.w * m(7)
  t.z = v.x * m(8)  + v.y * m(9)  + v.z * m(10) + v.w * m(11)
  t.w = v.x * m(12) + v.y * m(13) + v.z * m(14) + v.w * m(15)

  v = t
END SUB

'
' Set identity
'
SUB mIdentity (m() AS SINGLE)
  m(0)  = 1: m(1)  = 0: m(2)  = 0: m(3)  = 0
  m(4)  = 0: m(5)  = 1: m(6)  = 0: m(7)  = 0
  m(8)  = 0: m(9)  = 0: m(10) = 1: m(11) = 0
  m(12) = 0: m(13) = 0: m(14) = 0: m(15) = 1
END SUB

'
' Translation matrix
'
SUB mOffset (m() AS SINGLE, tX AS SINGLE, tY AS SINGLE, tZ AS SINGLE)
  DIM mT(15) AS SINGLE

  mT(0)  = 1: mT(1)  = 0: mT(2)  = 0: mT(3)  = tX
  mT(4)  = 0: mT(5)  = 1: mT(6)  = 0: mT(7)  = tY
  mT(8)  = 0: mT(9)  = 0: mT(10) = 1: mT(11) = tZ
  mT(12) = 0: mT(13) = 0: mT(14) = 0: mT(15) = 1

  mProduct m(), mT()
END SUB

'
' Combine matrix m2 on top of matrix m1
'
SUB mProduct (m1() AS SINGLE, m2() AS SINGLE)
  DIM t(15) AS SINGLE

  t(0) = m2(0) * m1(0) + m2(1) * m1(4) + m2(2) * m1(8) + m2(3) * m1(12)
  t(1) = m2(0) * m1(1) + m2(1) * m1(5) + m2(2) * m1(9) + m2(3) * m1(13)
  t(2) = m2(0) * m1(2) + m2(1) * m1(6) + m2(2) * m1(10) + m2(3) * m1(14)
  t(3) = m2(0) * m1(3) + m2(1) * m1(7) + m2(2) * m1(11) + m2(3) * m1(15)
  t(4) = m2(4) * m1(0) + m2(5) * m1(4) + m2(6) * m1(8) + m2(7) * m1(12)
  t(5) = m2(4) * m1(1) + m2(5) * m1(5) + m2(6) * m1(9) + m2(7) * m1(13)
  t(6) = m2(4) * m1(2) + m2(5) * m1(6) + m2(6) * m1(10) + m2(7) * m1(14)
  t(7) = m2(4) * m1(3) + m2(5) * m1(7) + m2(6) * m1(11) + m2(7) * m1(15)
  t(8) = m2(8) * m1(0) + m2(9) * m1(4) + m2(10) * m1(8) + m2(11) * m1(12)
  t(9) = m2(8) * m1(1) + m2(9) * m1(5) + m2(10) * m1(9) + m2(11) * m1(13)
  t(10) = m2(8) * m1(2) + m2(9) * m1(6) + m2(10) * m1(10) + m2(11) * m1(14)
  t(11) = m2(8) * m1(3) + m2(9) * m1(7) + m2(10) * m1(11) + m2(11) * m1(15)
  t(12) = m2(12) * m1(0) + m2(13) * m1(4) + m2(14) * m1(8) + m2(15) * m1(12)
  t(13) = m2(12) * m1(1) + m2(13) * m1(5) + m2(14) * m1(9) + m2(15) * m1(13)
  t(14) = m2(12) * m1(2) + m2(13) * m1(6) + m2(14) * m1(10) + m2(15) * m1(14)
  t(15) = m2(12) * m1(3) + m2(13) * m1(7) + m2(14) * m1(11) + m2(15) * m1(15)

  ' Replace target matrix
  FOR i% = 0 TO 15
    m1(i%) = t(i%)
  NEXT i%
END SUB

'
' Rotation matrix
'
SUB mRotate (m() AS SINGLE, x AS INTEGER, d AS SINGLE)
  DIM mT(15) AS SINGLE, r AS SINGLE, s AS SINGLE, c AS SINGLE

  r = d * .01745329#     ' Convert degrees to radians
  s = SIN(r): c = COS(r) ' Get sine and cosine

  SELECT CASE x
    CASE 1 ' X axis
      mT(0)  = 1: mT(1)  = 0: mT(2)  = 0: mT(3)  = 0
      mT(4)  = 0: mT(5)  = c: mT(6) = -s: mT(7)  = 0
      mT(8)  = 0: mT(9)  = s: mT(10) = c: mT(11) = 0
      mT(12) = 0: mT(13) = 0: mT(14) = 0: mT(15) = 1
    CASE 2 ' Y axis
      mT(0)  = c: mT(1)  = 0: mT(2)  = s: mT(3)  = 0
      mT(4)  = 0: mT(5)  = 1: mT(6)  = 0: mT(7)  = 0
      mT(8) = -s: mT(9)  = 0: mT(10) = c: mT(11) = 0
      mT(12) = 0: mT(13) = 0: mT(14) = 0: mT(15) = 1
    CASE 3 ' Z axis
      mT(0)  = c: mT(1) = -s: mT(2)  = 0: mT(3)  = 0
      mT(4)  = s: mT(5)  = c: mT(6)  = 0: mT(7)  = 0
      mT(8)  = 0: mT(9)  = 0: mT(10) = 1: mT(11) = 0
      mT(12) = 0: mT(13) = 0: mT(14) = 0: mT(15) = 1
    CASE ELSE
      EXIT SUB
  END SELECT

  mProduct m(), mT()
END SUB

'
' Scaling matrix
'
SUB mScale (m() AS SINGLE, sX AS SINGLE, sY AS SINGLE, sZ AS SINGLE)
  DIM mT(15) AS SINGLE

  mT(0) = sX: mT(1)  = 0: mT(2)   = 0: mT(3)  = 0
  mT(4)  = 0: mT(5) = sY: mT(6)   = 0: mT(7)  = 0
  mT(8)  = 0: mT(9)  = 0: mT(10) = sZ: mT(11) = 0
  mT(12) = 0: mT(13) = 0: mT(14)  = 0: mT(15) = 1

  mProduct m(), mT()
END SUB

'
' Set point
'
SUB vSet (v AS VEC3, x AS SINGLE, y AS SINGLE, z AS SINGLE)
  v.x = x: v.y = y: v.z = z: v.w = 1
END SUB