My approach is to calculate two tangent vectors parallel to axis X and Y respectively. Then calculate the cross product to find the normal vector.
The tangent vector is given by the line that crosses the middle point on the two nearest segments as is shown in the following picture.
I was wondering whether there is a more direct calculation, or less expensive in terms of CPU cycles.
You can actually calculate it without a cross product, by using the "finite difference method" (or at least I think it is called in this way).
Actually it is fast enough that I use it to calculate the normals on the fly in a vertex shader.
// # P.xy store the position for which we want to calculate the normals
// # height() here is a function that return the height at a point in the terrain
// read neightbor heights using an arbitrary small offset
vec3 off = vec3(1.0, 1.0, 0.0);
float hL = height(P.xy - off.xz);
float hR = height(P.xy + off.xz);
float hD = height(P.xy - off.zy);
float hU = height(P.xy + off.zy);
// deduce terrain normal
N.x = hL - hR;
N.y = hD - hU;
N.z = 2.0;
N = normalize(N);