86

I want to find out the clockwise angle between 2 vectors(2D, 3D).

The clasic 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.

Do you know a direct way of computing clockwise angle?

Mircea Ispas
  • 18,498
  • 26
  • 109
  • 202
  • 6
    Why not use `std::atan2()`? –  Dec 28 '12 at 08:53
  • 2
    How do you define "clockwise angle" for vectors in 3D? – Martin R Dec 28 '12 at 09:13
  • @H2CO3 This seems the best solution for 2D angles. – Mircea Ispas Dec 28 '12 at 09:14
  • @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
  • 4
    @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

9 Answers9

210

2D case

Just like the dot product is proportional to the cosine of the angle, the determinant is proprortional 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 counter-clockwise 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))

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.

MvG
  • 51,562
  • 13
  • 126
  • 251
  • 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
  • 11
    @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
  • Or that. :D Thanks. More elegant your way. – rbaleksandar Dec 13 '16 at 21: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
5

To compute angle you just need to call atan2(v1.s_cross(v2), v1.dot(v2)) for 2D case. Where s_cross is scalar analogue of cross production (signed area of parallelogram). For 2D case that would be wedge production. For 3D case you need to define clockwise rotation because from one side of plane clockwise is one direction, from other side of plane is another direction =)

Edit: this is counter clockwise angle, clockwise angle is just opposite

kassak
  • 3,631
  • 23
  • 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
  • 1
    @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
4

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). I'm posting it on the off chance that others find it helpful.

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.

Community
  • 1
  • 1
sircolinton
  • 5,486
  • 3
  • 16
  • 17
2

Scalar (dot) product of two vectors lets you get the cosinus 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 z coordinate) is angle is clockwise or not (i.e. should you extract it from 360 degrees or not).

Nickolay Olshevsky
  • 12,356
  • 1
  • 28
  • 44
1

For a 2D 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;
}
nichole
  • 111
  • 1
  • 6
0

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 loosing some precision in angle discretization and you are ok with loosing 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 linux kernel):

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

Since one of the simplest and most elegant solutions is hidden in one 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 3D 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)
Carlos Borau
  • 1,334
  • 1
  • 22
  • 29
-1

A formula for clockwise angle,2D case, between 2 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)))
-3

just copy & 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);

you're welcome ;-)

Red
  • 1
  • 1
  • 5
    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