Deriving the Perspective Projection Matrix

This is an updated repost from my old, deactivated blog or from notes of my past studies. The date follows the original time of writing.

Hi there, how are your molecules doing?

Today we boldly go into one of the most important transformations in 3D Computer Graphics, namely the Perspective Projection. It doesn’t matter how deep your knowledge of 3D graphics is: whether you wrote your own rasterizer from scratch and implemented the transformation or used it wrapped in some library with your beloved graphics API without knowing its foundations, you’ve used this matrix!

So, in this post we’re going to see how we can come up with the classic perspective projection matrix from the foundations of geometry and linear algebra. After that, we’ll derive the D3D matrix and leave the GL version as an exercise for the reader.

The Geometry of Projection

Let’s start! The first thing I want you to pay attention to is figure 1; all the explanation I’m gonna give you in the next few minutes is based on this picture, so activate your chameleon eyes and follow me. We know for a fact that most objects in a scene are composed of vertices (in other words, the object is a discrete group of vertices), thus we “loop” over all the vertices and transform the object as we like. This fact is great, as it allows us to perform parallel work in a scene, working on a lot of vertices at the same time. Therefore, to project an object we need to project all its vertices, in which case the order doesn’t matter, as no spatial relation between vertices is needed for this operation.

Well, so how do we project? It’s very simple: looking at figure 1, we can decompose $v$ into two components, $v_z$ and $v_y$ ; then by similarity of triangles we get:

v'(v'_x,v'_y,v'_z) = v'(\frac{v_x}{v_z}d, \frac{v_y}{v_z}d, d) \tag{1}

If you’re wondering what this $d$ is and how I got $v'_x$ , keep reading; if not, skip to the next section! Well, $v'_x$ you can get by doing the same procedure but with the $z$ and $x$ axes. What about $d$ ? $d$ is actually the distance of the projection plane (the orange line in figure 1 is a slice of the projection plane) from the center of projection $C$ , so in this case $d$ is the bottom line segment of the smaller triangle; that’s how it ended up in the formula. An interesting point to note is that d is also actually the zoom: the greater d is, the bigger the objects appear!

Now that we know how to get $v'$ , we need a way to represent the transformation with a matrix. Enter the Homogeneous Coordinates!

Homogeneous Coordinates

When working with linear transformations (rotation, scaling, etc.) in 3D, we need a 3×3 matrix, but to work with affine transformations like translation we need one more dimension, thus we work with 4×4 matrices. The perspective projection also isn’t a linear transformation, therefore we also use 4×4 matrices!

This is where homogeneous coordinates come in handy; the notion of perspective projection and the homogeneous coordinates representation are a very good match. Take a look at figure 2:

We have a similar situation to figure 1, but now we are in 4D space (Oh yeah!), as p is a 3D vector with one more dimension w (or up, if you like!); we also have $C$ , the center of projection, but here it is called the point at infinity, where $w = 0$ ; you’ll see why shortly. When $w = 1$ , we consider the point $p$ to be in homogeneous form and we call it $p_w$ , thus the name. So all the points in the direction between $C$ and $p$ are the same point $p_w$ after the homogenization process, in more mathematical terms:

\forall p = (x,y,z,w)

after homogenization we get:

p_w = (\frac{x}{w},\frac{y}{w},\frac{z}{w},1)

Now you can see that if $w \rightarrow 0 \Rightarrow p \rightarrow C$ , in other words, if $w$ tends to zero, $p$ tends to the center of projection or the point at infinity.

So far so good, but how does this help us exactly? Well, we just saw how to represent our perspective-projected vertex. If you are reading this post, then surely you have read or heard something about the hardware z-divide; this divide is actually a homogenization if we store the $z$ value in the $w$ coordinate:

1 -> Before the transformation we have:

v = (x,y,z,1)

2 -> After the projection we have:

v_p = (x,y,z,\frac{z}{d})

3 -> And if we homogenize the vertex we get:

v_w = (d\frac{x}{z},d\frac{y}{z},d,1)

Let’s analyze these 3 steps more carefully. The first interesting thing to note is that in 1 the vertex v is homogenized and also in 3D space, while in 2 it isn’t, as it is only in projective space (I believe this is the right term!); but after the homogenization in 3 the vertex returns to 3D space, though not in the same location, as it was projected. Another point worth mentioning is that the actual projection happens in 3, during the homogenization process; also, 3 normally happens in dedicated hardware on GPUs.

So our mission is to find a matrix $P$ which transforms a vertex $v$ in format 1 and returns $v_p$ in format 2:

(x,y,z,\frac{z}{d}) = (x,y,z,1)P(d)

Transforming to matrix form and using the OpenGL column major form, we get:

\begin{pmatrix} x \\ y \\ z \\ \frac{z}{d} \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & \frac{1}{d} & 0 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix}

Therefore:

P(d) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & \frac{1}{d} & 0 \end{pmatrix}

Note that if you’re using right-hand coordinates like OpenGL does, you’ll have a $-\frac{1}{d}$ instead of a $\frac{1}{d}$

The matrix P we got is the perspective projection matrix [1], but the matrices we normally see in the D3D/GL docs are quite different from this one, so what is missing? Those matrices actually do more than projection: they transform the frustum-shaped space into a cuboid-shaped one, thus making clipping easy to do. Let’s see how Direct3D does it!

From Frustum to Canonical Cuboid

We are almost there! Now we have to find a way to convert our frustum space into a cuboid one; we have to maintain proportion in relation to size and distance between our objects; in other words, we need to do a mapping between these two spaces. Look at figure 3:

Obs: $p1$ should be $(r,t,n)$ not $(r,t,f)$ , my mistake!

We have our frustum space and 2 interesting points, $p_0$ and $p_1$ ; these 2 points are the extrema of the near plane, its min and max. As we are interested in placing all the points in this plane in canonical space we must map them from the projective plane (near plane) to the canonical plane. I wanna make a personal comment about this process: all the books and sources I’ve read so far (including Real-Time Rendering Third Edition) don’t explain in much detail this final part, in which you select the near plane and assume all the points are inside it. I believe that this is intentional, to make you think, and what I’m gonna present to you now is how I think about it. I’m gonna explain the x-axis case; you can infer the y one by symmetry. The condition a point must satisfy to be visible is that it MUST be inside the view frustum (ignoring all relations between objects, therefore ignoring hidden-surface removal algorithms for now!). This is the basic principle of the view frustum itself, thus if you draw a ray between the center of projection and the point itself, it MUST intersect the near plane; otherwise the point is outside the view frustum (Draw and see for yourself!).

Now that we understand why (hopefully! If not, draw again!) we map from the near plane, it is time to do some calculations: We have all the points in the following format (forget about the z for now; think of a slice of the frustum in the x direction, if you like!):

\forall p(x,y) \{ p(x,y) | l < x < r \land b < y < t \}

and we wanna map them to the following format:

\forall p(x,y) \{ p(x,y) | -1 < x < 1 \land -1 < y < 1 \}

Let’s do some calibrations for the x case (Sorry! I mean calculations, I’m not Garrus!):

- 1 < 2\frac{x - l}{r - l} - 1 < 1

Therefore, x and y (it’s all about symmetry, my friend!) in canonical space are:

x_c = 2\frac{x - l}{r - l} - 1 \land y_c = 2\frac{y - b}{t - b} - 1

We now work out the $z$ case. There is a point I have to make, though: D3D actually maps $z$ from 0 to 1, not -1 to 1 as GL does. We now have 2 new variables, $n$ and $f$ : $n$ is the near plane (the blue one in figure 3) and $f$ the far plane (the green one in figure 3). Working out the math we have:

0 < \frac{z - n}{f - n} < 1

Therefore our vertex in canonical space is:

v_c = (2\frac{x - l}{r - l} - 1, 2\frac{y - b}{t - b} - 1, \frac{z - n}{f - n})

Now we must transform our vertex v into its canonical space equivalent v_c using a matrix transformation:

v(x,y,z,1)T = v_c

Using matrix notation and the D3D row major notation, we have (if you are not comfortable with the idea of translating operations from algebra to matrix form, well, now would be a good moment to start training!):

\begin{pmatrix} 2\frac{x - l}{r - l} - 1 & 2\frac{y - b}{t - b} - 1 & \frac{z - n}{f - n} & 1 \end{pmatrix} = \begin{pmatrix} x & y & z & 1 \end{pmatrix} \begin{pmatrix} \frac{2}{r - l} & 0 & 0 & 0 \\ 0 & \frac{2}{t - b} & 0 & 0 \\ 0 & 0 & \frac{1}{f - n} & 0 \\ -\frac{r + l}{r - l} & -\frac{t + b}{t - b} & -\frac{n}{f - n} & 1 \end{pmatrix}

We now have our mapping matrix, which is a concatenation of both a scaling and a translation transformation; think about it! It makes total sense, as we scaled the frustum and then translated it:

S(s)T(t) = \begin{pmatrix} \frac{2}{r - l} & 0 & 0 & 0 \\ 0 & \frac{2}{t - b} & 0 & 0 \\ 0 & 0 & \frac{1}{f - n} & 0 \\ -\frac{r + l}{r - l} & -\frac{t + b}{t - b} & -\frac{n}{f - n} & 1 \end{pmatrix}

If you have some familiarity with orthogonal projection, you’ve probably noticed that this isn’t a perspective projection; you can see that there is no perspective division in the transformation, therefore we must account for the perspective projection.

Let’s analyze the x case (as y is symmetrical, I leave it to you!). We know that x is in the following format when using orthogonal projection:

x_p = \frac{2x}{r - l} - \frac{r + l}{r - l}

But we also know that a perspective-projected x coordinate must be scaled by d and divided by z, thus:

x_p = \frac{2}{r - l}(\frac{dx}{z}) - \frac{r + l}{r - l}

By symmetry we also know $y_p$ :

y_p = \frac{2}{t - b}(\frac{dy}{z}) - \frac{t + b}{t - b}

Now things tend toward a more engineering side than a mathematical one. As we saw earlier, the z divide is done by the hardware, so we have to store the $z$ for further homogenization, but if you look at both $x_p$ and $y_p$ , you’ll see that they have already been divided by $z$ , thus making our math wrong; we must change that! Another point to also consider is the translation factor, which will be divided by $z$ (this isn’t what we want by any means!). Therefore we multiply both $x_p$ and $y_p$ by $z$ , and we get the coordinates in the correct format:

$x_{p}z = \frac{2d}{r - l}x - z\frac{r + l}{r - l}$ $y_{p}z = \frac{2d}{t - b}y - z\frac{t + b}{t - b}$

We now look at $z_p$ . Following the same parameters as before, we wanna write $z_p$ in the following way (scaling + translation):

z_{p}z = pz + q

But why do we do this with $z_p$ instead of just setting it to $d$ ? Well, we saw an explanation of why our mapping from a frustum to a cuboid works (and we made sure it was right by assuming a perspective-projected point); now we must remember that our vertex may be out of the frustum vision range, therefore we must map the depth of our frustum to the $[0,1]$ range for further clipping. Now note that $z=n \rightarrow z_{p}=0 \land z=f \rightarrow z_{p}=1$ .

Working out the $z = n$ case, we get:

q = -pn

And substituting $q = -pn$ for the $z = f$ case:

p = \frac{f}{f - n}

Therefore:

z_{p}z = \frac{f}{f - n}z -\frac{nf}{f - n}

I just wanna point out that the best explanation I found of the $z_{p}z$ derivation was given by Joe Farrell in an article he wrote for Code Guru back in 2005, thus my explanation in this part was based on his article.

And finally, transforming to matrix form, we have the Direct3D Perspective Projection Matrix [2]; we also set our projective plane $d = n$ :

P_{[0,1]} = \begin{pmatrix} \frac{2n}{r - l} & 0 & 0 & 0 \\ 0 & \frac{2n}{t - b} & 0 & 0 \\ -\frac{r + l}{r - l} & -\frac{t + b}{t - b} & \frac{f}{f - n} & 0 \\ 0 & 0 & -\frac{fn}{f - n} & 1 \end{pmatrix}

Working with a symmetrical projective plane, we have $|t| = |b|$ and $|r| = |l|$ , thus:

P_{[0,1]} = \begin{pmatrix} \frac{2n}{r - l} & 0 & 0 & 0 \\ 0 & \frac{2n}{t - b} & 0 & 0 \\ 0 & 0 & \frac{f}{f - n} & 0 \\ 0 & 0 & -\frac{fn}{f - n} & 1 \end{pmatrix}

That is the matrix that the old D3DX function D3DXMatrixPerspectiveLH implements.

And to end with a final blast, we look at the perspective projection as a function of the fov, the width and the height of the projective plane! Just pay attention: instead of $b - t$ and $r - l$ , we now use $h$ and $w$ , which is the same for the symmetric case.

Look at the figure above, simple stuff, right? Just remember your trig!

tan(\frac{\theta}{2}) = \frac{w}{2n}

Inverting we have:

\frac{1}{tan(\frac{\theta}{2})} = cot(\frac{\theta}{2}) = \frac{2n}{w}

Well, this result looks familiar, doesn’t it? Sure it does, as it is the scaling factor of the symmetrical matrix we got above, with the respective $x$ and $y$ versions, therefore:

P_{[0,1]} = \begin{pmatrix} cot(\frac{\theta_{w}}{2})& 0 & 0 & 0 \\ 0 & cot(\frac{\theta_{h}}{2}) & 0 & 0 \\ 0 & 0 & \frac{f}{f - n} & 0 \\ 0 & 0 & -\frac{fn}{f - n} & 1 \end{pmatrix}

In case you just wanna work with the height fov, we can use the height fov plus the aspect ratio $\alpha$ to get the width fov:

\alpha = \frac{w}{h}

\frac{cot(\frac{\theta_{h}}{2})}{\alpha} = \frac{\frac{2n}{h}}{\frac{w}{h}} = \frac{2n}{w}

And we have:

P_{[0,1]} = \begin{pmatrix} \frac{cot(\frac{\theta_{h}}{2})}{\alpha}& 0 & 0 & 0 \\ 0 & cot(\frac{\theta_{h}}{2}) & 0 & 0 \\ 0 & 0 & \frac{f}{f - n} & 0 \\ 0 & 0 & -\frac{fn}{f - n} & 1 \end{pmatrix}

This is exactly the matrix D3DXMatrixPerspectiveFovLH implements.

Well, my friend, our exploration is finished for today! I really hope that this post was of some assistance to you during your studies!

If you find any typos or mistakes, just let me know and I’ll fix them as fast as possible!

LIVE LONG & PROSPER

References

[1] Foley, van Dam, Feiner, Hughes, Computer Graphics Principles and Practice 2nd, Addison-Wesley, 1993
[2] Akenine-Möller, Haines, Hoffman, Real-Time Rendering 3rd, CRC, 2008