Let me try to make my comment into a complete answer.
The general idea is to build a linear system using those 6 point pairs and solve for the desired 12 unknowns. You may find this paper[1] useful if you want to learn the theoretical background, and you can find some more details in my blog post[2].
Here is my solution using Mathematica.
In the first step I will build a symbolic array which contains all variables.
v = Table[Symbol["m" <> ToString[i]], {i, 0, 15}];
M = ArrayReshape[v, {4, 4}]
Next we write down the coordinates for all vertices.
k = f/ n;
X = {
{l, t, -n, 1},
{l, b, -n, 1},
{r, b, -n, 1},
{r, t, -n, 1},
{k l, k t, -f, 1},
{k l, k b, -f, 1},
{k r, k b, -f, 1},
{k r, k t, -f, 1}};
U = {
{-1, 1, -1},
{-1, -1, -1},
{1, -1, -1},
{1, 1, -1},
{-1, 1, 1},
{-1, -1, 1},
{1, -1, 1},
{1, 1, 1}};
where l
,r
,b
,t
,n
,and f
denote left, right, bottom, top, near, and far planes respectively,
and 'X' is the coordinates in world space, U
is the coordinates in NDC.
Then next we implement a function Generate
to generate an equation for each pair of coordinates:
Generate[{i_, j_}] := M[[j]].X[[i]] - U[[i, j]] M[[4]].X[[i]] == 0;
Because we do not care about the third row,
we only solve for the unknowns in partialv
:
partialv = {m0, m1, m2, m3, m4, m5, m6, m7, m12, m13, m14, m15};
Now we can build our linear system and solve it:
solve[subset_] := Solve[
Map[Generate, Table[{i, j}, {i, subset}, {j, 1, 2}], {2}] // Flatten,
partialv] // Simplify
Let me explain it a little bit.
Because only 6 pairs of points are given, we select a subset
of 6 elements from the 8 pairs of points.
Then we call Generate
function on this subset
and the first two rows of matrix(we do not care about the third row. And the fourth row is the "homogeneous" part, which is also ignored).
Finally we enumerate all subsets from these 8 pairs of points, and try to find a solution:
solutions = Map[solve, Subsets[{1, 2, 3, 4, 5, 6, 7, 8}, {6}]];
We can view one solution by:
MatrixForm[ArrayReshape[v /. solutions[[1]], {4, 4}]]
$$
\left(
\begin{array}{cccc}
\text{m0} & 0 & \frac{\text{m0} (l+r)}{2 n} & 0 \\
0 & \frac{\text{m0} (l-r)}{b-t} & \frac{\text{m0} (l-r) (b+t)}{2 n (b-t)} & 0 \\
\text{m8} & \text{m9} & \text{m10} & \text{m11} \\
0 & 0 & \frac{\text{m0} (l-r)}{2 n} & 0 \\
\end{array}
\right)
$$
The m0
is the "unknown scalar factor" which is mentioned in problem statement.
In fact if we let m0
equals 2 n/(r - l)
,
we would get the conventional projection matrix used in OpenGL:
% /. m0 -> 2 n/(r - l)
% // Simplify // MatrixForm
$$
\left(
\begin{array}{cccc}
\frac{2 n}{r-l} & 0 & \frac{l+r}{r-l} & 0 \\
0 & -\frac{2 n}{b-t} & -\frac{b+t}{b-t} & 0 \\
\text{m8} & \text{m9} & \text{m10} & \text{m11} \\
0 & 0 & -1 & 0 \\
\end{array}
\right)
$$
[1] Heckbert, Paul S. Fundamentals of texture mapping and image warping. MS thesis. University of California, Berkeley, 1989.
[2] Inferring projective mappings Inferring projective mappings | linearconstraints.net