In a software I am developing, I use a triangle defined in 3D by three points $A$, $B$, and $C$. In order to display this triangle on the screen, I use a simple projection as follow (this is basically a projection with a bit of perspective):
$$
x^{\prime}_i = -Sx_i\frac{z_i+N}{F-N}
$$
$$
y^{\prime}_i = -Sy_i\frac{z_i+N}{F-N}
$$
$x_i$, $y_i$ and $z_i$ are the coordinates of any points inside the triangle $ABC$. $x^{\prime}_i$ and $y^{\prime}_i$ are the corresponding coordinates on the screen after projection. $S$, $N$ and $F$ are three constants. I know the coordinates of $A$, $B$ and $C$.
Now, I have a projected point $P^{\prime}$ with coordinates $(x^{\prime}_P, y^{\prime}_P)$. I would like to find $P$, the unprojected version of $P^{\prime}$ with coordinates $(x_P, y_P, z_P)$. I know that $P$ is inside the triangle $ABC$, and I know the coordinates of $P^{\prime}$. What I don't know, is the value of $x_P$, $y_P$ and $z_P$. Is there a way to find them?
What I have tried so far, is to use the dot product of $\vec{AP}$ and $\vec{n}$, the normal of the triangle. Because they are perpendicular, the dot product is equal to zero:
$$
\vec{AP}\cdot\vec{n} = 0
$$
which gives me this after decomposition:
$$
n_{x}\left(\frac{x^{\prime}_{P}(F-N)}{S(-z_{P}-N)}-x_{A}\right) + n_{y}\left(\frac{y^{\prime}_{P}(F-N)}{S(-z_{P}-N)}-y_{A}\right) + n_{z}(z_{P}-z_{A})=0
$$
Then I obtain a quadratic equation: $az^{2}_{P} + bz_P + c$ with:
$$
a=-Sn_{z}
$$
$$
b=S(n_{x}x_{A}+n_{y}y_{A}+n_{z}z_{A}-n_{z}N)
$$
$$
c=n_{x}(SNx_{A}+x^{\prime}_{P}(F-N)) + n_{y}(SNy_{A}+y^{\prime}_{P}(F-N)) + SNn_{z}z_{A}
$$
I am a bit stuck here, because most of the time, I have two solutions. Considering the nature of the problem, I was not expecting two solutions, and I am not sure why I get a quadratic equation. Also, I was expecting to get the same result whether I use $A$, $B$ or $C$ as a reference point for the dot product, but I don't. I suppose I made a mistake somewhere, but I just can't find a good way to get the coordinates of $P$ knowing $P^{\prime}$, $A$, $B$, $C$, $A^{\prime}$, $B^{\prime}$, $C^{\prime}$.
Best Answer
As kindly suggested by @Claude in comments, I checked the OpenGL documentation. They use barycentric coordinates to find the value I am looking for.
With the projected triangle $A^{\prime}B^{\prime}C^{\prime}$, and the point $P^{\prime}$ located inside the triangle, we want to find the coordinate $z_P$ that corresponds to the $z$ coordinate of the point $P^{\prime}$ in the unprojected space.
First, we compute the barycentric coordinates in the projected space (page 478 of the documentation): $$ \alpha = \frac{\mathrm{area}(PBC)}{\mathrm{area}(ABC)},\space\space\space \beta = \frac{\mathrm{area}(PAC)}{\mathrm{area}(ABC)},\space\space\space \gamma = \frac{\mathrm{area}(PBA)}{\mathrm{area}(ABC)}\tag{1} $$
Second, the documentation gives the following equation that interpolates any parameter $f$ in the triangle while taking into account the perspective. In the manual, the equation has ref $(14.9)$ page 479. $$ f_P = \frac{\alpha f_A\mathbin{/}w_A+\beta f_B\mathbin{/}w_B+\gamma f_C\mathbin{/}w_C}{\alpha\mathbin{/}w_A+\beta\mathbin{/}w_B+\gamma\mathbin{/}w_C} \tag{2} $$
In my case, $f_A$, $f_B$, $f_C$ and $f_P$ are the $z$ coordinates of the points $A$, $B$, $C$ and $P$ respectively. The values $w_A$, $w_B$ and $w_C$ are the clip coordinates. To get these coordinates, I had to change multiple things in the projection method. Using the information from here, I reconstructed a projection matrix like this:
$$ M_{proj} = \begin{pmatrix} \frac{N}{S} & 0 & 0 & 0 \\ 0 & \frac{N}{S} & 0 & 0 \\ 0 & 0 & -\frac{F + N}{F - N} & -2\frac{FN}{F-N} \\ 0 & 0 & -1 & 0 \end{pmatrix} \tag{3}$$
with $F$ the far plan, $N$ the near plan, $S$ the tangent of the field of view (in my case 45º but many people use 60º). The projected coordinates is obtained by multiplying the matrix and the unprojected coordinates of $A$, $B$ and $C$ and using $1$ for the 4th element of the matrix: $$ A^{\prime}=\begin{pmatrix} x^{\prime}_A \\ y^{\prime}_A \\ z^{\prime}_A \\ w_A \end{pmatrix}=M_{proj} \begin{pmatrix} x_A \\ y_A \\ z_A \\ 1 \end{pmatrix} \tag{4}$$
Now we can compute $z_P$ using the equation $(2)$, the barycentric weights found in $(1)$ and the clip coordinates found in $(4)$.