Because this is a projection from 3D down to 2D, it's not invertible. So we can't recover the triangle \$\triangle PMN\$ that was used to create the original projection. But we can find a member of a family of triangles that fits that bill.

For simplicity, let's choose the member of that family that lies in the same plane as our mesh triangle \$\triangle ABC\$. That means our points \$PMN\$ are linear combinations of the points \$ABC\$.

We can start with point \$P\$, which sits at the point in the plane where \$u = v = 0\$.

Let's define...

$$u_1 = U_b - U_a \quad u_2 = U_c - U_a\\ v_1 = V_b - V_a \quad v_2 = V_c - V_a$$

Now we'll solve...

$$\begin{align} \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix} \begin{bmatrix}s\\t\end{bmatrix} &=\begin{bmatrix}-U_a\\-V_a\end{bmatrix}\\ \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix}^{-1} \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix} \begin{bmatrix}s\\t\end{bmatrix} &=\begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix}^{-1}\begin{bmatrix}-U_a\\-V_a\end{bmatrix}\\ \begin{bmatrix}s\\t\end{bmatrix} &=\frac 1 {u_1v_2 - u_2v_1} \begin{bmatrix}v_2 & -u_2\\-v_1 & u_1\end{bmatrix} \begin{bmatrix}-U_a\\-V_a\end{bmatrix} \end{align}$$

Now we have

$$P = A + s(B-A) + t(C-A)$$

Next we want to find the point where \$(u, v) = (1, 0)\$, which will give us point \$M\$, and \$(u, v) = (0, 1)\$, which gives us point \$N\$, completing our texture space basis. We can proceed the same way...

$$\begin{bmatrix}s\\t\end{bmatrix} =\frac 1 {u_1v_2 - u_2v_1} \begin{bmatrix}v_2 & -u_2\\-v_1 & u_1\end{bmatrix} \begin{bmatrix}1 -U_a\\-V_a\end{bmatrix}\\ M = A + s(B - A) + t (C-A)$$

$$\begin{bmatrix}s\\t\end{bmatrix} =\frac 1 {u_1v_2 - u_2v_1} \begin{bmatrix}v_2 & -u_2\\-v_1 & u_1\end{bmatrix} \begin{bmatrix}-U_a\\1 - V_a\end{bmatrix}\\ N = A + s(B - A) + t (C-A)$$

In my tests with random triangles, using single-precision floats throughout, this correctly reproduces the UV mapping to within 0.001, which I'd attribute to rounding errors along the way. If you have more sensible triangles, and use double-precision intermediates, you could probably reduce the deviation even further.

Because this is a projection from 3D down to 2D, it's not invertible. So we can't recover the triangle \$\triangle PMN\$ that was used to create the original projection. But we can find a member of a family of triangles that fits that bill.

For simplicity, let's choose the member of that family that lies in the same plane as our mesh triangle \$\triangle ABC\$. That means our points \$PMN\$ are linear combinations of the points \$ABC\$. That is, they'll be of the form:

$$A + s (B - A) + t (C - A)$$

...for some real numbers \$s\$ and \$t\$. Taking a look at the UV coordinates of our points in \$\triangle ABC\$, we can define:

$$u_1 = U_b - U_a \quad u_2 = U_c - U_a\\ v_1 = V_b - V_a \quad v_2 = V_c - V_a$$

And with that, we can take any point \$(s, t)\$ and get its uv coordinates under this mapping with a matrix multiplication:

$$\begin{bmatrix}u\\v\end{bmatrix} = \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix} \begin{bmatrix}s\\t\end{bmatrix} + \begin{bmatrix}U_a\\V_a\end{bmatrix}$$

You can verify that for the points A, B, and C, this gives exactly the UVs we expect, and so it gives correct UVs for any linear combination of them too. And since we decided to place our \$\triangle PMN\$ in the same plane, this formula also gives the corresponding UVs for \$P\$, \$M\$, and \$N\$.

Now we can work backwards from their 2D UV coordinates to find their 3D positions on the plane.

By construction, the plane projection mapping algorithm always assigns \$P\$ the coordinates \$(u, v) = (0, 0)\$ (we always subtract \$P\$ from the input point, getting zero if \$P\$ was the input, and the subsequent multiplications keep it at zero)

Now we'll solve...

$$\begin{align} \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix} \begin{bmatrix}s_P\\t_P\end{bmatrix} + \begin{bmatrix}U_a\\V_a\end{bmatrix} &=\begin{bmatrix}0\\0\end{bmatrix}\\ \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix} \begin{bmatrix}s_P\\t_P\end{bmatrix} &=\begin{bmatrix}-U_a\\-V_a\end{bmatrix}\\ \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix}^{-1} \begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix} \begin{bmatrix}s_P\\t_P\end{bmatrix} &=\begin{bmatrix}u_1 & u_2\\v_1 & v_2\end{bmatrix}^{-1}\begin{bmatrix}-U_a\\-V_a\end{bmatrix}\\ \begin{bmatrix}s_P\\t_P\end{bmatrix} &=\frac 1 {u_1v_2 - u_2v_1} \begin{bmatrix}v_2 & -u_2\\-v_1 & u_1\end{bmatrix} \begin{bmatrix}-U_a\\-V_a\end{bmatrix} \end{align}$$

Now we have

$$P = A + s_P(B-A) + t_P(C-A)$$

There's a similar trick for finding \$M\$ and \$N\$.

First, note that we get our \$u\$ coordinate for a point \$X\$ as \$U_X = \frac {\vec U \cdot (X - P)} {\vec U \cdot (M - P)}\$. So if our input point was \$M\$ itself, then this can only give us the number one!

The v coordinate is a bit harder to see, but note that \$\vec V = (M - P) \times \vec {PMPN}\$. So it must be orthogonal to \$(M - P)\$, so the dot product \$\vec V \cdot (M - P)\$ must be zero.

The same argument applies in turn to \$N\$, and so we get a clearer picture of what these three points signify: \$P\$, \$M\$, and \$N\$ define three corners of our UV-mapping square:

$$P \rightarrow (0, 0)\\ M \rightarrow (1, 0)\\ N \rightarrow (0, 1)$$

So, we can substitute those UV values, and solve for \$M\$ and \$N\$ in exactly the same way:

$$\begin{bmatrix}s_M\\t_M\end{bmatrix} =\frac 1 {u_1v_2 - u_2v_1} \begin{bmatrix}v_2 & -u_2\\-v_1 & u_1\end{bmatrix} \begin{bmatrix}1 -U_a\\-V_a\end{bmatrix}\\ M = A + s_M(B - A) + t_M (C-A)$$

$$\begin{bmatrix}s_N\\t_N\end{bmatrix} =\frac 1 {u_1v_2 - u_2v_1} \begin{bmatrix}v_2 & -u_2\\-v_1 & u_1\end{bmatrix} \begin{bmatrix}-U_a\\1 - V_a\end{bmatrix}\\ N = A + s_N(B - A) + t_N(C-A)$$

In my tests with random triangles, using single-precision floats throughout, this correctly reproduces the UV mapping to within 0.001, which I'd attribute to rounding errors along the way. If you have more sensible triangles, and use double-precision intermediates, you could probably reduce the deviation even further.