WebGL Engine from Scratch 11: Normal Transforms and Linear Algebra Functions
Tags: linearalgebra, vanillajs, webgl
Published: 2022-01-26 05:48:37

So I made a bit of mistake before by thinking that I could just apply the model space transformation to normals to get the correct result.  Turns out that's not the case.  It is for just rotation and uniform scaling (which is what I tested) but it's not actually true for non-uniform scaling.

You can see what I mean:

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/rr8dbjda9a93wid5i85i.png)

This is a sphere that has been compressed in Y and lengthened in X.  The shader shows the normals.  We can see that it has a fair bit of magnitude going in X (red) but not a lot in Y (green).  But we expect the normals to be pointed upward, so it should be more green. 

To fix this we need to multiply by the transpose of the inverse of the model view matrix.  Why exactly is better explained by a mathematician as I'm not going to derive it.  But what this does mean is that we now have to deal with more of the seedy parts of linear algebra.

I'm not going to pretend any of this is optimized, it's not and there are almost certainly more efficient ways to do it but this is translating the hand math into algorithms.

## Getting an inverse

So if you search this should be something you can do in GLSL using the `inverse` function.  As far as I can tell it's valid for OpenGL ES 3.0 which is what WebGL 2.0 should use.  However, trying it doesn't seem to work:

```glsl
mat4 test = inverse(mat4(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0));
```
```
 '=' : cannot convert from 'const mediump float' to 'mediump 4X4 matrix of float'
```
This mean the equal failed because inverse returned a float, not a matrix, so whatever `inverse` does in WebGL's GLSL, it's not what we want.

We need to compute it on the CPU side.  This has benefits, namely that the inverse of a matrix is a costly operation and we wouldn't want to do it once per vertex anyway.

To get the inverse of a matrix we start with the formula:

```js
export function scaleMatrix(matrix, s){
	const result = [];
	for (let row = 0; row < matrix.length; row++) {
		const newRow = [];
		for (let col = 0; col < matrix[row].length; col++) {
			newRow.push(matrix[row][col] * s);
		}
		result.push(newRow);
	}
	return result;
}

export function getInverse(matrix){
	return scaleMatrix(getAdjugate(matrix), 1 / getDeterminant(matrix));
}
```
`scaleMatrix` is just an analog for `scaleVector` applied to matrices.  We're multiplying the adjugate with 1 over the determinant.

### Determinant

The determinant is formed by iterating through the top row of a matrix.  Just a note that we'll always be dealing with square matrices, but these won't work if they aren't.  We take each element and then multiply it by the determinant of the sub matrix formed by deleting the row and column that the element appeared in.  We alternate adding and subtracting each term.  So for a 3x3:

```
[
 [1, 2, 3],
 [4, 5, 6],
 [7, 8, 9]
]
```
The determinant is:

```
determinant = 1 * determinant([[5,6],[8,9]]) - 2 * determinant([[4,6],[7,9]]) + 3 * determinant([[4,5],[7,8]])
```
If the matrix is size 2 then the determinant is just:

```
//2-length
determinant = matrix[0][0] * matrix[1][1]) - (matrix[0][1] * matrix[1][0]
```
Using these we can build an algorithm to find the determinant.

```js
export function getDeterminantSubmatrix(matrix, row, col){
	const result = [];
	for(let i = 0; i < matrix.length; i++){
		const newRow = [];
		if(i === row) continue;
		for(let j = 0; j < matrix[i].length; j++){
			if(j === col) continue;
			newRow.push(matrix[i][j]);
		}
		result.push(newRow);
	}
	return result;
}

export function getDeterminant(matrix){
	let result = 0;

        if(matrix.length === 1) return matrix[0][0]; //[edit] I forgot this case previously
	if (matrix.length === 2 && matrix[0].length === 2) return (matrix[0][0] * matrix[1][1]) - (matrix[0][1] * matrix[1][0]);

	for(let i = 0; i < matrix[0].length; i++){
		if(i % 2 === 0){
			result += matrix[0][i] * getDeterminant(getDeterminantSubmatrix(matrix, 0, i));
		} else {
			result -= matrix[0][i] * getDeterminant(getDeterminantSubmatrix(matrix, 0, i));
		}
	}

	return result;
}
```
This is a direct translation.  If the matrix is 2x2 we can directly find the determinant, otherwise, we break it down and recursively find the determinants using `getDeterminantSubmatrix` to find the submatricies.

Next we need to get the adjugate.  The adjugate is the transpose of a cofactor matrix.  So what's a cofactor matrix?  It's a matrix where each term is a cofactor of the original element.  So what's a cofactor?  A cofactor is the same as the determinant except it will be negated at every other value.  Well, almost every other value, it's actually closer to a checkerboard pattern where each row is offset by an extra 1.

```js
export function getCofactor(matrix, row, col){
	const determinant = getDeterminant(getDeterminantSubmatrix(matrix, row, col));
	return (row + col) % 2 === 1
		? -determinant
		: determinant;
}
```
So the cofactor matrix is a matrix where we've gotten the cofactor for each element, easy enough:

```js
export function getCofactorMatrix(matrix){
	const result = [];
	for(let row = 0; row < matrix.length; row++){
		const newRow = [];
		for (let col = 0; col < matrix[row].length; col++) {
			newRow.push(getCofactor(matrix, row, col));
		}
		result.push(newRow);
	}
	return result;
}
```
Finally we can get the transpose which is the matrix flipped 90 degrees.  We already had this, but it didn't work for arbitrary matrix sizes.

```js
export function transpose(matrix) {
	const result = [];
	for(let row = 0; row < matrix.length; row++){
		const newRow = [];
		for(let col = 0; col < matrix[row].length; col++){
			newRow.push(matrix[col][row]);
		}
		result.push(newRow);
	}
	return result;
}
```
And when we put it all together we get the inverse.

Next we need to multiply the view and model matrices.  We have a function of that already but it makes sense to make it work for arbitrary sizes.  To do this for each row of the first term and row in the second term we need to get the dot product:

```js
export function getColumn(matrix, col){
	const result = [];
	for(let row = 0; row < matrix.length; row++){
		result.push(matrix[row][col]);
	}
	return result;
}

//A's rows must equal B's columns, no check is given
export function multiplyMatrix(a, b) {
	const result = [];
	for (let row = 0; row < a.length; row++) {
		const newRow = [];
		for (let col = 0; col < b[row].length; col++) {
			newRow.push(dotVector(a[row], getColumn(b, col)));
		}
		result.push(newRow);
	}

	return result;
}
```
The comment above about the size is that we cannot multiply matrices that do not match up in rows and columns.  But assuming they do we just iterate through each row by each column and do the dot product.

One thing to note is that our old dot product function won't quite work, we need to expand it to arbitrary size:

```js
//vectors must have same length!
export function dotVector(a, b){
	let result = 0;
	for(let i = 0; i < a.length; i++){
		result += a[i] * b[i];
	}
	return result;
}
```
One final tool we'll need is the take a 4x4 matrix and convert to a 3x3.  This is as simple as removing the last row and column:

```js
export function trimMatrix(matrix, height, width) {
	const result = [];
	for (let row = 0; row < height; row++) {
		const newRow = [];
		for (let col = 0; col < width; col++) {
			newRow.push(matrix[row][col]);
		}
		result.push(newRow);
	}
	return result;
}
```

We can wire this up by adding another uniform:

```js
bindUniforms(modelMatrix){
	const program = this.context.getParameter(this.context.CURRENT_PROGRAM)
	const modelMatrixLocation = this.context.getUniformLocation(program, "uModelMatrix");
	this.context.uniformMatrix4fv(modelMatrixLocation, false, modelMatrix);
	const viewMatrix = this.cameras.default.getViewMatrix();
	const normalMatrixLocation = this.context.getUniformLocation(program, "uNormalMatrix");
	const viewModelInverseTranspose = transpose(getInverse(trimMatrix(multiplyMatrix(asMatrix(modelMatrix, 4, 
4), asMatrix(viewMatrix, 4, 4)), 3, 3)));
	this.context.uniformMatrix3fv(normalMatrixLocation, false, viewModelInverseTranspose.flat());
}
```
This gets a bit complicated because the class gives back a flattened array but we're using shaped arrays to represent matrices.  We could rewrite the code, but I don't want to do that right now, instead we'll create another helper to shape the flat array back into a matrix.

```js
export function asMatrix(array, height, width){
	const result = [];
	for (let row = 0; row < height; row++) {
		const newRow = [];
		for (let col = 0; col < width; col++) {
			newRow.push(array[row * width + col]);
		}
		result.push(newRow);
	}
	return result;
}
```
In the vertex shader we can now use this:

```js
//vertex shader
precision mediump float;

uniform mat4 uProjectionMatrix;
uniform mat4 uModelMatrix;
uniform mat4 uViewMatrix;
uniform mat3 uNormalMatrix;

attribute vec3 aVertexPosition;
attribute vec3 aVertexColor;
attribute vec2 aVertexUV;
attribute vec3 aVertexNormal;

varying vec4 vColor;
varying vec2 vUV;
varying vec3 vNormal;
varying vec3 vPosition;

void main() {
	gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix * vec4(aVertexPosition, 1.0);
	vUV = aVertexUV;
	vColor = vec4(aVertexColor, 1.0);
	vNormal = uNormalMatrix * aVertexNormal;
	vPosition = vec3(uModelMatrix * vec4(aVertexPosition, 1.0));
}
```
This gives us the correct normals:

![Image description](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/waru1cpdbpbgquv5gzjm.png)
 
Note that the magnitude of the normals is not normalized in the pixel shader so it's actually a little more green than you might expect.

## Generalizing other functions

Since we've started to re-write most of the `vector.js` library, we can apply all the operations to arbitrary sized matrices and vectors.  This will be slightly less efficient as there will be overhead for function calls and boundary conditions.  One strategy is to make specific functions for certain sizes like 2x2, 3x3, 4x4 since you'll nearly always be dealing with a fixed size square matrix.  I didn't bother with that.  Another idea could be to represent the matrices as flat buffers which could make certain operations a bit faster since we need to convert them to buffers anyway (see issue above). I went with a straightforward approach, just implementing the raw algorithms in an elegant way, no concerns about optimization.

```js
//vectors
export function addVector(a, b) {
	return a.map((x, i) => x + b[i]);
}

export function subtractVector(a, b) {
	return a.map((x, i) => x - b[i]);
}

export function scaleVector(vec, s) {
	return vec.map(x => x * s);
}

export function divideVector(vec, s) {
	return vec.map(x => x / s);
}
```
For matrices we can use a helper function `mapMatrix` to make traversal easier and the code more compact, but it has a slight performance penalty.  Still smaller code is usually better on the web if it's not extreme slow down.  In the same spirit you don't need to add the functions you don't use (scale and divide, add and subtract are technically the same thing for example), I'm just listing them here as a reference.

```js
export function mapMatrix(matrix, func){
	const result = [];
	for (let row = 0; row < matrix.length; row++) {
		const newRow = [];
		for (let col = 0; col < matrix[row].length; col++) {
			newRow.push(func(matrix[row][col]));
		}
		result.push(newRow);
	}
	return result;
}
export function addMatrix(a, b){
	return mapMatrix(a, (x, r, c) => x + b[r][c]);
}

export function subtractMatrix(a, b) {
	return mapMatrix(a, (x, r, c) => x - b[r][c]);
}

export function scaleMatrix(matrix, s) {
	return mapMatrix(matrix, (x, r, c) => x * s);
}

export function divideMatrix(matrix, s) {
	return mapMatrix(matrix, (x, r, c) => x / s);
}
```

What about the cross product?  It doesn't usefully generalize.  The general form of it is the "exterior product" or "wedge product" which gets into a bunch of stuff that's currently unnecessary.

Anyway this was a bit of a strange and less exciting detour into the world of linear algebra.  But now we have more tools in our toolbox.