How do you find the Y position of a point between four vertices? HLSL

Let's say there is a grid terrain for a game composed of tiles made of two triangles - made from four vertices. How would we find the Y (up) position of a point between the four vertices?

Terrain

I have tried this:

float diffZ1    = lerp(heights[0], heights[2], zOffset);
float diffZ2    = lerp(heights[1], heights[3], zOffset);
float yPosition = lerp(diffZ1, diffZ2, xOffset);

Where z/yOffset is the z/y offset from the first vertex of the tile in percent / 100. This works for flat surfaces but not so well on bumpy terrain.

I expect this has something to do with the terrain being made from triangles where the above may work on flat planes. I'm not sure, but does anybody know what's going wrong?

This may better explain what's going on here:

enter image description here

In the code above "heights[]" is an array of the Y coordinate of surrounding vertices v0-3. Triangle 1 is made of vertex 0, 2 and 1. Triangle 2 is made of vertex 1, 2 and 3.

I wish to find coordinate Y of p1 when its x,y coordinates lay between v0-3.

So I have tried determining which triangle the point is between through this function:

bool PointInTriangle(float3 pt, float3 pa, float3 pb, float3 pc)
{
  // Compute vectors        
  float2 v0 = pc.xz - pa.xz;
  float2 v1 = pb.xz - pa.xz;
  float2 v2 = pt.xz - pa.xz;

  // Compute dot products
  float dot00 = dot(v0, v0);
  float dot01 = dot(v0, v1);
  float dot02 = dot(v0, v2);
  float dot11 = dot(v1, v1);
  float dot12 = dot(v1, v2);

  // Compute barycentric coordinates
  float invDenom = 1.0f / (dot00 * dot11 - dot01 * dot01);
  float u = (dot11 * dot02 - dot01 * dot12) * invDenom;
  float v = (dot00 * dot12 - dot01 * dot02) * invDenom;

  // Check if point is in triangle
  return (u >= 0.0f) && (v >= 0.0f) && (u + v <= 1.0f);
}

This isn't giving me the results I expected

I am then trying to find the y coordinate of point p1 inside each triangle:

// Position of point p1
float3 pos = input[0].PosI;

// Calculate point and normal for triangles
float3 p1 = tile[0];
float3 n1 = (tile[2] - p1) * (tile[1] - p1); // <-- Error, cross needed
       // = cross(tile[2] - p1, tile[1] - p1);
float3 p2 = tile[3];
float3 n2 = (tile[2] - p2) * (tile[1] - p2);  // <-- Error
       // = cross(tile[2] - p2, tile[1] - p2); 
float newY = 0.0f;

// Determine triangle & get y coordinate inside correct triangle
if(PointInTriangle(pos, tile[0], tile[1], tile[2]))
{
    newY = p1.y - ((pos.x - p1.x) * n1.x + (pos.z - p1.z) * n1.z) / n1.y;
}
else if(PointInTriangle(input[0].PosI, tile[3], tile[2], tile[1]))
{
    newY = p2.y - ((pos.x - p2.x) * n2.x + (pos.z - p2.z) * n2.z) / n2.y;
}

Using the following to find the correct triangle:

if((1.0f - xOffset) <= zOffset)
    inTri1 = true;

And correcting the code above to use the correct cross function seems to have solved the problem.

enter image description here

Solution

Because your 4 vertices may not be on a plane, you should consider each triangle separately. First find the triangle that the point resides in, and then use the following StackOverflow discussion to solve for the Z value (note the different naming of the axes). I personally like DanielKO's answer much better, but the accepted answer should work too:

Linear interpolation of three 3D points in 3D space

EDIT: For the 2nd part of your problem (finding the triangle that the point is in): Because the projection of your tiles onto the xz plane (as you define your coordinates) are perfect squares, finding the triangle that the point resides in is a very simple operation. Here I'll use the terms left-right to refer to the x axis (from lower to higher values of x) and bottom-top to refer to the z axis (from lower to higher values of z).

Each tile can only be split in one of two ways. Either (A) via a diagonal line from the bottom-left corner to the top-right corner, or (B) via a diagonal line from the bottom-right corner to the top-left corner.

For any tile that's split as A: Check if x' > z', where x' is the distance from the left edge of the tile to the point, and z' is the distance from the bottom edge of the tile to the point. If x' > z' then your point is in the bottom-right triangle; otherwise it's in the upper-left triangle.
For any tile that's split as B: Check if x" > z', where x" is the distance from the right edge of your tile to the point, and z' is the distance from the bottom edge of the tile to the point. If x" > z' then your point is in the lower-left triangle; otherwise it's in the upper-right triangle.

(Minor note: Above I assume your tiles aren't rotated in the xz plane; i.e. that they are aligned with the axes. If that's not correct, simply rotate them to align them with the axes before doing the above checks.)