# 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!

• Not sure what you mean by "local tilt angle" ... are we talking about a two-dimensional surface or a one-dimensional curve? – Zubin Mukerjee Nov 29 '14 at 19:38
• I have a two-dimentional height field, representing a three-dimentional terrain. By "local tilt angle" I mean the angle between a vector pointing straight up and the normal of the terrain in a point (x,y). – The Oddler Nov 29 '14 at 19:42
• Okay. Makes sense now, thanks :) – Zubin Mukerjee Nov 29 '14 at 19:42

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}}$$

• I love how you call it 'our terrain', like we're really working on it together :D Which, even though we're probably thousands of miles appart, is actually true. Thanks for this very clear explanation! – The Oddler Nov 29 '14 at 20:12
• Whoops, I mean... ahem... I guess you can keep it.... Anyway, in response to your other comment, we have $$\sin(\arccos(x)) = \sqrt{1 - x^2}$$ (and you're welcome) – Omnomnomnom Nov 29 '14 at 20:16
• Slight error before, see my edit. – Omnomnomnom Nov 29 '14 at 20:18
• Thanks a lot mate! That is indeed a much better way to calculate $\sin(\arccos(x))$. – The Oddler Nov 29 '14 at 20:23

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}$.

• Option 1 seems perfect. However, since the code is done literally millions of times per seconds I would like not to use $\sin(\arccos(...))$. Is there a better way to get $\sin\theta$ from this? – The Oddler Nov 29 '14 at 20:10
• $\sin^2\theta=1-\cos^2\theta=(a^2+b^2)/(a^2+b^2+1)$. Are you sure it is the sine that you need ? – Yves Daoust Nov 30 '14 at 9:11