Calculating Directional Shadow Map using Camera Frustum

I'm trying to calculate the 8 corners of the view frustum so that I can use them to calculate the ortho projection and view matrix needed to calculate shadows based on the camera's position. Currently, I'm not sure how to convert the frustum corners from local space into world space. Currently, I have calculated the frustum corners in local space as follows (correct me if I'm wrong):

float tan = 2.0 * std::tan(m_Camera->FOV * 0.5);
float nearHeight = tan * m_Camera->Near;
float nearWidth = nearHeight * m_Camera->Aspect;
float farHeight = tan * m_Camera->Far;
float farWidth = farHeight * m_Camera->Aspect;

Vec3 nearCenter = m_Camera->Position + m_Camera->Forward * m_Camera->Near;
Vec3 farCenter = m_Camera->Position + m_Camera->Forward * m_Camera->Far;

Vec3 frustumCorners = {
nearCenter - m_Camera->Up * nearHeight - m_Camera->Right * nearWidth, // Near bottom left
nearCenter + m_Camera->Up * nearHeight - m_Camera->Right * nearWidth, // Near top left
nearCenter + m_Camera->Up * nearHeight + m_Camera->Right * nearWidth, // Near top right
nearCenter - m_Camera->Up * nearHeight + m_Camera->Right * nearWidth, // Near bottom right

farCenter - m_Camera->Up * farHeight - m_Camera->Right * nearWidth, // Far bottom left
farCenter + m_Camera->Up * farHeight - m_Camera->Right * nearWidth, // Far top left
farCenter + m_Camera->Up * farHeight + m_Camera->Right * nearWidth, // Far top right
farCenter - m_Camera->Up * farHeight + m_Camera->Right * nearWidth, // Far bottom right
};

How do I move these corners into world space?

Update:

I'm still not sure if what I'm doing is right. I've also attempted to build the ortho projection by looping through the frustum corners and getting the min and max x,y,z coordinates. Then simply setting the values of the projection as:

left = minX
right = maxX
top = maxY
botton = minY
near = - maxZ
far = -minZ

I've searched on the internet but all the tutorials use hard coded values so the shadow maps aren't applicable to an open world but a restricted portion of the scene. Any help? Pseudocode is preferred as my linear algebra skills (and reading skills) aren't that great

Not sure why you are doing this "by hand", but it should be pretty straight forward if you use matrices for that.

Assume your models are all already transformed to world space, then you usually do the following to get them to clip space:

$$P_c = M_{vc} \cdot M_{wv}\cdot P_w$$

$$\P_c\$$ is the projected point in clip space. $$\M_{vc}\$$ is the projection matrix that transforms from view space to clip space, $$\M_{wv}\$$ the world to view space transformation matrix and $$\P_w\$$ is a point in world space. Now what you want to do is the opposite transformation which simply is:

$$P_w = M_{wv}^{-1} \cdot M_{vc}^{-1} \cdot P_c$$

The 8 corners of your view frustum are just the edge points of a cube in clip space (the space after projection - also see this link and the learnOpenGL coordinate system tutorial). For OpenGL, the 8 points have x-,y- and z- values that are either -1 or 1.

So all you need are the inverse matrices and the 8 corner points of your clip space. In case you don't use an existing linear algebra library and don't know how to invert a 4x4 matrix, check this StackOverflow question. You can actually simplify this formula for the projection matrix $$\M_{vc}\$$ since it contains a lot of zeros.

Important Update: The points in clip space need to be 4d vectors with $$\w=1\$$ and after the application of the inverse projection, all points have to be divided by the resulting $$\w\$$ component

Additionally, you can avoid calculating the inverse with the general formula for the matrix $$\M_{wv}\$$ since it usually is just a composition of a rotation and a subsequent translation:

$$M_{wv} = T \cdot R$$

So

$$M_{wv}^{-1} = R^{-1} \cdot T^{-1}$$

The inverse of a translation is just subtracting the translation from all points. Don't create the inverse matrix of $$\T\$$, just subtract the part that corresponds to the translations from your points. For a rotation, it's inverse is equal to its transposed matrix:

$$R^{-1} = R^T$$

EDIT:

Even though clip space ranges from -1 to 1 in all directions, you might need to use 0 as the lower or upper limit for the z-direction of the view frustum points in clip space since if I remember correctly, the actual "screen position" is at z=0.

UPDATE

As you asked for some code in the comments, I wrote a small Python script that shows you everything about the solution I proposed. The full script is at the end of this answer and includes basically all calculations. Only for the inverse and matrix multiplication, I used NumPy. Regarding the calculation of the inverse matrix, check the link in the previous part of my answer.

You can copy the script and play around with the different camera parameters at the beginning of the script to see how they affect the result. It should run in any Python 3 environment with NumPy.

Now let's go to the important substeps. It is actually not that much. First, we define the view frustum in clip space. As I said before, it is just a cube with coordinate ranges -1 to 1. The important part here is, that you use 4d coordinates, where the w component is 1:

points_clip = np.array(
[
[-1, -1, -1, 1],
[ 1, -1, -1, 1],
[-1,  1, -1, 1],
[ 1,  1, -1, 1],
[-1, -1,  1, 1],
[ 1, -1,  1, 1],
[-1,  1,  1, 1],
[ 1,  1,  1, 1],
],
dtype=float,
)

The corresponding plot looks like this: Now we calculate the perspective matrix and the world to view space matrix:

M_wv = get_world_to_view_matrix(camera_pitch, camera_yaw, camera_position)
M_vc = get_perspective_mat(field_of_view, z_near_plane, z_far_plane, aspect_ratio)

If you use OpenGl 3.3 or higher, you should have them already, since you need them in your shaders. If not, look at the corresponding function definitions in the full script (further references: learnOpenGL - Transformations, OGLdev - Perspective Projection).

Now we calculate the inverse matrices that we will need. Optionally, we can multiply them to get a single transformation matrix:

M_vw = np.linalg.inv(M_wv)
M_cv = np.linalg.inv(M_vc)
# M_cw = np.matmul(M_vw, M_cv) # direct clip to world transformation

Note that the indices of the result matrices are switched since they transform in the opposite direction. Now, all we need to do is to multiply each point with the corresponding transformation matrix, and very important, divide by the resulting w component afterward. I forgot to mention it in my original answer (I just found that out myself while I wrote the script ;) ). I transformed the points to world space and also to view space so that you can see the intermediate results:

points_view = []
points_world = []
for i in range(8):
points_view.append(np.matmul(M_cv, points_clip[i]))
points_view[i] = points_view[i] / points_view[i]
points_world.append(np.matmul(M_vw, points_view[i]))
# alternative
# points_world.append(np.matmul(M_cw, points_clip[i]))
# points_world[i] = points_world[i] / points_world[i]

Here is a plot of the frustum after the transformation to view space: The green dot is the actual camera position. Here is the final result in world space: I hope that helps you to understand the approach. If not, copy the script and run it. There are also some alternative calculations included, that you can use (In case you get an error, you missed probably to include some other comment). If there remain any questions, don't hesitate to ask, but consider creating a new question referring to this one for more complex questions to avoid extensive chatting in the comments.

How do I move these corners into world space?

In case you want to use your own calculated points and want to transform them from view space to world space, all you need to do is to multiply them with the inverse world-to-view matrix. This is the same as this part of the script:

for i in range(8):
...
points_world.append(np.matmul(M_vw, points_view[i]))

Just use your own points as points_view.

Some further remarks

• If you run the script, don't be confused if some frustum plots look distorted. This is because the axes are distorted. Resize the window manually, until the axes are identically scaled.
• In case you look into the "render function" don't get confused, that the y and z-axis are swapped. This is just to get the z-axis into the horizontal plane.
• Be aware, that you might not be able to copy/translate the code directly to use it with your program. In case you defined some coordinate directions differently, you might need to add or modify some transformations.

Full Python script

import numpy as np
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
import matplotlib.pyplot as plt

# setup --------------------------------------------------------------------------------

camera_position = [3, 0, 1]
camera_yaw = 20
camera_pitch = 30

field_of_view = 70
z_near_plane = 0.5
z_far_plane = 3
aspect_ratio = 16 / 9

# functions ----------------------------------------------------------------------------

def render_frustum(points, camera_pos, ax):
line_indices = [
[0, 1],
[0, 2],
[0, 4],
[1, 3],
[1, 5],
[2, 3],
[2, 6],
[3, 7],
[4, 5],
[4, 6],
[5, 7],
[6, 7],
]
for idx_pair in line_indices:
line = np.transpose([points[idx_pair], points[idx_pair]])
ax.plot(line, line, line, "r")
ax.set_xlim([-5, 5])
ax.set_ylim([-5, 5])
ax.set_zlim([-5, 5])
ax.set_xlabel("x")
ax.set_ylabel("z")
ax.set_zlabel("y")
ax.plot([-5, 5], [0, 0], [0, 0], "k")
ax.plot([0, 0], [-5, 5], [0, 0], "k")
ax.plot([0, 0], [0, 0], [-5, 5], "k")
if camera_pos is not None:
ax.scatter(
camera_pos, camera_pos, camera_pos, marker="o", color="g", s=30
)

def get_perspective_mat(fov_deg, z_near, z_far, aspect_ratio):
fov_rad = fov_deg * np.pi / 180
f = 1 / np.tan(fov_rad / 2)

return np.array(
[
[f / aspect_ratio, 0, 0, 0],
[0, f, 0, 0],
[
0,
0,
(z_far + z_near) / (z_near - z_far),
2 * z_far * z_near / (z_near - z_far),
],
[0, 0, -1, 0],
]
)

return np.array(
[[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]], dtype=float
)

return np.array(
[[c, 0, s, 0], [0, 1, 0, 0], [-s, 0, c, 0], [0, 0, 0, 1]], dtype=float
)

def get_translation_mat(position):
return np.array(
[
[1, 0, 0, position],
[0, 1, 0, position],
[0, 0, 1, position],
[0, 0, 0, 1],
],
dtype=float,
)

def get_world_to_view_matrix(pitch_deg, yaw_deg, position):
pitch_rad = np.pi / 180 * pitch_deg
yaw_rad = np.pi / 180 * yaw_deg

orientation_mat = np.matmul(
)
translation_mat = get_translation_mat(-1 * np.array(position, dtype=float))
return np.matmul(orientation_mat, translation_mat)

# script -------------------------------------------------------------------------------

points_clip = np.array(
[
[-1, -1, -1, 1],
[ 1, -1, -1, 1],
[-1,  1, -1, 1],
[ 1,  1, -1, 1],
[-1, -1,  1, 1],
[ 1, -1,  1, 1],
[-1,  1,  1, 1],
[ 1,  1,  1, 1],
],
dtype=float,
)

M_wv = get_world_to_view_matrix(camera_pitch, camera_yaw, camera_position)
M_vc = get_perspective_mat(field_of_view, z_near_plane, z_far_plane, aspect_ratio)

M_vw = np.linalg.inv(M_wv)
M_cv = np.linalg.inv(M_vc)
# M_cw = np.matmul(M_vw, M_cv) # direct clip to world transformation

# alternative:
# M_wc = np.matmul(M_vc, M_wv)
# M_cw = np.linalg.inv(M_wc)

points_view = []
points_world = []
for i in range(8):
points_view.append(np.matmul(M_cv, points_clip[i]))
points_view[i] = points_view[i] / points_view[i]
points_world.append(np.matmul(M_vw, points_view[i]))
# alternative
# points_world.append(np.matmul(M_cw, points_clip[i]))
# points_world[i] = points_world[i] / points_world[i]

# plot everything ----------------------------------------------------------------------

plt.figure()
ax_clip_space = plt.gca(projection="3d")
render_frustum(points=points_clip, camera_pos=None, ax=ax_clip_space)
ax_clip_space.set_title("view frustum in clip space")

plt.figure()
ax_view = plt.gca(projection="3d")
render_frustum(points=points_view, camera_pos=[0, 0, 0], ax=ax_view)
ax_view.set_title("view frustum in view space")

plt.figure()
ax_world = plt.gca(projection="3d")
render_frustum(points=points_world, camera_pos=camera_position, ax=ax_world)
ax_world.set_title("view frustum in world space")

plt.show()
• I'm not sure what you mean by hand? My linear algebra is a bit weak (even though I've written my own math library for this project). I can only do it when I can mentally visualise things. That's why I calculated it that way. Is it possible for you to make that into some kind of pseudocode which is easier for me to understand? Jun 4 '20 at 17:43
• @Sammi3 With "by hand" I just meant, why you don't use matrices and vectors ;) --- I can try to give you a small example, but first I need to know what you are struggling with. If you would use a LinAlg lib like GLM, Eigen, or an own implementation, transforming those equations into code shouldn't be too hard. So let me know what you don't understand :) Jun 5 '20 at 6:07
• I use a custom math library but it is basically based off GLM. It doesn't have all GLM's feature but it has all the main vector and matrix methods. The main part that I don't understand is what vectors I should transform. So I'm assuming that the inverse view projection matrix is equivalent to doing: Inverse(view) * Inverse(projection). Do I then multiply that by the vectors that I created in the original question? Jun 6 '20 at 13:30
• @Sammi3 Updated my answer with a full python script. Look into the Update section. I hope that helps. Jun 7 '20 at 17:48
• Thank you! The visual aids were perfect btw, helped me see what was actually happening! Jun 7 '20 at 20:08