Compute the vertex normals of a triangle mesh generated from heightmap


I am currently dealing with the problem with vertex normals. I have a triangle mesh generated from a heightmap with uniform grid triangulation as shown in the figure(I am sorry it is drawn by hand, so not quite precise).

I know there are 2 ways to compute the vertex normals:

  1. Use a central difference filter directly on the 2D heightmap texture. This way is pretty easy to implement. (Generating a normal map from a height map?)

  2. Iterate all the triangles(faces). For each face, compute its face normal, and add it into each of the three vertex normals that the face contributes to. After then, do the normalization for all of them (cited from Help)

The first method can be parallelized with CUDA without a race condition. However, when I tried to implement the second method, I found that it is hard to avoid the problem. Because each vertex is normally shared by 6 triangles. So I would like to ask whether there is a solution that can avoid race condition but also preserve a good performance?