You can also use lower exponents for specular highligts, it will make interpolation artifacts less noticeable.

But maybe there is something to improve in your algorithm after all.

It looks like the 5th vertice, the one within a quad, gets a wrong normal. How do you calculate its position and normal ?

I would do something like this :

- at regular grid points, z = f(x, y)
- within each quad, add a 5th vertex, lets call it S, z = average of all four neighbors.
- each grid vertex has its normal defined from the average of plane normals from each of the 8 adjacents triangles. Don’t forget to renormalize after average !
- each S vertex has its normal defined from average of 4 adjacent vertex. Don’t forget to renormalize after average !

This is will bring the result as near as possible as the per-fragment lighting. Increase mesh subdivision if needed.

If the S vertex is defined also with z = f(x, y), then its normal should also be average of 4 adjacent triangles. With also renormalization after.