How to split quads of a height-map terrain?

Recently I needed to import some basic terrain height maps for a prototype project. It quickly became obvious that in order to get "smooth" geometry (whether you want it smooth for rendering, collision detection, or both) it's important to pay attention to how the terrain quads are triangulated. The problem is that even on a smooth heightmap most quads will not be planar, so two different ways to split the quad into two triangles will result in different geometry. Most of the time one of these is much worse than the other (sometimes neither is ideal, but at least you pick the "less bad" one). I'd also like to note that adding a vertex in the middle of the quad to get four triangles instead of two actually tends to make things worse on top of doubling the triangle count.

After searching the web for an algorithm and finding out that some big-name engines (eg UE3 docs about the issue just to name one example) just seem to let you set the orders manually, I ended up rolling my own algorithm. I'm not going to claim this is the best possible algorithm, nor am I going to claim it's novel or anything, but since it's very simple, works almost perfectly as long as the source data isn't awfully noisy (in which case any order tends to be bad), and I haven't found any other algorithms published I thought I'd share it here:

So, first do any non-uniform scaling of the heightmap terrain data such that we have a 2d array of 3d coordinates of the terrain vertices. Now for each vertex, calculate diagonal (one for each diagonal direction, on-axis doesn't work nearly as well because the problems occur in diagonal directions) two-sided finite-difference (eg [-1,0,1] kernel) "tangents" for each point (in 3d), then do a cross-product of the two "tangents" to get a "normal" and normalize the result. These "normals" are not that great for shading (you want to use something more sophisticated to calculate final normals), but take a dot-products between these "normals" at diagonally opposite corners of each quad, compare the results from the two diagonal pairs, and split the quad by adding an edge between the pair that gives the larger dot-product (ie have more "similar" normals) and you have an algorithm that will choose the better split order with "surprisingly large probability", even when neither of them is ideal.

That's not very complicated, nor did it take very long for me to come up with the algorithm. I'm sure more intelligent people could come up with much better algorithms, but if you just need "something" (so you don't need to rely on manual ordering) the above should provide you with a decent starting point.