13
$\begingroup$

When applying the classical formula for the angle between two vectors:

$$\alpha = \arccos \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{\|\mathbf{v_1}\| \|\mathbf{v_2}\|}$$

one finds that, for very small/acute angles, there is a loss of precision and the result is not accurate. As explained in this Stack Overflow answer, one solution is to use the arctangent instead:

$$\alpha = \arctan2 \left(\|\mathbf{v_1} \times \mathbf{v_2}\|, \mathbf{v_1} \cdot \mathbf{v_2} \right) $$

And this indeed gives better results. However, I wonder if this would give bad results for angles very close to $\pi / 2$. Is it the case? If so, is there any formula to accurately compute angles without checking for a tolerance inside an if branch?

$\endgroup$
  • 1
    $\begingroup$ This is going to depend on the implementation of the two parameter inverse tangent function. The slow, stable versions switch conditionally between working with x/y and y/x to maintain precision, while the fast ones just stick things in the right quadrant and thus aren't any more precise than the one parameter version. $\endgroup$ – origimbo Aug 23 '17 at 10:51
  • $\begingroup$ You should define "loss of precision": suppose the right answer is $\alpha$ and you get instead $\alpha + \Delta$. Do you need $\Delta \ll \alpha$ or $\Delta \ll \pi$ is sufficient? $\endgroup$ – Stefano M Aug 23 '17 at 20:54
  • $\begingroup$ In this case, the right answer was $\alpha$ and I got $\alpha \cdot 10^8$, both $\ll 1$. $\endgroup$ – astrojuanlu Aug 24 '17 at 7:40
18
$\begingroup$

(I have tested this approach before, and I remember it worked correctly, but I haven't tested it specifically for this question.)

As far as I can tell, both $\|\mathbf{v}_1\times \mathbf{v}_2\|$ and $\mathbf{v}_1\cdot \mathbf{v}_2$ can suffer from catastrophic cancellation if they are almost parallel/perpendicular—atan2 can't give you good accuracy if either input is off.

Start by reformulating the problem as finding the angle of a triangle with side lengths $a=|\mathbf{v}_1|$, $b=|\mathbf{v}_2|$ and $c=|\mathbf{v}_1-\mathbf{v}_2|$ (these are all accurately computed in floating point arithmetic). There is a well-known variant of Heron's formula due to Kahan (Miscalculating Area and Angles of a Needle-like Triangle), which allows you to compute the area and angle (between $a$ and $b$) of a triangle specified by its side lengths, and do so numerically stably. Because the reduction to this subproblem is accurate as well, this approach should work for arbitrary inputs.

Quoting from that paper (see p.3), assuming $a\geq b$, $$ \mu = \begin{cases} c-(a-b),&\text{if }b\geq c\geq 0,\\ b-(a-c),&\text{if }c>b\geq 0,\\ \text{invalid triangle},&\text{otherwise} \end{cases}$$ $$ \mathrm{angle} = 2\arctan\left( \sqrt{\frac{((a-b)+c)\mu}{(a+(b+c))((a-c)+b)}} \right)$$ All the parentheses here are placed carefully, and they matter; if you find yourself taking the square root of a negative number, the input side lengths are not the side lengths of a triangle.

There is an explanation of how this works, including examples of values for which other formulas fail, in Kahan's paper. Your first formula for $\alpha$ is $C''$ on page 4.

The main reason I suggest Kahan's Heron's formula is because it makes a very nice primitive—lots of potentially tricky planar geometry questions can be reduced to finding the area/angle of an arbitrary triangle, so if you can reduce your problem to that, there is a nice stable formula for it, and there's no need to come up with something on your own.

Edit Following Stefano's comment, I made a plot of relative error for $v_1=(1,0)$, $v_2=(\cos\theta, \sin\theta)$ (code). The two lines are the relative errors for $\theta=\epsilon$ and $\theta=\pi/2-\epsilon$, $\epsilon$ going along the horizontal axis. It seems that it works. enter image description here

$\endgroup$
  • $\begingroup$ Thanks for the link and the answer! Unfortunately the second formula I wrote does not appear in the article. On the other hand, this method can get a bit complex, as it requires projection in 2D. $\endgroup$ – astrojuanlu Aug 23 '17 at 17:18
  • 2
    $\begingroup$ @astrojuanlu There is no projection to 2d here: whatever the two 3d vectors are, they define a single (planar) triangle between them—you only need to know its side lengths. $\endgroup$ – Kirill Aug 23 '17 at 17:22
  • $\begingroup$ You are right, my comment makes no sense. I was thinking in coordinates instead of lengths. Thanks again! $\endgroup$ – astrojuanlu Aug 23 '17 at 18:13
  • 2
    $\begingroup$ @astrojuanlu One more thing I want to note: it seems that there is a formal proof that the area formula is accurate in How to Compute the Area of a Triangle: a Formal Revisit, Sylvie Boldo, using Flocq. $\endgroup$ – Kirill Aug 23 '17 at 21:51
  • $\begingroup$ Excellent answer, but I dispute that you can always accurately compute $c$ in floating point arithmetic. In fact if $c < \epsilon \cdot \min(a,b)$ then catastrophic cancellations occurs in computing the components of $(\mathbf{v}_1-\mathbf{v}_2)$. $\endgroup$ – Stefano M Aug 23 '17 at 21:51
7
$\begingroup$

The efficient answer to this question is, not too surprisingly, in another note by Velvel Kahan:

$$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|,\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}-\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|\right)$$

where I use $\arctan(x,y)$ as the angle made by $(x,y)$ with the horizontal axis. (You might have to flip the order of the arguments in some languages.)

(I gave a Mathematica demonstration of Kahan's formula here.)

$\endgroup$
  • $\begingroup$ You mean $\arctan2$? $\endgroup$ – astrojuanlu Sep 1 '17 at 14:19
  • 1
    $\begingroup$ I'm used to just depicting the two-argument arctangent as $\arctan(x,y)$, yes. In a language like FORTRAN, the equivalent would be ATAN2(Y, X). $\endgroup$ – J. M. Sep 1 '17 at 14:42

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.