Having only a perspective matrix, how can I accurately obtain the frustrum's origin point?

Question

Considering as a perspective matrix the result of multiplying a view matrix (in opengl GL_MODELVIEW) by a projection matrix (GL_PROJECTION) how can I get the coordenates of the origin (the point of intersection of the lateral fustrum planes, or top of the pyramid)?

persp_mat_4x4 = view_mat_4x4 × proj_mat_4x4

If I had the view_matrix available, I could calculate the origin simply by calculating the inverse of that matrix and reading the values in the 4th column.

orig = (viewmat_inv_4x4[0][3], viewmat_inv_4x4[1][3], viewmat_inv_4x4[2][3])

One option would be to calculate the intersection point of 3 planes of the fustrum. (homogeneous 4d vector representations of planes):

plane_orig_v4 = persp_mat_4x4[3].xyzw
plane_left_v4 = plane_orig_v4 + persp_mat_4x4[0].xyzw
plane_bott_v4 = plane_orig_v4 + persp_mat_4x4[1].xyzw

orig = intersect_plane_plane_plane(plane_orig_v4, plane_left_v4, plane_bott_v4)

But this solution is neither accurate nor efficient. I wanted a simple algorithm to calculate that orig. I hope I have been clear and that you can help me :)

Answer 1

After much analysis, I have found that the best solution is to obtain the intersection of the 3 planes v4 represented by lines 0, 1 and 3 of the matrix.

Here is a pseudo code:

float determinant_m3(float a1, float a2, float a3,
                     float b1, float b2, float b3,
                     float c1, float c2, float c3)
{
    float ans;
    ans = (a1 * (b2 * c3 - b3 * c2) -
           b1 * (a2 * c3 - a3 * c2) +
           c1 * (a2 * b3 - a3 * b2));

    return ans;
}

bool isect_plane_plane_plane_v3(
        vec4 plane_a[4], vec4 plane_b[4], vec4 plane_c[4],
        vec3 r_isect_co[3])
{
    float det;
    det = determinant_m3(UNPACK3(plane_a), UNPACK3(plane_b), UNPACK3(plane_c));

    if (det != 0.0) {
        r_isect_co =
        vec3(plane_b.xyz.cross(plane_c.xyz) * -plane_a[3] +
             plane_c.xyz.cross(plane_a.xyz) * -plane_b[3] +
             plane_a.xyz.cross(plane_b.xyz) * -plane_c[3]) / det;

        return true;
    }
    else {
        return false;
    }
}

bool get_origin_from_projmat(mat4 projmat, vec3 r_origin)
{
    return isect_plane_plane_plane_v3(
                projmat[0], projmat[1], projmat[3], r_origin);
}