TRTC – Matrices – Chapter 3

Remembering the times back when I was in the high school. I wish that some one would have had the wisdom to motivate me by telling that matrices are crucial when doing 3D stuff on computers. It might have helped to keep the motivation up what comes to studying math 🙂

So I’m not going to focus on how to implement the matrix datastrcutures. Like in the previous post I will use the simd framework which have all the building blocks that are needed. There is one twist between the book and the simd. The book is talking matrices as row based. So when indexing matrix in the book the number of rows comes first and the number of columns comes second.

Well the simd matrices works opposite. First comes the columns and then the rows. So to create matrix like is:

/*
|  1   |   2   |   3   |   4   |
|  5.5 |   6.5 |   7.5 |   8.5 |
|  9   |  10   |  11   |  12   |
| 13.5 |  14.5 |  15.5 |  16.5 |
*/

the book uses notation like M21 to index the value at third row on a second column which is 10. (Rows and columns indexing starts at zero). In my code the the same is equal to M12 where the column comes first and then the row. I’m not actually sure if I can use these simd data structures like this just by taking this difference in account.

By the way there is over 17 test cases in this chapter so I’m not going to put them all here. They can be easily found in the source code. I concentrate on those points that I think were challenging for me or where I think I found some easy or maybe clean way to implement some part. So lets dive in.

Creating a Matrix

As mentioned above the simd matrices are column based. This next test case is an example of creating and testing 2×2 matrix but the same principles are used with 3×3 and 4×4 matices.

/**
 Chapter 3.
 2nd test case.
 2x2 matrix.
 */
func testRGMatrix2x2Creation() {
    // Given
    // This is the matrix we are trying to create and test that the creation sucseeds.
    // |  -3  |   5  |
    // |   1  |  -2  |

    // When
    // Since our matrix is column based we create 4 RGTuple variables using values from the columns above.
    let m2x2 = RGMatrix2x2(simd_float2(x: -3, y: 1), simd_float2(x: 5, y: -2))

    
    // Then
    XCTAssertEqual(m2x2.columns.0.x == -3, true)
    XCTAssertEqual(m2x2.columns.0.y ==  1, true)
    XCTAssertEqual(m2x2.columns.1.x ==  5, true)
    XCTAssertEqual(m2x2.columns.1.y == -2, true)
}

So basically in the case above what is needed is two tuples with 2 values in them. For matrices of sizes 2×2 and 3×3 we use simd types float2 and float3 to create the columns in our case. For 4×4 matrices we can use the RGTupe as well. It’s the alias for float4.

Multiply Matrices

Multiplication is one of the key elements when dealing with 3D graphics in computers. By using multiplication we could “join” together different transformations. There are transformations like scaling, rotation, etc. and each of them could be represented by a matrix.

Basically when we have a some vector we could multiply that with a matrix that represents for example the scaling matrix. After the multiplication we will have the scaled vector as our result. Then we could multiply that with another matrix representing maybe the rotation and again we get a new vector that is also a result of that rotation.

By doing things this way works but there is easier solution. We could multiply all those transformation matrices together and then multiply our vector with that one matrix to represent the same effect as doing all the steps separately. This is more efficient and shows the usefulness of matrix multiplication. I think this is really amazing thing with matrices. It’s like grouping multiple answers together and then answer to all questions an once 🙂

For multiplication I will use the building functionality. Just create two matrices and to multiply use * like with scalars. The code is very simple.

In this test case we also need to test if the result of the multiplication is same as the matrix given in the book. For the testing I use the simd simd_equal() function. In this test case this works but there will be case a where that function is NOT OK. We will come to that later.

The Identity Matrix

The identity matrix is like the number 1, but for matrices.

Jamis Buck

To generate identity matrix I can use the build in functionality form simd or to make it a bit more convenient I created a extension to the RGMatrix4x4. Here are both techniques. The latter uses internal the first technique mentioned.

// To use the buildin init method to create an identity matrix.
let M = RGMatrix4x4.init(diagonal: RGTuple(repeating: 1))

// OR to create and use the extension to make it even simpler.

public extension RGMatrix4x4 {
    // Static method to create identity matrix.
    internal static func identity() -> RGMatrix4x4 {
        return RGMatrix4x4.init(diagonal: RGTuple(repeating: 1))
    }
}

// And to use the identity matrix.
let M = RGMatrix4x4.identity()

Transposing Matrices

For transposing matrix I will use also the simd build in functionality. The test case is easy to follow just check the code. What transposing is is something I am not going to explain here. You should read the book. In Wikipedia site there is actually a quite good animation that visualises what happens in traposing.

Invertin Matrices

What inverting means is basically that when multiplying matrices like A with B we get C and to go back from C to A we multiply C with the inverse of B. So that’s why inverse is the keyword here. So B:s inverse matrix is what is needed.

Determinant

First thing we need to get inverse matrix is a determinant so that must be tested. I will continue to use the build in functionality here. The test case for calculating determinant from 2×2 matrix is like this:

/**
 Chapter 3.
 9th test case. Inverting matrices.
 */
func testInvertingMatrices() {
    // Given
    let m = RGMatrix2x2(float2(x: 1, y: -3), float2(x: 5, y: 2))
    
    // Then
    XCTAssertEqual(m.determinant == 17, true)
}

Submatrix

The goal is to invert matrix and next thing needed is submatrix. This is something that is left when single column and row is deleted form the matrix. So the size of the matrix is also changed. When deleting a row and column from 4×4 matrix the result is 3×3 matrix and so on.

There’s no magic there, and not even really any math.

Jamis Buck

So this must be easy…well not that easy. Not for me anyway. When I tried to think how to construct the new matrix after deleting row and column I had to take a break, breaks. Well two breaks…no three…not sure if there were more of them. Lots of coffee was needed any way. So the solution at the end is probably no the best but it works 🙂

For the submatrix functionality there is no direct method in simd. So I needed to create new method and I will implement it as an extension. Extensions are really great feature in Swift. You can easily add functionality to existing datatypes without messing up the original datatype.

I wanted to implement this functionality so that the method is easy to use. In the books submatrix function interface submatrix(matrix, row, column) you could easily think that row and column arguments are of integer type. But what if a call is made using for example a row value that is greater than the actual number of rows? Error checking is needed or then make the implementation to return nil if input is not valid or something else.

I decided to create custom Swift enum to limit the scope of the input. In Swift language the enums are quite capable for all sorts of things. Here is the enums for this purpose.

/**
 RGPosition is an enum type used when index is needed.
 It's generic type to be used when ever index is limit
 from 0 to 2. Ex. in 3x3 matrix or float3 vector.
 */
public enum RGPosition3: Int {
    case First
    case Second
    case Third
}

/**
 RGPosition is type used when index is needed.
 It's generic type to be used when ever index is limit
 from 0 to 3. Ex. in 4x4 matrix or float4 vector.
 */
public enum RGPosition4: Int {
    case First
    case Second
    case Third
    case Fourth
}

In the test case I can then use these enums in the function call (function implementation will be explained later). Test case for 3×3 to 2×2 matrix is as follows:

/**
 Chapter 3.
 11th test case.
 Submatrices from 3x3 to 2x2 matrix.
 */
func testSubmatrices3x3to2x2() {
    // When
    // |  1 | 5 |  0 |
    // | -3 | 2 |  7 |
    // |  0 | 6 | -3 |
    let M = RGMatrix3x3(
        float3(x: 1, y: -3, z: 0),
        float3(x: 5, y: 2, z: 6),
        float3(x: 0, y: 7, z: -3)
    )
    
    let sub = RGMatrix2x2(float2(x: -3, y: 0), float2(x: 2, y: 6))
    
    // Then
    let subM = M.submatrix(c: .Third, r: .First)
    // | -3 | 2 |
    // |  0 | 6 |
    XCTAssertEqual(simd_equal(sub, subM) , true)
}

In my case the vectors are the columns. This means that the values in the vector as rows. Same for all the data types from float2 to float3 and to RGTuple which is same as float4. First thing to implement is a method that returns float2 from float3 after removing the selected row. Implementation is here:

/**
 In simd we handle vectors as 1 dimensional column matrix.
 So the values x, y and z could be handeld as rows.
 x
 y
 z
 */
func removeRow(_ row: RGPosition3) -> float2 {
    switch row {
    case .First:
        return float2(x: self.y, y: self.z)
    case .Second:
        return float2(x: self.x, y: self.z)
    case .Third:
        return float2(x: self.x, y: self.y)
        }
}

Equal approach is used for other vector datatypes and this is a place where I can use the RGPosition3 enum. I think this is a quite simple implementation.

Next thing is to implement the actual submatrix method. Well the same kind of technique could be used with matrix. I show you the code first so it’s easier to explain.

func submatrix(c columnToDelete: RGPosition3, r rowToDelete: RGPosition3) -> RGMatrix2x2 {
    switch columnToDelete {
    case .First:
        return RGMatrix2x2(self.columns.1.removeRow(rowToDelete), self.columns.2.removeRow(rowToDelete))
    case .Second:
        return RGMatrix2x2(self.columns.0.removeRow(rowToDelete), self.columns.2.removeRow(rowToDelete))
    case .Third:
        return RGMatrix2x2(self.columns.0.removeRow(rowToDelete), self.columns.1.removeRow(rowToDelete))
    }
}

Same RGPosition3 could be used here because the dimensional limits are the same for rows and for columns. The inputs for the method is column and row I want to get rid off.

Swift switch case block is used to check first what column is to be deleted. The resulting RGMatrix2x2 is constructed from the remaining two columns. Then the method created earlier for the vector types is used to remove the specific row.The exactly same approach is used for the RGMatrix4x4 too.

Minors

The minor of an element at row i and column j is the determinant of the submatrix at (i,j).

Jamis Buck

I have already the functionality to create submatrix and the determinant so this is easy part. Just put those two thing to gather in a method of it’s own. So there nothing fancy here. The code is easy to follow. To mention here the result of a minor is a scalar value.

Cofactor

In short the cofactor is based on the knowledge that the sum of the row and column of a minor is either odd or even. If it’s odd the cofactor is negated minor and if it’s even the cofactor is same as minor.

The test case you find in the code but here is the implementation. It’s implemented using the good old extension. Internaly it uses the minor method and then just uses basic if else bock to either negate the minor or return the original.

func cofactor(c columnToDelete: RGPosition3, r rowToDelete: RGPosition3) -> Float {
    let minor = self.minor(c: columnToDelete, r: rowToDelete)
    if (columnToDelete.rawValue + rowToDelete.rawValue) % 2 == 0 {
        return minor
        }
    else {
        return -minor
        }
}

At this point the book says to test the determinant of the larger matrices. I found that even when there is no methods for submatrix or cofactor in simd there is however method (actually a calculated property) for determinant of float4x4 matrix. I bet the simd implementation is good to go with.

What the determinant is needed for is to find if the matrix have inverse matrix at all. If determinant is 0 then there is no such thing and if it is anything else then there is.

The ultimate goal was to somehow calculate the inverse matrix. And guess what. For that purpose there is also build in stored property in simd. So I am not going to make this any harder at this point. I will count on that. And here is the test:

/**
 Chapter 3.
 16th test case.
 Test if matrix is invertible.
 */
func testInverseMatrix() {
    // Given
    // | -5 |  2 |  6 | -8 |
    // |  1 | -5 |  1 |  8 |
    // |  7 |  7 | -6 | -7 |
    // |  1 | -3 |  7 |  4 |
    
    // |  0.21805 |  0.45113 |  0.24060 | -0.04511 |
    // | -0.80827 | -1.45677 | -0.44361 |  0.52068 |
    // | -0.07895 | -0.22368 | -0.05263 |  0.19737 |
    // | -0.52256 | -0.81391 | -0.30075 |  0.30639 |
    
    // When
    let M = RGMatrix4x4(RGTuple(x: -5, y: 1, z: 7, w: 1),
                        RGTuple(x: 2, y: -5, z: 7, w: -3),
                        RGTuple(x: 6, y: 1, z: -6, w: 7),
                        RGTuple(x: -8, y: 8, z: -7, w: 4)
    )
    
    let inv = RGMatrix4x4(RGTuple(x: 0.21805, y: -0.80827, z: -0.07895, w: -0.52256),
                          RGTuple(x: 0.45113, y: -1.45677, z: -0.22368, w: -0.81391),
                          RGTuple(x: 0.24060, y: -0.44361, z: -0.05263, w: -0.30075),
                          RGTuple(x: -0.04511, y: 0.52068, z: 0.19737, w: 0.30639)
    )
    
    // These are here so we could see the differences in the values.
    print(M.inverse)
    print(inv)
    // Then
    XCTAssertEqual(M.inverse.isEqualTo(inv), true)
}

Aha! Even now when I use the simd:s build in inverse functionality I have to do the comparison of the two matrices. Well this was the thing I said earlier that I will “explain later”. Previously when comparing matrices the values had been precise so the comparison had worked well using the build in simd_equla() function for matrices. But at the case above the floating point values are not exactly equal to what are given and what the inverse returns. That why I have to implement the isEqaulTo method for the matrices. This is how it looks:

func isEqualTo(_ m: RGMatrix4x4) -> Bool {
    if self.columns.0.isEqualTo(m.columns.0) && self.columns.1.isEqualTo(m.columns.1) && self.columns.2.isEqualTo(m.columns.2) && self.columns.3.isEqualTo(m.columns.3) {
        return true
    }
    else {
        return false
    }
}

It looks quite the same for the 3×3 matrix as well.

Last thing is to test the original problem. The A * B = C and C * inverse(B) = A. So that is here and it’s the final piece.

/**
 Chapter 3.
 17th test case.
 Test original iverse case A * B = C so C * ivnerse(B) = A
 */
func testInverseMultiplicationByProduct() {
    // Given the 4x4 matric A:
    // |  3 | -9 |  7 |  3 |
    // |  3 | -8 |  2 | -9 |
    // | -4 |  4 |  4 |  1 |
    // | -6 |  5 | -1 |  1 |

    // And the following 4x4 matrix B:
    // |  8 |  2 |  2 |  2 |
    // |  3 | -1 |  7 |  0 |
    // |  7 |  0 |  5 |  4 |
    // |  6 | -2 |  0 |  5 |

    let A = RGMatrix4x4(RGTuple(x: 3, y: 3, z: -4, w: -6),
                        RGTuple(x: -9, y: -8, z: 4, w: 5),
                        RGTuple(x: 7, y: 2, z: 4, w: -1),
                        RGTuple(x: 3, y: -9, z: 1, w: 1)
    )
    
    let B = RGMatrix4x4(RGTuple(x: 8, y:  3, z: 7, w: 6),
                        RGTuple(x: 2, y: -1, z: 0, w:-2),
                        RGTuple(x: 2, y:  7, z: 5, w: 0),
                        RGTuple(x: 2, y:  0, z: 4, w: 5)
    )
    
    // When
    let C = A * B
    let iB = B.inverse
    // Then
    XCTAssertEqual((C * iB).isEqualTo(A), true)
    // To note the iverse(B) * C is not same as C * inverse(B)
    // To try out commen out line belowe and test.
    // XCTAssertEqual((iB * C).isEqualTo(A), true)
}

Whoah. Lots of stuff. Even without implementing all the functionality by my self there was lots of things to get familiar with. For now I will omit the “Putting it together” section. I will add it here later. You should check the source code. There on place in the code where I used the row based init method to create a matrix to follow the books approach. I leave it to you to find it 🙂 Until next time.