I have a data for a 3D function plot (equally distributed values of f(x, y) in a rectangle). I've built a 3d tetrahedral(a.k.a pyramidal) mesh for these values.

Now I have to calculate normals for this mesh, but I want them to be smooth. I came up with 2 approximate algorithms, but both of them produce similar artifacts like this:

(look at the top right part of the plot)

Can anyone please give me a hint about a correct algorithm?