Direct way of computing the clockwise angle between two vectors

Question

I want to find out the clockwise angle between two vectors (two-dimensional or three-dimensional).

The classic way with the dot product gives me the inner angle (0-180 degrees), and I need to use some if statements to determine if the result is the angle I need or its complement.

Is there a direct way of computing the clockwise angle?

@MartinR "clockwise" is a generic term to say I want the angle in a specific "direction", not in the nearest "direction". Nickolay O. specified in his answer a way of describind this "direction" — Mircea Ispas, Dec 28 '12 at 09:17
if statement is one of the language primitives. You don't have much left if you can't use them. — SomeWittyUsername, Dec 28 '12 at 09:21
@Felics read my answer. Clockwise direction is not well-defined in 3d. It is planar term — kassak, Dec 28 '12 at 09:22
@icepack The problem is not the "if", it's the additional computatio to be able to use the "if" - like one possible unnecessary cross product — Mircea Ispas, Dec 28 '12 at 09:23
@Felics: "clockwise" is well-defined in 2D, but not in 3D. Checking the z-coordinate of the cross product (as in Nickolay O.'s answer) would mean in 3D: "clockwise for an observer looking from above on the x/y plane." — Martin R, Dec 28 '12 at 09:34
@Felics Also, I should note that you could not define 3D clockwise angle continuously because of hairy ball theorem http://en.wikipedia.org/wiki/Hairy_ball_theorem You would always have pair of vectors, epsilon-movement of one of which would lead to instant switch of clock-wisiness and as a result angle sign — kassak, Dec 28 '12 at 09:58

score 276 · Accepted Answer · edited Apr 18 '23 at 11:18

276

2D case

Just like the dot product is proportional to the cosine of the angle, the determinant is proportional to its sine. So you can compute the angle like this:

dot = x1*x2 + y1*y2      # Dot product between [x1, y1] and [x2, y2]
det = x1*y2 - y1*x2      # Determinant
angle = atan2(det, dot)  # atan2(y, x) or atan2(sin, cos)

The orientation of this angle matches that of the coordinate system. In a left-handed coordinate system, i.e. x pointing right and y down as is common for computer graphics, this will mean you get a positive sign for clockwise angles. If the orientation of the coordinate system is mathematical with y up, you get counterclockwise angles as is the convention in mathematics. Changing the order of the inputs will change the sign, so if you are unhappy with the signs just swap the inputs.

3D case

In 3D, two arbitrarily placed vectors define their own axis of rotation, perpendicular to both. That axis of rotation does not come with a fixed orientation, which means that you cannot uniquely fix the direction of the angle of rotation either. One common convention is to let angles be always positive, and to orient the axis in such a way that it fits a positive angle. In this case, the dot product of the normalized vectors is enough to compute angles.

dot = x1*x2 + y1*y2 + z1*z2    # Between [x1, y1, z1] and [x2, y2, z2]
lenSq1 = x1*x1 + y1*y1 + z1*z1
lenSq2 = x2*x2 + y2*y2 + z2*z2
angle = acos(dot/sqrt(lenSq1 * lenSq2))

Note that some comments and alternate answers advise against the use of acos for numeric reasons, in particular if the angles to be measured are small.

Plane embedded in 3D

One special case is the case where your vectors are not placed arbitrarily, but lie within a plane with a known normal vector n. Then the axis of rotation will be in direction n as well, and the orientation of n will fix an orientation for that axis. In this case, you can adapt the 2D computation above, including n into the determinant to make its size 3×3.

dot = x1*x2 + y1*y2 + z1*z2
det = x1*y2*zn + x2*yn*z1 + xn*y1*z2 - z1*y2*xn - z2*yn*x1 - zn*y1*x2
angle = atan2(det, dot)

One condition for this to work is that the normal vector n has unit length. If not, you'll have to normalize it.

As triple product

This determinant could also be expressed as the triple product, as @Excrubulent pointed out in a suggested edit.

det = n · (v1 × v2)

This might be easier to implement in some APIs, and gives a different perspective on what's going on here: The cross product is proportional to the sine of the angle, and will lie perpendicular to the plane, hence be a multiple of n. The dot product will therefore basically measure the length of that vector, but with the correct sign attached to it.

edited Apr 18 '23 at 11:18

Peter Mortensen

30,738
21
105
131

answered May 14 '13 at 13:19

MvG

57,380
22
148
276

4

Have an upvote - I can't be bothered figuring out if the other answers are correct or not, yours is the clearest and most readable, so it's the one that helped me. – Excrubulent Jul 18 '13 at 15:22
2

For the 2D I'm getting (0,180) and (-180,0). One can check when the result is negative and add 360 in order to get a nice clockwise angle (for example if it's -180 adding 360 results in 180, for -90 adding 360 results in 270 etc.). Don't know if it's just my calculation or the implementation of the `qAtan2(y, x)` (from the Qt framework) but if someone has the same problem as me this might help. – rbaleksandar Dec 13 '16 at 18:35
17

@rbaleksandar: `atan2` usually is in the range [-180°,180°]. To get [0°,360°] without a case distinction, one can replace `atan2(y,x)` with `atan2(-y,-x) + 180°`. – MvG Dec 13 '16 at 18:54
MvG's 2D answer gives the counter closewise angle instead of a closewise one, doesn't it? – sofname Dec 15 '17 at 03:32
1

@jianz: The angle is the positive angle with respect to the coordinate system. If *x* is right and *y* is up, then the angle is counter-clockwise. If *y* is down, it's clockwise. Most computer graphics environments use the latter. If you want to reverse the orientation, simply change the order of the inputs, which will flip the sign of the determinant. – MvG Dec 15 '17 at 10:35
@Mvg: Thanks for the clarification. In my case, I was working with a geographic coordinate reference system that x is right and y is up. Thanks again. – sofname Dec 18 '17 at 02:36
@HariKaramSingh: I'm sure I didn't write down a derivation, or need one to write this answer. What would you like to see derived? The first sentence in the 2d section is knowlegde I've been using so often I don't remember where I first heard it, but Wikipedia knows about it, too \[[1](https://en.wikipedia.org/wiki/Dot_product#Geometric_definition), [2](https://en.wikipedia.org/wiki/Determinant#2_%C3%97_2_matrices)\]. From that to the `atan2` invocation is essentially a single step. The other sections are essentially variants of this theme. – MvG Dec 21 '17 at 21:04
Is it just me, or is the triple product method by fast the most performant one? – Tara Apr 24 '18 at 05:06
1

@Tara: That will depend a lot on the environment. Mainly CPU vs. GPU and compiler optimizations. Naively the triple product would be 9 multiplications and 5 additions/subtractions, while my component-wise formula is 12 multiplications and 5 additions/subtractions. I guess there are compilers which would apply the [distributive law](https://en.wikipedia.org/wiki/Distributive_property) as part of the optimization, reducing the number of multiplications to match what the triple product does out of the box. I wouldn't be surprised if GPUs had highly optimized implementations for vector operations. – MvG Apr 30 '18 at 07:48
@MvG: I think I misread your initial answer... I thought you got rid of the call to `atan2()`. However the code for the triple product just replaces the determinant, leaving the rest of the code the same. That wasn't immediately obvious to me. But thank you anyway for the details. – Tara May 07 '18 at 05:20
3

Noooooo don't ever take acos of a dot product! That's mathematically correct but horribly inaccurate in practice. You could replace your 3d method with another atan2(det,dot); in this case det would be the length of the cross product. – Don Hatch Dec 13 '18 at 09:17
It appears that your answer for 2D is incorrect, please you elaborate how this worked? – Hamish Gibson Oct 28 '20 at 00:14
I don't understand this part please help. in 3d case you said "One special case is the case where your vectors are not placed arbitrarily, but lie within a plane with a known normal vector n". if my vectors lie with in a plane then normal would be v1xv2, right? but then triple product become zero because of `n.(v1xv2)` – M.kazem Akhgary May 17 '21 at 06:08
@M.kazemAkhgary: If you have n := v₁ ×v ₂ then n∙n is not zero but the squared norm of n. But in most circumstances that n will not have the unit length I required in my post, so you'd have to normalize first. The whole point of case 3 however is handle the situation where you already know the normal vector so that it can define the orientation of the plane independent of the other two vectors. In that situation, using the cross product is of little help. The comment of Don Hatch above suggests the cross product would still be useful for numeric stability, but my answer doesn't cover that. – MvG May 17 '21 at 19:32
@DonHatch could you elaborate or link to a separate discussion? – N4ppeL Sep 01 '21 at 12:52
I compared and benchmarked this using numpy: the max difference comparing 100k random vectors was 2.5097014310498933e-13, which seems small enough to be accurate. However, the atan2 version is a little faster since there is no need to check for handling div by zero – N4ppeL Sep 01 '21 at 14:32
2

@N4ppeL For more about poor behavior of of acos of dot product, try this (question asked on 2 different sites, with different answers and references): https://math.stackexchange.com/questions/1143354/numerically-stable-method-for-angle-between-3d-vectors/1782769 https://scicomp.stackexchange.com/questions/27689/numerically-stable-way-of-computing-angles-between-vectors – Don Hatch Sep 11 '21 at 11:36
Just wondering if there's any python package that covers all these and other related use cases, as I write these things over and over again every time I hit them, usually screwing them up three times before getting them right for the particular use case following paper and pencil work – matanster Jun 01 '22 at 13:50

score 6 · Answer 2 · edited Apr 18 '23 at 11:21

This answer is the same as MvG's, but explains it differently (it's the result of my efforts in trying to understand why MvG's solution works).

The anti-clockwise angle theta from x to y, with respect to the viewpoint of their given normal n (||n|| = 1), is given by

atan2( dot(n, cross(x,y)), dot(x,y) )

(1) = atan2( ||x|| ||y|| sin(theta), ||x|| ||y|| cos(theta) )

(2) = atan2( sin(theta), cos(theta) )

(3) = anti-clockwise angle between x axis and the vector (cos(theta), sin(theta))

(4) = theta

where ||x|| denotes the magnitude of x.

Step (1) follows by noting that

cross(x,y) = ||x|| ||y|| sin(theta) n,

and so

dot(n, cross(x,y))

= dot(n, ||x|| ||y|| sin(theta) n)

= ||x|| ||y|| sin(theta) dot(n, n)

which equals

||x|| ||y|| sin(theta)

if ||n|| = 1.

Step (2) follows from the definition of atan2, noting that atan2(cy, cx) = atan2(y,x), where c is a scalar. Step (3) follows from the definition of atan2. Step (4) follows from the geometric definitions of cos and sin.

score 5 · Answer 3 · edited Apr 18 '23 at 10:58

5

To compute the angle you just need to call atan2(v1.s_cross(v2), v1.dot(v2)) for the two-dimensional case. Where s_cross is the scalar analogue of cross production (signed area of parallelogram).

For the two-dimensional case that would be wedge production.

For the three-dimensional case you need to define the clockwise rotation, because from one side of plane clockwise is one direction, from other side of plane is another direction =)

This is the counterclockwise angle, and the clockwise angle is just opposite.

edited Apr 18 '23 at 10:58

Peter Mortensen

30,738
21
105
131

answered Dec 28 '12 at 09:14

kassak

3,974
1
25
36

v1.cross(v2) is a vector, not a scalar and can't be used like this. Nickolay O. describes in his answer how to find out 'direction' of the angle. One way to get 2D angle is: angle = atan2f(v2.x, v2.y) - atan2f(v1.x, v1.y) – Mircea Ispas Dec 28 '12 at 09:27
2

@Felics In 2D cross production often means wedge production http://en.wikipedia.org/wiki/Wedge_product That is signed area of parallelogram. For 2D case that formula is absolutely correct as it dot = |v1||v2|*cos and cross = |v1||v2|sin. That is why atan2 gives correct angle in whole circle range. And as I said for 3d case you need to make some assumptions to have some extension of clockwise orientation – kassak Dec 28 '12 at 09:35
1

@Felics: Note that `atan2f` has the y-coordinate as first argument, so it should be `angle = atan2f(v2.y, v2.x) - atan2f(v1.y, v1.x)`. – Martin R Dec 28 '12 at 09:38
1

@kassak: You could replace `cross` and `dot` by the explicit formula in the 2D case, that would remove all doubts about `cross` returning a 3D vector (but that is only a suggestion, which you can ignore). - Otherwise I like this solution, because it requires only one `atan2f` function call. – Martin R Dec 28 '12 at 09:45
@Martin R thanks for good advice. I've made some corrections to make meaning of formula clearer – kassak Dec 28 '12 at 09:50
your formulas doesn't work for example for vectors (-1, 0) and (0, 1) – user2083364 Oct 17 '13 at 08:47

score 4 · Answer 4 · edited Apr 18 '23 at 11:34

4

Since one of the simplest and most elegant solutions is hidden in one of the comments, I think it might be useful to post it as a separate answer.

acos can cause inaccuracies for very small angles, so atan2 is usually preferred. For the three-dimensional case:

dot = x1*x2 + y1*y2 + z1*z2
cross_x = (y1*z2 – z1*y2)
cross_y = (z1*x2 – x1*z2) 
cross_z = (x1*y2 – y1*x2)
det = sqrt(cross_x*cross_x + cross_y*cross_y + cross_z*cross_z)
angle = atan2(det, dot)

edited Apr 18 '23 at 11:34

Peter Mortensen

30,738
21
105
131

answered May 27 '21 at 09:17

Carlos Borau

1,433
1
24
34

2

thx. Linking to the respective answer or at least citing the user would be a nice thing to add – N4ppeL Sep 01 '21 at 12:54
This will cause a compilation error due to the EN DASHes. Something like "`someFile.c:44:1: error: stray ‘\342’ in program. someFile.c:44:1: error: stray ‘\200’ in program. someFile.c:44:1: error: stray ‘\223’ in program`". The [EN DASH](https://www.utf8-chartable.de/unicode-utf8-table.pl?start=8192&number=128) (Unicode code point U+2013) is can be searched for by `\x{2013}` in regular expression mode in most text editors – Peter Mortensen Apr 18 '23 at 11:51
cont' - The three numbers are the octal numbers for the [UTF-8](https://en.wikipedia.org/wiki/UTF-8) byte sequence for U+2013 (hexadecimal 0xE2 0x80 0x93) – Peter Mortensen Apr 18 '23 at 11:58
The canonical for those kind of errors is *[Compilation error: stray ‘\302’ in program, etc](https://stackoverflow.com/questions/19198332/compilation-error-stray-302-in-program-etc)*. – Peter Mortensen Apr 18 '23 at 12:09
Where was this copied from? Some web page? – Peter Mortensen Apr 18 '23 at 12:09

score 2 · Answer 5 · edited Apr 18 '23 at 10:37

2

The scalar (dot) product of two vectors lets you get the cosine of the angle between them.

To get the 'direction' of the angle, you should also calculate the cross product. It will let you check (via the z coordinate) of the angle is clockwise or not (i.e., should you extract it from 360 degrees or not).

edited Apr 18 '23 at 10:37

Peter Mortensen

30,738
21
105
131

answered Dec 28 '12 at 09:04

Nickolay Olshevsky

13,706
1
34
48

1

Even this is correct it is what I want to avoid - to compute some value and determine if the computed value represents my angle or my angle's complement. – Mircea Ispas Dec 28 '12 at 09:13
I want to know if this is possible:) Why to use some ineficient way of doing things if there is (maybe!) a better way. If there is no better way I will the "standard" thing, but it's always good to ask for better! – Mircea Ispas Dec 28 '12 at 13:37
Actually, standard ways not always efficient ) – Nickolay Olshevsky Dec 28 '12 at 19:52
@NickolayOlshevsky What do you mean exactly by *check via z coordinate*, how can I go about doing this? – Ogen Apr 12 '14 at 14:36
You should check the sign of z coordinate, as far as I remember. – Nickolay Olshevsky Apr 13 '14 at 15:36

score 2 · Answer 6 · edited Apr 18 '23 at 11:19

For a two-dimensional method, you could use the law of cosines and the "direction" method.

To calculate the angle of segment P3:P1 sweeping clockwise to segment P3:P2.

    P1     P2

        P3

    double d = direction(x3, y3, x2, y2, x1, y1);

    // c
    int d1d3 = distanceSqEucl(x1, y1, x3, y3);

    // b
    int d2d3 = distanceSqEucl(x2, y2, x3, y3);

    // a
    int d1d2 = distanceSqEucl(x1, y1, x2, y2);

    //cosine A = (b^2 + c^2 - a^2)/2bc
    double cosA = (d1d3 + d2d3 - d1d2)
        / (2 * Math.sqrt(d1d3 * d2d3));

    double angleA = Math.acos(cosA);

    if (d > 0) {
        angleA = 2.*Math.PI - angleA;
    }

This has the same number of transcendental operations as suggestions above and only one more or so floating point operation.

The methods it uses are:

 public int distanceSqEucl(int x1, int y1,
    int x2, int y2) {

    int diffX = x1 - x2;
    int diffY = y1 - y2;
    return (diffX * diffX + diffY * diffY);
}

public int direction(int x1, int y1, int x2, int y2,
    int x3, int y3) {

    int d = ((x2 - x1)*(y3 - y1)) - ((y2 - y1)*(x3 - x1));

    return d;
}

score 0 · Answer 7 · edited Apr 18 '23 at 11:00

If by "direct way" you mean avoiding the if statement, then I don't think there is a really general solution.

However, if your specific problem would allow losing some precision in angle discretization and you are ok with losing some time in type conversions, you can map the [-pi,pi] allowed range of phi angle onto the allowed range of some signed integer type. Then you would get the complementarity for free. However, I didn't really use this trick in practice. Most likely, the expense of float-to-integer and integer-to-float conversions would outweigh any benefit of the directness. It's better to set your priorities on writing autovectorizable or parallelizable code when this angle computation is done a lot.

Also, if your problem details are such that there is a definite more likely outcome for the angle direction, then you can use compilers' builtin functions to supply this information to the compiler, so it can optimize the branching more efficiently. E.g., in case of GCC, that's __builtin_expect function. It's somewhat more handy to use when you wrap it into such likely and unlikely macros (like in the Linux kernel):

#define likely(x)      __builtin_expect(!!(x), 1)
#define unlikely(x)    __builtin_expect(!!(x), 0)

score 0 · Answer 8 · edited Apr 18 '23 at 12:04

For the two-dimensional case, atan2 can easily calculate the angle between a (1, 0) vector (the x-axis) and one of your vectors.

The formula is:

Atan2(y, x)

So you can easily calculate the difference of the two angles relative to the x-axis:

angle = -(atan2(y2, x2) - atan2(y1, x1))

Why is it not used as default solution? atan2 is not efficient enough. The solution from the top answer is better. Tests on C# showed that this method has 19.6% less performance (100 000 000 iterations).

It's not critical, but unpleasant.

So, other information that can be useful:

The smallest angle between outer and inner in degrees:

abs(angle * 180 / PI)

Full angle in degrees:

angle = angle * 180 / PI
angle = angle > 0 ? angle : 360 - angle

or

angle = angle * 180 / PI
if (angle < 0)
    angle = 360 - angle;

score -1 · Answer 9 · edited Apr 18 '23 at 11:25

-1

A formula for clockwise angle, two-dimensional case, between two vectors, (xa,ya) and (xb,yb).

Angle(vec.a-vec,b) =
  pi()/2*((1 + sign(ya))*
  (1 - sign(xa^2)) - (1 + sign(yb))*
  (1 - sign(xb^2))) + pi()/4*
  ((2 + sign(ya))*sign(xa) - (2 + sign(yb))*
  sign(xb)) + sign(xa*ya)*
  atan((abs(ya) - abs(xa))/(abs(ya) + abs(xa))) - sign(xb*yb)*
  atan((abs(yb) - abs(xb))/(abs(yb) + abs(xb)))

edited Apr 18 '23 at 11:25

Peter Mortensen

30,738
21
105
131

answered Nov 23 '18 at 06:44

theodore panagos

77
1

1

what the language is it? – gdbdable Jun 17 '21 at 13:23
Why *"`vec.a-vec,b`"* (seems strangely asymmetrical)? Do you mean *"`vec.a-vec.b`"*? – Peter Mortensen Apr 18 '23 at 11:27
OK, the OP has [left the building](https://www.explainxkcd.com/wiki/index.php/979:_Wisdom_of_the_Ancients): *"Last seen more than 3 years ago"* – Peter Mortensen Apr 18 '23 at 11:30

score -4 · Answer 10 · edited Apr 18 '23 at 11:31

-4

Just copy and paste this:

angle = (acos((v1.x * v2.x + v1.y * v2.y)/((sqrt(v1.x*v1.x + v1.y*v1.y) * sqrt(v2.x*v2.x + v2.y*v2.y))))/pi*180);

edited Apr 18 '23 at 11:31

Peter Mortensen

30,738
21
105
131

answered Jun 29 '19 at 13:58

Red

1
1

7

Although this code snippet may answer the question, including an explanation of *why* and *how* it helps solve the problem improves the quality and longevity of your answer, especially regarding older questions like this. See ["How do I write a good answer?"](https://stackoverflow.com/help/how-to-answer). – slothiful Jun 29 '19 at 14:34

Direct way of computing the clockwise angle between two vectors

10 Answers10

2D case

3D case

Plane embedded in 3D

As triple product

Linked