Height at 2D coordinate on a 3D rectangular surface The Problem:
How can I obtain every 3D coordinate on a rectangular surface given x and z? For those who are visual, picture looking down on the surface, and finding the height at where the x and z grid lines intersect under the surface. Please note that the surface is not infinite.
Background Information:
The coordinate system I'm using has the z axis pointing forward, the x axis pointing to the right, and the y axis pointing up. I have a rectangular surface rotated and offset in the 3D space. I know the coordinates for the four corners of the rectangle. This is for a personal programming project. The algorithm will be executed once on tens of thousands of surfaces. The surfaces are bounded by the four points.
Research:
I know that you can find the scalar equation for the plane the rectangle rests on and use that to find the y at a (x, ?, z) coordinate, however that won't tell me if the coordinate at (x, y, z) lies on the rectangle or not, only that it lies on the same plane as the rectangle.
Notation:
Lets call the four corners of the rectangular surface P, and the four corners thereof P1, P2, P3, and P4.
Note: Please use layman terms. I'm 17 in college algebra.
 A: Note that the normal to your plane (rectangle) is given by:
$${\bf n} = ( {\bf p_2} - {\bf p_1} ) \times ( {\bf p_3} - {\bf p_1} )$$
Where $\times$ is the cross product. 
Note how only three of the four corners are needed to describe the plane. Now use the fact that any point on your plane $\bf p$ is described by:
$${\bf n}\cdot({\bf p_1} - {\bf p})=0$$
Or, in other terms:
$$y = \frac{n_x(x_1 - x)+n_z(z_1-z)+n_y(y_1)}{n_y}$$
A: OK. I'm going to say that $P_1$ through $P_4$ have coordinates $(x_1, z_1)$ through $(x_4, z_4)$, and the associated heights are $y_1$ through $y_4$. And you've got a point $(x, z)$ at which you'd like the know the $y$-coordinate. Here are a few steps:
Compute $v = P_2 - P_1$, termwise;
Compute $w = P_3 - P_1$, termwise;
Compute $n = v \times w$, a cross-product. [You may have to look up the definition]. 
$n$ is a list of three coordinates, let's say $n = (A, B, C)$. 
The equation of your plane is then 
$$
Ax + By + Cz + D = 0.
$$
But we need to know what $D$ is. To find it, plug in the coordinates of the point $P_1$ to get
$$
D = -(Ax_1 + By_1 + Cz_1). 
$$
OK, so now you have $A, B, C,$ and $D$. 
Now you've got some point $(x, z)$ in the $xz$-plane, and want to find its y-coordinate. You compute this:
$$
y = -\frac{Ax + Cz + D}{B}.
$$
And that's it!
Let me make this concrete. Suppose that $P_1 = (1, 3, 4)$ (i.e., you've got the point $(1, 4)$, and the surface is 3 units above it), $P_2 = (2, 5, 2)$, and $P_3 = (1, 0, 7)$. Then I get
$$
v = (2-1, 5-3, 2-4) = (1, 2, -2)
$$
and
$$
w = (1-1, 0-3, 7-4) = (0, -3, 3).
$$
Now $n$ is the cross product of these, which is [0, -3, -3].
So $A = 0, B = -3, C = -3$. To find $D$, we plug in $P_1 = (1, 3, 4)$ into $Ax + By + Cz + D = 0$, i.e.
$$
0 \cdot 1 + (-3) \cdot 3 + (-3) \cdot 4 + D = 0
$$
which we can rewrite as
$$
-21 + D = 0,
$$
so $D = 21$. 
Given a point $(x, z)$, we compute
$$
y = -\frac{Ax + Cz + D}{B} = -\frac{0x -3z + 21}{-3}
$$
and the point $(x, y, z)$ will be on the surface. 
Of course, this requires that $(x, 0, z)$ lie within the projection of the rectangle onto the $xz$ plane. Let me describe how to do that. 
I realize that above, I shifted a bit, and said that each of $P_1$ through $P_4$ had three coordinates, but I want to go back to the two-coordinates version of things. So for the remainder of this explanation, $P_1 = (x_1, z_1)$, etc., and we've got a point $(x, y)$ that might or might not be in the rectangle bounded by $P_1$ through $P_4$, and we want a yes or no answer to that question. 
For this, I want to work with vectors, which, for this explanation, are just lists of 2 numbers, except that I'll use brackets instead of parentheses to distinguish them. I need, also, the notion of dot-product. That's defined like this: 
$$
\begin{bmatrix}a & b\end{bmatrix} \cdot \begin{bmatrix}c & d\end{bmatrix} = ac + bd.
$$
If I have two points, say $Q_1 = (3, 6)$ and $Q_2 = (4, 1)$, then their difference (term by term) gets written as a vector:
$$
Q_2 - Q_1 = (4, 1) - (3, 6) = \begin{bmatrix}1 & -5\end{bmatrix}.
$$
Intuition: the vector $\begin{bmatrix}1 & -5\end{bmatrix}$, meaning "one unit to the right, 5 units down" is how you have to move from $Q_1$ to get to $Q_2$. 
I want to define one more thing -- an operation on vectors: 
$$
R(\begin{bmatrix}a & b\end{bmatrix}) = \begin{bmatrix}-b & a\end{bmatrix}.
$$
 You can think of $R$ as "rotate this thing 90 degrees counterclockwise." So if $\begin{bmatrix}a & b\end{bmatrix}$ is a displacement in the northeast direction, then $R(\begin{bmatrix}a & b\end{bmatrix})$ is a displacement in the northwest direction. 
From the points $P_1$ and $P_2$, you can form the equation of the line in the plane that contains them. Indeed, you can probably do this with the algebra that you know. But I'm going to do it a little differently, and say that the equation is this:
$$
( (x, z) - P_1 ) \cdot  R(P_2 - P_1) = 0. 
$$
A point $(x, y)$ is on the line only if it makes the equation above be true. In fact, we can write this:
$$
h = ( (x, z) - P_1 ) \cdot  R(P_2 - P_1).
$$
Then if $h$ is zero, $(x, y)$ is on the line; if $h < 0$, then $(x, y)$ is to the left of the line (as you go from $P_1$ toward $P_2$, and if $h > 0$, then $(x, y)$ is on the other side of the line. 
So here's the deal: to test whether $P = (x, z)$ is in the rectangle, we're going to test whether it's on opposite sides of the edges from $P_1$ to $P_4$ and $P_2$ to $P_3$, and also on opposite sides of the edges from $P_1$ to $P_2$ and $P_4$ to $P_3$. The order's very important there.
Let me write some pseudocode:
input: four points (x1, z1), (x2, z2), (x3, z3), (x4, z4), that form a non-self-intersecting quadrangle, and a test point (x, z)
Output: true if (x, z) is in the quadrangle defined by the other four points. 
q1 = side(x1, z1, x4, z4, x, z) * side(x2, z2, x3, z3, x, z);
q2 = side(x1, z1, x2, z2, x, z) * side(x4, z4, x3, z3, x, z);
if (q1 < 0) and (q2 < 0) return true
else return false
function side(x1, y1, x2, y2, x, y)
a = x2 - x1;  // [a b] is vector between the two points
b = y2 - y1;
u = -b;       // [u v] is rotated vector
v = a;
h = (x - x1) * u + (y - y1) * v;  // dot product in disguise
if (h < 0) return -1;
if (h = 0) return 0;
if (h > 0) return 1; 
And that pretty much does the job. 
