Model-View-Projection Matrices Derivation

A vertex shader processes vertices, simple enough, but no matter what we do, at the end of the day, we're just turning knobs on a complicated state machine that rasterizes stuff on a screen. The ultimate rasterization entrails are not controllable, but how we map points in three-dimensions onto a two-dimensional plane for them to be rasterized is up to us. The wizard beards of yore tackled these issues and set up a cannonical solution, the MVP matrix (and a lot of other stuff, as always)

Say we have a collection of points (vertices) that are ordered intelligently to allow for rasterization (a mesh). These points all have positions that are defined with respect to a local origin. A model matrix is simply a 4x4 matrix (so as discussed before, it contains the information for position, rotation and scale) that places that mesh in the shared vector space of all the meshes (world space).

As an example, taken from Marco Alamia here are three different meshes (instances of same mesh) with three different world positions. I think it's slightly confusing to even have a seperate "model space". Intuitively, I think it makes more sense to imagine all of these superimposed on each other at the world's origin and then simply, scaled, rotated and translated as desired

Now for the more interesting part, let's write a simple model matrix to check our understanding in our graphics API. The "Hello, World!" for me, a mere mortal, is the screen quad, so two triangles that cover the screen. However you want to process them, arrayed or index drawing, there are 4 different vertices

vertexPositions = [-1, +1, 0,
                   -1, -1, 0,
                   +1, -1, 0,
                   +1, +1, 0
                   ];

I'll color it with shadertoy's default vec3 col to ease the monotony of using a primary color.

vec3 col = 0.5 + 0.5*cos(t +uv.xyx + vec3(0,2,4));

Our Hello Quad

As a simple test, let's scale the \(xy\) plane down to a quarter, translate circularly in the plane and spin around the \(z\) axis. I think the circular translation should be pretty obvious, but the rotation matrix is slightly more msyterious and because we're all about first principles, here is a derivation of it (necessary trig derivation included). (Note, this is only a rotation around a principal axis, namely the \(z\) axis, this is analogous to any other rotation around either \(x\) or \(y\))

Math is nice, but how to write a matrix that our shader program can understand? A real "gotcha" for wanting to test this is out is the confusion between how matrices are represented mathematically and how they're laid out in an array to be processed by the API. It's doubly confusing when seeing different resources that are geared toward a specific API that will use a different matrix writing convention, for example OpenGL or D3D/XNA (Column-major notation vs. Row-major notation respectively).

The Column-major notation matrix notation used in OpenGL documentation does not describe in-memory layout for OpenGL matrices OpenGL matrices have the same memory layout as DirectX matrices. OpenGL matrices are 16-value arrays with base vectors laid out contiguously in memory. You can read about this here

model =[ 0.25 * cos(  t ),   -0.25 * sin(t), 0,    0,
         0.25 * sin(  t ),    0.25 * cos(t), 0,    0,
                        0,                0, 1,    0,
         0.50 * sin(0.5t),  0.5 * cos(0.5t), 0,    1
       ];

//... in vertex shader ...
gl_Position =  model * vec4(vertexPos, 1.0);

Here is our Model Demo

And that's pretty much it, really. If you're comfortable with using 4x4 matrices to encode any rigid transformation, there's nothing more to it.

This was of course a very simple example, but any arbitrary position and orientation could be expressed with a model matrix. We're limited right now in what we can test by our ability to see it. The rendering context will only show whatever is in the normalized device coordinates. This bring us to our next object, how to match our scene geometry to this rendering constraint.

This is really where things become less obvious and most resources fail to sufficiently explain what's going on. See this for example

This brings me to what I think is the most unecessary and confusing aspect of using these different transformations.

Textbooks and tutorials are constantly talking about different spaces and transforming from one space to another and while this has some organizational/ pedagogical merit in it, there's really only one coordinate space.

The now often invoked authority of the wizard beards of yore deigned this to be a good solution:

1. Take everything that exists in the coordinate space bring it in front of a another thing in the coordinate space called the camera.
2. Project it all onto a plane relative to that camera.
3. Map it to a unit volume.
4. Rasterize.

That's an over simplification, but that's what's happening.

The previous link says "When we're talking about camera/view space we're talking about all the vertex coordinates as seen from the camera's perspective as the origin of the scene"

And then they go on to talk about the Gram-Schmidt process and give a magical thing called a lookAt matrix, explaining "A great thing about matrices is that if you define a coordinate space using 3 perpendicular (or non-linear) axes you can create a matrix with those 3 axes plus a translation vector and you can transform any vector to that coordinate space by multiplying it with this matrix.".

This is technically true, but what's really happening is not all clear and why this is true is not explained.

Very clearly, all we're trying to do is transform the virtual camera/ pinhole camera placeholder such that it's centered at the origin and pointing down an axis. (by convention: \(-z\) in OpenGL, \(+z\) in D3D/XNA but it really could be whatever if the API were constructed differently).

That's it.

Of course if we do that, everything else needs to correspondingly transform to keep the congruence of space, because that's how reality works. The relative position of something to the camera does not change.

So how and why are we doing this?

Here is our View Demo

As you can see, we have identical results to our model test as expected, that's definitely encouraging :)

Note that I've rounded the \(z\) value to its closest integer (one of the vertices' \(z\) value is like -1.005 or something and just barely gets clipped , anything outside the cannonical view volume gets clipped or culled, see picture below.

A visual summary of what we're going to be doing to help with intuition:

Perspective projection just employs the mathematical object of a frustum, the portion of a solid that lies between one or two (in our case, "near" and "far" planes) parallel planes cutting it.

All that previous work was to get it along the \(z\) axis to aide in projecting and clipping along these frustum planes

Deriving the \(x\) and \(y\) coordinates of our projected point

Note that both expressions are identical up to a change of variable name, so you can directly see that both \(x\) and \(y\) have the same mapping.

Let's pause and think about this for a minute.

The derived expressions both have \(z\) dependence and we need to ask ourselves, does this makes sense?

Let's consider two points in the frustrum with identical \(x\) values, but different \(z\) values

\((x_1, y_1, z_1)\) & \((x_2, y_2, z_2)\) such that \(x_1 = x_2; z_1 < z_2\)

We need to check the limiting cases of either one and a general comparison of the two. Remember we did all that work to situate the frustrum along the negative z axis. So \(z\) can range from the near plane value, \(n\), to far plane value, \(f\).

For a given point:

As z approaches n, the x value goes to just its x value, so that makes perfect sense.

As z approaches f, the x value goes to some fraction of x, \(n/f\), anything more and it's clipped. This is kind of inconclusive at first blush for me, so let's compare.

Comparison:

The point with greater z value will always be smaller, closer to z axis for an identical x value as n is constant. So closer points seem "larger", spread further from the camera axis (z axis), and further points seem "smaller" closer to the camera axis (z axis). Again The expression for the projected y value is identical and an identical argument would be made.

Quick recap:

We know from the beginning of the conversation that we can represent any arbitrary linear transformation as a matrix. By augmenting the 3D position vector of our vertices to 4D we can hit it with a 4x4 matrix that simultaneously scales, rotates and translates our vector. In turn, for the sake of expediting and standardizing the projection and clipping processes, we transformed this now 4x1 vector (the product of the model and position vector) with another 4x4 transformation, the view transformation. This resulted in another 4x1 vector, now said to be in camera/ eye space. We did the math by hand and showed how to go from this point to the cannonical view volume where such things like the Sutherland-Hodgman clipping algorithm are executed.

How to represent this as a matrix so we can concatenate it with our view and model transformation?

Let's see what we got once more for convenience:

\[x_{unit-cube} = \frac {2nx} {(r - l)z} - \frac {r + l} {r - l}\] \[y_{unit-cube} = \frac {2nx} {(t - b)z} - \frac {t + b} {t - b}\]

Working from first principles like this, the expressions seem inextricabley mixed or indeterminate, but we have a way out, namely we need to take in consideration what the API is doing more fully.

After the vertex shader is run, it's final output is recorded (called transform feedback) and a number of other transformations and fixed function operations are executed in vertex post processing one of them being perspective division.

Whatever the final 4x1 vector after the concatenated transformations is: \[\begin{bmatrix} x \\ y \\ z \\ w \\ \end{bmatrix}\]

the vectors values are divided by its \(w\) value: \[\begin{bmatrix} \frac x w \\ \frac y w \\ \frac z w \\ 1 \\ \end{bmatrix}\]

Trying to square the geometrically derived results with a projection matrix as seen in a textbook was really confusing until I read about this. This problem is very representative of working with a complicated API in general. It has all sorts of small details that are there, silently working away, always for good or necessary reasons, but fairly opaque unless you're an initiate. I can only offer my commiserations to you and the wasted time we'll both spend trying to learn similar things in the future.

It's kind of circular, I'm sure the wizard beards who came before me made this Vertex Post-Processing step because of the hariness of the projected expressions; but for our sake trying to work backwards, to get the correct result, if we're going to divide by \(w\) no matter what, then we need to make sure whatever is in the divisor of \(x\) expression (\(z\)) is the value that \(w\) will take on after the transformation. Since \(w\) up until this point will be \(1\), we need only set the \(z\) column of the \(w\) row to be \(1\)

That is to say:

\[\begin{bmatrix} & & & \\ & & & \\ & & & \\ 0 & 0 & 1 & 0 \\ \end{bmatrix} \begin{bmatrix} x \\ y \\ z \\ 1 \\ \end{bmatrix}\]

\[=\begin{bmatrix} \\ \\ \\ z \\ \end{bmatrix}\]

Where I've left the other elements blank, because there is no way to reason about them before establishing this point.

Equipped with that bit of knowledge, we can see how the expressions for \(x\) and \(y\) will fit in the matrix.

Let's see them again for our mind's eye's convenience:

\[x_{unit-cube} = \frac {2nx} {(r - l)z} - \frac {r + l} {r - l}\]

\[y_{unit-cube} = \frac {2ny} {(t - b)z} - \frac {t + b} {t - b}\]

\[\begin{bmatrix} \frac {2n} {(r - l)} & 0 & -\frac {r + l} {r - l} & 0 \\ & \frac {2n} {(t - b)} & -\frac {t + b} {t - b} & 0 \\ & & & \\ 0 & 0 & 1 & 0 \\ \end{bmatrix}\begin{bmatrix} x \\ y \\ z \\ 1 \\ \end{bmatrix}\] \[=\begin{bmatrix} \frac {2nx} {(r - l)} -\frac {(r + l)z} {r - l} \\ \frac {2ny} {(t - b)} -\frac {(t + b)z} {t - b} \\ \\ z \\ \end{bmatrix}\]

And it's clear that with the vertex post-processing \(w\) division, we get back the projected expressions for our 4x1 resultant vector's \(x\) and \(y\) values \[=\begin{bmatrix} \frac {2nx} {(r - l)z} -\frac {(r + l)} {r - l} \\ \frac {2ny} {(t - b)z} -\frac {(t + b)} {t - b} \\ \\ 1 \\ \end{bmatrix}\]

Good progress, but what about an expression for \(z\)?

We can play the same game of linearly mapping from the view frustrum to the cnanonical vieww volume, but we need to remember that we'll be dividing by the \(w\) value eventually and that the \(w\) value takes on whatever the original \(z\) value was.

A linear mapping is of the form \(z^{\prime} = mz + b\)

This is really all we want, but whatever our mapped \(z\) value is will eventually be divided by the \(w\) value, that is to say the \(z\) value , so for the sake of constructing a matrix, we must offset this by changing the expression to be \(z^{\prime}z = mz + b\)

Our finished Perspective matrix (This agrees with figure 4.68 in Real Time Rendering 3rd Ed.): \[\begin{bmatrix} \frac {2n} {(r - l)} & 0 & -\frac {r + l} {r - l} & 0 \\ 0 & \frac {2n} {(t - b)} & -\frac {t + b} {t - b} & 0 \\ 0 & 0 & \frac {f + n} {f - n} & - \frac {2fn} {f-n} \\ 0 & 0 & 1 & 0 \\ \end{bmatrix}\]

Checking ourselves: \[\begin{bmatrix} \frac {2n} {(r - l)} & 0 & -\frac {r + l} {r - l} & 0 \\ 0 & \frac {2n} {(t - b)} & -\frac {t + b} {t - b} & 0 \\ 0 & 0 & \frac {f + n} {f - n} & - \frac {2fn} {f-n} \\ 0 & 0 & 1 & 0 \\ \end{bmatrix}\begin{bmatrix} x \\ y \\ z \\ 1 \\ \end{bmatrix}\] \[=\begin{bmatrix} \frac {2nx} {(r - l)} -\frac {(r + l)z} {r - l} \\ \frac {2ny} {(t - b)} -\frac {(t + b)z} {t - b} \\ \frac {(f + n)z} {f - n} - \frac {2fn} {f-n} \\ 1 \\ \end{bmatrix}\]

And after the \(w\) division, our final 4x1 vector:

\[=\begin{bmatrix} \frac {2nx} {(r - l)z} -\frac {(r + l)} {r - l} \\ \frac {2ny} {(t - b)z} -\frac {(t + b)} {t - b} \\ \frac {(f + n)} {f - n} - \frac {2fn} {(f-n)z} \\ z \\ \end{bmatrix}\]

It’s as simple as that; now, let's try to implement it.

Let's keep things nice and simple. Let's choose our near plane to be -1 and our far plane to be -10. Let's also choose a square near plane to project to so the subtractive terms in \(x\) and \(y\) goes to zero, say -1 and 1 for right and left and top and bottom respectively.

Plugging in our chosen parameters, we get our projection matrix

// ... in main passed as uniforms
projection = [
-1,  0,    0, 0,
0, -1,    0, 0,
0,  0, 11/9, 1,
0,  0, 2.22, 0
]);

Excited, we immedediately hook it up to our vertex shader and…. nothing happens

We've already discussed the inherit difficulties of working with a complicated API, but if you naively try this in your vertex shader, you won't get the result back you expect. I ran into a real snag trying to naively put in our derived matrix into my vertex shader. There is yet more post-processing I was unaware of, see the OpenGL matrix in Real Time Rendering for a further discussion (equation 4.69 in the third edition)

There is always a solution. Debugging, we can manually make our vertex shader coords what we mathematically think they should be and get our perspective that way.

//... in vertex shader ...
vec4 pos =  round(view * vec4(vertexPos, 1.0));
pos.z -= 5. * abs(sin(0.2 * time)) + 1.;
vec4 projPos = projection * pos;
projPos /= projPos.w;
gl_Position = projPos;

And we have perspective, I've never been more excited about a rectangle in my life.

Our first projection test

Just to really prove it, let's spin it around the \(y\) axis while we do the same translation

// Please note this is pseudo-ish code.
translationX = 0.5 * sin(0.5 * t);
translationY = 0.5 * cos(0.5 * t);
model =
[
cos(t),            0,  sin(t), 0,
0,            1,       0, 0,
sin(t),            0,  cos(t), 0,
translationX, translationY,       0, 1
];

//... in vertex shader ...
vec4 pos =   model * round(view * vec4(vertexPos, 1.0));
pos.z -= 5. * abs(sin(0.2 * time)) + 1.;
vec4 projPos = projection * pos;
projPos /= projPos.w;
gl_Position = projPos;

Our second projection test