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?

Sprage

Why not use std::atan2()? – Ayers

How do you define "clockwise angle" for vectors in 3D? – Germinal

@H2CO3 This seems the best solution for 2D angles. – Sprage

@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" – Sprage

if statement is one of the language primitives. You don't have much left if you can't use them. – Craniometer

@Felics read my answer. Clockwise direction is not well-defined in 3d. It is planar term – Swap

@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 – Sprage

@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." – Germinal

@Felics Also, I should note that you could not define 3D clockwise angle continuously because of hairy ball theorem 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 – Swap

Answer 1

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.

Range 0 – 360°

Most atan2 implementations will return an angle from [-π, π] in radians, which is [-180°, 180°] in degrees. If you need positive angles [0, 2π] or [0°, 360°] instead, you can just add 2π to any negative result you get. Or you can avoid the case distinction and use atan2(-det, -dot) + π unconditionally. If you are in a rare setup where you need the opposite correction, i.e. atan2 returns non-negative [0, 2π] and you need signed angles from [-π, π] instead, use atan2(-det, -dot) - π. This trick is actually not specific to this question here, but can be applied in most cases where atan2 gets used. Remember to check whether your atan2 deals in degrees or radians, and convert between these as needed.

Meister

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. – Dvinsk

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. – Console

@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°. – Meister

MvG's 2D answer gives the counter closewise angle instead of a closewise one, doesn't it? – Storekeeper

@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. – Meister

@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. – Storekeeper

@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, 2]. From that to the atan2 invocation is essentially a single step. The other sections are essentially variants of this theme. – Meister

Is it just me, or is the triple product method by fast the most performant one? – Hibbler

@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 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. – Meister

@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. – Hibbler

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. – Roa

It appears that your answer for 2D is incorrect, please you elaborate how this worked? – Psychophysics

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) – Eakins

@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. – Meister

@DonHatch could you elaborate or link to a separate discussion? – Adnopoz

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 – Adnopoz

@Adnopoz For more about poor behavior of of acos of dot product, try this (question asked on 2 different sites, with different answers and references): math.stackexchange.com/questions/1143354/… scicomp.stackexchange.com/questions/27689/… – Roa

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 – Worried

Answer 2

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.

Fructose

Answer 3

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.

Swap

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) – Sprage

@Felics In 2D cross production often means wedge production 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 – Swap

@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). – Germinal

@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. – Germinal

@Martin R thanks for good advice. I've made some corrections to make meaning of formula clearer – Swap

your formulas doesn't work for example for vectors (-1, 0) and (0, 1) – Schopenhauerism

Answer 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)

Holarctic

thx. Linking to the respective answer or at least citing the user would be a nice thing to add – Adnopoz

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 (Unicode code point U+2013) is can be searched for by \x{2013} in regular expression mode in most text editors – Punt

cont' - The three numbers are the octal numbers for the UTF-8 byte sequence for U+2013 (hexadecimal 0xE2 0x80 0x93) – Punt

The canonical for those kind of errors is Compilation error: stray ‘\302’ in program, etc. – Punt

Where was this copied from? Some web page? – Punt

Answer 5

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).

Vitascope

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. – Sprage

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! – Sprage

Actually, standard ways not always efficient ) – Vitascope

@NickolayOlshevsky What do you mean exactly by check via z coordinate, how can I go about doing this? – Melise

You should check the sign of z coordinate, as far as I remember. – Vitascope

Answer 6

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

Karakoram

Answer 7

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)

Raggedy

Answer 8

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;

Cooperage

Answer 9

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)))

Zadazadack

what the language is it? – Tait

Why "vec.a-vec,b" (seems strangely asymmetrical)? Do you mean "vec.a-vec.b"? – Punt

OK, the OP has left the building: "Last seen more than 3 years ago" – Punt

Answer 10

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);

Underscore

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?". – Jun