Generating a normal map from a height map?

Example GLSL code from my water surface rendering shader:\n\n#version 130\nuniform sampler2D unit_wave\nnoperspective in vec2 tex_coord;\nconst vec2 size = vec2(2.0,0.0);\nconst ivec3 off = ivec3(-1,0,1);\n\n vec4 wave = texture(unit_wave, tex_coord);\n float s11 = wave.x;\n float s01 = textureOffset(unit_wave, tex_coord, off.xy).x;\n float s21 = textureOffset(unit_wave, tex_coord, off.zy).x;\n float s10 = textureOffset(unit_wave, tex_coord, off.yx).x;\n float s12 = textureOffset(unit_wave, tex_coord, off.yz).x;\n vec3 va = normalize(vec3(size.xy,s21-s01));\n vec3 vb = normalize(vec3(size.yx,s12-s10));\n vec4 bump = vec4( cross(va,vb), s11 );\n\n\nThe result is a bump vector: xyz=normal, a=height",

\nMy thinking is that each pixel of the image represents a lattice coordinate in a 256 x 256 grid (hence, why there are 257 x 257 heights). That would mean that the normal at coordinate (i, j) is determined by the heights at (i, j), (i, j + 1), (i + 1, j), and (i + 1, j + 1) (call those A, B, C, and D, respectively).\n\nNo. Each pixel of the image represents a vertex of the grid, so intuitively, from symmetry, its normal is determined by heights of neighboring pixels (i-1,j), (i+1,j), (i,j-1), (i,j+1).\nGiven a function f : ℝ2 → ℝ that describes a surface in ℝ3, a unit normal at (x,y) is given by\nv = (−∂f/∂x, −∂f/∂y, 1) and n = v/|v|.\nIt can be proven that the best approximation to ∂f/∂x by two samples is archived by:\n∂f/∂x(x,y) = (f(x+ε,y) − f(x−ε,y))/(2ε)\nTo get a better approximation you need to use at least four points, thus adding a third point (i.e. (x,y)) doesn't improve the result.\nYour hightmap is a sampling of some function f on a regular grid. Taking ε=1 you get:\n2v = (f(x−1,y) − f(x+1,y), f(x,y−1) − f(x,y+1), 2)\nPutting it into code would look like:\n// sample the height map:\nfloat fx0 = f(x-1,y), fx1 = f(x+1,y);\nfloat fy0 = f(x,y-1), fy1 = f(x,y+1);\n\n// the spacing of the grid in same units as the height map\nfloat eps = ... ;\n\n// plug into the formulae above:\nvec3 n = normalize(vec3((fx0 - fx1)/(2*eps), (fy0 - fy1)/(2*eps), 1));\n",

A common method is using a Sobel filter for a weighted/smooth derivative in each direction.\n\nStart by sampling a 3x3 area of heights around each texel (here, [4] is the pixel we want the normal for).\n\n[6][7][8]\n[3][4][5]\n[0][1][2]\n\n\nThen,\n\n//float s[9] contains above samples\nvec3 n;\nn.x = scale * -(s[2]-s[0]+2*(s[5]-s[3])+s[8]-s[6]);\nn.y = scale * -(s[6]-s[0]+2*(s[7]-s[1])+s[8]-s[2]);\nn.z = 1.0;\nn = normalize(n);\n\n\nWhere scale can be adjusted to match the heightmap real world depth relative to its size.",