Local tilt angle based on height-field I'm trying to implement a mathematical formula in a program I'm making, but while the programming is no problem I'm having trouble with some of the math. I need to calculate $\sin(\alpha(x,y))$ with $\alpha(x,y)$ the local tilt angle in $(x,y)$.
I have $2$-dimentional square grid, with at each point the height, representing a $3$-dimentional terrain. To find the tilt in a point I can use the height of its direct neighbors. So $h(x+1,y)$ can be used, however $h(x+2,y)$ not. I also know the distance between two neighboring points ($dx$). By tilt I mean the angle between the normal at a point on the terrain and a vector pointing straight up.
This seems like a not too hard problem, but I can't seem to figure out how to do it. Anyone got a good way to do this?
Thanks!
 A: A helpful construct here would be the normal vector to our terrain.
Our terrain is modeled by the equation
$$
z = h(x,y)
$$
Or equivalently, 
$$
z - h(x,y) = 0
$$
We can define $g(x,y,z) = z - h(x,y)$.  It turns out that the vector normal to this level set is given by
$$
\operatorname{grad}(g) =
\newcommand{\pwrt}[2]{\frac{\partial #1}{\partial #2}}
\left\langle \pwrt{g}{x},\pwrt gy, \pwrt gz \right \rangle = 
\left\langle -\pwrt{h}{x},-\pwrt hy, 1 \right \rangle := v(x,y)
$$
We can calculate the angle between this normal and the vertical $\hat u = \langle 0,0,1 \rangle$ using the formula
$$
\cos \theta = \frac{u \cdot v}{\|u\| \|v\|}
$$
in particular, we find that
$$
\cos \theta = \frac{\hat u \cdot v}{\|\hat u\| \|v\|} = 
\frac{1}{\sqrt{1 + \left( \pwrt hx \right)^2 + \left( \pwrt hy \right)^2}}
$$
We may approximate
$$
\pwrt hx(x,y) \approx \frac{h(x+dx,y) - h(x-dx,y)}{2(dx)}\\
\pwrt hy(x,y) \approx \frac{h(x,y+dy) - h(x,y-dy)}{2(dy)}
$$

Note: since you have to calculate $\sin \theta$, you find
$$
\sin \theta = \sqrt{1 - \cos^2 \theta} = 
\frac{\sqrt{\left( \pwrt hx \right)^2 + \left( \pwrt hy \right)^2}}{\sqrt{1 + \left( \pwrt hx \right)^2 + \left( \pwrt hy \right)^2}}
$$
A: Option 1: estimate the partial derivatives using the finite difference scheme $\frac{h(x+1,y)-h(x-1,y)}{2\Delta x}$, $\frac{h(x,y+1)-h(x,y-1)}{2\Delta x}$ and use the normalvector $(h_x',h_y',-1)\to\cos\theta=1/\sqrt{h_x^2+h_y^2+1}$.
Option 2: fit a least squares plane $z=ax+by+c$ in the 3x3 neighborhood and use $(a,b,-1)\to\cos\theta=1/\sqrt{a^2+b^2+1}$.
