I only tried this very briefly, and seems to get correct results. However, I'm not sure and I may only have tried it trivial cases so it may not work in some more exotic (but still common) cases. I'll give it a shot anyway, but keep that in mind. Note that you need 4-dimensional homogeneous coordinates for this to work, and I'll assume column vectors.

First, let's focus on the problem in the first image only; what is the coordinate of P in the reference coordinates system spanned by the edges of the quad?

Construct the matrix T = [D-A, B-A, 0, A]. That is, the first column is the vector D-A, the second column is the vector B-A, the third column is all zeros and the fourth column is the vector A. The position of P in T's coordinate frame can now be expressed as P' = T[sup]+[/sup] * P, where T[sup]+[/sup] = (T[sup]T[/sup] * T)[sup]-1[/sup] is the pseudo-inverse of T. Now you have the relative position of P within the quad, so that P=(0,0) is A, P = (1,1) is C, and so on.

To map that relative coordinate to the coordinate system panned by the UV-coordinates, you all you should have to do is multiply P' by the matrix T[sub]u[/sub] = [D[sub]u[/sub]-A[sub]u[/sub], B[sub]u[/sub]-A[sub]u[/sub], 0, 0]. That is P[sub]u[/sub] = T[sub]u[/sub] * P'.

You have to keep one thing in mind here; the quad has to be planar in all its attributes. That is, if you can express the coordinate C as A + x*(D-A) + y*(B-A), a linear combination of the basis vectors and an offset, then you must be able to express C[sub]u[/sub] as A[sub]u[/sub] + x*(D[sub]u[/sub]-A[sub]u[/sub]) + y*(B[sub]u[/sub]-A[sub]u[/sub]) using the same x and y. If that is not the case, the quad is not planar. As shown in your image, the quad is not planar unless there is perspective that "distorts" the top quad (that is, it is tilted with B closer and D further away from the view point). That is why you need 4-dimensional homogeneous for this to work; you need the W-component to control the perspective if the quad isn't planar in the first three dimensions.

You may have noticed that my solution does not use all four coordinates, but only three of them. You can split the quad into triangles and use any combination of three points. Just replace A, B and D in my solution with whichever points you select. My solution is still valid, so you don't have to look for another solution. Unless, of course, my solution is totally wrong to begin with, but that's another story...

The C-coordinate is effectively part of the solution in that C and C[sub]u[/sub] are constrained to form a planar quad. If you split the quad into triangles which is always planar, then the constraints on the virtual C and C[sub]u[/sub] are satisfied and are no longer needed for the actual numerical solution.

A different idea: Compute the cross-ratios (AB,AD;AC,AP) and (DA,DC;DB,DP), then find the lines AuPu and DuPu using the fact that those cross-ratios are preserved by projective transformations, then intersect AuPu and DuPu.

NOTE: The cross-ratio of four lines that go through a point can be computed by intersecting those lines with an arbitrary line that doesn't pass through the common point and computing the cross-ratio of the resulting points on that line. I don't know if the Wikipedia article is clear enough about that.

Yet another idea: Compute the projective matrix M that maps ([0 0 1], [1 0 1], [0 1 1], [1 1 1]) to (A, B, C, D) and the matrix Mu that maps those same points to (Au, Bu, Cu, Du). Then compute Mu*inverse(M). This matrix is a 3x3 matrix that corresponds to a projective transformation that maps (A, B, C, D) to (Au, Bu, Cu, Du), and therefore it will map P to Pu (in homogenous coordinates, so plug in a "1" as a third coordinate for P and remember to divide the result by the third coordinate of Pu).

The only tricky part here is how to compute the matrix that maps ([0 0 1], [1 0 1], [0 1 1], [1 1 1]) to any given 4 points. But fear not, I did it for you:
Matrix3D compute_projective_matrix(Point2D A, Point2D B, Point2D C,Point2D D) {
double aux_a = B.x-D.x;
double aux_b = C.x-D.x;
double aux_c = D.x+A.x-B.x-C.x;
double aux_d = B.y-D.y;
double aux_e = C.y-D.y;
double aux_f = D.y+A.y-B.y-C.y;
double h = (aux_a*aux_f-aux_d*aux_c)/(aux_a*aux_e-aux_d*aux_b);
double g = (aux_c-aux_b*h)/aux_a;

Matrix3D result;
result.m11 = B.x*(g+1.0)-A.x;
result.m12 = C.x*(h+1.0)-A.x;
result.m13 = A.x;
result.m21 = B.y*(g+1.0)-A.y;
result.m22 = C.y*(h+1.0)-A.y;
result.m23 = A.y;
result.m31 = g;
result.m32 = h;
result.m33 = 1.0;
return result;
}


EDIT: Now that I look at my code again, it looks like it would divide by 0 if B.x == D.x, which is something I didn't expect. I might have made a mistake somewhere... I think the idea for a solution is sound. I'll try to fix the code when I find time.

Thank you Alvaro. So all I'd need to do to get Pu is to multiply P by the Matrix returned by the above function?

Not quite that simple. I explained it above the code, but I'll try to be more explicit: You need to call this function twice:
// Pseudocode follows
Matrix3D M = compute_projective_matrix(A,B,C,D);
Matrix3D Mu = compute_projective_matrix(Au,Bu,Cu,Du);
Matrix3D T = Mu*inverse(M);

You then multiply T by a column vector (P.x, P.y, 1.0) and if the result is (R.x, R.y, R.z), then Pu is (R.x/R.z, R.y/R.z).