2
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ Not sure what you mean by "local tilt angle" ... are we talking about a two-dimensional surface or a one-dimensional curve? $\endgroup$ – Zubin Mukerjee Nov 29 '14 at 19:38
  • $\begingroup$ 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). $\endgroup$ – The Oddler Nov 29 '14 at 19:42
  • $\begingroup$ Okay. Makes sense now, thanks :) $\endgroup$ – Zubin Mukerjee Nov 29 '14 at 19:42
2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – The Oddler Nov 29 '14 at 20:12
  • $\begingroup$ 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) $\endgroup$ – Omnomnomnom Nov 29 '14 at 20:16
  • $\begingroup$ Slight error before, see my edit. $\endgroup$ – Omnomnomnom Nov 29 '14 at 20:18
  • $\begingroup$ Thanks a lot mate! That is indeed a much better way to calculate $\sin(\arccos(x))$. $\endgroup$ – The Oddler Nov 29 '14 at 20:23
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – The Oddler Nov 29 '14 at 20:10
  • $\begingroup$ $\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 ? $\endgroup$ – Yves Daoust Nov 30 '14 at 9:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.