What is the projection matrix of reverse (Byzantine) perspective? I would like to construct a projection matrix for reverse perspective.
I'm using OpenGL and tried to modify concepts from this excelent tutorial.
I came up with:
$$
\begin{bmatrix}
2\frac{(\text{near}-M)}{\text{right}-\text{left}} &  0  & \frac{\text{right} + \text{left}}{\text{right}-\text{left}} & M\frac{\text{right} + \text{left}}{\text{right}-\text{left}}\\
0  &  2\frac{(\text{near}-M)}{\text{top}-\text{bottom}} & \frac{\text{top} + \text{bottom}}{\text{top}-\text{bottom}} & M\frac{\text{top} + \text{bottom}}{\text{top}-\text{bottom}}\\
0 & 0 & \frac{2M-\text{near}-\text{far}}{\text{far}-\text{near}} & \frac{2(M-\text{far}) \cdot \text{near}}{\text{far}-\text{near}} + M\\
0 & 0 & -1 & -M
\end{bmatrix}
$$
$$ \text{Variables define viewing frustum. }M \text{ is the tip; point } (0,0,-M) \text{ in eye space.}$$
To avoid ambiguity near point M I put far plane closer to the camera. To debug this, I inverted points from NDC [-1,1] cube and they seem to fit the desired frustum. Sadly my models don't render correctly when using this matrix. They don't show up when they clearly are in the viewing frustum. I verify this by rendering the scene from other cameras using regular perspective or orthographic projection.
To clarify, what kind of perspective I'm talking about, here are links:


*

*https://en.wikipedia.org/wiki/Reverse_perspective

*True Reverse Perspective
Code in c++ using glm:
template <typename T>
auto frustum(T left, T right, T bottom, T top, T near, T M, T far) {
    glm::tmat4x4<T, glm::defaultp> result(0);
    result[0][0] = (static_cast<T>(2) * (near - M)) / (right - left);
    result[2][0] = (right + left) / (right - left);
    result[3][0] = (right + left) / (right - left) * M;

    result[1][1] = (static_cast<T>(2) * (near - M)) / (top - bottom);
    result[2][1] = (top + bottom) / (top - bottom);
    result[3][1] = (top + bottom) / (top - bottom) * M;

    result[2][2] = ((static_cast<T>(2) * M) - near - far) / (far - near);
    result[3][2] = static_cast<T>(2) * (M - far) * near / (far - near) + M;

    result[2][3] = static_cast<T>(-1);
    result[3][3] = static_cast<T>(-M);
    return result;
}

 A: Instead of just plopping some random transformation matrix here, let's investigate how one constructs these odd projection matrices instead.
We have point $(x_i, y_i, z_i)$ in world coordinates, and we wish to project it to a 2D plane, say at $(\chi_i, \gamma_i, -1)$, with the projective point at $(0, 0, z_p)$.  Since these three points are always collinear, on the same line, we have
$$\left\lbrace ~ \begin{aligned}
\chi_i &= 0 + \lambda (x_i - 0) \\
\gamma_i &= 0 + \lambda (y_i - 0) \\
-1 &= z_p + \lambda (z_i - z_p) \\
\end{aligned} \right . \quad \iff \quad \left\lbrace ~ \begin{aligned}
\chi_i &= \displaystyle \frac{z_p + 1}{z_p - z_i} x_i \\
\gamma_i &= \displaystyle \frac{z_p + 1}{z_p - z_i} y_i \\
\end{aligned} \right .$$
OpenGL projection works thus:
$$\left[\begin{matrix} \chi_i \\ \gamma_i \\ \zeta_i \end{matrix}\right] =
\left[\begin{matrix} X_i / W_i \\ Y_i / W_i \\ Z_i / W_i \end{matrix}\right], \quad
\left[\begin{matrix} X_i \\ Y_i \\ Z_i \\ W_i \end{matrix}\right] =
\left[\begin{matrix} x_x & y_x & z_x & t_x \\ x_y & y_y & z_y & t_y \\ x_z & y_z & z_z & t_z \\ x_w & y_w & z_w & t_w \end{matrix}\right]
\left[\begin{matrix} x_i \\ y_i \\ z_i \\ 1 \end{matrix}\right]$$
where the visible region is $-1 \le X_i \le +1$, $-1 \le Y_i \le +1$, $-1 \le Z_i \le +1$.  (The last one is the reason I chose the projective plane to be $z = -1$.)  These expand to
$$\left\lbrace ~ \begin{aligned}
\chi_i &= \displaystyle \frac{ x_x x_i + y_x y_i + z_x z_i + t_x }{x_w x_i + y_w y_i + z_w z_i + t_w } \\
\gamma_i &= \displaystyle \frac{ x_y x_i + y_y y_i + z_y z_i + t_y }{x_w x_i + y_w y_i + z_w z_i + t_w } \\
\zeta_i &= \displaystyle \frac{ x_z x_i + y_z y_i + z_z z_i + t_z }{x_w x_i + y_w y_i + z_w z_i + t_w } \\
\end{aligned} \right.$$
To match the projection we want, remembering that $\zeta_i = -1$, we need
$$\left\lbrace ~ \begin{aligned}
X_i &= (z_p + 1) x_i \\
Y_i &= (z_p + 1) y_i \\
Z_i &= z_i - z_p \\
W_i &= z_p - z_i \\
\end{aligned} \right .$$
although we can apply any nonzero multiplier (say, $\beta$) to all four without affecting the result. One possible solution is, listing only nonzero coefficients for the projection matrix,
$$\begin{aligned}
x_x &= z_p + 1, \\
y_y &= z_p + 1, \\
z_z &= 1, \\
t_z &= -z_p, \\
z_w &= -1, \\
t_w &= z_p \\
\end{aligned} \quad \iff \quad \left[\begin{matrix}
z_p + 1 & 0 & 0 & 0 \\
0 & z_p + 1 & 0 & 0 \\
0 & 0 & 1 & -z_p \\
0 & 0 & -1 & z_p \\
\end{matrix} \right]$$
Before applying perspective projection, it is useful to translate and scale the coordinates.  This is done thus:
$$\left[\begin{matrix} x_{new} \\ y_{new} \\ z_{new} \\ 1 \end{matrix} \right] =
\left[\begin{matrix} S_x & 0 & 0 & T_x \\ 0 & S_y & 0 & T_y \\ 0 & 0 & S_z & T_z \\ 0 & 0 & 0 & 1 \end{matrix}\right] \left[\begin{matrix} x_{old} \\ y_{old} \\ z_{old} \\ 1 \end{matrix}\right]$$
To combine transformations, we multiply the matrices together, oldest rightmost, newest leftmost.  Combining these two matrices yields the projection matrix
$$\left[\begin{matrix}
S_x (z_p + 1) & 0 & 0 & T_x (z_p + 1) \\
0 & S_y (z_p + 1) & 0 & T_y (z_p + 1) \\
0 & 0 & S_z & T_z - z_p \\
0 & 0 & -S_z & z_p - T_z \\
\end{matrix}\right]$$
where $S_x$, $S_y$, and $S_z$ are the coordinate axis scale factors, $z_p+1$ is the distance from projection plane to projective point prior to scaling, and $(T_x, T_y, T_z)$ is the translation vector (movement) in scaled coordinates.
With scaling but no translation, and using $d = z_p + 1$ for the distance from projection plane to projective point prior to scaling, the projection matrix is
$$\left[\begin{matrix}
S_x d & 0 & 0 & 0 \\
0 & S_y d & 0 & 0 \\
0 & 0 & S_z & 1 - d \\
0 & 0 & -S_z & d - 1 \\
\end{matrix}\right]$$
Remembering that OpenGL only displays the projected coordinates within the axis-aligned cube $(-1,-1,-1)-(+1,+1,+1)$, we can find out where the eight vertices of the view frustum (truncated pyramid with a rectangular base) in world coordinates are.  Like OP shows, these coordinates are useful for us humans to define the projection parameters.  That last projection matrix implements
$$\left\lbrace ~ \begin{aligned}
X &= d S_x x \\
Y &= d S_y y \\
Z &= S_z z - d + 1 \\
\end{aligned}\right . , \quad \begin{aligned}
-1 \le & X \le +1 \\
-1 \le & Y \le +1 \\
-1 \le & Z \le +1 \\
\end{aligned}$$
so if we use $\text{width}$ for the range of visible $x$ world coordinates ($\text{left} = -\text{width}/2$, $\text{right} = \text{width}/2$), $\text{height}$ for the range of visible $y$ world coordinates ($\text{top} = -\text{height}/2$, $\text{bottom} = \text{height}/2$), $\text{near}$ for the minimum visible $z$ coordinates, and $\text{far}$ for the maximum visible $z$ coordinates, we can solve for the four constants:
$$\left\lbrace ~ \begin{aligned}
d &= \displaystyle \frac{2 \text{ far}}{\text{far} - \text{near}} \\
S_x &= \displaystyle \frac{\text{far} - \text{near}}{\text{far} \text{ width}} \\
S_y &= \displaystyle \frac{\text{far} - \text{near}}{\text{far} \text{ height}} \\
S_z &= \displaystyle \frac{2}{\text{far} - \text{near}} \\
\end{aligned} \right .$$
which means the projection matrix becomes
$$\left[\begin{matrix}
\displaystyle \frac{2}{\text{ width}} & 0 & 0 & 0 \\
0 & \displaystyle \frac{2}{\text{ height}} & 0 & 0 \\
0 & 0 & \displaystyle \frac{2}{\text{far}-\text{near}} & \displaystyle \frac{\text{near} + \text{far}}{\text{near} - \text{far}} \\
0 & 0 & \displaystyle \frac{2}{\text{near}-\text{far}} & \displaystyle \frac{\text{far} + \text{near}}{\text{far} - \text{near}} \\
\end{matrix}\right]$$
This same procedure works, if you know how to map an arbitrary point $(x, y, z)$ to a projection plane $(\chi, \gamma, -1)$.  This mapping is usually a line through some fixed projection point $(x_p, y_p, z_p)$, making all three points collinear.
A: I'm not sure that I can help with the reasoning, but I use the following scheme for a smooth transition from the classical perspective to the reverse perspective.
For Ws = We = 1, we obtain the canonical matrix orthogonal to the projection. With Ws = near and We = far we get the classical perspective projection matrix. With Ws = far and We = near, we obtain a matrix for the inverse perspective.
I used such an algorithm to implement a camera with a reverse perspective effect for Blender 3D.

Further, unicam is a parameter that varies from 0 to 1. If unicam = 0 then the perspective effect. If unicam = 1 then the effect of reverse perspective.

$
width = right-left \\
height = top-bottom \\
depth = far-near \\
Ws = near+depth*unicam \\
We = far-depth*unicam \\
A = \frac{-(We+Ws)}{depth}\\
B = We+A*far\\
C = \frac{-(We-Ws)}{depth} \\
D = We+C*far \\
Px = 2*Ws/width \\
Pz = -C*(right + left)/width \\
Pw = D*(right + left)/width \\
Qy = 2*Ws/height \\
Qz = -C*(top + bottom)/height \\
Qw = D*(top + bottom)/height \\
$
$
\begin{bmatrix}
Px& 0& Pz& Pw \\
0& Qy& Qz& Qw \\
0& 0& A& B  \\
0& 0& C& D \\
\end{bmatrix}
$

Experiments in Reverse Perspective, Written by Paul Bourke
Reverse perspective rendering
Video examples

Code in c++
bool reverse_perspective(float m[4][4], const float left, const float right, const float bottom, const float top, const float near, const float far, const float unicam)
{
    float A, B, C, D, Ws, We, width, height, depth;
    width = right-left;
    height = top-bottom;
    depth = far-near;

    Ws = near+depth*unicam;
    We = far-depth*unicam;

    C = -(We-Ws)/depth;
    D = We+C*far;
    A = -(We+Ws)/depth;
    B = We+A*far;

    m[0][0] = 2*Ws/width;
    m[1][0] = 0.0f;
    m[2][0] = -C*(right + left)/width;
    m[3][0] = D*(right + left)/width;

    m[0][1] = 0.0f;
    m[1][1] = 2*Ws/height;
    m[2][1] = -C*(top + bottom)/height;
    m[3][1] = D*(top + bottom)/height;

    m[0][2] = 0.0f;
    m[1][2] = 0.0f;
    m[2][2] = A;
    m[3][2] = B;

    m[0][3] = 0.0f;
    m[1][3] = 0.0f;
    m[2][3] = C;
    m[3][3] = D;

    return true;
}

