The question is a slight extension of the question given here: Check whether a point is within a 3D Triangle
There is an elegant solution to this given by W. Heidrich, Journal of Graphics, GPU, and Game Tools,Volume 10, Issue 3, 2005.
Let $\vec{u}=P_2-P_1$, $\vec{v}=P_3-P_1$, $\vec{n}=\vec{u}\times\vec{v}$,
$\vec{w}=P-P_1$. We then have directly the barycentric coordinates of the projection $P'$ of $P$ onto $T$ as
- $\gamma=\left[(\vec{u}\times\vec{w})\cdot\vec{n}\right]/\vec{n}^2$
- $\beta=\left[(\vec{w}\times\vec{v})\cdot\vec{n}\right]/\vec{n}^2$
- $\alpha=1-\gamma-\beta$
The coordinates of the projected point is
- $P'=\alpha P_1+\beta P_2 +\gamma P_3$
The point $P'$ lies inside $T$ if
- $0\leq\alpha\leq 1$,
- $0\leq\beta\leq 1$, and
- $0\leq\gamma\leq 1$.
I'll start by looking about how you can describe the projection of an arbitrary point $Q$ onto your plane, then derive a matrix notation from that.
Projecting a single point
The point $Q$ and its projection form a line which is perpendicular to the plane. All these perpendicular lines meet in a point $F=(A,B,C,0)^T$, which is the point at infinity orthogonal to your plane. A generic point on the line connecting $Q$ and $F$ can be described as a linear combination $Q+\lambda F$. (This inhomogenous formulation using a single parameter excludes $F$ itself, which is all right since $F$ will not be the projected point in any case.) Now you are looking for the point where this line connecting $F$ and $Q$ intersects the plane, i.e. $\left<Q+\lambda F,h\right>=0$, where $h=(A,B,C,D)^T$ is the vector describing your plane. This you can solve for $\lambda$:
$$\left<Q+\lambda F,h\right>=\left<Q,h\right> + \lambda\,\left<F,h\right>=0\\
\lambda=-\frac{\left<Q,h\right>}{\left<F,h\right>}$$
To avoid the division, you can also use a multiple of that point:
$$Q+\lambda F\;\sim\;
\left<F,h\right>\,Q - \left<Q,h\right>F
$$
This is the projection of a generic point $Q$.
Finding the matrix
If you assume the coordinates of $Q$ to be variables, then the coordinates of the result will be linear in these variables, so you can interpret the whole operation as a linear map and therefore write it as a matrix. You can also find the formula of that matrix like this:
$$
\left<F,h\right>\,Q - \left<Q,h\right>F =
\left<F,h\right>\mathbb 1\cdot Q - Fh^T\cdot Q
= \left(\left<F,h\right>\mathbb 1 - Fh^T\right)Q \\
= \begin{pmatrix}
B^{2} + C^{2} & -A B & -A C & -A D \\
-A B & A^{2} + C^{2} & -B C & -B D \\
-A C & -B C & A^{2} + B^{2} & -C D \\
0 & 0 & 0 & A^{2} + B^{2} + C^{2}
\end{pmatrix}\cdot Q
$$
If anyone reading this post should be interested in central instead of orthogonal projection, simply use that center of projection as $F$, and the above formula will help you compute the matrix for that as well.
Best Answer
Let $u = \vec{OA}, w = \vec{OB}$. Let $b = \vec {OC}$, and $x = \vec{OP}$.
Let $S = [u \quad w \quad v]$. We note that $S^{-1}x$ is simply the projection of $S^{-1}b$ onto the $xy$-plane.
However, projecting onto the $xy$-plane is easy. So, all together, we can simply compute $$ x = P_{uw}b = S P_{xy} S^{-1}b = S \pmatrix{1&0&0\\0&1&0\\0&0&0} S^{-1}b $$