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][3]
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][3]
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[0]], points[idx_pair[1]]])
ax.plot(line[0], line[2], line[1], "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[0], camera_pos[2], camera_pos[1], 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],
]
)
def get_rotation_mat_x(angle_rad):
s = np.sin(angle_rad)
c = np.cos(angle_rad)
return np.array(
[[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]], dtype=float
)
def get_rotation_mat_y(angle_rad):
s = np.sin(angle_rad)
c = np.cos(angle_rad)
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]],
[0, 1, 0, position[1]],
[0, 0, 1, position[2]],
[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(
get_rotation_mat_x(-pitch_rad), get_rotation_mat_y(-yaw_rad)
)
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][3]
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][3]
# 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()
