Understanding the math behind perspective matrix in OpenGL

Question

I've been trying to figure out the math behind perspective matrix for 2 weeks now but I'm failing badly. I understand the theory behind the perspective matrix but I am not sure how the math works.

The code:

def perspective(fov: Float, aspect: Float, zNear: Float, zFar: Float) = {
  val h = Math.tan(fov * 0.5f).toFloat
  val c00 = 1.0f / (h * aspect)
  val c11 = 1.0f / h
  val c22 = (zFar + zNear) / (zNear - zFar)
  val c23 = (zFar + zFar) * zNear / (zNear - zFar)

  Matrix4(
    c00,    0,      0,     0,
    0,    c11,      0,     0,
    0,      0,    c22,   c23,
    0,      0,     -1,     0
  )
}

What I understand

I've seen video tutorials by Jorge Rodriguez & Arpan Pathak but I cannot fully relate it to the perspective matrix in the function above.

• Congruent triangles: Following Arpan's video I understand that to project a point P=(x, y, z) onto the a 2D XY plane I need to create a frustum and which is then used to find it's projection P^'=(x^', y^', z^') using congruent triangles relationship. Arpan's video and the final matrix makes sense but I do not see how it relates to the perspective function above? Here is my attempt.

• 

• Following the above I switched the frustum in opposite direction to visualise how it would apply in OpenGL using the perspective function above with zFar and zNear parameters but it's no way close to the perspective matrix in the code.

• 

Answer 1

Your h at the bottom is missing the left/right and top/bottom parameters; cf. to the code at the top where they are expressed as angle fov and aspect ratio.

Thanks to the 1 in the w-row of the matrix you get the z value as 4th coordinate, because this is the other parameter for the new x and y values. A high x (same for y) value at large distance/depth shall be displayed near the center.

The idea of the perspective projection is to choose a frustum and turn it into a cube ranging from -1 to 1 in all 3 dimensions - the NDC coordinates.

So for x and y the relation is:

Xndc = h * X / Z

This divide by Z cannot be expressed in a matrix and is done separately. h stands for the width and heigth of the near plane.

The depth Zndc is normalized between the near and far plane. If it is done linearly you get A*z + B = 0, with A*near + B = -1 and A*far + B = 1.

But depth wants to be stored non-linearly, with more precision at the near side. So the general equation is (A*z+B)/z, or A + B/z. And again you can define the two planes:

-1 = A + B/near
 1 = A + B/far

Now A and B have to be expressed in terms of near and far.

A = 1 - B/f into first equation -->

-1 = 1 - B/f + B/n           minus one

-2 = -B/f + B/n              multiply by f*n

-2fn = -Bn + Bf              factor out B

B(f-n) = -2fn                divide by f-n    

B = -2fn / (f-n)

And this B is exactly your c23, with signs and subtraction reversed. A is even easier to calculate.

For the matrix, A + B/z is not possible, so we have to switch to the (A*z+B)/z form with extra divide by z as for x and y.

Putting it together gives

Zndc = ( Z*(f+n)/(f-n) - 2fn/(f-n) )  /  Z

where the linear part Z*A + B is done by the matrix and W is assumed to be one so B is really a constant.

This brings me back to the lone "1", or rather "-1" ??, on the bottom row of the matrix. But here is my Frustum matrix, in glsl column-major (sigh) layout:

/* Frustum with aspect ratio */
double fside = 1 / u.frhtan);
double pldiff = u.frf - u.frn;
mat4 F = {
    vec4(fside, 0,        0, 0),
    vec4(0, fside*u.aspr, 0, 0),
    vec4(0, 0,           (u.frf + u.frn) / pldiff, 1),
    vec4(0, 0,           -2*u.frf*u.frn  / pldiff, 0)
};

What I changed is 1) diff is positive ie. far - near 2) "-1" is "1" 3) "2fn" is "-2fn".

This avoids the ominous Z-flip and seems to work without needing glDepthRange or glDepthFunc. It uses directly A and B as derived above (B at least I showed).

So the "lone 1" is the matrix field that stores the original z as w-value, because x, y and z all need that division to become NDC.

Where you have "0 0 -1 0" for z it should be "0 0 A B" in terms of near and far plane. Z needs to be scaled, just like X and Y, not flipped.