TerrainY: Calculating terrain height in a given point


Introduction and Requirements

This page (fully written by me, Mr Brain, between January 10 and 12 of 2008, for the 5th Age) explains progressively how to build a function that returns the height of the surface of the terrain in any given location of a world with terrain (both the location and the height in centimetres). Why is this useful? Because with it, you can make bots that move on top of the terrain instead of going through it. You can also detect when your users are underground, no matter how complex your terrain is. You can automatically place particle emitters that make it rain or snow all the way to the ground level of your complex terrain without actually filling the underground portion of your world with stray particles. And there are many other useful applications.

Intellectual Property is for losers, so this text is free for reproduction, retransmission and usage, commercial or not; I only ask that you always credit me whenever you use this text or the images and code it contains.

You should be able to follow this procedure reasonably well as long as you have basic knowledge of the C language and high-school level trigonometry and geometry. It shouldn't take a rocket scientist to convert the code into a different language if needed. Remember that floorf and sqrtf require libm (linux); Make sure you #include <math.h> and compile with -lm .

WARNING: Values returned by this function are based on the original/flat terrain level without having in account the TERRAIN ELEVATION OFFSET world attribute. If yours is different than zero, you may need to add its value to the return value of this function for *some* practical applications.

The code was debugged, and I think I fixed all the syntax errors I originally made, but please let me know if I skipped any of them when fixing them here. My e-mail address is myshelter )at( gmail.com .


Algorithm and code

Terrain is a matrix of triangles. Each cell contains two triangles next to each other so that, if the terrain is flat, both triangles form a square that covers the entire cell. Let's give the name 'diagonal' to the junction between the two triangles that goes across the middle of one cell. After some observation, I have noticed that the diagonal does not connect the same points for each cell; In some cells, it connects the NE point to the SW point, while in others it connects the NW to the SE. If a cell has one of these diagonal directions, then the four adjacent cells will have the opposite directions, so that each 2x2 cells form an X-shape as represented below. The ground zero cell has its diagonal in the STRAIGHT direction, SE to NW. The blue lines are the diagonals and the red lines are the cartesian X and Z axes of the world.

The first thing we have to do is convert the standard format coordinates into SDK format. This is very easy - coordinates with N or W in front of them stay the same, all we have to do is to remove the N or W. On the other hand, S and E coordinates become negative, and the character also disappears, for instance, 2S becomes -2.

Now, if you look at my drawing, you should see that there is a very easy method to distinguish cells with straight diagonals and cells with inverted diagonals: If you sum the coordinates and the result is ODD, the diagonal is INVERTED, otherwise (if it's EVEN) the diagonal is STRAIGHT. Therefore, we get:

int isdiagonalstraight (int x, int z) {
  x = (x < 0 ? x * -1 : x);
  z = (z < 0 ?z * -1 : z);
  return !((x + z) % 2);
}

That out of the way, we can get started on the actual algorithm. The SDK treats object positions as centimetres, with the origin (position 0,0) being where the red axes cross in the above drawing. Each cell has 10m, so 1000cm. All positions in AW are in cm and no smaller unit is permitted; Positions are rounded to the nearest cm when necessary. So we should develop this algorithm as cm and convert to and from terrain coordinates if necessary.

What lies beneath...

As I see it, there are four distinct cases here: Either our position is in a CORNER (junction between 4 cells), in an EDGE (junction between 2 cells), in a DIAGONAL (junction between the 2 triangles in a cell) or in the middle of a triangle. We can treat these on a case by case basis. For that purpose, we should devise a function that can tell us what case we're dealing with.

Let x and z be the position of the object, in centimetres (different from the x and z in the function above). If the remainders of the integer division (also called modulus) between x and 1000 and z and 1000 are both zero, it means we're right on top of a corner, because the number 1000 fits an exact amount of times in both of our coordinates. If only ONE of them is zero, but the other isn't, then we're in an edge between two cells.

int terrainy (int x, int z) {
  if (x % 1000 == 0 && z % 1000 == 0) {
    /* Corner */
  } else if (x % 1000 == 0) {
    /* Vertical edge */
  } else if (z % 1000 == 0) {
    /* Horizontal edge */
  } else {
    /* Middle of the cell */
  }
}

If none of them is zero, then we are in the middle of a cell. If we're in the middle of a cell, we must find out if we're on the diagonal or on a triangle, and which triangle would also help.

For that purpose, we must find the orientation of the diagonal using the first function. We can also SHIFT our coordinate system so the southeast corner of the cell we're working in becomes position (0,0) (temporarily). We perform both of these steps by obtaining the cell coordinates of the cell we're working in, which is easy: First you must divide your absolute position values by 1000, which will return a decimal value. Then you obtain the floor of that value, which is the greatest integer value no greater than the floating point value. For instance, if you look at the above drawing, you'll see that it's obvious that -600, -700 is in the 1s 1e cell. -600/1000=-0.6 and -700/1000=-0.7. Then the floor of those values is -1 and -1 (because -0.6 and -0.7 are greater than -1), which is the expected value.

We then feed the cell coordinates, which we will call cx and cz, to the isdiagonalstraight function, and also add -1000*cx to x and -1000*cz to z, which will shift the coordinate system as intended. Therefore:

int terrainy (int x, int z) { int cx = 0, cz = 0; if (x % 1000 == 0 && z % 1000 == 0) { /* Corner */ } else if (x % 1000 == 0) { /* Vertical edge */ } else if (z % 1000 == 0) { /* Horizontal edge */ } else { cx = (int) floorf ((float) x / 1000); cz = (int) floorf ((float) z / 1000); x -= 1000*cx; z -= 1000*cz; if (isdiagonalstraight (cx, cz)) { /* Find position */ } else { /* Find position */ } } }

Now that we shifted our coordinate system, we're always working as if we were in the ground zero cell. This means our shifted x and z coordinates only range between 0 and 999. Now, considering the cell has a STRAIGHT diagonal, you can look at the drawing to see that there's a very easy way to determine where our point is within the cell: If x = z, we're right on top of the diagonal. If z > x, we're in the upper triangle. If z < x, we're in the lower triangle. Conversely, if the diagonal is INVERTED, then all we have to do is see which of the above conditions would apply in a cell with FLIPPED coordinates. However, we can't flip BOTH of them, because we'd end up with an inverted diagonal again, so in our case we'll just use 1000-x and z. If z = 1000-x, we're on top of our inverted diagonal. If z > 1000-x, we're in the upper triangle, and if it's < we're in the lower triangle. Here's the revised code:

int terrainy (int x, int z) {
  int cx = 0, cz = 0;
  if (x % 1000 == 0 && z % 1000 == 0) {
    /* Corner */
  } else if (x % 1000 == 0) {
    /* Vertical edge */
  } else if (z % 1000 == 0) {
    /* Horizontal edge */
  } else {
    cx = (int) floorf ((float) x / 1000);
    cz = (int) floorf ((float) z / 1000);
    x -= 1000*cx;
    z -= 1000*cz;
    if (isdiagonalstraight (cx, cz)) {
      if (z == x) {
        /* Diagonal */
      } else if (z > x) {
        /* Upper */
      } else {
        /* Lower */
      }
    } else {
      if (z == 1000 - x) {
        /* Diagonal */
      } else if (z > 1000 - x) {
        /* Upper */
      } else {
        /* Lower */
      }
    }
  }
}

With this code, we know EXACTLY what lies under the point where we want to determine the altitude of the terrain: If it's a surface, an edge or a vertex.

Terrain in corners

Now, we have to make a small break in our train of thought to introduce a new abstract function:

int terrainheight (int x, int z);

This function returns the actual terrain height value for a CELL; This is obtained from the SDK through the use of query functions. Because arrays can't have negative indices, we can't just stick our data into an array and use the array directly in this algorithm; A function that shifts positions based on the amount of pages retrieved must be used, but since that procedure isn't important for this algorithm, I won't include it here. Just assume that terrainheight returns the terrain cell height values :P You will have to make your own terrainheight.

As some of you may know, the height of a terrain cell is actually the height of it's SOUTHEAST corner (I'm tired of telling this to Max). Therefore, we can fill the corners part of our algorithm very easily:

int terrainy (int x, int z) {
  int cx = 0, cz = 0;
  cx = (int) floorf ((float) x / 1000);
  cz = (int) floorf ((float) z / 1000);
  if (x % 1000 == 0 && z % 1000 == 0) {
    return terrainheight (cx, cz);
  } else if (x % 1000 == 0) {
    /* Vertical edge */
  } else if (z % 1000 == 0) {
    /* Horizontal edge */
  } else {
    x -= 1000*cx;
    z -= 1000*cz;
    if (isdiagonalstraight (cx, cz)) {
      if (z == x) {
        /* Diagonal */
      } else if (z > x) {
        /* Upper */
      } else {
        /* Lower */
      }
    } else {
      if (z == 1000 - x) {
        /* Diagonal */
      } else if (z > 1000 - x) {
        /* Upper */
      } else {
        /* Lower */
      }
    }
  }
}

Terrain in borders

Wasn't that was easy? Now let's look at the edges between two cells. These edges are straight lines between two cell corners. So to calculate the height of the terrain in an edge, we must use some good old trigonometry. Consider that the two corners of a cell are always aligned to an axis, X or Z. So they both have the same X or Z coordinate (but not both, of course). With the remaining coordinate and the height of the terrain, we can build a triangle. For example, imagine that Z is fixed (Horizontal edge):

Since Z is always the same, we're only working with X and Y coordinates. We know the X coordinates of the two corners that are linked by the EDGE we are on, and we can obtain their Y coordinates (Cy and Cy+Ay in the picture) using the terrainheight function. Now, there are two things that may happen here: Either the two corners are at the same height (Ay = 0), or they are not. If they are at the same height, and knowing that they are connected by a straight line, it stands to reason that any point on the line is also at that height. If they are not, then one of them must be higher than the other, and we have a situation similar to the one in the above picture. Cy is the height of the lowest of the two corners, and Ay is the difference between both heights.

We all went to school, so we know how to solve this problem, right? Both triangles (the one that ends in the green line, which is the position whose height we want, and the one that ends in the black line) superimpose, so we can use either the SINE or the TANGENT of the angle alpha to obtain the portion of the green line that stands above the dX line (the portion below also has height Cy, so we only have to add that value afterwards). Looking at the smaller triangle, tg alpha = Y/B, and looking at the bigger triangle, tg alpha = Ay/dX. As we know, the distance between two cell corners is always 1000cm, so dX=1000 => tg alpha = Ay/1000. As such, we can obtain the equation Y/B = Ay/1000 <=> Y = Ay*B/1000 . B is eiher x - (cx*1000) or (cx+1)*1000 - x, depending on what corner has the lowest height. So we only have to calculate Ay and then stick this formula into the program (no, no actual tangents are required, whew). The value we want will be Cy+Y. Also, Vertical edges will work just like horizontal edges but we use the Z coordinate (which is the one that varies along the edge) instead of the X one. So we can now fill the next two sections of the function:

int terrainy (int x, int z) {
  double cx = 0, cz = 0, b = 0, ay = 0;
  cx = (int) floorf ((float) x / 1000);
  cz = (int) floorf ((float) z / 1000);
  if (x % 1000 == 0 && z % 1000 == 0) { /* Corner */
    return terrainheight (cx, cz);
  } else if (x % 1000 == 0) { /* Vertical edge */
    if (terrainheight (cx, cz) == terrainheight (cx, cz + 1))
      return terrainheight (cx, cz);
    else if (terrainheight (cx, cz) < terrainheight (cx, cz + 1)) {
      ay = terrainheight (cx, cz + 1) - terrainheight (cx, cz);
      b = z - (cz * 1000);
      return ay * b / 1000 + terrainheight (cx, cz);
    } else {
      ay = terrainheight (cx, cz) - terrainheight (cx, cz + 1);
      b = (cz + 1) * 1000 - z;
      return ay * b / 1000 + terrainheight (cx, cz + 1);
    }
  } else if (z % 1000 == 0) { /* Horizontal edge */
    if (terrainheight (cx, cz) == terrainheight (cx + 1, cz))
      return terrainheight (cx, cz);
    else if (terrainheight (cx, cz) < terrainheight (cx + 1, cz)) {
      ay = terrainheight (cx + 1, cz) - terrainheight (cx, cz);
      b = x - (cx * 1000);
      return ay * b / 1000 + terrainheight (cx, cz);
    } else {
      ay = terrainheight (cx, cz) - terrainheight (cx + 1, cz);
      b = (cx + 1) * 1000 - x;
      return ay * b / 1000 + terrainheight (cx + 1, cz);
    }
  } else {
    x -= 1000*cx;
    z -= 1000*cz;
    if (isdiagonalstraight (cx, cz)) {
      if (z == x) {
        /* Diagonal */
      } else if (z > x) {
        /* Upper */
      } else {
        /* Lower */
      }
    } else {
      if (z == 1000 - x) {
        /* Diagonal */
      } else if (z > 1000 - x) {
        /* Upper */
      } else {
        /* Lower */
      }
    }
  }
}

Terrain in diagonals

Next up are diagonals. Regardless of them being straight or inverted, diagonals are always straight lines between two corners. Therefore, they are also EDGES; The only difference is that instead of cx to cx+1 or cz to cz+1, they connect opposite corners of a cell. Otherwise, they also form a triangle similar to the one in the drawing above, so the method is exactly the same and we can reuse most of our code.

Let's look first at the variables along the Y axis, Cy and Ay. Despite the corners now being in opposite sides of the cell, these variables are still the same and can be retrieved in the same way. But the horizontal variables, dX and B, are now different. dX is no longer always 1000; However, it is the hypothenusa of a triangle it makes with two sides of the cells, whose lengths are always 1000, as we know. So it's still a constant and it can be calculated through the pythagorean theorem, dX = sqrt (2 * 1000^2) = 1414.21 (approximately).

Finally, we have to obtain B. We will use the STRAIGHT diagonal cells as an example and then convert the results to the inverted diagonal cells. B is not hard to obtain as long as we have offset our coordinate system so we're always (artificially) working with the GZ cell. We know that, inside the blocks marked /* Diagonal */ in the above code, x and z (which have been modified above) essentially represent the target position WITHIN the cell. This position is above the diagonal, so we can use the x or z coordinates with the pythagorean theorem to obtain B (one of the coordinates will suffice because we know they are the same). B = sqrt (2 * z^2) . Once we have B, we can apply the same formula as we did for edges, Ay * B / dX + Cy, to obtain the height of the terrain point. In INVERTED diagonal cells, since we flipped around the Z axis, we will use the Z coordinate to obtain B (if we'd flipped around X, we'd use X; we can also use 1000 - x but that would be unneccessary calculations which would make the function heavier).

int terrainy (int x, int z) {
  double cx = 0, cz = 0, b = 0, ay = 0;
  cx = (int) floorf ((float) x / 1000);
  cz = (int) floorf ((float) z / 1000);
  if (x % 1000 == 0 && z % 1000 == 0) { /* Corner */
    return terrainheight (cx, cz);
  } else if (x % 1000 == 0) { /* Vertical edge */
    if (terrainheight (cx, cz) == terrainheight (cx, cz + 1))
      return terrainheight (cx, cz);
    else if (terrainheight (cx, cz) < terrainheight (cx, cz + 1)) {
      ay = terrainheight (cx, cz + 1) - terrainheight (cx, cz);
      b = z - (cz * 1000);
      return ay * b / 1000 + terrainheight (cx, cz);
    } else {
      ay = terrainheight (cx, cz) - terrainheight (cx, cz + 1);
      b = (cz + 1) * 1000 - z;
      return ay * b / 1000 + terrainheight (cx, cz + 1);
    }
  } else if (z % 1000 == 0) { /* Horizontal edge */
    if (terrainheight (cx, cz) == terrainheight (cx + 1, cz))
      return terrainheight (cx, cz);
    else if (terrainheight (cx, cz) < terrainheight (cx + 1, cz)) {
      ay = terrainheight (cx + 1, cz) - terrainheight (cx, cz);
      b = x - (cx * 1000);
      return ay * b / 1000 + terrainheight (cx, cz);
    } else {
      ay = terrainheight (cx, cz) - terrainheight (cx + 1, cz);
      b = (cx + 1) * 1000 - x;
      return ay * b / 1000 + terrainheight (cx + 1, cz);
    }
  } else {
    x -= 1000*cx;
    z -= 1000*cz;
    if (isdiagonalstraight (cx, cz)) {
      if (z == x) { /* Diagonal */
        if (terrainheight (cx, cz) == terrainheight (cx + 1, cz + 1))
          return terrainheight (cx, cz);
        else if (terrainheight (cx, cz) < terrainheight (cx + 1, cz + 1)) {
          ay = terrainheight (cx + 1, cz + 1) - terrainheight (cx, cz);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx, cz);
        } else {
          ay = terrainheight (cx, cz) - terrainheight (cx + 1, cz + 1);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx + 1, cz + 1);
        }
      } else if (z > x) {
        /* Upper */
      } else {
        /* Lower */
      }
    } else {
      if (z == 1000 - x) { /* Diagonal */
        if (terrainheight (cx + 1, cz) == terrainheight (cx, cz + 1))
          return terrainheight (cx + 1, cz);
        else if (terrainheight (cx + 1, cz) < terrainheight (cx, cz + 1)) {
          ay = terrainheight (cx, cz + 1) - terrainheight (cx + 1, cz);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx + 1, cz);
        } else {
          ay = terrainheight (cx + 1, cz) - terrainheight (cx, cz + 1);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx, cz + 1);
        }
      } else if (z > 1000 - x) {
        /* Upper */
      } else {
        /* Lower */
      }
    }
  }
}

Alright, that does it for all the vertices and edges! All that's left now are the triangle surfaces. And this is where things get a little more complicated...

Terrain in the triangle surfaces

Let us think of the upper triangle of a cell with a straight diagonal first. As you can see, there is no simple trigonometry method for finding the height in the middle of the triangle surface. We could use multiple triangles or obtain lines and intersect them, but I think that would be an unnecessary hassle, so what I'll do is consider the vertical line (parallel to the Y axis) that crosses the triangle in the point where we want to calculate the height and then obtain its intersection with the surface. The Y coordinate of this intersection is the height of the terrain in that point.

First of all, we must obtain the equation of the plane defined by the triangle. We will write this equation in the form N dot (P - A) = 0 (dot is the dot product), where N is the normal to the plane, P is a point in the plane - in this case our desired point - and A is another point in the plane, I'm using A but you can also use B or C. The important step of this operation is to obtain a vector that is normal (perpendicular) to the plane. This can be done using the cross product, as the cross product between two vectors results in a perpendicular vector. Using the three corners of the triangle, we can easily obtain two vectors, which we will call U and V: U = B - A and V = C - A. I'm not going to explain the nuances of the cross product here, but considering the cross product U x V = (Wx, Wy, Wz) the formulas for each coordinate are Wx = UyVz-VyUz, Wy = -(UxVz-UzVx), Wz = UxVy-UyVx.

The second step is to obtain the equation of the line. This equation is in the format P = P0 + tV, where P is an arbitrary point in the line, P0 is a known point in the line, V is a vector that defines the line and t is a number that modifies the said vector. We can obtain a vector by subtracting two known points in the line. As a matter of fact, we know an infinite amount of points in the line - any point (x, Y, z), where x and z are the coordinates passed to our function, is in the line. so let's use (x, 0, z) as our P0, and (0, 1, 0) as our vector.

Here are the two equations:
(Wx, Wy, Wz) dot (P - A) = 0
P = (x, 0, z) + t(0, 1, 0)

The intersection between the line and the plane is where both conditions are met, so we're doing to insert the equation of the line in the equation of the plane:

(Wx, Wy, Wz) dot ((x, 0, z) + t(0, 1, 0) - A) = 0

Solving this for t:

t = (Wx, Wy, Wz) dot (A - (x, 0, z)) / (Wx, Wy, Wz) dot (0, 1, 0)

And inserting it back in the equation of the line:

P = (x, 0, z) + ((Wx, Wy, Wz) dot (A - (x, 0, z)) / (Wx, Wy, Wz) dot (0, 1, 0))(0, 1, 0)

At this point, we can replace A with (Ax, Ay, Az):

P = (x, 0, z) + ((Wx, Wy, Wz) dot ((Ax, Ay, Az) - (x, 0, z)) / (Wx, Wy, Wz) dot (0, 1, 0))(0, 1, 0)

Now let's solve this:

P = (x, 0, z) + ((Ax-x)Wx + AyWy - (Az-z)Wz) / Wy * (0, 1, 0)

Looking at the equation we can see that P will be a point (x, Y, z) where Y is the result of the expression to the right of the +. Since we only care about the Y (we already know where x and z are), we can discard that part now:

Y = ((Ax-x)Wx + AyWy + (Az-z)Wz) / Wy

Remember that our coodinate system is offset so our cell is always the GZ cell. As such, when we obtain A from the SE corner, we can replace it with (0, Ay, 0) - this happens in 3 out of our 4 triangles:

Yse = (-xWx + AyWy - zWz) / Wy

For the upper triangle of the cells with the inverted diagonal, we're using the NW corner as our A, so the formula is:

Ynw = ((1000-x)Wx + AyWy + (1000-z)Wz) / Wy

And those are our formulas to get Y. In the lower triangle of cells with a straight diagonal we can use the same method, merely using the other corner (SW) of the cell as our C. For cells with an inverted diagonal, the lower triangle can still use this method if B and C are the SW and NE corners of the cell, but the upper triangle no longer contains our SE point, so A will have to be the NW point (cx+1, cz+1). After making these assertions, we can optimize our calculation of the normal vector.

Considering straight diagonal, upper triangle, A=SE, B=NW and C=NE. Then U = B - A = (Bx-Ax, By-Ay, Bz-Az) = (1000-0, By-Ay, 1000-0) = (1000, By - Ay, 1000). Likewise, V = C - A = (Cx-Ax, Cy-Ay, Cz-Az) = (0-0, Cy-Ay, 1000-0) = (0, Cy - Ay, 1000).
For the lower triangle, A=SE, B=NW and C=SW. Then U = (1000, By - Ay, 1000) (still the same), V = C - A = (1000, Cy - Ay, 0).
For the inverted diagonal, upper triangle, A=NW, B=SW, C=NE, so U = B - A = (0, By - Ay, -1000), V = C - A = (-1000, Cy - Ay. 0).
For the lower triangle, A=SE, B=SW, C=NE mean that U = B - A = (1000, By - Ay, 0), V = C - A = (0, Cy - Ay, 1000).

We can still perform the calculations for the Y coordinate of the cross product in advance. Wy = -(UxVz-UzVx), so:
For the straight diagonal, upper triangle: Wy = -(1000*1000 - 1000*0) = -1000000
Lower triangle: Wy = -(1000*0 - 1000*1000) = 1000000
Inverted diagonal, upper triangle: Wy = -(0*0 - -1000*-1000) = 1000000
Lower triangle: Wy = -(1000*1000 - 0*0) = -1000000

Let's put this into code. As with the edges, I'm also going to add an optimization for immediately returning the height of one corner of the triangle if the three corners are at the same height (this will do wonders for the speed of the function in flat areas):

int terrainy (int x, int z) {
  double cx = 0, cz = 0, b = 0, ay = 0, wx = 0, wy = 0, wz = 0;
  cx = (int) floorf ((float) x / 1000);
  cz = (int) floorf ((float) z / 1000);
  if (x % 1000 == 0 && z % 1000 == 0) { /* Corner */
    return terrainheight (cx, cz);
  } else if (x % 1000 == 0) { /* Vertical edge */
    if (terrainheight (cx, cz) == terrainheight (cx, cz + 1))
      return terrainheight (cx, cz);
    else if (terrainheight (cx, cz) < terrainheight (cx, cz + 1)) {
      ay = terrainheight (cx, cz + 1) - terrainheight (cx, cz);
      b = z - (cz * 1000);
      return ay * b / 1000 + terrainheight (cx, cz);
    } else {
      ay = terrainheight (cx, cz) - terrainheight (cx, cz + 1);
      b = (cz + 1) * 1000 - z;
      return ay * b / 1000 + terrainheight (cx, cz + 1);
    }
  } else if (z % 1000 == 0) { /* Horizontal edge */
    if (terrainheight (cx, cz) == terrainheight (cx + 1, cz))
      return terrainheight (cx, cz);
    else if (terrainheight (cx, cz) < terrainheight (cx + 1, cz)) {
      ay = terrainheight (cx + 1, cz) - terrainheight (cx, cz);
      b = x - (cx * 1000);
      return ay * b / 1000 + terrainheight (cx, cz);
    } else {
      ay = terrainheight (cx, cz) - terrainheight (cx + 1, cz);
      b = (cx + 1) * 1000 - x;
      return ay * b / 1000 + terrainheight (cx + 1, cz);
    }
  } else {
    x -= 1000*cx;
    z -= 1000*cz;
    if (isdiagonalstraight (cx, cz)) {
      if (z == x) { /* Diagonal */
        if (terrainheight (cx, cz) == terrainheight (cx + 1, cz + 1))
          return terrainheight (cx, cz);
        else if (terrainheight (cx, cz) < terrainheight (cx + 1, cz + 1)) {
          ay = terrainheight (cx + 1, cz + 1) - terrainheight (cx, cz);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx, cz);
        } else {
          ay = terrainheight (cx, cz) - terrainheight (cx + 1, cz + 1);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx + 1, cz + 1);
        }
      } else if (z > x) { /* Upper */
        if (terrainheight (cx, cz) == terrainheight (cx, cz + 1) && terrainheight (cx, cz) == terrainheight (cx + 1, cz + 1))
          return terrainheight (cx, cz);
        else {
          wx = (terrainheight (cx + 1, cz + 1) - terrainheight (cx, cz)) * 1000 - (terrainheight (cx, cz + 1) - terrainheight (cx, cz)) * 1000;
          wy = -1000000;
          wz = (terrainheight (cx, cz + 1) - terrainheight (cx, cz)) * 1000;
          return (-x * wx + terrainheight (cx, cz) * wy - z * wz) / wy;
        }
      } else { /* Lower */
        if (terrainheight (cx, cz) == terrainheight (cx + 1, cz) && terrainheight (cx, cz) == terrainheight (cx + 1, cz + 1))
          return terrainheight (cx, cz);
        else {
          wx = (terrainheight (cx + 1, cz) - terrainheight (cx, cz)) * -1000;
          wy = 1000000;
          wz = (terrainheight (cx + 1, cz) - terrainheight (cx, cz)) * 1000 - (terrainheight (cx + 1, cz + 1) - terrainheight (cx, cz)) * 1000;
          return (-x * wx + terrainheight (cx, cz) * wy - z * wz) / wy;
        }
      }
    } else {
      if (z == 1000 - x) { /* Diagonal */
        if (terrainheight (cx + 1, cz) == terrainheight (cx, cz + 1))
          return terrainheight (cx + 1, cz);
        else if (terrainheight (cx + 1, cz) < terrainheight (cx, cz + 1)) {
          ay = terrainheight (cx, cz + 1) - terrainheight (cx + 1, cz);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx + 1, cz);
        } else {
          ay = terrainheight (cx + 1, cz) - terrainheight (cx, cz + 1);
          b = (int) sqrtf (2 * z^2);
          return ay * b / 1414.21 + terrainheight (cx, cz + 1);
        }
      } else if (z > 1000 - x) { /* Upper */
        if (terrainheight (cx + 1, cz + 1) == terrainheight (cx + 1, cz) && terrainheight (cx + 1, cz + 1) == terrainheight (cx, cz + 1))
          return terrainheight (cx + 1, cz + 1);
        else {
          wx = (terrainheight (cx, cz + 1) - terrainheight (cx + 1, cz + 1)) * 1000;
          wy = 1000000;
          wz = (terrainheight (cx + 1, cz) - terrainheight (cx + 1, cz + 1)) * 1000;
          return ((1000-x) * wx + terrainheight (cx + 1, cz + 1) * wy + (1000-z) * wz) / wy;
        }
      } else { /* Lower */
        if (terrainheight (cx, cz) == terrainheight (cx + 1, cz) && terrainheight (cx, cz) == terrainheight (cx, cz + 1))
          return terrainheight (cx, cz);
        else {
          wx = (terrainheight (cx + 1, cz) - terrainheight (cx, cz)) * 1000;
          wy = -1000000;
          wz = (terrainheight (cx, cz + 1) - terrainheight (cx, cz)) * 1000;
          return (-x * wx + terrainheight (cx, cz) * wy - z * wz) / wy;
        }
      }
    }
  }
}

And that's our completed function. I removed the members of the normal vector coordinate formulas that were being multiplied by 0, since it made no sense to keep them in the code; I hope that doesn't confuse anyone.

By now, you may be wondering why we went through all the trouble of treating points on a case by case basis if the triangle formulas would be good enough to solve the problem for any point. The answer is simple: Optimization. The triangle formulas require 10 calculations (12 for the upper inverted triangle), while the diagonal formulas only use 8 calculations, the edge formulas 6 calculations and of course when the object is in a corner or in a flat edge/diagonal/triangle no calculations are required (note: cx and cz offset calculations are not included). This allows the function to run faster for a great number of objects, saving precious CPU cycles. Besides, case-by-case analysis allowed us to tackle the problem progressively, which is also good for learning.

If you'd like to use this function with BERRAIN dump files, try this example file. You'll need a BERRAIN dump file (change the name in main) and the loadbtr function, which can be downloaded here.